Perfect transmission at oblique incidence by trigonal warping in graphene P-N junctions
Shu-Hui Zhang Wen Yang
cond-mat.mes-hall



shuhuizhang@mail.buct.edu.cn

wenyang@csrc.ac.cn
^1College of Science, Beijing University of Chemical Technology, Beijing,
100029, China
^2Beijing Computational Science Research Center, Beijing 100193, China


We develop an analytical mode-matching technique for the tight-binding model
to describe electron transport across graphene P-N junctions. This method
shares the simplicity of the conventional mode-matching technique for the
low-energy continuum model and the accuracy of the tight-binding model over a
wide range of energies. It further reveals an interesting phenomenon on a
sharp P-N junction: the disappearance of the well-known Klein tunneling (i.e.,
perfect transmission) at normal incidence and the appearance of perfect
transmission at oblique incidence due to trigonal warping at energies beyond
the linear Dirac regime. We show that this phenomenon arises from the
conservation of a generalized pseudospin in the tight-binding model. We expect
this effect to be experimentally observable in graphene and other Dirac
fermions systems, such as the surface of three-dimensional topological insulators.



72.10.Bg, 73.40.Lq, 72.80.Vp, 73.23.Ad
Perfect transmission at oblique incidence by trigonal warping in graphene P-N junctions
    Wen Yang^2
    November 8, 2016
=======================================================================================






§ INTRODUCTION




Klein tunneling <cit.>, the unimpeded penetration of relativistic
particles regardless of the height and width of potential barriers, is an exotic
effect compared with the exponential-decaying transmission of nonrelativistic
particles <cit.>. In 2006, the seminal theoretical work of
Katsnelson et al. <cit.> brought about the
possibility of demonstrating Klein tunneling across electrostatic junctions in
graphene. This proposal has stimulated widespread theoretical
<cit.>
and experimental
<cit.>
interest and shed light on possible electronic applications
<cit.>.


Graphene is a 2D layer of carbon atoms on a honeycomb lattice [Fig.
<ref>(a)]. The conduction band and the valence band touch each
other at six Dirac points, but only two of them are inequivalent, as denoted
by 𝐊 and 𝐊^' in the right panel of Fig.
<ref>(a). According to group theory <cit.>,
the D_6h point group symmetry of the honeycomb lattice determines
the global hexagonal symmetry of the graphene energy band over the
first Brillouin zone, while the local D_3h symmetry of each
valley determines the local triangular symmetry of the energy band in this
valley. Close to the Dirac point, say 𝐊, the energy dispersion as a
function of the reduced wave vector 𝐪≡𝐤
-𝐊 is linear and isotropic. Away from the Dirac point, however, the
D_3h local symmetry becomes important and the constant energy
contour approaches a regular triangle with a side length 2π/3 [see Fig.
<ref>(b)]. This is commonly referred to as trigonal
warping, which increases rapidly with energy [Fig. <ref>(c)]. Note
that throughout this work we use the carbon-carbon bond length
a=0.142 nm as the unit of length and the nearest-neighbor hopping energy
t_0=2.7 eV as the unit of energy <cit.>.



Trigonal warping is a well-known feature of the graphene energy band beyond
the linear regime <cit.> and its influence on the electron
transport is receiving growing interest, including the observation of broken
chirality by trigonal warping in transport measurements
<cit.>, the influence
<cit.> of trigonal warping on the famous Veselago
focusing across graphene P-N junctions
<cit.>, and potential
applications of trigonal warping for producing valley-polarized electrons in
N-P-N junction <cit.>, double barriers <cit.>,
and other junctions <cit.>. To incorporate trigonal warping,
the simplest method <cit.> is to introduce
nonlinear corrections to the widely used low-energy linear continuum model.
The nearest-neighbor tight-binding model captures the hexagonal lattice
symmetry and hence is commonly used
<cit.> to study the
trigonal warping effect over a wide energy range. At energies below 1 eV, this model gives very accurate energy bands, but its deviation from the ab initio calculation becomes significant at energies approaching the Van Hove singularity at ∼3 eV <cit.>.
However, transport calculations in the tight-binding model are usually based on
the recursive Green's function method
<cit.> and hence are
numerical. Moreover, previous studies on Klein tunneling mainly focus on the
low-energy linear regime of graphene, leaving the effect of trigonal warping
largely unexplored. In particular, it is well known that the Klein tunneling
(i.e., perfect transmission) at normal incidence originates from the
conservation of the pseudospin <cit.>.
However, this simple physical picture is based on the linear continuum model
and hence is limited to the low-energy linear Dirac regime. It is not clear
whether a similar physical picture exists in the tight-binding model and over
a wide range of energies.

In this paper, we develop an analytical mode-matching technique in the
tight-binding model for electron transport across graphene P-N junctions. The
key is to introduce a titled coordinate system, reduce the 2D junction into a
1D chain, and then perform mode-matching at the P-N interface, in a similar
way to the mode matching in the continuum model
<cit.>. As a result, this method shares
the simplicity of the mode-matching method in the low-energy continuum model
and the accuracy of the tight-binding model over a wide range of energies. To
focus on trigonal warping, we consider the electron transmission across a
sharp P-N interface along the zigzag direction in the energy range
E∈[-1,+1], where the intervalley scattering is absent
<cit.>; but trigonal warping is significant and can be
described reasonably by the tight-binding model. The results show that beyond
the linear regime, the Klein tunneling at normal incidence becomes imperfect;
i.e., finite backscattering occurs. Interestingly, we find that perfect
transmission is still possible, but the critical incident angle for perfect
transmission deviates from zero, and the deviation increases with increasing
trigonal warping. We introduce the concept of a generalized pseudospin in the
tight-binding model and show that its conservation across the P-N interface is
responsible for the perfect transmission at oblique incidence. This
generalizes the well-known pseudospin picture for perfect transmission,
previously limited to the linear Dirac regime, to a wide energy range
E∈[-1,+1]. The continuum model with a second-order nonlinear correction
fails to describe this phenomenon quantitatively. We expect this phenomenon to
be experimentally observable.

The rest of this paper is organized as follows. In Sec. II, we introduce the
tilted coordinates, present the analytical mode-matching method, and introduce
the generalized pseudospin picture for perfect transmission. In Sec. III, we
discuss the perfect transmission at oblique incidence and the failure of the
continuum model. Finally, we give a brief conclusion in Sec. IV.



§ MODE-MATCHING TECHNIQUE IN TILTED COORDINATES




For specificity, we consider a graphene P-N junction with a sharp, zigzag
interface separating the left region and the right region [Fig. <ref>
(a)], leaving the more general cases to the end of this section. The
tight-binding Hamiltonian of the junction is the sum of the
(nearest-neighbor) tight-binding Hamiltonian for uniform graphene
<cit.> and the on-site junction potential, which takes a
constant value V_L (V_R) for all the carbon sites in
the left (right) region [see the inset of Fig. <ref>(a)]. The Fermi
level E_F, which can be tuned by electric gating, determines the doping
level in the left (right) region as ε_L≡
E_F-V_L (ε_R≡ E_F-V_R),
where a positive (negative) doping level means electron or N (hole or P)
doping. Here we consider a P-N junction with ε_L>0 and
ε_R<0; i.e., the left region is N doped and the right
region is P doped.

