Quantum transport across van der Waals domain walls in bilayer graphene
Hasan M. Abdullah, B. Van Duppen, M. Zarenia, H. Bahlouli, F. M. Peeters
QQuantum transport across van der Waals domain walls in bilayer graphene
Hasan M. Abdullah,
1, 2, 3, ∗ B. Van Duppen, † M. Zarenia, H. Bahlouli,
1, 2 and F. M. Peeters Department of Physics, King Fahd University of Petroleum and Minerals, 31261 Dhahran, Saudi Arabia Saudi Center for Theoretical Physics, 31261 Dhahran, Saudi Arabia Department of Physics, University of Antwerp, Groenenborgerlaan 171, B-2020 Antwerp, Belgium (Dated: July 2, 2018)Bilayer graphene can exhibit deformations such that the two graphene sheets are locally detachedfrom each other resulting in a structure consisting of domains with different inter-layer coupling.Here we investigate how the presence of these domains affect the transport properties of bilayergraphene. We derive analytical expressions for the transmission probability, and the correspond-ing conductance, across walls separating different inter-layer coupling domain. We find that thetransmission can exhibit a valley-dependent layer asymmetry and that the domain walls have aconsiderable effect on the chiral tunnelling properties of the charge carriers. We show that trans-port measurements allow one to obtain the strength with which the two layers are coupled. Weperformed numerical calculations for systems with two domain walls and find that the availability ofmultiple transport channels in bilayer graphene modifies significantly the conductance dependenceon inter-layer potential asymmetry.
PACS numbers: 73.20.Mf, 71.45.GM, 71.10.-w
I. INTRODUCTION
A decade ago, researchers started investigatinggraphene and its associated multilayers for use as a ba-sis for next generation of fast and smart electronic logicgates. The absence of a band gap leads to different pro-posals for gap generation . For example, by changingthe size of the graphene flakes into nanoribbons or quan-tum dots, one can control the energy gap through sizequantization . Important experimental advances wereachieved in recent years which enabled the fabrication ofgraphene based electronic devices at the nano scale .The increasing control over the structure of grapheneflakes allowed for new devices that could constitute thebuilding blocks for a fully integrated carbon based elec-tronics. An example of this is deformed bilayer graphene,where the two layers are not aligned due to a mismatchin orientation or stacking order resulting in e.g. twistedbilayer graphene. Its electronic structure is strongly dif-ferent from normal bilayer graphene and exhibits verypeculiar properties such as the appearance of additionalDirac cones .Recent experiments have shown that epitaxialgraphene can form step-like bilayer/single layer (SL/BL)interfaces or that it is possible to create bilayer grapheneflakes that are connected to single layer grapheneregions . The appearance of these structures fueledtheoretical and experimental investigations on the be-havior of massless and massive particles in such junc-tions. For example, few works have investigated differentdomain walls that separate, for instance, different typeof stacking or even different number of layers .These theoretical investigations showed that the trans-mission probabilities through SL/BL interfaces exhibitsa valley-dependent asymmetry which could be used forvalley-based electronic applications . Other theoreti-cal and experimental works focused on the emergence of Landau levels, edge state properties and peculiar trans-port properties in such systems . Bilayer grapheneflake sandwiched between two single zigzag or armchairnanoribbons was also investigated and it was foundthat the conductance exhibits oscillations for energieslarger than the inter-layer coupling.Most of these recent theoretical works considered do-main walls separating patches of bilayer graphene withdifferent stacking type or where only a single layer wasconnected to a bilayer graphene sheet. Very recently,however, a number of new bilayer graphene platformshave been synthesized. These consist of regions where thecoupling between the two graphene layers is changed. Forexample in the case of folded graphene a part of thefold forms a coupled bilayer structure, while another partof it is uncoupled . One has also observed systemswith domain walls separating regions of different Bernalstacking . In general, these systems can be modelledas being composed of two single layers of graphene (2SL)which are locally bound by van der Waals interaction intoan AA- or AB-stacked bilayer structure.Here, we present a systematic study of electrical trans-port across domain walls separating regions of differentinter-layer coupling. We discuss the dependance on thecoupling between the graphene layers, on the distancebetween subsequent domain walls and on local electro-static gating. For completeness, we also present all pos-sible combinations of locally detached bilayer systems.Analytical expressions for the transport across a singledomain wall are also obtained. These results can serveas a guide for future experiments.From a theoretical point of view, one can wonder howcharge carriers will respond to transitions between sys-tems that have completely different transport properties.For example, single layer graphene and AA-stacked bi-layer graphene are known to feature Klein tunnellingat normal incidence while AB-stacked bilayer graphene a r X i v : . [ c ond - m a t . m e s - h a ll ] O c t shows anti-Klein tunnelling . It is, therefore, inter-esting to investigate under which conditions these pe-culiar chirally-assisted tunnelling properties pertain incombined systems, as well as to investigate how the pres-ence of multiple transport channels changes the transportproperties.From our study we obtain useful analytical expressionsfor the transmission probability across a single domainwall. These results also show that the effect of localgating is to break the symmetry between the two lay-ers and to introduce a valley-dependent angular asym-metry, which could be used for a layer-dependent valley-filtering device. We show that the inter-layer couplingstrength and stacking has a characteristic effect on theconductance across a domain wall which can be used tomeasure structural deformations in bilayer graphene. Wefind that the presence of multiple conductance channelsin bilayer graphene can modify the dependance of theconductance on an applied inter-layer potential differ-ence from constructive to destructive. Finally, we showthat transitions in-between AA-stacked and AB-stackedbilayer graphene systems largely conserve the parity ofthe transport channel.The paper is organized as follows. In Sec. II, we dis-cuss the formalism, explain the geometry of the inves-tigated domain walls, and define the possible scatteringprocesses between the different transport modes. In Sec.III, we give analytical expressions for the transmissionprobabilities through one domain wall and analyze howthe symmetry between the graphene layers can be brokenby electrostatic potentials. An overview of the numeri-cal results for more complex set-ups consisting of multipledomain walls and gates is presented in Sec. IV. Finally, inSec. V we briefly summarize the main points of this paperand comment on possible experimental signatures of thepresence of coupling domain walls in bilayer graphene. II. MODEL
Single layer graphene consists of two inequivalent sub-lattices, denoted as α and β , with interatomic distance a = 0 .
142 nm and that are coupled in the tight binding(TB) formalism by γ = 3 eV . It has a gapless en-ergy spectrum with band crossings at the so-called Diracpoints K and K (cid:48) that are located at the corners of theBrillouin zone. The energy dispersion around one of thesepoints is depicted in Fig. 1(a).Bilayer graphene consists of two single layers ofgraphene which can be stacked in two stable configu-rations: AB-stacked bilayer graphene (AB-BL) or AA-stacked bilayer graphene (AA-BL). In AB-BL, atom α is placed directly above atom β with inter-layer coupling γ ≈ . as shown in Fig. 1(b). It has a parabolicdispersion relation with four bands. Two of them touchat zero energy, whereas the other two bands are splitaway by an energy γ . The skew hopping parameters γ and γ between the other two sublattices are negligible since they have insignificant effect on the transmissionprobabilities and band structure at high energies .In AA-BL two single layers of graphene are placedexactly on top of each other such that the structurebecomes mirror-symmetric. Atoms α and β in thetop layer are located directly above atoms α and β in the bottom layer, with direct inter-layer coupling γ ≈ . , see Fig. 1(c). AA-BL has a linear en-ergy spectrum with two Dirac cones shifted in energy byan amount of ± γ as depicted in Fig. 1(c) by the fullcurves. A. Geometries
We consider four different junctions that can be madeof from the building blocks depicted in Fig. 1: monolayer,AA- stacked and AB-stacked bilayer graphene. Withoutloss of generality, we assume that the charge carriers arealways propagating from the left to the right hand side.Then we consider three different configurations: ( I ) astructure where the leads on the left ( x <
0) and on theright hand side ( x > d ) consist of two decoupled singlelayers while in between they are connected into an AB-BL (AA-BL) configuration. This is depicted in Fig. 2(a).We will refer to such a structure as 2SL-AB-2SL (2SL-AA-2SL). ( II ) A structure where the middle region ismade up of two decoupled monolayers and whose leadsare AB (AA) stacked bilayer graphene. This is depictedin Figs. 2(b, d). Such a configuration henceforth will berefereed to as AB-2SL-AB (AA-2SL-AA). ( III ) A struc-ture where a domain wall separates an AB (AA) stackedstructure from two decoupled single layers. We will as-sign the abbreviation 2SL-AB (2SL-AA) to this structureif the charge carriers are incident on one of the two sep-arated layers or AB-2SL (AA-2SL) if the coupled bilayerstructure is connected to the source. This is depicted inFig. 2(c). (IV) left and right leads are bilayer graphenewith AA- and AB stacking, respectively, separated by adomain where the two layers are completely decoupled(AA-2SL-AB), see Fig. 2(e). To describe transport inthe above mentioned structures, we allow for scatteringbetween the layers as well as between the different prop-agating modes in an AB-BL or between the two Diraccones in AA-BL. In the next section, we describe thetransport modes in 2SL and BL and how charge carrierscan be scattered in-between them.
B. Scattering definitions
In this section we define the model Hamiltonian thatdescribes the different structures. For this purpose we usea suitable basis defined by Ψ = (Ψ α , Ψ β , Ψ α , Ψ β ) T , whose elements refer to the sublattices in each layer. The (a) (b) ( c) AB-BL
AA-BLSL
SLG E k γ AB - BL kE AA - BL2 γ E k
FIG. 1: (Colour online) Lattice structure with their cor-responding energy spectrum of (a) Monolayer graphene,(b) AB-stacked bilayer graphene, (c) AA-stacked bilayergraphene. The dashed curves correspond to the spectrumin the presence of a finite bias. general form of the Hamiltonian near the K-point reads H = V v F π † τ γ v F π V ζγ τ γ τ γ ζγ V v F π † τ γ v F π V . (1)The coupling between the two graphene layers is con-trolled by the parameters τ and ζ through which we can“ switch on ” or “ switch off ” the inter-layer hopping be-tween specific sublattices. This allows to model differentstackings by assigning different values to these parame-ters. For τ = ζ = 0, the two layers are decoupled andthe Hamiltonian reduces to two independent SL sheets.To achieve AA-stacking we select τ = 1 and ζ = 0 whilefor AB-stacking we need τ = 0 and ζ = 1. In Eq. (1) v F ≈ m/s is the Fermi velocity of charge carriersin each graphene layer, π = p x + ip y denotes the mo-mentum, V and V are the potentials on layers 1 and2. In the present study, we only apply these potentialsin the intermediate region. We assume that the domainwall is oriented in the y -direction and of infinite length.Therefore, the system is translational invariant and themomentum p y is conserved. This enables us to write thewave function as Ψ ( x, y ) = e ik y y Φ ( x ).
1. Decoupled graphene layers
The eigenfunctions of the 2SL Hamiltonian are thoseof the isolated graphene sheet, Φ = (cid:18) φ φ (cid:19) , φ j = (cid:18) µ − j − µ + j (cid:19) (cid:18) e ik j x e − ik j x (cid:19) , (2)where j = 1 , k j = (cid:113) ( (cid:15) + s j δ ) − k y with s j =sgn( j − . µ ± j = ( k j ± ik y ) / ( (cid:15) + s j δ ), (cid:15) = E − v , δ = ( V − V ) / , v = ( V + V ) /
2. Introducing tb tb
TopBottom2SL 2SLBL(a) tb AB-BL (d)2SL d x x d x d d x x d FIG. 2: (Colour online) Different geometries for bilayer andtwo decoupled graphene layer interfaces with schematic rep-resentation of the transmission probabilities. (a) AA or ABstacking bilayer graphene sandwiched between two SL gara-phene layers (2SL-AA(AB)-2SL), (b) AB-BL leads with 2SLas intermediate region (AB-2SL-AB), (c) two single gaphenelayers connected to AB-BL(2SL-AB), and (d) similar to (b)but now with AA-BL as the leads with two upper (red)-lower(blue) shifted Dirac cones (AA-2SL-AA). (e) left and rightleads are bilayer graphene with different stacking connectedto the two decoupled graphene sheets (AA-2SL-AB). The pos-sible transmission processes between the different conductionchannel are indicated above the respective junctions. the length scale l = (cid:126) v F /γ , which represents the inter-layer coupling length, allows us to define the followingdimensionless quantities: (cid:15) → (cid:15)γ , v → v γ , δ → δγ , k y → lk y , and (cid:126)r → (cid:126)rl . (3)Notice that for the two stacking configurations, γ wasfound to be different. For the AB-BL the value is γ ≈ . γ ≈ . .In order to discuss the different scattering modes, weintroduce the notation A outgoingincoming , where A can stand fortransmission ( T ) or reflection ( R ) probabilities and theindexes denote the mode by which the particles are in-coming or outgoing. Fig. 2 depicts all possible transitionsthat are considered in the present work. Fig. 2(a) showsall possible transmission processes in a 2SL-BL-2SL sys-tem where t denotes the top layer on either side and b thebottom layer. For example, T bt denotes a particle comingthrough the top layer and exiting on the bottom layer.
2. AB-stacking
For AB-BL there are two branches corresponding topropagating modes. These branches correspond to thewave vector k ± given by k ± = (cid:104) − k y + (cid:15) + δ ± (cid:112) (cid:15) (1 + 4 δ ) − δ (cid:105) / , (4)The modes presented in Eq. (4) labeled by “ k + ” corre-spond to eigenstates that are odd under layer inversion,while the “ k − ”modes are even. These modes are shown,respectively, in blue and red in Fig. 2(b). This meansthat there are two available channels for transmission ata given energy, and an additional two for the reflectionprobabilities. Note that for energies 0 < E < γ , there isonly one propagating mode and one transmission and re-flection channel. Similarly, the wave function of AB-BLcan be written as Ψ( x, y ) = GM ( x ) Ce ik y y , (5)where M ( x ) corresponds to a 4 × C depend on the propagating region, and G is given by G = ξ + − − ξ ++ ξ −− − ξ − + ρ + ρ + ρ − ρ − ζ ++ − ζ + − ζ − + − ζ −− , (6)where ξ ±± = ( k ± ± ik y ) /E − δ, ρ ± = ( (cid:15) − δ ) (cid:2) − (( k ± ) + k y ) / ( (cid:15) − δ ) (cid:3) and ζ ±± = ( (cid:15) − δ ) ρ ± ξ ±± / ( (cid:15) + δ ).The use of the matrix notation will prove to be veryuseful to construct the transfer matrix as outlined below.
3. AA-stacking
In the case of an AA-BL, the corresponding wave func-tion can be written similar to Eq. (5) but now with thematrix G given by G = ξ − + ξ ++ ξ −− ξ + − ζ − + ζ ++ ζ −− ζ + − ρ + ρ + ρ − ρ − , (7) where ρ ± = (cid:15) (cid:2) − ( k y + ( k ± ) ) + ( (cid:15) − δ ) + 1) (cid:3) , ξ ±± =( ρ ± + δ + (cid:15) )( ik y ± k ± ) / ( δ − (cid:15) + 1) and ζ ± s = ( ξ ±± − ρ ± ( ik y ± k ± ) / ( (cid:15) + δ ). To investigate when scatteringbetween the Dirac cones of AA-BL is allowed or forbid-den, one can apply a unitary transformation that formssymmetric and anti-symmetric combinations of the topand bottom layer. This yields a Hamiltonian in the basis Ψ = 2 − / (Ψ α +Ψ α , Ψ β +Ψ β , Ψ α − Ψ α , Ψ β − Ψ β ) T of the form: H AA = γ + v v F π † − δ v F π γ + v − δ − δ − γ + v v F π † − δ v F π − γ + v . (8)For δ = 0, this Hamiltonian is block-diagonal and rep-resents two Dirac cones as shown in Fig. 1(c). The twocones correspond to modes with wave vector k ± given by k ± = (cid:20) − k y + (cid:16) (cid:15) ± (cid:112) (1 + δ ) (cid:17) (cid:21) / . (9)In Fig. 1(c) the blue bands corresond to the odd k + modes and red bands, denoting the even modes, are givenby the k − wavevector. In these equations, v denotes theenergy shift of the whole spectrum. This shift can bechosen zero by assigning the same magnitude but differ-ent signs to the electrostatic potentials on both layers V = − V . Eq. (8) shows that for zero electric field( δ = 0) both cones are decoupled and the scattering be-tween them is strictly forbidden. This was used before inRef. [53] to propose AA-BL as a potential candidate for “cone-tronics” based devices. However, this protectedcone transport is broken for finite bias ( δ (cid:54) = 0) and hencescattering between the cones is allowed. Furthermore,one might wonder if the charge carriers stay within theircone transport through a domain consisting of two de-coupled layers.
4. Scattering probability
In order to calculate the scattering probability in thereflection and transmission channel, we use the transfermatrix method together with boundary conditions thatrequire the eigenfunctions in each domain to be contin-ious for each sublattice . To conserve probability cur-rent we normalize transmission probabilities T and reflec-tion probabilities R such that (cid:88) i,j (cid:16) T ji + R ji (cid:17) = 1 , (10)where, the index i refers to the incoming mode whilethe index j denotes the outgoing mode. For a coupledbilayer the different modes are labelled by “ − ” for themodes that are even under in-plane inversion and by “+”for odd modes. For a decoupled 2SL system, we employthe notation t for the top layer and b for the bottom - AB R, RR, II, RR, R I, I ( a ) v v + γ v - γ k E R, R I, IR, I AB - ( b ) v γ k E FIG. 3: (Colour online) Schematic diagrams, for one domainwall separating 2SL and AB-BL, showing the regions wherethe modes ( k + , k − ) in AB-BL are either real (propagating)or imaginary (evanescent). (a) shows the bands of pristine2SL and gated AB-BL and vice versa in (b). In the yellowregion both modes are real (R, R), while one of them is realand the other is imaginary as in the green (I, R) and pink (R,I) regions. In the gray region both modes are imaginary (I,I). Blue, red and dashed black bands correspond to k + , k − and 2SL modes, respectively. layer. For example, for the system 2SL-AB-2SL and foran incident particle in the top layer of 2SL gives T tt + T bt + R tt + R bt = 1. In Fig. 2 all possible transition probabilitiesare shown schematically.
5. Conductance
To obtain measureable quantities, we finally calculatethe zero temperature conductance that can be obtainedfrom the Landauer-B¨uttiker formula where we have tosum over all the transmission channels, G ji ( E ) = G L y π (cid:90) + ∞−∞ dk y T ji ( E, k y ) , (11)with L y the length of the sample in the y -direction and G = 4 e /h . The factor 4 comes from the valley and spindegeneracy in graphene. The total conductance of anyconfiguration is the sum of all available channels G T = (cid:80) i,j G ji . III. TRANSMISSION ACROSS A SINGLEDOMAIN WALL
Here we will present analytical expressions for thetransmission probabilities of transport across a single do-main wall. These analytical expressions will shed light onthe requirements for transport across a domain wall andhow local electrostatic gating can affect these transportproperties. By doing so, we encounter that curiously,electrostatic gates can break the symmetry between thelayers in the transmission probability if there are evanes-cent modes in the system. The breaking of the layersymmetry results in an asymmetric angular distributionof the transmission probability as will be shown further. - ° - ° - ° ° ° R bt T b + R bb T b - ( a ) - AA v = - ° - ° - ° ° ° R bt R bb T b + T b - ( b ) - AA v = γ - ° - ° - ° ° ° R -- T - b T + b R ++ ( c ) AA - v = - ° - ° - ° ° ° R ++ R -- T - b T + b ( d ) AA - v = γ FIG. 4: (Colour online) The angle-dependent transmissionand reflection probabilities through (a, b) 2SL-AA and (c, d)AA-2SL systems. The systems in (b, d) are the same as in (a,c), respectively, but where now the right side of the junction issubjected to an electrostatic potential of strength v = 1 . γ . In the system 2SL-AA R b ( t ) b = R t ( b ) t and T ± b = T ± t while R − + = R + − = 0 and T b ± = T t ± in AA-2SL system. In all panels E = 1 . γ . We consider a situation where two propagating modesexist in the AB-BL or AA-BL. This requires some cau-tion in defining the incident angle in the calculation ofthe transmission probabilities. Failing to do so may re-sult in erroneous results such as transmission exceedingunity or unexpected symmetry features . Consideringone domain wall, the simplest configuration, separating2SL and either AA or AB-BL allows to obtain analyticexpressions for the transmission probabilities. The in-cident angle for each propagating mode depends on thetype of layer stacking in the incident region. Hence, forcharge carriers incident from 2SL we define k j = E cos φ, k y = E sin φ. (12)On the other hand, when charge carriers are incidentfrom AB-BL we need to define incident angle for eachmode separately such that k ± = (cid:112) E ± E cos φ, k y = (cid:112) E ± E sin φ. (13)Finally, if charge carriers incident from AA-BL the asso-ciated angle is defined as k ± = ( E ±
1) cos φ, k y = ( E ±
1) sin φ. (14)A straightforward calculation results in the transmis-sion probability for charge carriers incident from 2SL andimpinging on AA-BL T ± j = 2( (cid:15) + v )( ± (cid:15) )Re( k ± ) k j (cid:104) ( ± (cid:15) + k ± sec φ ) + ( ∓ v ) tan φ (cid:105) , (15)while for the reverse configuration (AA-2SL) it is givenby T j ± = 2 (cid:15) Re( k j )cos φ (cid:104) ( (cid:15) + k j sec φ ) + ( ∓ v ) tan φ (cid:105) . (16)Similar as performed for the AA-BL Hamiltonian, alsothe AB-BL Hamiltonian can be expressed in terms ofsymmetric and anti-symmetric combinations of the twolayers. This manipulation allows to determine a closed-form expression for the transmission probability of the2SL-AB structure. The derivation is outlined in Ap-pendix A and results in T ± j = 4Re( k ± ) η (cid:104) η + (Im( k ∓ ) + κ j v sin φ ) (cid:105) C + (cid:80) m =1 C m cos( mφ ) , (17)with η = (cid:15) cos φ and κ j = +1( −
1) for j = b ( t ). For thereverse configuration (AB-2SL) the transmission proba-bilities are T j ± = 4Re( k j ) k ± λ [ µ ± + κ j v sin φ Im( k ∓ )] | Q ± | , (18)where λ, C m , µ ± , and Q ± are functions defined in Ap-pendix A.For a domain wall separating 2SL and AA-BL, thetransmission probabilities are always symmetric with re-spect to normal incidence as indicated in Eqs. (15,16).In other words, for the 2SL-AA T ± b ( φ ) = T ± t ( φ ) and sim-ilarly T b ± ( φ ) = T t ± ( φ ) for AA-2SL configuration, and thissymmetry still holds when the right side of the junction isgated ( v (cid:54) = 0). We will refer to this symmetry as “ layersymmetry” since it is a consequence of the equivalence of2SL layers and the symmetric coupling of the AA-BL.Notice that Klein tunnelling for normal incidence inSL and AA-BL is also conserved in the combined struc-ture. For example, in 2SL-AA and for normal incidence( φ = 0), the modes become k j = (cid:15) + v , k ± = ± (cid:15) and hence Eq. (15) reads T ± j = 1 / . Then, for chargecarriers propagating in the bottom (top) layer it may be transmitted into k + or k − states and thus the total prob-ability is T + b ( t ) + T − b ( t ) = 1 / / R b ( t ) b = R t ( b ) t = 0.In an analogous manner it can be shown that for normalincidence Eq. (16) gives T j ± = 1 / . Turning now to the 2SL-AB/AB-2SL case, one can in-fer from Eqs. (17,18) that for v = 0 the layer symmetryholds since the only term carrying asymmetric features isproportional to v . However, for v (cid:54) = 0 it is striking thatdespite the fact that a homogeneous electrostatic poten-tial does not break any in-plane symmetry in the system,layer symmetry is broken. This leads to an angular asym-metry in the transmission channel, i.e. T ± b ( φ ) = T ± t ( − φ )for 2SL-AB and T b ± ( φ ) = T t ± ( − φ ) for AB-2SL. Upon fur-ther analysis of Eqs. (17,18), one notices that this asym-metric feature is present in regions in the ( E, k y ) planewhere one of the two modes is propagating while theother is evanescent. In Figs. 3(a,b) we show a diagramfor these different regions associated with 2SL-AB andAB-2SL, respectively. The layer symmetry is broken inthe green and pink regions while in the yellow regionslayer symmetry holds.The mechanism for breaking the layer symmetry inconfigurations consisting of AB-BL is attributed only tothe evanescent modes. For example, in 2SL-AB (see Fig.3) the transmission probability for charge carriers to betransmitted into k + from either bottom or top layers of2SL is T + j = 4Re( k + ) η (cid:104) η + (Im( k − ) + κ j v sin φ ) (cid:105) C + (cid:80) m =1 C m cos( mφ ) , (19)where κ b ( t ) = 1( − T + b ( φ ) = T + t ( − φ ) , only when v (cid:54) = 0 and Im( k − ) (cid:54) = 0 which is satisfied in the pink andgray regions in Fig. 3(a). However in the gray regionthere are no k + propagating states and consequently thetransmission probabilities T + j are zero. The same anal-ysis applies also to T − j where the asymmetric feature ispreserved only when Im( k + ) (cid:54) = 0 as shown by the greenregion in Fig. 3(a). For AB-2SL configuration, the layerasymmetry is only reflected in the T j + , see Eq. (18),since Im( k − ) (cid:54) = 0 corresponds to the pink region in Fig.3(b). While for T j − , the k − propagating states are onlyavailable for E > γ (yellow region in Fig. 3(b)) whichcoincides with Im( k + ) = 0. Thus, the layer symmetryis always conserved in T j − as it can be seen in Eq. (18).Now it is clear why layer symmetry is not broken in theAA-BL configuration; because there are always two prop-agating modes associated with any energy value.The breaking of angular symmetry in this situation isqualitatively similar to that obtained in AB-BL sub-ject to an inter-layer bias. One can connect this layerasymmetry in the vicinity of the two valleys K and K (cid:48) through time-reversal symmetry. The Hamiltonian H K (cid:48) can be related to the Hamiltonian H K through the trans-formation H K (cid:48) ( k ) = Θ H K ( − k )Θ − , (20)where Θ is the time-reversal symmetry operator. Thisimplies, for example in the T + b ( t ) channel, that chargecarriers moving from right to left and scattered from thebottom layer to k + in K valley are equivalent to thosescattered from top layer to k + but moving in the oppositedirection in the vicinity of K (cid:48) . If layer symmetry holdsin the vicinity of one of the valleys, then the transmissionprobabilities of charge carriers moving in the opposite di-rections must be the same. It is worth pointing out herethat the layer asymmetry in the K valley is reversed inthe K (cid:48) valley and hence the overall symmetry of the sys-tem is restored. Therefore, the macroscopic time reversalsymmetry is preserved. IV. NUMERICAL RESULTS
We first present the results for transmission, and re-flection probabilities and for the conductance in the caseof domain walls separating 2SL and AA-BL structures.The different regions as defined in Fig. 3 are superim-posed as dashed black and white curves. Moreover, incalculating the transport properties we considered differ-ent magnitudes for the electrostatic potential v and bias δ applied to the drain structure. A. AA-Stacking
1. 2SL-AA/AA-2SL
We consider charge carriers tunnelling through 2SL-AA and AA-2SL systems. In Fig. 4(a) we show thetransmission and reflection probabilities for charge carri-ers impinging on pristine AA-BL as a function of incidentangle φ . As a result of the layer symmetry, charge car-riers incident from bottom/top layer of 2SL and trans-mitted into the lower Dirac cone ( k + ) in the AA-BL willhave the same transmission probability T + b = T + t . Simi-larly, for those charge carriers transmitted into the uppercone, they will also have the same probability T − b = T − t regardless which layer they are incident from.This symmetry stems from the fact that the wavefunc-tion in the 2SL are a superposition of two spinors cor-responding to the two sublattices while in AA-BL it isa superposition of four. For this reason, charge carriersincident from top or bottom layer of 2SL have the samedynamics and hence share their transmission probabil-ity. A partial reflection into the same layer, R bb = R tt is shown in Fig. 4(a), which corresponds to evanescentmodes associated with the upper Dirac cone ( k − ). Asin transmission, charge carriers can be back scatteredbetween the layers. However, the absence of the electro-static potential results in a small scattering current as FIG. 5: (Colour online) Density plot of the transmission andreflection probabilities through 2SL-AA-2SL as a function ofFermi energy and transverse wave vector k y with v = δ = 0and width of the AA-BL d = 25 nm. depicted in Fig. 4(a). In addition, scattering back fromtop to bottom layer or vice versa occurs also with thesame reflection probabilities R tb = R bt .Because of chiral decoupling of oppositely propagatingwaves in AA-BL and in SL, back-scattering is forbid-den for normal incidence ( φ = 0) and thus the reflectionprobabilities for each channel are zero, i.e. R b ( t ) b (0) = R t ( b ) t (0) = 0. This is associated with perfect tunnelling T + b (0) + T − b (0) = T + t (0) + T − t (0) = 1. The effect holdsfor all forthcoming structures composed of AA-BL and2SL.Fig. 4(b) shows the numerical results of the same sys-tem, 2SL-AA, but now in the AA region, the potentialis increased to v = 1 . γ . This shifts the two Diraccones in energy to v ± γ . As a result of the presence ofthe electrostatic potential, a strong scattered reflection R tb /R bt takes place when there are no propagating modesin the AA section.In Figs. 4(c,d), we show the reversed configuration,i.e. an AA-2SL system. The transmission and reflectionprobabilities for zero ( v = 0) and with nonzero ( v =1 . γ ) electrostatic potentials applied to 2SL are reportedin panels (c) and (d) respectively. Similar to the 2SL-AAsystem, we can note that layer symmetry still holds suchthat T b + = T t + and T b − = T t − . Furthermore, we find strongnon-scattered reflection in the R ++ and R −− channels thatis associated with evanescent modes on both sides of AA-BL and 2SL whereas the scattered reflection channels R + − and R − + are always zero due to the protected conetransport discussed earlier. FIG. 6: (Colour online) The same as in Fig. 5, but now with v = 1 . γ . Red and white dashed curves correspond to thelower and upper Dirac cones in AA-BL, respectively, whilethe black dashed curves are the bands of 2SL.
2. 2SL-AA-2SL
In this Section, we show the results of transport acrosstwo domain walls forming a system with three regions;where AA-BL is sandwiched between two regions of 2SL,see Fig. 2(a). Such a system can exhibit a strong layerselectivity when current flows through the intermediateregion , i.e. AA-BL. This behaviour has already beeninvestigated in Ref. [59]. Here, however, we go in muchmore detail to show how the different transmission andreflection channels are affected by the electrostatic po-tential or finite bias applied to the intermediate region.In Figs. 5 and 6 we show the scattered andnon-scattered channels for transmission and reflectionfor pristine AA-BL and with electrostatic potential ofstrength v = 1 . γ , respectively. Layer symmetry ispreserved in both reflection and transmission channelsas clarified in Figs. 5. and 6 also show strong scat-tered transmission, especially for normal incidence whichcan be altered depending on the width of the AA-BL.When an electrostatic potential is applied to the middledomain, resonances appear in the transmission probabil-ities for v + γ > E > v − γ as shown in Fig. 6. Thisis a consequence of the finite size of the AA-BL and thepresence of charge carriers with different chirality in thementioned range of energies . Introducing a finite bias δ = 0 . γ on AA-BL breaks the layer symmetry of thesystem. As a result, T bb (cid:54) = T tt and R bb (cid:54) = R tt . However, itis still preserved in the scattered channels T bt = T tb and R bt = R tb ( see Fig. 7).It is worth mentioning here that the finite bias does FIG. 7: (Colour online) The same as in Fig. 5, but now with v = 1 . γ and δ = 0 . γ . not break the angular symmetry with respect to normalincidence in the transmission and reflection probabilitiesas it does for normal AB-BL . This is a manifestationof the symmetric inter-layer coupling in AA-BL.
3. AA-2SL-AA
In this system we interchange the AA-BL and 2SL asshown in Fig. 2(d). In this case, scattering is definedbetween the two cones in the AA-BL regions. In Figs.8 and 9 we show the transmission and reflection proba-bilities between the two Dirac cones through the pristine2SL and in the presence of an electrostatic potential, re-spectively. The first and the last rows of Figs. (8) and(9) show the non-scattered transmission and reflectionprobabilities corresponding to the lower and upper Diraccones, respectively. We notice that Klein tunnelling ispreserved at normal incidence. This shows that Kleintunnelling in AA-stacked bilayer graphene is a robustfeature that is insensitive to local changes in the inter-layer coupling. On the other hand we see that scatteringbetween two different Dirac cones remains strictly for-
FIG. 8: (Colour online) Density plot of the transmission andreflection probabilities through AA-2SL-AA as a function ofFermi energy and transverse wave vector k y with v = δ = 0and width of the 2SL d = 25 nm. bidden even with a local decoupling of the two layers.Therefore, these devices could be used for conetronics.As a result, in the second row of Figs. 8 and 9 thescattered transmission and reflection channels are zero T − + = T + − = R − + = R + − = 0.In Fig. 10 we plot the transmission and reflection prob-abilities for a potential strength v = 1 . γ and inter-layer bias δ = 0 . γ . The shift in the bands of the top(white) and bottom (red) layer of 2SL is due to the inter-layer bias which couples the two Dirac cones as shownin Eq. (8). Therefore, the suppression of the scatter-ing transmission and reflection probabilities due to theprotected cone transport does not hold anymore. It is,therefore, possible that scattering between different conestakes place as clarified in the second row of Fig. 10.
4. Conductance
The conductance of two and three-block systems isshown in Figs. 11 and 12, respectively. For the two
FIG. 9: (Colour online) The same as in Fig. 8, but now with v = 1 . γ . systems 2SL-AA and AA-2SL with pristine AA-BL and2SL, the conductance for different channels is shown inFigs. 11(a, b). It shows that the conductance of thesetwo systems are identical. Referring to Figs. 4(a, c)we notice that the transmission probabilities for pristine2SL-AA and AA-2SL are quite different. However, thecorresponding conductances (see Fig. 11) exhibit timereversal symmetry in spite of the fact that the domainwall separates two different systems. This is a strongpoint which can be verified experimentally even in thecase of zero electrostatic potential.Adding an electrostatic potential to one of the twosides leads to different behavior in the conductance ofthe above mentioned two systems as depicted in Figs.11(c,d). In Fig. 11(c) the charge carriers incident from2SL and impinging on AA-BL whose bands are shiftedby v . Each conductance channel gives zero at E = 0 dueto the absence of propagating states in the 2SL at thisenergy, even though there are propagating states avail-able in AA-BL corresponding to two cones. We notealso that G ± b = G ± t are almost zero at upper and lowercones v ± γ as a result of the absence of states at thesepoints as seen in Fig. 11(c). In Fig. 11(d) we see that0 FIG. 10: (Colour online) The same as in Fig. 8, but nowwith v = 1 . γ and δ = 0 . γ . Red and white dashed curvescorrespond to the bands of bottom and top layers of 2SL,respectively, while the black dashed curves are the AA-BLbands. the conductance of different channels is not zero in con-trast to the previous case because here at E = 0 thereare propagating states available in both AA-BL and 2SL.Furthermore, all channels have one minimum, due to thelack of states, at E = v which corresponds to the Diraccone in 2SL shifted by v while G t/b − has also anotherminimum at the upper cone E = γ as shown in Fig.11(d). Finally, for comparison we add in Figs. 11(e, f)the conductance that will be measured in the absence ofa domain wall for 2SL-2SL and AA-AA junctions with v = 0 (blue curves). Our results indicate that domainwalls are experimentally identifiable channels even in theabsence of a gate. As a reference we also calculate thetotal conductance in the presence of an electrostatic po-tential ( v = 1 . γ ) as shown with black curves in Figs.11(e, f) which corresponds, in this case, to the usual p-njunctions in single-layer graphene and AA-BL, respec-tively.The conductance of three-block systems is shown inFig. 12 where left and right panels correspond to AA- 2SL-AA and 2SL-AA-2SL structure, respectively. Pro-tected cone transport leads to zero conductance in thescattered channels G + − = G − + = 0 as shown in Fig. 12(a).A close inspection also reveals that G −− = G ++ at E = 0with finite and non-zero values, regardless of the factthat in the 2SL region there are no available propagat-ing states. This is attributed to the evanescent modes in2SL at E = 0 which are responsible for ballistic transportin graphene . We thus also expect that G −− (red curvein Fig. 12(a)) should be exactly zero at the Dirac cone E = γ as a result of the absence of propagating statesin the leads at this energy.By shifting the bands of 2SL using a local potentialwith strength v = 1 . γ , a local minimum appears inthe conductance G T at E = v which corresponds to theposition of the charge-neutrality point in 2SL as shownin Fig. 12(c). This minimum can be obtained by align-ing the upper cone in AA-BL and the Dirac cone in 2SLsuch that they are located at the same energy, this canbe achieved by choosing v = γ . The main differenceintroduced by applying an inter-layer bias is the brokenprotected cone transport where now G − + = G + − (cid:54) = 0 asdepicted in Fig. 12(e). For completeness, we performedsimilar calculations but now with 2SL as the leads (2SL-AA-2SL) and the results for the conductance with pris-tine, gated and biased AA-BL are shown in Figs. 12(b, d,f), respectively. Here, all conductance channels are zeroat E = 0 such that G tt = G bb and G bt = G tb as shown inFigs. 12(b, d). Similarly, the main features in Fig. 12(f)are in qualitative agreement with those shown in Figs.12(b, d) but now the tunnelling equivalence through thesame channel is broken so that G tt (cid:54) = G bb . This is a di-rect consequence of the perpendicular electric field whichleads to the breaking of the inter-layer sublattice equiva-lence. The peaks appearing in the total conductance aredue to the finite size of the AA-BL region. B. AB-Stacking
1. 2SL-AB/AB-2SL
In this section, we evaluate how the stacking of theconnected region changes the transport properties acrossa domain wall. The angle-dependent transmission andreflection probabilities for pristine systems 2SL-AB areplotted in Fig. 13(a). The charge carriers can be inci-dent from the two layers in the 2SL structure and im-pinge on AB-BL where, depending on their energy, theycan access only one propagating mode k + or two k ± ifthe energy is large enough. Scattering from the top orbottom layer of 2SL into one of these modes is equiv-alent T ± t = T ± b as well as backscattering R t ( b ) t = R b ( t ) b and hence, as before, layer symmetry is preserved (seeFig. 13(a)). In Fig. 13(b) we show results with theAB-BL region subjected to an electrostatic potential ofstrength v = 1 . γ . Surprisingly, we see that the layersymmetry is broken and an asymmetric feature with re-1 - AA G b + = G t + G b - = G t - ( a ) G T / γ G l / G L y AA - G + t = G + b G - t = G - b ( b ) G T / γ G l / G L y - AA ( c ) G b + = G t + G b - = G t - G T / γ G l / G L y G T AA - G + t = G + b G - t = G - b ( d ) / γ G l / G L y ( e ) - / γ G l / G L y AA - AA ( f ) / γ G l / G L y FIG. 11: (Colour online) Conductance of two-block system fordifferent magnitudes of the applied gate: (a, b) v = δ = 0,(c, d) v = 3 γ / δ = 0 . G T is the total conductanceobtained by summation of all possible channels, (e, f) the totalconductance for 2SL-2SL and AA-AA junctions, respectively,with v = 0 (blue curves) and v = 1 . γ (black curves). spect to normal incidence shows up in the transmis-sion and non-scattered reflection probabilities, see Ap-pendix A, such that [ T /R ]( φ ) = [ T /R ]( − φ ) . For exam-ple, T ± b ( φ ) = T ± t ( − φ ) as well as the non-scattered reflec-tion channels R bb ( φ ) = R tt ( − φ ) as discussed in Sec. III.This asymmetric feature can be understood by resortingto the bands on both sides of the junction, where dueto the electrostatic potential the band alignment of 2SLand AB-BL is altered. In this case, the center of theAB-BL band is shifted upwards in energy with respectto the crossing of the 2SL energy bands. The originof such asymmetry is a direct consequence of the asym-metric coupling in AB-BL which leads to shifting of thebands by γ . Therefore, at low energy | E − v | < γ there is only one propagating mode k + (i.e one type ofcharge carrier ) and consequently only T + b ( t ) is available.For larger energy, on the other hand, there are two modesavailable giving rise to four channels T ± b ( t ) .The angular asymmetry feature is present only in theregion in the ( E, k y )-plane where there is only one prop-agating mode. This can be also understood as a mani-festation of the asymmetric amplitude of the wave func-tion in the AB-BL side due to the evanescent modes inthis region . The theory of tunnelling through an in-terface of monolayer and bilayer was presented earlier and such asymmetry was noticed as well. Moreover, in AA - - AA G -+ = G +- G T G ++ G -- ( a ) / γ G l / G L y G T - AA - G bt = G tb G bb = G tt ( b ) / γ G l / G L y AA - - AA G T G ++ G -- G -+ = G +- ( c ) / γ G l / G L y - AA - G T G bt = G tb G bb = G tt ( d ) / γ G l / G L y AA - - AA G tb = G bt G T G ++ G -- ( e ) / γ G l / G L y - AA - G T G tb = G bt G bb ( f ) G tt / γ G l / G L y FIG. 12: (Colour online) Conductance of three-block systemwith different magnitudes of the applied gate: (a, b) v = δ =0, (c, d) v = 3 γ / δ = 0 and (e, f) v = 3 γ / δ = 0 . γ . G T is the total conductance obtained by summation of allpossible channels. our case there are two single layer graphene sheets con-nected to the bottom and top layers of the bilayer sys-tem but the asymmetric feature in Ref. [26] will be re-covered when considering only one propagation channel.For instance, the transmission probabilities T ± t and T ± b presented in Fig. 13(b) show the same asymmetric fea-tures discussed in Ref. [26]. This asymmetry feature isreversed in the other valley, so that the total transmis-sion or reflection averaged over both layers is symmetricas can be seen from Fig. 13(b). However, this valley-dependent angular asymmetry could also be used for thebasis of a layer-dependent valley-filtering device as pro-posed in other works .The above analogy, which is discriminating betweenthe presence of one or two modes, applies also to thenon-scattered reflection probabilities R bb and R tt . Thesenon-scattered currents are carried by the states localizedon the disconnected sublattices α and β , as seen in Fig.1. In that case, there is one traveling mode and thus,inherently, a layer asymmetric feature will be present. Incontrast, for the scattered channels R tb and R bt the chargecarriers must jump between the layers of AB-BL. This oc-curs through the localized states on the connected sub-lattices α and β where there are two travelling modesand, hence, these probabilities exhibit layer symmetry asshown in Fig. 13(b). In the AB-2SL configuration, wherecharge carriers incident from the AB-BL impinge on the2 - ° - ° - ° ° ° R bb R bt T b + T b - ( a ) - AB v = - ° - ° - ° ° ° T t - R tt T t + ( b ) - AB v = γ - ° - ° - ° ° ° R -- T - b T + b R ++ ( c ) AB - v = - ° - ° - ° ° ° R +- T + t ( d ) AB - v = γ FIG. 13: (Colour online) The angle-dependent transmissionand reflection probabilities through (a, b) 2SL-AB and (c, d)AB-2SL junctions. The systems in (b, d) are the same as in(a, c), respectively, but where the right side of the junction issubjected to an electrostatic potential of strength v = 1 . γ . In (a) E = 1 . γ for all channels while in (b) E = 1 . γ for T + b ( t ) and E = 0 . γ for the rest of the channels and in (c,d) E = (0 . , . γ for R ++ /T b ( t )+ and R −− /T b ( t ) − , respectively.We choose energy values in (b, d) such that they correspondto only one propagating mode in the AB-BL region. T t ± = T b ± with partial reflection associatedwith the non-scattered channels R −− and R ++ . While forthe scattered channels R − + and R − + are almost zero. Thisis due to efficient transmission resulting from the absenceof the electrostatic potential in the 2SL. An electrostaticpotential of strength v = 1 . γ induces a scattering be- FIG. 14: (Colour online) Density plot of the transmission andreflection probabilities through 2SL-AB-2SL as a function ofFermi energy and transverse wave vector k y with v = δ = 0 . tween the two modes in the reflection channels so thatnow R + − = R − + (cid:54) = 0 as depicted in Fig. 13(d). In addition,it breaks the band alignment and gives rise to the layerasymmetry feature in the transmission probabilities T b ( t )+ where only one travelling mode exists i.e. E < γ . Thus, T b ( t ) − always preserves layer symmetry in this case, seeFig. 13(d), because the mode k − exists for E > γ wherealso the mode k + is available as discussed above. Thisis also the same reason that configurations consisting ofAA-BL always preserve layer symmetry. Indeed, AA-BLdoes not have a region in the ( E, k y )-plane with only onepropagating mode, and there are always two travellingmodes for all energies.
2. 2SL-AB-2SL
Different configurations have been proposed to connecta single layer to the AB-stacked bilayer graphene .Now, two SL are connected to the AB-stacked bilayer,see Fig. 2(a). In Fig. 14 we show the dependenceof the transmission and reflection probabilities on thetransverse wave vector k y and the Fermi energy. It ap-pears that all channels are symmetric with respect to nor-mal incidence since the Dirac cones of AB and 2SL arealigned. It also implies that scattered and non-scatteredchannels of the transmission and reflection are equivalentsuch that ( T /R ) tb = ( T /R ) bt and ( T /R ) tt = ( T /R ) bb (seeFig. 14).Another interesting feature of this configuration is thatfor E < γ the scattered and non-scattered transmissions3 FIG. 15: (Colour online) The same as in Fig. 14, but nowwith v = 3 γ / are equal T ji = T ii . In this energy regime such device canbe used as an electronic beam splitter .Fig. 15 displays the same plot as in Fig. 14 but withan electrostatic potential on the AB-BL region. Thereis an important difference as compared to the pristineAB-BL case, the layer symmetry is broken such that T bt ( k y ) = T tb ( − k y ) as clarified in Fig. 15. This can bealso understood by pointing out that charge carriers scat-tered from top to bottom when moving from left to rightin the K valley are equivalent to charge carriers scatter-ing from bottom to top when moving oppositely in thesecond valley K (cid:48) .Introducing a finite bias ( δ >
0) to the AB-BL regionalong with an electrostatic potential ( v >
0) will shiftthe bands and opens a gap in the spectrum. As a resultof the presence of a strong electric field, the transmissionchannels are completely suppressed inside the gap due tothe absence of traveling modes as seen in Fig. 16. More-over, non-zero asymmetric reflection appears in the gapas well as a violation of the equivalence of non-scatteredtransmission channels. This is a result of the breakingof inter-layer sublattice equivalence . In addition, somelocalized states appear inside the “Mexican hat” of the FIG. 16: (Colour online) The same as in Fig. 14, but nowwith v = 3 γ / δ = 0 . γ . New localized states appear insidethe “Mexican hat” shape of the low energy bands of AB-BLdue to the strong gate potential. low energy bands where they are pushed by the strongelectric field ( δ = 0 . γ ), see Fig. 16.There is a link between the transmission probabili-ties of our system 2SL-AB-2SL and those investigatedby Gonz´alez et al. . The channels T bb and T tb are qual-itatively equivalent to those obtained in Ref. [38]. Forexample, T tb shows electron-hole ( e − h ) and δ → − δ sym-metry whereas T bb exhibits another symmetry which canbe obtained under the exchange ( e, δ ) ↔ ( h, − δ ). Theresults in Fig. 17 are in good agreement with those ofRef. [38] where we fix v = 0 and d = 25 nm.4 FIG. 17: (Colour online) Transmission probabilities as func-tion of Fermi energy and bias for normal incidence.FIG. 18: (Colour online) Density plot of the transmission andreflection probabilities through AB-2SL-AB as a function ofFermi energy and transverse wave vector k y with v = δ = 0and d = 25nm.
3. AB-2SL-AB
For leads composed of AB-BL where the intermediateregion is pristine 2SL, we show the results in Fig. 18 forthe transmission and reflection probabilities. Now chargecarriers will scatter between the different modes of the
FIG. 19: (Colour online) The same as in Fig. 18, but herewith v = 3 γ / AB-BL on the left and right leads as shown in Fig. 2(b).As expected, all channels are symmetric and as a resultof the finite size of the 2SL region, resonances appearin T as shown in Fig. 18. These so-called Fabry-P´erotresonances appear at quantized energy levels E nSL ( k y ) = (cid:114) k y + (cid:16) nπd (cid:17) . (21)This is the dispersion relation for modes confined in the2SL region with width d .The results presented in Fig. 18 reveal no scatteringbetween the two modes k + and k − and charge carriers areonly transmitted or reflected through the same channelfrom which they come from. Unexpectedly, introducingan electrostatic potential induces a strong scattering inthe reflection channels ( R − + = R + − (cid:54) = 0) and very weakscattering in the transmission channels ( T − + = T + − (cid:54) = 0),as seen in Fig. 19. When the 2SL are biased, the Diraccones at bottom and top layers will be shifted up (whitedashed lines) and down (red dashed lines) in energy, re-spectively (see Fig. 20). This bias will strengthen thecoupling between the two modes resulting in a strong5 FIG. 20: (Colour online) The same as in Fig. 18, but herewith v = 3 γ / δ = 0 . γ . Red and white dashed curvescorrespond to the bands of bottom and top layers of 2SLwhile the black dashed curves are the AB-BL bands. scattering between them. In addition, the inversion sym-metry is broken due to the bias leading to an asymmetrywith respect to normal incidence.
4. Conductance
The conductance of the two-block system consisting of2SL and BA-BL is shown in Fig. 21 for different valuesof the applied gate voltage. Figs. 21(a,b) reveal thatthe system where charge carriers are incident from the - AB G b + = G t + G b - = G t - ( a ) G T / γ G l / G L y AB - G + t = G + b G - t = G - b ( b ) G T / γ G l / G L y - AB G b - = G t - G b + = G t + G T ( c ) / γ G l / G L y AB - G T G + t = G + b G - t = G - b ( d ) / γ G l / G L y ( e ) - / γ G l / G L y AB - AB ( f ) / γ G l / G L y FIG. 21: (Colour online) Conductance of different junctionsfor different magnitudes of the applied gate: (a, b) v = δ = 0,(c, d) v = 3 γ / δ = 0, (e, f) the total conductance for 2SL-2SL and AB-AB junctions, respectively, with v = 0(bluecurves) and v = 1 . γ (black curves). G t ( b )+ = G + t ( b ) are contributing to the total conductance G T startingfrom E = 0 where the k + mode exists. On the contrary, G t ( b ) − = G − t ( b ) only contributes when E > γ where k − states are available and this appears as a sharp increasein G T at E = γ . On the other hand, considering anapplied electrostatic potential on the right side of thetwo-block system will break this equivalence as seen inFigs. 21(c,d). In addition, as a result of the shift ofthe Dirac cone in AB-BL (see Fig. 21(c)) or 2SL (seeFig. 21(d)) due to the electrostatic potential, all con-ductance channels are zero at E = v . Similar to theAA-BL case, the conductances of the pristine systems2SL-AB/AB-2SL (see Figs. 21(a, b)) clearly preservethe time reversal symmetry. Even though, both systemshave different transmission probabilities as can be seenfrom Figs. 13(a, c). We also show in Figs. 21(e, f) thetotal conductance in the absence of domain wall in 2SL-SL and AB-AB systems, respectively, for v = 0 (bluecurves) and v = 1 . γ (black curves). This shows thattransport channels in the presence of domain walls areexperimentally recognisable.In Fig. 22 we show the conductance in a 2SL-AB sys-tem as a function of the bias for transport using a singleFig. 22(a) or a double Fig. 22(b) mode. The results6 G t + G T G b + - - δ / γ G l / G L y ( a ) G t - G b - G b + G t + G T - - δ / γ G l / G L y ( b ) FIG. 22: (Colour online) Conductance across the 2SL-ABsystem as a function of the bias on the AB-BL with v = 0 . (a) and (b) correspond to the single and double modes regimewith E = 0 . γ and E = 1 . γ , respectively. With G ± T = G ± t + G ± b . G +- = G -+ G -- G ++ G T AB - - AB ( a ) / γ G l / G L y G bt = G tb G bb = G tt - AB - G T ( b ) / γ G l / G L y G T AB - - AB G -- G ++ G +- = G -+ ( c ) / γ G l / G L y - AB - G T G bb = G tt G tb = G bt ( d ) / γ G l / G L y G +- = G -+ G -- G ++ G T ( e ) AB - - AB / γ G l / G L y G tt G tb = G bt G bb G T - AB - ( f ) / γ G l / G L y FIG. 23: (Colour online) Conductance of different junctionsfor different magnitudes of the applied gate: (a, b) v = δ =0, (c, d) v = 3 γ / δ = 0 and (e, f) v = 3 γ / δ = 0 . γ . show that the contribution from the top and bottom lay-ers to the conductances have opposite behaviours as afunction of the inter-layer bias. The total conductance G T , however, has a convex form, increasing with the ap-plication of an inter-layer bias. From Fig. 22(b), on theother hand, we see that when a second mode is available,four channels contribute to the conductance and the to-tal conductance assumes a concave form, i.e. decreasingwith increasing inter-layer bias. This is a characteristicexperimental feature that can signal the presence of asecond mode of propagation.For the three-block system we show the conductanceof the configuration AB-2SL-AB and 2SL-AB-2SL in theleft and right columns of Fig. 23, respectively. The re- sulting conductance of the first configuration shows onlytwo non-zero channels G ++ and G −− , while the scatteredones G − + = G + − = 0 since T − + = T + − = 0 (see Fig. 23(a)).Furthermore, for low energy G T = G ++ since the mode k − is not available in this regime but it starts conduct-ing when E > γ . The applied electrostatic potentialon the 2SL keeps the scattered conductance channels atzero and a minimum in the conductance appears aroundthe shifted Dirac cone E = v of the 2SL as depicted inFig. 23(c). As pointed out before, if the Fermi energyapproaches the strength of the electrostatic potential, anon-zero minimum is present in the conductance becausecharge carriers can be transmitted through a width d of2SL via evanescent modes . In Fig. 23(f) this minimumdisappears and the conductance dramatically increasesat E = γ . This is because the bias will couple the twomodes and two additional scattered channels G − + and G + − start conducting. The resonant peaks resulting in theconductance, see Figs. 23(a,c,e), are due to the finitesize of the intermediate region and hence strictly dependon its width d .On the other hand, the conductance of the configura-tion 2SL-AB-2SL has different features. In Fig. 23(b)the four channels, in contrast to the previous configura-tion, start conducting from E = 0. This possess layersymmetry such that G tt = G bb and G tb = G bt . Of partic-ular importance is the equivalence of the four channelsfor E < γ while for E > γ charge carriers stronglyscatter between the layers ( i.e. G ji > G ii ) as shown inFig. 23(b). This equivalence of the four channels in theregime E < γ vanishes when an electrostatic potentialis applied ( v >
0) to the intermediate region as seenin Fig. 23(d). However, the scattered and non-scatteredconducting channels are still equivalent in this case where G t ( b ) t = G b ( t ) b with G ji > G ii for all energy ranges, see Fig.23(d).As discussed before, the most characteristic feature ofthe inter-layer bias in the AB-BL is the opening of a gapin the energy spectrum between v ± δ which is reflectedin the conductance as seen in Fig. 23(f). The resonantsharp peaks in the conductance near the edges of thegap result from the localized states inside the Mexicanhat of the low energy bands. Another consequence of theinter-layer bias is the breaking of the equivalence in thenon-scattered conducting channels where now G tt (cid:54) = G bb as seen in Fig. 23(f). C. AA-2SL-AB
Here we consider the case where the leads consist ofBL with different stackings separated by two uncoupledgraphene sheets. Such a structure can be formed if in thedecoupled region one of the graphene sheets has largerlattice constant, e.g. due to strain, leading to an inter-layer shift when the two layers couple.Notice that the inter-layer coupling strength γ differsfor the two bilayer structures. Their ratio is γ AA /γ AB ≈ FIG. 24: (Colour online)Density plot of the transmissionand reflection probabilities through AA-2SL-AB junction as afunction of Fermi energy and transverse wave vector k y with v = 1 . γ , δ = 0 and d = 25nm. The superimposed dashedcurves represent the bands of AB-BL(black), AA-BL(green)and 2SL (white), with γ being the inter-layer coupling ofAB-BL. / . To account for this difference the energy isnormalized to γ AB such that the upper Dirac cone ofpristine AA-BL is now located at E = 1 / E = 1 as in the previous sections. In the junction AA-2SL-AB the charge carriers incident from AA-BL andtransmitted through 2SL into AB-BL. The results for thetransmission and reflection probabilities of this junctionare shown in Fig. 24 for v = 1 . γ , δ = 0 and d = 25 nm.The carriers incident from lower( k + )/upper( k − ) Dirac ( a ) T -- R -+ T ++ R +- / γ T , R T -- T ++ T AB T AA ( b ) / l T AA - - AB ( c ) G T G -+ G +- G ++ G -- / γ G l / G L y G -+ G +- G -- G ++ G T AA - - AB ( d ) / γ G l / G L y AA - - AB ( e ) G +- G -+ G ++ G -- G T / γ G l / G L y FIG. 25: (Colour online) (a) Transmission and reflection prob-abilities for normal incidence for v = 3 γ / , δ = 0. (b)Transmission probabilities with normal incidence for AA-BL(AB-BL) n-p-n junction, green (black) curves. Blue (red)curves are the non-zero channels T ++ ( T −− ) in AA-2SL-AB.All energies are considered to be less than the electrostaticpotential strength. Conductance of AA-2SL-AB junction fordifferent magnitudes of the applied gate: (c) v = δ = 0, (d) v = 3 γ / δ = 0, (e) v = 3 γ / δ = 0 . γ , with γ beingthe inter-layer coupling of AB-BL. cones in AA-BL can be transmitted into one of the modes( k + or k − ) in the AB-BL, see Fig. 2(e). On the otherhand, the reflection process occurs between the intra- orinter-cone in the AA-BL.Remarkably, Fig. 24 shows that the scattered trans-mission probabilities are very small and that almost alltransmission is carried by the non-scattered channels.This is not immediately expected since a priori the k + -mode in AA-BL is not related to the k + -mode in AB-BL. However, both modes have the same parity underin-plane inversion, showing that this feature is robustagainst variations in the inter-layer coupling.In contrast to the AA-2SL-AA junction where the scat-tering between lower and upper cones is forbidden in caseof zero bias, here the two cones are coupled even withoutbias. This results in non-zero reflection in the scatteredchannels R − + and R + − .For normal incidence, the scattered transmission ( T − + and T + − ) and the non-scattered reflection ( R ++ and R −− )channels are zero (see Fig. 24) because in that case boththe AA and AB Hamiltonian are block diagonal in theeven and odd modes basis. Now, we can investigate Kleintunnelling when transitioning in-between the two types8of stacking. For this, we show the non-zero channels oftransmission and reflection for normal incidence in Fig.25(a). We find that in contrast with the AA-2SL-AAcase, perfect Klein tunnelling does not occur in the junc-tion AA-2SL-AB. However, as shown in Fig. 25(b), wedo find that the transmission probability does not de-pend on the length or even presence of the 2SL region,in contrast to the previous cases with two domain walls.For δ (cid:54) = 0 the coupling between the different modes isstrengthened and, hence, strong scattering in the trans-mission and reflection channels occurs. Furthermore, thesymmetry with respect to normal incidence in the reflec-tion and transmission channels is broken.The conductance for the discussed structure is shownin Figs. 25(c, d, e) for ( v = δ = 0) , ( v = 1 . γ , δ = 0)and ( v = 1 . γ , δ = 0 . γ ), respectively. For pristine2SL, the dominant channels are G ++ and G −− . Noticethat the latter one starts conducting only when E > γ and this shows up as a rapid increase in the total con-ductance G T at E = γ . The scattered channels G − + and G + − are only weakly contributing to the total con-ductance as a result of weak coupling of the modes. Incontrast to the junctions AA(AB)-2SL-AA(AB), in thiscase the scattered channels of the conductance are notequivalent G − + (cid:54) = G + − , see Fig. 25(c,d). This is becausethe scattering occurs between modes in bilayer grapheneof different stackings. The electrostatic potential intro-duces a minimum at E = v in the total conductance dueto the absence of propagating states at this energy in the2SL, see Fig. 25(d). Biasing the intermediate region(2SL) of the junction AA-2SL-AB provides propagatingstates at E = v , and hence removing the minima in G T as shown in Fig. 25(e). In addition, the contribution ofthe scattered channels G − + and G + − becomes more pro-nounced as a result of the strong coupling between themodes induced by the bias.Finally, notice that the counterpart junction AB-2SL-AA, represents the time-reversal case of the system dis-cussed above. We have verified that the transmissionchannels are equivalent in the absence of a bias. In thepresence of a bias, the angular symmetry is broken and,consequently, the reversed junction features the oppositeangular asymmetry, preserving time-reversal invariance. V. SUMMARY AND CONCLUSION
Using the four-band model we obtained the conduc-tance, transmission and reflection probabilities throughsingle and double domain walls separating two single lay-ers and AA/AB-stacked bilayer graphene. We discussedin detail the scattering mechanism from detached lay-ers to bilayer graphene and presented compact analyt-ical formulae for the transmission probabilities. Theseresults showed that one can find the inter-layer couplingstrength solely through measuring the conductance.We found that an electrostatic potential applied to AB-BL, in an 2SL-AB junction, breaks the layer symmetry in the single-valley transmission probability channels. Suchasymmetry originates from the asymmetric coupling inAB-BL and arises as a consequence of the mismatch inenergy between the 2SL and AB-BL Dirac cones causedby the electrostatic potentials applied to the AB-BL re-gion. Layer asymmetry exists when only one propagatingmode is present and hence is not seen in configurationsconsisting of AA-BL where the entire energy range is as-sociated with two transport channels.We have also evaluated the robustness of chirality-induced properties, such as Klein tunnelling and anti-Klein tunnelling, to scattering on domains without inter-layer coupling. We found that in domain walls separat-ing 2SL and AA-BL, Klein tunnelling is still preserved.On the other hand, for domain walls separating 2SL andAB-BL, the well known anti-Klein tunnelling in AB-BL isnot preserved any more, but neither is Klein tunnellingitself. Moreover, in two domain walls separating threeregions whose interlayer coupling is all different, i.e. theAA-2SL-AB case, we find that although perfect Kleintunnelling does not hold, the tunnelling does not dependon the thickness of the 2SL region either. This remark-able effect is attributed to a conservation of parity of themodes.Furthermore, we have found that a strong gate poten-tial difference allows some states to be localized insidethe Mexican hat of the low energy bands in the AB-BL.Those states contribute to the conductance and appearas sharp peaks at the two edges of the gap. We showedthat scattering between these modes, in the transmissionchannels, is not allowed in the configuration (AA/AB)-2SL-(AA/AB). However, such scattering can be inducedby applying an inter-layer bias on the 2SL which in addi-tion to shifting the bands of the top and bottom layers of2SL, also couples the modes. In contrast, we showed thatthe two modes of AA-BL are coupled even without bias-ing the system in the junction AA-2SL-AB and revealedthat the latter junction is equivalent to the AB-2SL-AA.In order to limit the number of parameters, throughthis article we only considered abrupt domain walls, how-ever, the results are robust against smoothness of thedomain walls. Our study reveals that the presence of the local domainwall in bilayer graphene samples change the transportproperties significantly. Our results may shed light on thedesign of electronic devices based on bilayer graphene.Finally, we showed that for a given sample with unknownsizes of local stacking domains, the average inter-layercoupling can be estimated through quantum transportmeasurements.
Acknowledgments
HMA and HB acknowledge the Saudi Center for The-oretical Physics (SCTP) for their generous support andthe support of KFUPM under physics research groupprojects RG1502-1 and RG1502-2. This work is sup-9ported by the Flemish Science Foundation (FWO-Vl) bya post-doctoral fellowship (BVD).
Appendix A: Functions definitions
The transmission probabilities are calculated by ap-plying appropriate boundary conditions at the 2SL-BLinterfaces together with the transfer matrix. After somecumbersome algebra, we obtain for 2SL-AB T ± j = 4Re( k ± ) η (cid:104) η + (Im( k ∓ ) + κ j v sin φ ) (cid:105) C + (cid:80) m =1 C m cos( mφ ) , (A1)where C = 2 (Im( k ∓ )Re( k ± )) + (cid:15) (cid:0) Im ( k ∓ ) + Re ( k ± ) (cid:1) + Γ ,Γ = 2 v − v E + 5 v E − v E + E , C = − (cid:15) Re( k ± ) (cid:2) (cid:0) v + Im ( k ∓ ) (cid:1) − v E + 3 E (cid:3) ,C = (cid:15) (cid:0) Im ( k ∓ ) + Re ( k ± ) (cid:1) + Γ ,Γ = E (cid:0) − v + 6 v E − v E + E (cid:1) , C = Re( k ± ) E (cid:0) v − v E + E (cid:1) , C = E ( E − v ) .Similarly, the transmission probabilities for the AB-2SL system are obtained as T j ± = 4Re( k j ) k ± λ [ µ ± + κ j v sin φ Im( k ∓ )] | Q ± | , (A2) µ ± = (cid:15) (cid:16) Im ( k ∓ ) + E (cid:17) − E ( ± E )( E + v ) sin φ (cid:112) E ( ± E ) ,λ = E (cid:112) E ( ± E ) ,Q ± = [ z − z ( k ± + i Im( k ∓ )) + z k ± Im( k ∓ )] , with z = 2 i [ v α − ik j E ] [ α ( − ik j + α ) + (cid:15)E ], z = E (cid:104) ( ik j + α ) − (cid:15) (cid:105) ,and finally z = 2 (cid:15) [ ik j + α ],where α = √ E ± E sin φ. ∗ Electronic address: [email protected] † Electronic address: [email protected] T. Ohta, A. Bostwick, T. Seyller, K. Horn, and E. Roten-berg, Science , 951 (2006). E. V. Castro, K. S. Novoselov, S. V. Morozov, N. M. R.Peres, J. M. B. L. dos Santos, J. Nilsson, F. Guinea, A. K.Geim, and A. H. C. Neto, J. Phys.: Condens. Matter ,175503 (2010). Zhang, Yuanbo, T.-Ta Tang, C. Girit, Z. Hao, M. C. Mar-tin, A. Zettl, M. F. Crommie, Y. R. Shen, and F. Wang,Nature (London) , 820 (2009). D. R. da Costa, A. Chaves, S. H. R. Sena, G. A. Farias,and F. M. Peeters, Phys. Rev. B , 045417 (2015). J. W. Gonz´alez, H. Santos, M. Pacheco, L. Chico, and L.Brey, Phys. Rev. B , 195406 (2010). M. Zarenia, A. Chaves, G. A. Farias, and F. M. Peeters,Phys. Rev. B , 245403 (2011). I. Silvestre, A. W. Barnard, S. P. Roberts, P. L. McEuen,and R. G. Lacerda, Appl. Phys. Lett. , 153105 (2015). X. Zhao, B. Zheng, T. Huang and C. Gao, Nanoscale ,9399 (2015). S.-H. Ji, J. B. Hannon, R. M. Tromp, V. Perebeinos, J.Tersoff, and F. M. Ross, Nat. Mater. , 114 (2012). J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H.Castro Neto, Phys. Rev. Lett. , 256802 (2007). E. J. Mele, Phys. Rev. B , 161405(R)(2010). J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H.Castro Neto, Phys. Rev. B , 155449 (2012). A. O. Sboychakov, A. L. Rakhmanov, A. V. Rozhkov, andF. Nori, Phys. Rev. B , 075402 (2015). M. Van der Donck, F.M. Peeters, and B. Van Duppen,Phys. Rev. B , 247401 (2016). M. Van der Donck, C. De Beule, B. Partoens, F.M. Peeters,and B. Van Duppen, 2D Mater. , 35015 (2016). L.-Jing, S. Li, J. Qiao, W. Zuo, and L. He,arXiv:1509.04405 (2015). W. Yan, S.Y. Li, L.J. Yin, J. B. Qiao, J.C. Nie, and L. He,Phys. Rev. B , 195408 (2016). K. W. Clark, X.-G. Zhang, G. Gu, J. Park, G. He, R. M.Feenstra, and A.-P. Li, Phys. Rev. X , 011021 (2014). M. Pelc, W. Jask´olski, A. Ayuela, and L. Chico, Phys. Rev.B , 085433 (2015). L. Ju, Z. Shi, N. Nair, Y. Lv, C. Jin, J. Velasco, C. O.-Aristizabal, H.A. Bechtel, M.C. Martin, A. Zettl, J. Ana-lytis, and F. Wang, Nature (London) , 650 (2015). S. Li, H. Jiang, J. Zhou, H. Liu, F. Zhang, and L. He,arXiv:1609.03313 (2016) F. Giannazzo, I. Deretzis, A. La Magna, F. Roccaforte,and R. Yakimova, Phys. Rev. B , 235422 (2012). A. Nourbakhsh, M. Cantoro, A. V. Klekachev, G. Pourtois,T. Vosch, J. Hofkens, M. H. van der Veen, M. M. Heyns,S. De Gendt, and B. F. Sels, J. Phys. Chem. C , 16619(2011). B. Z. Rameshti, M. Zareyan, and A. G. Moghaddam, Phys.Rev. B , 085403 (2015). T. Nakanishi, M. Koshino, and T. Ando, J. Phys.: Conf.Ser. , 012021 (2011). T. Nakanishi, M. Koshino, and T. Ando, Phys. Rev. B ,125428 (2010). J. Nilsson, A. H. Castro Neto, F. Guinea, and N. M. R. Peres, Phys. Rev. B , 165416 (2007). C. P. Puls, N. E. Staley, and Y. Liu, Phys. Rev. B ,235415 (2009). M. Koshino, T. Nakanishi, and T. Ando, Phys. Rev. B ,205436 (2010). J. Tian, Y. Jiang, I. Childres, H. Cao, J. Hu, and Y. P.Chen, Phys. Rev. B , 125410 (2013). D. Yin, W. Liu, X. Li, L. Geng, X.Wang, and P. Huai,Appl. Phys. Lett. , 173519 (2013). M. Berahman, M. Sanaee, and R. Ghayour, Carbon ,411 (2014). D. Dragoman, J. Appl. Phys. , 214312 (2013). Y. Hasegawa and M. Kohmoto, Phys. Rev. B , 125430(2012). H. Li, H. Li, Y. Zheng, and J. Niu, Physica B , 1385(2011). Z. X. Hu and W. Ding, Phys. Lett. A , 610 (2012). Y. Wang, J. Appl. Phys. , 164317 (2014). J. W. Gonz´alez, H. Santos, E. Prada, L. Brey, and L.Chico, Phys. Rev. B , 205402 (2011). B. Wang, M. Huang, N. Y. Kim, B. V. Cunning, Y. Huang,D. Qu, X. Chen, S. Jin, M. Biswal, X. Zhang, S. H. Lee,H. Lim, W. J. Yoo, Z. Lee, and R. S. Ruoff, Nano Lett. , 1467 (2017). J. C. Rode, C. Belke, H. Schmidt, and R. J. Haug, arXiv:1608.08133v1 (2016). M. Schmitz, S. Engels, L. Banszerus, K. Watanabe,T. Taniguchi, C. Stampfer, and B. Beschoten, p.arXiv:1704.04352 (2017). W. Yan, S.-Y. Li, L.-J. Yin, J.-B. Qiao, J.-C. Nie, andL. He, Phys. Rev. B , 195408 (2016). Y. Hao, L. Wang, Y. Liu, H. Chen, X. Wang, C. Tan,S. Nie, J. W. Suk, T. Jiang, T. Liang, J. Xiao, W. Ye, C. R.Dean, B. I. Yakobson, K. F. McCarty, P. Kim, J. Hone,L. Colombo, and R. S. Ruoff, Nat. Nanotechnol. , 426(2016). L.-J. Yin, W.-X. Wang, Y. Zhang, Y.-Y. Ou, H.-T. Zhang,C.-Y. Shen, and L. He, Phys. Rev. B , 081402(R) (2017). L.-J. Yin, H. Jiang, J.-B. Qiao, and L. He, Nat. Commun. , 11760 (2016). M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Nat. Phys. , 620 (2006). N. Stander, B. Huard, and D. Goldhaber-Gordon, Phys.Rev. Lett. , 026807 (2009). A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov,and A. K. Geim, Rev. Mod. Phys. , 109 (2009). Z. Q. Li, E. A. Henriksen, Z. Jiang, Z. Hao, M. C. Martin,P. Kim, H. L. Stormer, and D. N. Basov, Phys. Rev. Lett. , 037403 (2009). B. Van Duppen and F. M. Peeters, Phys. Rev. B ,205427 (2013). I. Lobato and B. Partoens, Phys. Rev. B , 165429 (2011). Y. Xu, X. Li, L. Reining, and J. Dong, Nanotechnology , 065711 (2010). M. Sanderson, Y. S. Ang, and C. Zhang, Phys. Rev. B ,245404 (2013). M. Barbier, P. Vasilopoulos, and F. M. Peeters, Phys. Rev.B. , 235408 (2010). B. Van Duppen, S. H. R. Sena, and F. M. Peeters, Phys.Rev. B , 195439 (2013). Ya. M. Blanter and M. Buttiker, Phys. Rep. , 1 (2000). S. B. Kumar and J. Guo, Appl. Phys. Lett. , 163102(2012). B. Van Duppen and F. M. Peeters, Appl. Phys. Lett. ,226101 (2012). H. M. Abdullah, M. Zarenia, H. Bahlouli, F. M. Peetersand B. Van Duppen, EPL , 1 (2016). I. Snyman and C. W. J. Beenakker, Phys. Rev. B ,045322 (2007). D. R. da Costa, A. Chaves, S. H. R. Sena, G. A. Farias,and F. M. Peeters, Phys. Rev. B , 045417 (2015). A. Rycerz, J. Tworzyd(cid:32)lo, and C. W. J. Beenakker, Nat.Phys. , 172 (2007). L. R. F. Lima, A. R. Hern´andez, F. A. Pinheiro and C.Lewenkopf, J. Phys.: Condens. Matter , 505303 (2016). P. Brandimarte, M. Engelund, N. Papior, A. Garcia-Lekue, T. Frederiksen, and D. S´anchez-Portal,arXiv:1611.03337(2016). M. R. Masir, P. Vasilopoulos, and F. M. Peeters, PhysicalReview B82