In addition to the ortho-normalized basis vectors (𝐞_x
,𝐞_y) of the conventional Cartesian coordinate system, we
introduce a tilted coordinate system characterized by the nonorthogonal basis
vectors (𝐞_X,𝐞_Y) [see Fig. <ref>(b)]. The
tilted vectors 𝐞_X and 𝐞_Y are actually the two
primitive vectors for the Bravais lattice of uniform graphene and are
connected to (𝐞_x,𝐞_y) by

    {[ 𝐞_X=√(3)(√(3)2𝐞_x-1
    2𝐞_y),;               𝐞_Y=√(3)𝐞_y, ].  ⇔{[ 𝐞_x=2𝐞_X+𝐞_Y3,;   𝐞_y=𝐞_Y√(3). ].

The Cartesian components k_x≡𝐤·𝐞_x and
k_y≡𝐤·𝐞_y of a wave vector 𝐤 in
the reciprocal space are connected to the tilted components k_X
≡𝐤·𝐞_X and k_Y≡𝐤·𝐞_Y via

    {[ k_X=√(3)(√(3)2k_x-12k_y),;              k_Y=√(3)k_y, ].  ⇔{[ k_x=2k_X+k_Y3,;   k_y=k_Y√(3). ].

For a position vector 𝐑 in the real space, we define the tilted
components R_X,R_Y through the expansion

    𝐑≡ R_X𝐞_X+R_Y𝐞_Y⇒{[        R_X=2R_x3,; R_Y=R_x+√(3)R_y3. ].

These definitions lead to the convenient expression 𝐤·𝐑=k_xR_x+k_yR_y=k_XR_X+k_YR_Y, although
𝐤≠ k_X𝐞_X+k_Y𝐞_Y and R_X,Y
≠𝐑·𝐞_X,Y since (𝐞_X,𝐞_Y)
are not ortho-normalized. An arbitrary vector 𝐎 can be denoted by
𝐎=(O_x,O_y)=(O_X,O_Y)_T.



 §.§ Tight-binding model for uniform graphene


Here we consider uniform graphene with ε_L
=ε_R=E_F to illustrate the usage of the tilted
coordinates and establish the relevant notations and some important concepts.

As shown in Fig. <ref>(a) or <ref>(a), the honeycomb lattice
of graphene consists of two sublattices (denoted by A and B) and each unit
cell contains two carbon atoms, one on each sublattice. The unit cell (m,n)
locates at 𝐑_m,n=m𝐞_X+n𝐞_Y=(m,n)_T, and the relative displacements of the two carbon atoms inside this unit
cell are τ_A=0 and τ_B=(1/2,-√(3)/2). The tight-binding Hamiltonian of uniform graphene is

    Ĥ_0=-∑_m,n(|m+1,n,A⟩+|m,n-1,A⟩+|m,n,A⟩)⟨
    m,n,B|-h.c.,

where |m,n,λ⟩ (λ=A,B) denotes the carbon atom on the
sublattice λ. To utilize the translational invariance of graphene with
the primitive vector 𝐞_Y along the Y axis, we make a Fourier
transform from the on-site basis {|m,n,λ⟩} to the hybrid basis
(i.e., on-site basis along the X axis and Bloch basis along the Y axis)

    |m,k_Y,λ⟩=1/√(N_Y)∑_ne^ink_Y
    |m,n,λ⟩,

where N_Y is the normalization length along the Y axis and k_Y
∈0,2π]. Under this basis Ĥ_0=∑_k_Yĥ(k_Y),
where

    ĥ(k_Y)=-∑_m[(1+e^ik_Y)|m,k_Y,A⟩+|m+1,k_Y
    ,A⟩]⟨ m,k_Y,B|-h.c..

For a given k_Y, the Hamiltonian ĥ(k_Y) describes a 1D chain
along the X axis [see Fig. <ref>(b)], with two orbitals A and B in
each unit cell. The 2×2 on-site energy of each unit cell is

    𝐡(k_Y)=-[
    [         0  1+e^ik_Y; 1+e^-ik_Y         0 ]]

and the 2×2 hopping matrix from unit cell m+1 to m is

    𝐭≡⟨ m|ĥ(k_Y)|m+1⟩=
    [  0  0; -1  0 ]
    .

So the Fourier transform reduces the 2D graphene into decoupled 1D chains
[Fig. <ref>(b)] parametrized by different k_Y's. 

Next, to utilize the translational invariance along the X axis with the
primitive vector 𝐞_X, we make another Fourier transform along X
to get the full Bloch basis:

    |𝐤,λ⟩≡1/√(N_X)∑_me^imk_X
    |m,k_Y,λ⟩=1/√(N)∑_m,ne^i𝐑_m,n
    ·𝐤|m,n,λ⟩

parameterized by the wave vector 𝐤=(k_X,k_Y)_T. This
transformation further decouples the lattice degree of freedom along the X
axis and leads to

    Ĥ_0=∑_𝐤[|𝐤,A⟩,|𝐤,B⟩
    ]𝐇_0(𝐤)
    [ ⟨𝐤,A|; ⟨𝐤,B| ]
    ,

where

    𝐇_0(𝐤)≡
    [     0  f(𝐤); f(-𝐤)     0 ]

is the 2×2 tight-binding Hamiltonian and

    f(𝐤)=-(1+e^ik_Y
    +e^-ik_X).


The eigenenergies are E^(±)(𝐤)=±
|f(𝐤)|, where

    | f(𝐤)|     =√(3+2[cos k_X+cos k_Y
    +cos(k_X+k_Y)])
         =√(3+2cos(√(3)k_y)+4cos3k_x/2cos√(3)k_y
    /2).

The corresponding eigenstates are |Φ^(±)(𝐤)⟩=[1,±
f^∗(𝐤)/| f(𝐤)| ]^T. Note that we
normalize every spinor to √(2) throughout this work.

In the Cartesian coordinates, the first Brillouin zone is a regular hexagon,
as shown in Fig. <ref>(a). In the tilted coordinates, the first
Brillouin zone is a square k_X,k_Y∈-π,π], as shown in Fig.
<ref>(b). The two inequivalent Dirac points are

    𝐊   ≡(  2π/3,2π/3√(3))
    =(  2π/3,2π/3)  _T,
    𝐊^'   ≡(  2π/3,-2π/3√(3)
    )  =(  4π/3,-2π/3)  _T
    ⟶(  -2π/3,-2π/3)  _T,

where in the last step of the second line we have shifted 𝐊^' by a reciprocal vector along the X axis.

The continuum model Hamiltonian is obtained from 𝐇_0(𝐤)
by considering 𝐤 in a given valley, say 𝐊, defining the
reduced wave vector 𝐪≡𝐤-𝐊, and expanding
f(𝐤) into Taylor series with respect to 𝐪. For example,
expanding f(𝐤) up to the first order of 𝐪 gives the
widely used linear continuum model. For the 𝐊 valley, we have
f(𝐤)≈ v_F|𝐪|e^-i(π/6-θ_𝐪) and
hence the linear continuum model 
    𝐇_𝐊(𝐪)=v_F|𝐪|
    [             0 e^-i(π/6-θ_𝐪);  e^i(π/6-θ_𝐪)             0 ]
    ,

where v_F≡3/2 is the Fermi velocity and θ_𝐪 is the
azimuth angle of 𝐪 in the Cartesian coordinates. For the
𝐊^' valley, f(𝐤)≈-v_F|𝐪
|e^i(π/6-θ_𝐪); thus

    𝐇_𝐊^'(𝐪)=-v_F|𝐪|
    [             0  e^i(π/6-θ_𝐪); e^-i(π/6-θ_𝐪)             0 ]
    .

The linear continuum model for either valley gives the same isotropic, linear
dispersion E^(±)(𝐪)=± v_F|𝐪| that totally ignores
the trigonal warping. The corresponding eigenstates are |Φ_𝐊
^(±)(𝐪)⟩=[1,± e^i(π/6-θ_𝐪)]^T for
the 𝐊 valley and |Φ_𝐊^'^(±)
(𝐪)⟩=[1,∓ e^-i(π/6-θ_𝐪)]^T for the
𝐊^' valley. By expanding f(𝐤) into high orders of
𝐪, trigonal warping can be included into the continuum model.

Now we discuss a distinguishing feature of the tight-binding model compared
with all the continuum models (including those with high-order corrections for
trigonal warping): the existence of “abnormal” evanescent states. Let us consider a given k_y
(or q_y) and determine the k_x (or q_x) of all the eigenstates on
the Fermi level. For q_y<q_F≡|E_F|/v_F, the linear continuum
model of either valley gives two traveling eigenstates characterized
by the reduced wave vector(±(q_F^2-q_y^2)^1/2,q_y) in the
Cartesian coordinate. For the tight-binding model, k_x is determined by
E_F=sgn(E_F)|f(𝐤)|, where 𝐤≡
(k_x,k_y) and sgn(x)≡+1 for x>0 and -1 for x<0. The
solutions can be visualized as the two intersection points between the
horizontal line k_y and the Fermi contour E_F=sgn
(E_F)|f(𝐤)| in the (k_x,k_y) plane; e.g., for E_F=0.75
and k_y=0.70, we obtain two traveling eigenstates [blue dots in Fig.
<ref>(a)]. When k_y approaches the Dirac point, the eigenstates from
the tight-binding model approach those from the linear continuum model.
Similarly, in the tilted coordinate system, the eigenstates with a given
k_Y on the Fermi level can also be determined as the two intersection
points between the horizontal line k_Y and the Fermi contour
E_F=sgn(E_F)|f(𝐤)| in the (k_X,k_Y)_T
plane, as shown in Fig. <ref>(b) for E_F=0.75 and k_Y=1.21.

Surprisingly, in addition to these two “normal” eigenstates, the tight-binding model also has two
extra, “abnormal” eigenstates. To make this
clear, we follow Ref. ZhangPRB2017 and solve the 1D
Schrödinger equation

    -𝐭^†|Φ(m-1)⟩+(E-𝐡(k_Y))
    |Φ(m)⟩-𝐭|Φ(m+1)⟩=0

for the 1D chain along the X axis under the Bloch condition |Φ
(m)⟩=e^imk_X|Φ⟩, where k_X may be either real or
complex. Equation (<ref>) can be written as

    𝐇_0(𝐤)|Φ⟩=E|Φ⟩,

where 𝐤≡(k_X,k_Y)_T=(k_x,k_y) is the wave
vector and 𝐇_0(𝐤) is the 2×2 Hamiltonian of the
tight-binding model [see Eq. (<ref>)]. For a given, real
wave vector 𝐤, Eq. (<ref>) reproduces the energy dispersions
and eigenstates of uniform graphene. Here we need to find all the
eigenstates with a given k_Y on the Fermi level. For this purpose, we let
λ≡ e^ik_X be the propagation phase along the +X axis by one
unit cell, set E=E_F, and rewrite Eq. (<ref>) as

    -𝐭^†+λ(E_F-𝐡(k_Y))
    -λ^2𝐭]|Φ⟩=0,

from which we obtain four solutions for λ (and hence k_X) and
|Φ⟩. Two solutions correspond to the “normal” eigenstates and can be explicitly obtained by
rewriting Eq. (<ref>) as an explicit quadratic equation for λ:

    (1+e^ik_Y)λ^2+[(3+2cos k_Y)-E_F^2]λ+(1+e^-ik_Y
    )=0.

The other two solutions are “abnormal”
 evanescent eigenstates, including a right-going one

    λ_0=0,   |Φ_0⟩=
    [ √(2);    0 ]

and a left-going one

    λ_∞=∞,  |Φ_∞⟩=
    [    0; √(2) ]
    .



Finally we discuss the physical meaning of these “abnormal” evanescent eigenstates. The wave function of the
right-going one [Eq. (<ref>)] starting from a unit cell m_0 is
defined in the half plane on the right of m_0 only (i.e., m≥ m_0):

    |Φ_0(m|m_0)⟩=λ_0^m-m_0|Φ_0⟩={[       |Φ_0⟩    (m=m_0),;           0 (m≥ m_0+1). ].

The wave function of the left-going one [Eq. (<ref>] starting from a
unit cell m_0 is defined in the half plane on the left of m_0 only
(i.e., for m≤ m_0):

    |Φ_∞(m|m_0)⟩=λ_∞^m-m_0|Φ_∞
    ⟩={[       |Φ_∞⟩    (m=m_0),;           0 (m≤ m_0-1). ].

Using 𝐭^†|Φ_0⟩=𝐭|Φ_∞⟩
=0, we can readily verify that |Φ_0(m|m_0)⟩ (|Φ_∞(m|m_0)⟩) indeed satisfies the Schrödinger equation Eq.
(<ref>) with m≥ m_0+1 (m≤ m_0-1), so it is indeed an
eigenstate on the Fermi level of uniform graphene, although it only exists in
the half plane m≥ m_0 (m≤ m_0). These “abnormal” evanescent eigenstates originate from the
singularity of the 2×2 hopping matrix between neighboring unit cells.
For example, the wave function |Φ_0(m|m_0)⟩ is nonzero on the
A site of the unit cell m=m_0, but this site cannot hop to the
neighboring unit cell m_0+1 on its right (see Fig. <ref>), so
|Φ_0(m|m_0)⟩ vanishes for m≥ m_0+1. Similarly, the wave
function |Φ_∞(m|m_0)⟩ is nonzero on the B site of the
unit cell m=m_0, but this site cannot hop to the neighboring unit cell
m_0-1 on its left, so |Φ_∞(m|m_0)⟩ vanishes for m≤
m_0-1.

The “abnormal” evanescent eigenstates do
not exist in an infinite uniform graphene, but they do appear near the
interface of the graphene junction. When considering the transmission of an
incident traveling wave, it is important to include these “abnormal” eigenstates.



 §.§ Mode-matching across P-N junctions


Due to the translational invariance along the Y axis, the scattering of
incident states with different k_Y's is decoupled, so we need only
consider the 1D chain with a given k_Y∈0,2π]. As shown
in Fig. <ref>, each unit cell contains two orbitals A and B. The
2×2 on-site energy of unit cell m is V_m𝐈_2×
2+𝐡(k_Y), where V_m is the on-site potential: V_m
=V_L for m≤0 and V_m=V_R for m≥1 [see
Fig. <ref>(b)]. The 2×2 hopping from unit cell m+1 to m is
𝐭 [Eq. (<ref>)]. The scattering state on the Fermi level E_F
satisfies the Schrödinger equation for the 1D chain:

    -𝐭^†|Ψ(m-1)⟩+(E_F-𝐡(k_Y)-V_m
    )|Ψ(m)⟩-𝐭|Ψ(m+1)⟩=0.

Next we follow similar procedures to the usual mode-matching technique for
continuum models.

As the first step, we need to find all the left-going and right-going
eigenstates of each uniform region on the Fermi level. Here right-going
(left-going) eigenstates refer to both evanescent eigenstates that decay along
the +X (-X) axis and traveling eigenstates whose group velocity ∂
E(𝐤)/∂ k_X along the +X axis of the tilted coordinates,
or equivalently the group velocity ∂ E(𝐤)/∂ k_x
along the +x axis of the Cartesian coordinates, is positive (negative)
<cit.>, where we have used ∂/∂ k_X=(2/3)∂
/∂ k_x according to Eq. (<ref>). These eigenstates can be
obtained by exactly the same way as the previous subsection, with the only
difference being E_F→ε_L=E_F
-V_L (E_F→ε_R=E_F
-V_R) for the left (right) region. For the region s
(=L or R), we obtain one right-going traveling
eigenstate and one left-going traveling eigenstate as the two intersection
points between the given horizontal line k_Y and the Fermi contour
|ε_s|=|f(𝐤)| in the 𝐤=(k_X,k_Y
)_T plane. For each region, we also obtain two “abnormal” evanescent eigenstates: Eqs. (<ref>) and
(<ref>). For clarity, we use λ_α≡ e^ik_X,α
(α=i,r,t) for the propagation phases of different traveling eigenstates
over one unit cell along the X axis: α=i (α=r) for the
right-going (left-going) eigenstate in the left N region and α=t for
the right-going traveling eigenstate in the right P region. The condition that
𝐤_i,𝐤_r, and 𝐤_t lie on the Fermi contour
amounts to



    |ε_L|     =| f(𝐤
    _i)| =| f(𝐤_r)| ,
    |ε_R|     =| f(𝐤
    _t)|,

where f(𝐤) is defined in Eq. (<ref>). The wave vectors and spinors of the incident, reflection, and transmission
eigenstates are 𝐤_α≡(k_X,α,k_Y)_T
and |Φ_α⟩=[1,e^iφ_α]^T, where




    e^iφ_i   ≡f^∗(𝐤_i)/|
    f(𝐤_i)|,
    
    e^iφ_r   ≡f^∗(𝐤_r)/|
    f(𝐤_r)|,
    
    e^iφ_t   ≡-f^∗(𝐤_t)/|
    f(𝐤_t)|.



As the second step, we consider a right-going incident traveling wave in the
left N region


    |Φ_i(m)⟩=e^i(m-1)k_X,i|Φ_i⟩

and calculate the scattering state. By solving Eq. (<ref>) for
m=-∞,⋯,0, we obtain the scattering state in the left N region as
the sum of the incident wave and the reflection wave:

    |Ψ(m)⟩=|Φ_i(m)⟩+|Φ_r(m)⟩  (m≤1),

where

    |Φ_r(m)⟩=re^i(m-1)k_X,r|Φ_r⟩+r̃|Φ_∞(m|1)⟩

is a linear combination of the left-going traveling eigenstate and the
left-going “abnormal” evanescent eigenstate
in the left N region. By solving Eq. (<ref>) for m=1,2,⋯
,+∞, we obtain the scattering state in the right P region as the
transmission wave

    |Ψ(m)⟩=|Φ_t(m)⟩  (m≥0),

which is a linear combination of the right-going traveling eigenstates and the
right-going “abnormal” evanescent
eigenstate in the right P region:

    |Φ_t(m)⟩=te^imk_X,t|Φ_t⟩+t̃|Φ
    _0(m|0)⟩.

There are four unknown coefficients, including two reflection coefficients
r,r̃ and two transmission coefficients t,t̃.

As the final step, the overlap of Eqs. (<ref>) and (<ref>) at the
interface region m=0,1 gives the continuity condition |Φ_i
(m)⟩+|Φ_r(m)⟩=|Φ_t(m)⟩ for m=0,1 [see the
dashed square in Fig. <ref>(b)]. The equation at m=1 gives



    1+r       =e^ik_X,tt,
    
    e^iφ_i+e^iφ_rr+√(2)r̃     =e^ik_X,t
    e^iφ_tt.

The equation at m=0 gives




    e^-ik_X,i+e^-ik_X,rr       =t+√(2)t̃,
    
    e^-ik_X,ie^iφ_i+e^-ik_X,re^iφ_rr     
    =e^iφ_tt.

These four equations uniquely determine the four unknown coefficients
r,r̃,t,t̃. Remarkably, Eqs. (<ref>) and (<ref>)
form a closed set of equations for the traveling eigenstates alone,
from which we obtain




    r       =-e^i(φ_i-k_X,i)-e^i(φ_t-k_X,t)
    /e^i(φ_r-k_X,r)-e^i(φ_t-k_X,t),
    
    t       =e^-ik_X,t(1+r).

Moreover, the “abnormal” evanescent
eigenstates decay to zero over a single unit cell and hence do not contribute
to the reflection and transmission waves away from the interface; e.g.,
|Φ_r(m)⟩=re^i(m-1)k_X,r|Φ_r⟩ at m≤0 only
contains the left-going traveling eigenstate, while |Φ_t(m)⟩
=te^imk_X,t|Φ_t⟩ at m≥1 only contain the right-going
traveling eigenstates. Therefore, the presence of the “abnormal” evanescent eigenstates has no effect on
the scattering properties of a sharp interface along the zigzag axis. The
transmission probability of an incident traveling eigenstate with wave vector
k_Y on the Fermi level is


    T(k_Y)=1-|r(k_Y)|^2.

Perfect transmission corresponds to T(k_Y)=1, which amounts to
r(k_Y)=0 or equivalently

    e^i(φ_i-k_X,i)=e^i(φ_t-k_X,t)≠ e^i(φ
    _r-k_X,r).





Next we discuss the mode-matching condition at the interface m=0,1 (see Fig.
<ref>) in more detail to explain why the “abnormal” evanescent eigenstates do not affect the
transmission of the traveling eigenstate. The key is that the 2×1
spinor |Φ⟩=[a,b]^T means that the wave amplitude at the A site
is a, while that at the B site is b. Therefore, the matching of the
2×1 spinor wave function at m=0 and 1 amounts to the matching of
the wave amplitudes at the four sites in Fig. <ref>: Eq.
(<ref>) [Eq. (<ref>)] for the matching at the A (B) site of
unit cell m=1 and Eq. (<ref>) [Eq. (<ref>)] for the matching at
the A (B) site of m=0. However, the left-going ideal evanescent state
|Φ_∞(m|1)⟩ contained in |Φ_r(m)⟩ is nonzero
only on the B site of m=1, while the right-going ideal evanescent
eigenstate |Φ_0(m|0)⟩ contained in |Φ_t(m)⟩ is
nonzero only on the A site of m=0. Therefore, the matching conditions on
the A site of m=1 [Eq. (<ref>)] and on the B site of m=0 [Eq.
(<ref>)], as enclosed by the red box in Fig. <ref>, do not
contain any “abnormal” evanescent waves.
The matching of the traveling waves at these two sites uniquely determines the
reflection and transmission coefficients r and t for the traveling eigenstates.



 §.§ Generalized pseudospin in tight-binding model




In the linear continuum model, the perfect Klein tunneling at normal incidence
has a physically transparent interpretation as the conservation of the
pseudospin <cit.>. In the
tight-binding model, however, such a simple physical picture for the Klein
tunneling is still lacking, due to the broken chirality of the Dirac fermions
by the trigonal warping
<cit.>.

Here we demonstrate that in the tight-binding model, the perfect transmission
over the energy range E∈[-1,+1] can always be interpreted as the
conservation of a generalized pseudospin. The key is that the
mode-matching conditions for the traveling waves at the interface form a
closed set of equations [Eqs. (<ref>) and (<ref>)], which can be
put into the form

    |u_i⟩+r|u_r⟩=te^ik_X,t|u_t⟩,

where

    |u_α⟩≡
    [              1; e^i(φ_α-k_X,α) ]  (α=i,r,t).

The perfect transmission condition [Eq. (<ref>)] leads to |u_i
⟩=|u_t⟩≠|u_r⟩. Conversely, once |u_i
⟩=|u_t⟩≠|u_r⟩, we immediately obtain r=0 and hence
perfect transmission. Thus we arrive at the following necessary and
sufficient condition for perfect transmission:

    |u_i⟩=|u_t⟩≠|u_r⟩.

By regarding the A,B sites inside the red box of Fig. <ref> as a
special unit cell across the interface, |u_i⟩, r|u_r⟩, and
te^ik_X,t|u_t⟩ become the amplitudes of the incident,
reflection, and transmission waves inside this unit cell, as shown in Fig.
<ref>(b)-(d). Therefore, Eq. (<ref>) simply expresses the
continuity of the scattering wave function inside this special unit cell. Then
we can associate each spinor with a generalized pseudospin

    σ_α≡⟨ u_α|σ̂|u_α⟩/⟨ u_α|u_α⟩
    =(cos(φ_α-k_X,α),sin(φ_α-k_X,α)),

where σ̂=(σ̂_x,σ̂_y) and σ̂_x,y are Pauli matrices. Then the perfect
transmission condition amounts to

    σ_i=σ_t≠σ_r,

i.e., conservation of the generalized pseudospin across the graphene junctions.



 §.§ Generalizations to finite-width junctions




The mode-matching technique in the tilted coordinates can also deal with
finite-width junctions. Suppose the left N region, the right P region, and
the interface region consist of the unit cells (m,n) with m≤0, m≥
N+1, and 1≤ m≤ N, respectively. By using the Fourier transform Eq.
(<ref>) from the on-site basis into the hybrid basis |m,k_Y
,α⟩, we reduce the 2D junction into a 1D chain parameterized by a
given k_Y. Each unit cell contains two orbitals A and B. The
2×2 on-site energy of the unit cell m is V_m𝐈_2×
2+𝐡(k_Y), where the on-site potential V_m=V_L in
the left N region (m≤0), V_m=V_R in the right P region
(m≥ N+1), and V_m can be arbitrary in the interface region (1≤
m≤ N), as shown in Fig. <ref>(b). The 2×2 hopping
from unit cell m+1 to m is 𝐭 [Eq. (<ref>)].

For a right-going incident wave on the Fermi level from the left N region
|Φ_i(m)⟩=e^i(m-1)k_X,i|Φ_i⟩, the scattering state
is obtained by solving the 1D Schrödinger equation Eq. (<ref>).
Solving Eq. (<ref>) with m≤0 gives

    |Ψ(m)⟩=|Φ_i(m)⟩+re^i(m-1)k_X,r|Φ_r⟩
    +r̃λ_∞^(m-1)|Φ_∞⟩  (m≤1),

where the last two terms are reflection waves. Solving Eq. (<ref>)
with m≥ N+1 gives the transmission waves

    |Ψ(m)⟩=te^i(m-N)k_X,t|Φ_X,t⟩+t̃λ
    _0^(m-N)|Φ_0⟩ (m≥ N).

The unknown variablesr,r̃,t,t̃, and |Ψ(m)⟩ (2≤
m≤ N-1) can be uniquely determined by Eq. (<ref>) with
m=1,2,⋯,N. This
is reminiscent of the mode-matching method developed by Ando
<cit.>, albeit in the tilted coordinates. For large N, analytical solutions are no longer available, and numerical calculations are necessary.



§ NUMERICAL RESULTS AND DISCUSSIONS




The mirror symmetry of the P-N junction about the x axis guarantees the
transmission probability T(k_y) to be an even function of k_y. So we
limit our discussions to one valley, say 𝐊, and define the reduced
wave vector 𝐪≡𝐤-𝐊. For specificity we
calculate the transmission probability of electrons in the 𝐊 valley
across a symmetric graphene P-N junction with ε_L
=-ε_R=V_0; i.e., the electron doping in the left N
region is equal to the hole doping in the right P region. We present our
calculation results in the conventional Cartesian coordinates (see Fig.
<ref>) to make them accessible to readers that are more familiar
with the Cartesian coordinates. For a given q_y, the wave vectors of the
incident, reflection, and transmission waves correspond to the intersection
points of the horizontal line q_y and the Fermi contour in the
(q_x,q_y) plane, as shown in Fig. <ref>. The incident wave is
characterized by either q_y or the azimuth angle θ_i of
𝐪_i, so the transmission probability is denoted by T(q_y) or
T(θ_i). For clarity we define q_y,PT (θ
_PT) as the critical momentum (angle) leading to perfect
transmission: T(q_y,PT)=1 [T(θ_PT)=1]. For
reference, the conventional continuum model with a linear dispersion gives the
transmission probability T(θ_i)=cos^2θ_i (Ref.
CheianovPRB2006) and hence q_y,PT=θ
_PT=0.



 §.§ Numerical results




At the beginning, we plot in Fig. <ref> the exact results from the
tight-binding model for the transmission probability, as calculated from Eqs.
(<ref>) and (<ref>). At very low doping V_0=0.01, perfect
transmission occurs at q_y,PT=θ_PT=0, in
agreement with the well-known Klein tunneling in the widely used linear
continuum model <cit.>. However, with increasing doping
level, q_y,PT and θ_PT exhibit larger and
larger deviations from zero, indicating perfect transmission at oblique incidence.



In Fig. <ref>, we compare the exact results from the tight-binding
model to those from the continuum models. At very low doping V_0=0.01
[Fig. <ref>(a)], all the models agree with each other and give
θ_PT=0. However, with increasing doping level, the results
from the continuum models deviate more and more from the exact results.
Actually, the deviation is appreciable even at relatively low doping
V_0=0.1 [Fig. <ref>(b)], where the trigonal warping is small:
η=3% [Fig. <ref>(c)]. Interestingly, even the continuum model
with a quadratic correction (dashed orange line) to incorporate the trigonal
warping up to the lowest order (see the Appendix) does not agree with the exact results.



Next we illustrate the generalized pseudospin picture (see Sec.
<ref>) in the tight-binding model. As shown in Fig.
<ref>, for a given q_y, the wave vectors 𝐪
_i,𝐪_r,𝐪_t of the incident, reflection, and
transmission waves are obtained as the intersection points between a
horizontal line (see the orange dashed line for an example) corresponding to
this q_y and the Fermi contour. For low doping V_0=0.01, the condition
σ_i=σ_t≠σ_r
is satisfied at q_y=0, so perfect transmission occurs at normal incidence
θ_PT=0 [Fig. <ref>(a)]. For higher doping
V_0=0.5 [Fig. <ref>(b)], the texture of the generalized
pseudospins on the Fermi contour is twisted relative to those at low doping.
In this case, the condition σ_i=σ
_t≠σ_r occurs at θ_PT≈-14.6


^∘; i.e., perfect transmission occurs for oblique incidence.



 §.§ Experimental feasibility




Ever since the initial theoretical prediction <cit.>, many experimental efforts
<cit.>
have been devoted to Klein tunneling in graphene, culminating in the
experimental demonstration of the prominent angular dependence of the
transmission probability in graphene P-N junctions
<cit.>. In particular,
the recent fabrication of high-quality graphene P-N junctions with high doping
levels <cit.> makes the high-energy transmission across
graphene P-N junctions experimentally accessible; e.g., the perfect
transmission angle θ_PT may be extracted by measuring the
angular dependence of transmission probability. To serve future experiments,
we plot the perfect transmission momentum q_y,PT and the perfect
transmission angle θ_PT as functions of the doping level
V_0 in Fig. <ref>.



Up to now, we have focused on symmetric P-N junctions with a sharp interface
and equal doping level and trigonal warping on both sides. For asymmetric P-N
junctions, the trigonal warping in the N region is not equal to that in the P
region. This may weaken the effect of perfect transmission at oblique
incidence by shifting θ_PT towards zero. A smooth junction
potential usually suppresses the transmission as the incident angle increases
<cit.> and plays an important role in the electron optics in
the graphene P-N junction <cit.>. This may hinder the
experimental observation of perfect transmission at oblique incidence.
Fortunately, a very recent experiment <cit.> shows that
atomically sharp graphene P-N junctions can be fabricated on the copper
surface. The potential difference between the P and the N regions can reach
2V_0=660 meV, corresponding to a doping level V_0=0.33 eV and doping
density 10^13 cm^-2 [see the inset of Fig. <ref>(b)]. At this
doping level, the trigonal warping is η=4% and the perfect transmission
angle is θ_PT≈-3.5


^∘. Moreover, a previous theoretical study shows that the electron-electron
interaction <cit.> generally enhances the trigonal warping.
Therefore, we expect the perfect transmission at oblique incidence to be
experimentally accessible in the near future.

Finally, we note the intensive experimental activities in simulating the
honeycomb lattices by cold atoms and optical lattices <cit.>
. In particular, the possibility of Klein tunneling in these system has been
examined
<cit.> and
observed recently <cit.>. These simulated
graphene systems may provide an alternative platform for observing the perfect
transmission at oblique incidence. It would also be interesting to explore the
effect of trigonal warping on the transmission across bilayer graphene P-N
junctions <cit.>.





§ CONCLUSIONS


We have developed an analytical mode-matching technique to study the electron
transmission across graphene P-N junctions over a wide energy range. In
contrast to the well-known Klein tunneling at normal incidence for low
energies in the linear Dirac regime, we find that at energies beyond the
linear Dirac regime, the Klein tunneling at normal incidence becomes imperfect
and trigonal warping causes perfect transmission at oblique incidence. We show
that this phenomenon arises from the conservation of a generalized pseudospin
in the tight-binding model. This generalizes the well-known pseudospin picture
for perfect Klein tunneling, previously limited to low energies in the linear
Dirac regime, to all the energy ranges. The perfect transmission at oblique
incidence cannot be described by the continuum model, even after quadratic
corrections have been introduced to incorporate trigonal warping up to the
leading order. Our work may be relevant for the applications of the graphene P-N
junction in electronics and electron optics.



§ ACKNOWLEDGEMENTS


This work was supported by the National Key R&D Program of China (Grant
No. 2017YFA0303400), the NSFC (Grants No. 11504018, No. 11774021, No.
11274036, and No. 11322542), the MOST of China (Grant No. 2014CB848700), and the NSFC
program for “Scientific Research Center” (Grant No. U1530401). S.H.Z.
thanks the preliminary research program (YY1708) of the College of Science of
BUCT. We acknowledge the computational support from the Beijing Computational
Science Research Center (CSRC).






§ HAMILTONIAN OF UNIFORM GRAPHENE IN CARTESIAN COORDINATES




As shown in Fig. <ref>(a), the honeycomb lattice of graphene
consists of two sublattices (denoted by A and B) and each unit cell
contains two carbon atoms (or π_z-orbitals), one on each sublattice. The
tight-binding Hamiltonian of graphene is    
    Ĥ_0=-∑_⟨ i,j⟩|i,A⟩⟨
    j,B|-h.c.,

where ⟨ i,j⟩ sums over all the nearest-neighbor carbon pairs,
|i,λ⟩ (λ=A,B) is the π_z orbital on the sublattice
λ of unit cell i, and we have taken the nearest-neighbor hopping
energy as the unit of energy. We take the A site as the origin of each unit
cell, so the relative displacements of the two carbon atoms inside the unit
cell are τ_A=0 and τ_B=(1/2,-√(3)/2), where we have taken the carbon-carbon bond length as the unit of
length. In terms of the unit cell location 𝐑_i, the location of
the site λ in unit cell i is 𝐑_i,λ=𝐑
_i+τ_λ.

By Fourier transforming the real-space basis {|i,λ⟩} into the
momentum space basis {|𝐤,λ⟩}, the graphene
Hamiltonian can be transformed into the momentum space:

    Ĥ_0=∑_𝐤f(𝐤)|𝐤,A⟩⟨𝐤,B|+h.c.,

where f(𝐤)=-∑_i=1,2,3e^i𝐤·𝐝_i and
{𝐝_i} depend on the choices on the Fourier transform. Namely,
if we make the Fourier transform using the location {𝐑_iλ} of the carbon sites:

    |𝐤,λ⟩_I≡1/√(N)∑
    _ie^i𝐤·𝐑_iλ|i,λ⟩,

then 𝐝_1,𝐝_2,𝐝_3 denote the relative
displacements of the three nearest-neighbor B sites [empty circles
in Fig. <ref>(a)] with respect to the central A site
[blue filled circle in Fig. <ref>(a)]:

    𝐝_1^I=τ_B,  𝐝
    _2^I=(  1/2,√(3)/2)
    ,  𝐝_3^I=(  -1,0)  .

However, if we make the Fourier transform using the location {𝐑
_i} of the unit cells:


    |𝐤,λ⟩_II≡1/√(N)∑
    _ie^i𝐤·𝐑_i|i,λ⟩,

then 𝐝_1,𝐝_2,𝐝_3 denote the relative
displacements of the three nearest-neighbor unit cells [black filled
circles in Fig. <ref>(a)] with respect to the central unit
cell [blue filled circle in Fig. <ref>(a)]:

    𝐝_1^II=0,   𝐝_2^II=(0,√(3)),   𝐝_3^II=(  -3/2
    ,√(3)/2)  .

These two choices are connected by |𝐤,A⟩_II
=|𝐤,A⟩_I, |𝐤,B⟩_II
=e^-i𝐤·τ_B|𝐤,B⟩_I, 𝐝_i^II=𝐝_i^I-τ_B, and f_II(𝐤)= e^-i𝐤·τ_Bf_I(𝐤). Both cases satisfy
f^∗(-𝐤)=f(𝐤); thus under time-reversal operation
θ̂, which leaves |iλ⟩ invariant and brings
|𝐤,λ⟩ into θ̂|𝐤,λ⟩=|-𝐤,λ⟩, the graphene Hamiltonian remains
invariant:

    θ̂Ĥ_0θ̂^-1=∑_𝐤f^∗
    (𝐤)|-𝐤,A⟩⟨-𝐤,B|+h.c.=Ĥ_0.

Diagonalizing the Hamiltonian in the momentum space gives one conduction band
and one valence band E_±(𝐤)=±|f(𝐤)| or explicitly

    E_±(𝐤)=±√(3+2cos(√(3)k_y)+4cos3k_x/2
    cos√(3)k_y/2),

which touch each other (i.e., f(𝐤)=0) at six Dirac points at the
edge of the first Brillouin zone [Fig. <ref>(b)]. Since
𝐊_1,𝐊_3,𝐊_5 (𝐊_2,𝐊
_4,𝐊_6) only differ by a reciprocal vector, only two Dirac
points are inequivalent, e.g., 𝐊_1 and 𝐊_2.
Hereafter we denote 𝐊_1 as 𝐊 and 𝐊_2 as
𝐊^':

    𝐊=2π/3(  1,1/√(3))  ,𝐊
    ^'=2π/3(  1,-1/√(3))  .



The continuum model near the Dirac point 𝐊_j is obtained by
considering 𝐤≈𝐊_j and expanding f(𝐤)
into Taylor series of 𝐪≡𝐤-𝐊_j:

    Ĥ_j(𝐪)=∑_𝐪f_j(𝐪)|𝐊
    _j,A⟩⟨𝐊_j,B|+h.c.,

where f_j(𝐪)≡ f(𝐊_j+𝐪). Up to the second
order of 𝐪, we have

    f_1,I(𝐪)       =f_4,I^∗(-𝐪
    )=3/2qe^-iπ/6+iθ_𝐪-3/8q^2e^iπ
    /3e^-i2θ_𝐪,
    
    f_2,I(𝐪)       =f_5,I^∗(-𝐪
    )=3/2qe^-iπ/6-iθ_𝐪-3/8
    q^2e^iπ/3e^i2θ_𝐪,
    
    f_3,I(𝐪)       =f_6,I^∗(-𝐪
    )=3/2iqe^iθ_𝐪+3/8q^2e^-i2θ
    _𝐪,

for choice I, where θ_𝐪 is the azimuthal angle of
𝐪. This choice gives a concise form for the O(q^2) terms, so it
is usually used to study the trigonal warping effect of graphene
<cit.>. For choice II, the O(q^2) term is very
complicated (although different choices of the Fourier transform give the same
physics <cit.>), so we only give the Taylor expansion up to
O(q):

    f_1,II(𝐪)       =f_3,II(𝐪
    )=f_5,II(𝐪)=3/2qe^-iπ/6+iθ_𝐪
    ,
    
    f_2,II(𝐪)       =f_4,II(𝐪
    )=f_6,II(𝐪)=-3/2qe^iπ/6-iθ_𝐪.

For either case, we have the identity f_i+3(𝐪)=f_i^∗(-𝐪), which follows from the time-reversal symmetry f^∗(-𝐤)=f(𝐤).




74
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1url<#>1urlprefixURL 


[Klein(1929)]KleinZP1929
authorO. Klein,
  journalZeitschrift für Physik volume53,
  pages157 (year1929).

[Allain and Fuchs(2011)]AllainEPJB2011
authorP. E. Allain and
  authorJ. Fuchs,
  journalThe European Physical Journal B
  volume83, pages301 (year2011).

[Katsnelson et al.(2006)Katsnelson,
  Novoselov, and Geim]KatsnelsonNatPhys2006
authorM. I. Katsnelson,
  authorK. S. Novoselov,
  and authorA. K. Geim,
  journalNat. Phys. volume2,
  pages620 (year2006).

[Cheianov and Fal'ko(2006)]CheianovPRB2006
authorV. V. Cheianov and
  authorV. I. Fal'ko,
  journalPhys. Rev. B volume74,
  pages041403 (year2006).

[Pereira et al.(2006)Pereira, Mlinar,
  Peeters, and Vasilopoulos]PereiraPRB2006
authorJ. M. Pereira,
  authorV. Mlinar,
  authorF. M. Peeters,
  and
  authorP. Vasilopoulos,
  journalPhys. Rev. B volume74,
  pages045424 (year2006).

[Bai and Zhang(2007)]BaiPRB2007
authorC. Bai and
  authorX. Zhang,
  journalPhys. Rev. B volume76,
  pages075430 (year2007).

[Beenakker et al.(2008)Beenakker,
  Akhmerov, Recher, and Tworzydło]BeenakkerPRB2008
authorC. W. J. Beenakker,
  authorA. R. Akhmerov,
  authorP. Recher, and
  authorJ. Tworzydło,
  journalPhys. Rev. B volume77,
  pages075409 (year2008).

[Setare and Jahani(2010)]SetareJPCM2010
authorM. R. Setare and
  authorD. Jahani,
  journalJ. Phys. Condens. Matter volume22,
  pages245503 (year2010).

[Roslyak et al.(2010)Roslyak, Iurov,
  Gumbs, and Huang]RoslyakJPCM2010
authorO. Roslyak,
  authorA. Iurov,
  authorG. Gumbs, and
  authorD. Huang, journalJ.
  Phys. Condens. Matter volume22, pages165301
  (year2010).

[Zeb et al.(2008)Zeb, Sabeeh, and
  Tahir]ZebPRB2008
authorM. A. Zeb,
  authorK. Sabeeh, and
  authorM. Tahir,
  journalPhys. Rev. B volume78,
  pages165420 (year2008).

[Sonin(2009)]SoninPRB2009
authorE. B. Sonin,
  journalPhys. Rev. B volume79,
  pages195438 (year2009).

[Schelter et al.(2010)Schelter, Bohr, and
  Trauzettel]SchelterPRB2010
authorJ. Schelter,
  authorD. Bohr, and
  authorB. Trauzettel,
  journalPhys. Rev. B volume81,
  pages195441 (year2010).

[Yang et al.(2011)Yang, Huang, Lai, and
  Grebogi]YangPRB2011c
authorR. Yang,
  authorL. Huang,
  authorY.-C. Lai, and
  authorC. Grebogi,
  journalPhys. Rev. B volume84,
  pages035426 (year2011).

[Rozhkov et al.(2011)Rozhkov, Giavaras,
  Bliokh, Freilikher, and Nori]RozhkovPhysRep2011
authorA. Rozhkov,
  authorG. Giavaras,
  authorY. P. Bliokh,
  authorV. Freilikher,
  and authorF. Nori,
  journalPhys. Rep. volume503,
  pages77  (year2011).

[Liu et al.(2012)Liu, Bundesmann, and
  Richter]LiuPRB2012b
authorM.-H. Liu,
  authorJ. Bundesmann,
  and authorK. Richter,
  journalPhys. Rev. B volume85,
  pages085406 (year2012).

[Giavaras and Nori(2012)]GiavarasPRB2012
authorG. Giavaras and
  authorF. Nori,
  journalPhys. Rev. B volume85,
  pages165446 (year2012).

[Rodriguez-Vargas
  et al.(2012)Rodriguez-Vargas, Madrigal-Melchor, and
  Oubram]RodriguezJAP2012
authorI. Rodriguez-Vargas,
  authorJ. Madrigal-Melchor,
  and authorO. Oubram,
  journalJ. Appl. Phys. volume112,
  pages073711 (year2012).

[Popovici et al.(2012)Popovici, Oliveira,
  de Paula, and Frederico]PopoviciPRB2012
authorC. Popovici,
  authorO. Oliveira,
  authorW. de Paula, and
  authorT. Frederico,
  journalPhys. Rev. B volume85,
  pages235424 (year2012).

[Heinisch et al.(2013)Heinisch, Bronold,
  and Fehske]HeinischPRB2013
authorR. L. Heinisch,
  authorF. X. Bronold,
  and authorH. Fehske,
  journalPhys. Rev. B volume87,
  pages155409 (year2013).

[Logemann et al.(2015)Logemann,
  Reijnders, Tudorovskiy, Katsnelson, and Yuan]LogemannPRB2015
authorR. Logemann,
  authorK. J. A. Reijnders,
  authorT. Tudorovskiy,
  authorM. I. Katsnelson,
  and authorS. Yuan,
  journalPhys. Rev. B volume91,
  pages045420 (year2015).

[Oh et al.(2016)Oh, Coh, Son, and
  Cohen]OhPRL2016
authorH. Oh,
  authorS. Coh,
  authorY.-W. Son, and
  authorM. L. Cohen,
  journalPhys. Rev. Lett. volume117,
  pages016804 (year2016).

[Erementchouk et al.(2016)Erementchouk,
  Mazumder, Khan, and Leuenberger]ErementchoukJPCM2016
authorM. Erementchouk,
  authorP. Mazumder,
  authorM. A. Khan, and
  authorM. N. Leuenberger,
  journalJ. Phys. Condens. Matter volume28,
  pages115501 (year2016).

[Downing and Portnoi(2017)]DowningJPCM2017
authorC. A. Downing and
  authorM. E. Portnoi,
  journalJ. Phys. Condens. Matter volume29,
  pages315301 (year2017).

[Huard et al.(2007)Huard, Sulpizio,
  Stander, Todd, Yang, and Goldhaber-Gordon]HuardPRL2007
authorB. Huard,
  authorJ. A. Sulpizio,
  authorN. Stander,
  authorK. Todd,
  authorB. Yang, and
  authorD. Goldhaber-Gordon,
  journalPhys. Rev. Lett. volume98,
  pages236803 (year2007).

[Gorbachev et al.(2008)Gorbachev,
  Mayorov, Savchenko, Horsell, and Guinea]GorbachevNL2008
authorR. V. Gorbachev,
  authorA. S. Mayorov,
  authorA. K. Savchenko,
  authorD. W. Horsell,
  and authorF. Guinea,
  journalNano Lett. volume8,
  pages1995 (year2008).

[Stander et al.(2009)Stander, Huard, and
  Goldhaber-Gordon]StanderPRL2009
authorN. Stander,
  authorB. Huard, and
  authorD. Goldhaber-Gordon,
  journalPhys. Rev. Lett. volume102,
  pages026807 (year2009).

[Young and Kim(2009)]YoungNatPhys2009
authorA. F. Young and
  authorP. Kim, journalNat.
  Phys. volume5, pages222
  (year2009).

[Rossi et al.(2010)Rossi, Bardarson,
  Brouwer, and Das Sarma]RossiPRB2010
authorE. Rossi,
  authorJ. H. Bardarson,
  authorP. W. Brouwer,
  and
  authorS. Das Sarma,
  journalPhys. Rev. B volume81,
  pages121408 (year2010).

[Sajjad et al.(2012)Sajjad, Sutar, Lee,
  and Ghosh]SajjadPRB2012
authorR. N. Sajjad,
  authorS. Sutar,
  authorJ. U. Lee, and
  authorA. W. Ghosh,
  journalPhys. Rev. B volume86,
  pages155412 (year2012).

[Sutar et al.(2012)Sutar, Comfort, Liu,
  Taniguchi, Watanabe, and Lee]SutarNL2012
authorS. Sutar,
  authorE. S. Comfort,
  authorJ. Liu,
  authorT. Taniguchi,
  authorK. Watanabe, and
  authorJ. U. Lee,
  journalNano Lett. volume12,
  pages4460 (year2012).

[Rahman et al.(2015)Rahman, Guikema,
  Hassan, and Marković]RahmanAPL2015
authorA. Rahman,
  authorJ. W. Guikema,
  authorN. M. Hassan,
  and
  authorN. Marković,
  journalAppl. Phys. Lett. volume106,
  pages013112 (year2015).

[Gutirrez et al.(2016)Gutirrez,
  Brown, Kim, Park, and Pasupathy]GutierrezNatPhys2016
authorC. Gutirrez,
  authorL. Brown,
  authorC. J. Kim,
  authorJ. Park, and
  authorA. N. Pasupathy,
  journalNat. Phys. volume12,
  pages1069 (year2016).

[Chen et al.(2016)Chen, Han, Elahi,
  Habib, Wang, Wen, Gao, Taniguchi, Watanabe, Hone et al.]ChenScience2016
authorS. Chen,
  authorZ. Han,
  authorM. M. Elahi,
  authorK. M. M. Habib,
  authorL. Wang,
  authorB. Wen,
  authorY. Gao,
  authorT. Taniguchi,
  authorK. Watanabe,
  authorJ. Hone, et al.,
  journalScience volume353,
  pages1522 (year2016).

[Laitinen et al.(2016)Laitinen, Paraoanu,
  Oksanen, Craciun, Russo, Sonin, and Hakonen]LaitinenPRB2016
authorA. Laitinen,
  authorG. S. Paraoanu,
  authorM. Oksanen,
  authorM. F. Craciun,
  authorS. Russo,
  authorE. Sonin, and
  authorP. Hakonen,
  journalPhys. Rev. B volume93,
  pages115413 (year2016).

[Bai et al.(2017)Bai, Qiao, Jiang, Liu,
  and He]BaiPRB2017
authorK.-K. Bai,
  authorJ.-B. Qiao,
  authorH. Jiang,
  authorH. Liu, and
  authorL. He, journalPhys.
  Rev. B volume95, pages201406
  (year2017).

[Sajjad and Ghosh(2011)]SajjadAPL2011
authorR. N. Sajjad and
  authorA. W. Ghosh,
  journalAppl. Phys. Lett. volume99,
  pages123101 (year2011).

[Jang et al.(2013)Jang, Kim, Son,
  Atwater, and Goddard]JangPNAS2013
authorM. S. Jang,
  authorH. Kim,
  authorY.-W. Son,
  authorH. A. Atwater,
  and authorW. A.
  Goddard, journalProc. Natl. Acad. Sci.
  volume110, pages8786 (year2013).

[Wilmart et al.(2014)Wilmart, Berrada,
  Torrin, Nguyen, Fève, Berroir, Dollfus, and
  Plaçais]Wilmart2DM2014
authorQ. Wilmart,
  authorS. Berrada,
  authorD. Torrin,
  authorV. H. Nguyen,
  authorG. Fève,
  authorJ.-M. Berroir,
  authorP. Dollfus, and
  authorB. Plaçais,
  journal2D Materials volume1,
  pages011006 (year2014).

[Chen and Chang(2015)]ChenPRB2015
authorC.-C. Chen and
  authorY.-C. Chang,
  journalPhys. Rev. B volume92,
  pages245406 (year2015).

[de Sousa et al.(2017)de Sousa, Chaves,
  PereiraJr., and Farias]SousaJAP2017
authorD. J. P. de Sousa,
  authorA. Chaves,
  authorJ. M. PereiraJr.,
  and authorG. A.
  Farias, journalJ. Appl. Phys.
  volume121, pages024302
  (year2017).

[Tan et al.(2017)Tan, Elahi, Tsao, Habib,
  Barker, and Ghosh]Tan2017
authorY. Tan,
  authorM. M. Elahi,
  authorH.-Y. Tsao,
  authorK. M. M. Habib,
  authorN. S. Barker,
  and authorA. W. Ghosh,
  journalScientific Reports volume7,
  pages9714 (year2017).

[Bir and Pikus(1974)]PikusStrainBook1974
authorG. Bir and
  authorG. Pikus,
  titleSymmetry and Strain-induced Effects in Semiconductors
  (publisherWiley, New York, year1974).

[Castro Neto et al.(2009)Castro Neto,
  Guinea, Peres, Novoselov, and Geim]CastroRMP2009
authorA. H. Castro Neto,
  authorF. Guinea,
  authorN. M. R. Peres,
  authorK. S. Novoselov,
  and authorA. K. Geim,
  journalRev. Mod. Phys. volume81,
  pages109 (year2009).

[Wu et al.(2007)Wu, Li, Song, Berger, and
  de Heer]WuPRL2007a
authorX. Wu,
  authorX. Li,
  authorZ. Song,
  authorC. Berger, and
  authorW. A. de Heer,
  journalPhys. Rev. Lett. volume98,
  pages136801 (year2007).

[Dombrowski et al.(2017)Dombrowski,
  Jolie, Petrovi ćć, Runte, Craes, Klinkhammer,
  Kralj, Lazi ćć, Sela, and
  Busse]DombrowskiPRL2017
authorD. Dombrowski,
  authorW. Jolie,
  authorM. Petrovi ć
se ć, authorS. Runte,
  authorF. Craes,
  authorJ. Klinkhammer,
  authorM. Kralj,
  authorP. Lazi ćć, authorE. Sela,
  and authorC. Busse,
  journalPhys. Rev. Lett. volume118,
  pages116401 (year2017).

[Xing et al.(2010)Xing, Wang, and
  Sun]XingPRB2010
authorY. Xing,
  authorJ. Wang, and
  authorQ.-f. Sun,
  journalPhys. Rev. B volume81,
  pages165425 (year2010).

[Reijnders and Katsnelson(2017)]ReijndersPRB2017
authorK. J. A. Reijnders
  and authorM. I.
  Katsnelson, journalPhys. Rev. B
  volume96, pages045305
  (year2017).

[Cheianov et al.(2007)Cheianov, Fal'ko,
  and Altshuler]CheianovScience2007
authorV. V. Cheianov,
  authorV. Fal'ko, and
  authorB. L. Altshuler,
  journalScience volume315,
  pages1252 (year2007).

[Lee et al.(2015)Lee, Park, and
  Lee]LeeNatPhys2015
authorG.-H. Lee,
  authorG.-H. Park, and
  authorH.-J. Lee,
  journalNat. Phys. volume11,
  pages925 (year2015).

[Garcia-Pomar et al.(2008)Garcia-Pomar,
  Cortijo, and Nieto-Vesperinas]GarciaPomarPRL2008
authorJ. L. Garcia-Pomar,
  authorA. Cortijo, and
  authorM. Nieto-Vesperinas,
  journalPhys. Rev. Lett. volume100,
  pages236801 (year2008).

[Jr et al.(2009)Jr, Peeters, Filho, and
  Farias]JrJPCM2009
authorJ. M. P. Jr,
  authorF. M. Peeters,
  authorR. N. C. Filho,
  and authorG. A.
  Farias, journalJ. Phys. Condens. Matter
  volume21, pages045301
  (year2009).

[Reich et al.(2002)Reich, Maultzsch,
  Thomsen, and Ordejón]ReichPRB2002
authorS. Reich,
  authorJ. Maultzsch,
  authorC. Thomsen, and
  authorP. Ordejón,
  journalPhys. Rev. B volume66,
  pages035412 (year2002).

[Datta(1995)]DattaBook1995
authorS. Datta,
  titleElectronic Transport in Mesoscopic Systems
  (publisherCambridge University Press, Cambridge, England,
  year1995).

[Ferry and Goodnick(1997)]FerryBook1997
authorD. Ferry and
  authorS. Goodnick,
  titleTransport in Nanostructures
  (publisherCambridge University Press, Cambridge,
  year1997).

[Lewenkopf and Mucciolo(2013)]LewenkopfJCE2013
authorC. H. Lewenkopf
  and authorE. R.
  Mucciolo, journalJournal of Computational
  Electronics volume12, pages203
  (year2013).

[Zhang et al.(2017)Zhang, Yang, and
  Chang]ZhangPRB2017
authorS.-H. Zhang,
  authorW. Yang, and
  authorK. Chang,
  journalPhys. Rev. B volume95,
  pages075421 (year2017).

[Beenakker(2008)]BeenakkerRMP2008
authorC. W. J. Beenakker,
  journalRev. Mod. Phys. volume80,
  pages1337 (year2008).

[Ando(1991)]AndoPRB1991
authorT. Ando, journalPhys.
  Rev. B volume44, pages8017
  (year1991).

[Culcer and Winkler(2008)]CulcerPRB2008
authorD. Culcer and
  authorR. Winkler,
  journalPhys. Rev. B volume78,
  pages235417 (year2008).

[Ajiki and Ando(1996)]AjikiJPSJ1996
authorH. Ajiki and
  authorT. Ando, journalJ.
  Phys. Soc. Jpn. volume65, pages505
  (year1996).

[Ando et al.(1998)Ando, Nakanishi, and
  Saito]AndoJPSJ1998
authorT. Ando,
  authorT. Nakanishi,
  and authorR. Saito,
  journalJ. Phys. Soc. Jpn. volume67,
  pages2857 (year1998).

[Bai et al.(2017)Bai, Zhou, Wei,
  Qiao, Liu, Liu, Jiang, and He]BaiArxiv2017
authorK.-K. Bai,
  authorJ.-J. Zhou,
  authorY.-C. Wei,
  authorJ.-B. Qiao,
  authorY.-W. Liu,
  authorH.-W. Liu,
  authorH. Jiang, and
  authorL. He,
  journalarXiv:1705.10952 [cond-mat.mes-hall]
  (year2017), 1705.10952.

[Roldán et al.(2008)Roldán,
  López-Sancho, and Guinea]RoldanPRB2008
authorR. Roldán,
  authorM. P. López-Sancho,
  and authorF. Guinea,
  journalPhys. Rev. B volume77,
  pages115410 (year2008).

[Polini et al.(2013)Polini, Guinea,
  Lewenstein, Manoharan, and Pellegrini]PoliniNatNano2013
authorM. Polini,
  authorF. Guinea,
  authorM. Lewenstein,
  authorH. C. Manoharan,
  and
  authorV. Pellegrini,
  journalNat Nano volume8,
  pages625 (year2013).

[Longhi(2010)]LonghiPRB2010
authorS. Longhi,
  journalPhys. Rev. B volume81,
  pages075102 (year2010).

[Casanova et al.(2010)Casanova,
  García-Ripoll, Gerritsma, Roos, and Solano]CasanovaPRA2010
authorJ. Casanova,
  authorJ. J. García-Ripoll,
  authorR. Gerritsma,
  authorC. F. Roos, and
  authorE. Solano,
  journalPhys. Rev. A volume82,
  pages020101 (year2010).

[Zhang et al.(2012)Zhang, Xue, Yan, Wang,
  and Zhu]ZhangPRA2012
authorD.-W. Zhang,
  authorZ.-Y. Xue,
  authorH. Yan,
  authorZ. D. Wang, and
  authorS.-L. Zhu,
  journalPhys. Rev. A volume85,
  pages013628 (year2012).

[Fang et al.(2016)Fang, Zhang, Louie, and
  Chan]FangPRB2016
authorA. Fang,
  authorZ. Q. Zhang,
  authorS. G. Louie, and
  authorC. T. Chan,
  journalPhys. Rev. B volume93,
  pages035422 (year2016).

[Ozawa et al.(2017)Ozawa, Amo, Bloch, and
  Carusotto]OzawaPRA2017
authorT. Ozawa,
  authorA. Amo,
  authorJ. Bloch, and
  authorI. Carusotto,
  journalPhys. Rev. A volume96,
  pages013813 (year2017).

[Gerritsma et al.(2011)Gerritsma, Lanyon,
  Kirchmair, Zähringer, Hempel, Casanova, García-Ripoll, Solano, Blatt,
  and Roos]GerritsmaPRL2011
authorR. Gerritsma,
  authorB. P. Lanyon,
  authorG. Kirchmair,
  authorF. Zähringer,
  authorC. Hempel,
  authorJ. Casanova,
  authorJ. J. García-Ripoll,
  authorE. Solano,
  authorR. Blatt, and
  authorC. F. Roos,
  journalPhys. Rev. Lett. volume106,
  pages060503 (year2011).

[Salger et al.(2011)Salger, Grossert,
  Kling, and Weitz]SalgerPRL2011
authorT. Salger,
  authorC. Grossert,
  authorS. Kling, and
  authorM. Weitz,
  journalPhys. Rev. Lett. volume107,
  pages240401 (year2011).

[Tudorovskiy et al.(2012)Tudorovskiy,
  Reijnders, and Katsnelson]TudorovskiyPS2012
authorT. Tudorovskiy,
  authorK. J. A. Reijnders,
  and authorM. I.
  Katsnelson, journalPhys. Scr.
  volumeT146, pages014010
  (year2012).

[Culcer and Winkler(2009)]CulcerPRB2009a
authorD. Culcer and
  authorR. Winkler,
  journalPhys. Rev. B volume79,
  pages165422 (year2009).

[Bena and Montambaux(2009)]BenaNJP2009
authorC. Bena and
  authorG. Montambaux,
  journalNew J. Phys. volume11,
  pages095003 (year2009).




