Origin of nonlocal resistance in multiterminal graphene on hexagonal-boron-nitride: Fermi surface edge currents rather than Fermi sea topological valley currents

The recent observation [R. V. Gorbachev et al., Science {\bf 346}, 448 (2014)] of nonlocal resistance $R_\mathrm{NL}$ near the Dirac point (DP) of multiterminal graphene on aligned hexagonal boron nitride (G/hBN) has been interpreted as the consequence of topological valley Hall currents carried by the Fermi sea states just beneath the bulk gap $E_g$ induced by the inversion symmetry breaking. However, the valley Hall conductivity $\sigma^v_{xy}$, quantized inside $E_g$, is not directly measurable. Conversely, the Landauer-B\"{u}ttiker formula, as numerically exact approach to observable nonlocal transport quantities, yields $R_\mathrm{NL} \equiv 0$ for the same simplistic Hamiltonian of gapped graphene that generates $\sigma^v_{xy} \neq 0$. We combine ab initio with quantum transport calculations to demonstrate that G/hBN wires with zigzag edges host dispersive edge states near the DP that are absent in theories based on the simplistic Hamiltonian. Although such edge states exist also in isolated zigzag graphene wires, aligned hBN is required to modify their energy-momentum dispersion and generate $R_\mathrm{NL} \neq 0$ near the DP persisting in the presence of edge disorder. Concurrently, the edge states resolve the long-standing puzzle of why the highly insulating state of G/hBN is rarely observed. We conclude that the observed $R_\mathrm{NL}$ is unrelated to Fermi sea topological valley currents conjectured for gapped Dirac spectra.

The recent observation [R. V. Gorbachev et al., Science 346, 448 (2014)] of nonlocal resistance RNL near the Dirac point (DP) of multiterminal graphene on aligned hexagonal boron nitride (G/hBN) has been interpreted as the consequence of topological valley Hall currents carried by the Fermi sea states just beneath the bulk gap Eg induced by the inversion symmetry breaking. However, the valley Hall conductivity σ v xy , quantized inside Eg, is not directly measurable. Conversely, the Landauer-Büttiker formula, as numerically exact approach to observable nonlocal transport quantities, yields RNL ≡ 0 for the same simplistic Hamiltonian of gapped graphene that generates σ v xy = 0. We combine ab initio with quantum transport calculations to demonstrate that G/hBN wires with zigzag edges host dispersive edge states near the DP that are absent in theories based on the simplistic Hamiltonian. Although such edge states exist also in isolated zigzag graphene wires, aligned hBN is required to modify their energy-momentum dispersion and generate RNL = 0 near the DP persisting in the presence of edge disorder. Concurrently, the edge states resolve the long-standing puzzle of why the highly insulating state of G/hBN is rarely observed. We conclude that the observed RNL is unrelated to Fermi sea topological valley currents conjectured for gapped Dirac spectra.
The recent measurements [1] of a sharply peaked nonlocal resistance R NL in a narrow energy range near the Dirac point (DP) of multiterminal graphene on hexagonal boron nitride (G/hBN) heterostructures have been interpreted as the manifestation of the valley Hall effect (VHE) [2][3][4][5]. In this interpretation, injecting charge current I 3 between leads 3 and 4 of the device illustrated in Fig. 1 generates a VH current in the first crossbar flowing from lead 1 to lead 2, which traverses the channel of length L ( 4 µm in the experiments [1]), and it is finally converted into a nonlocal voltage V NL between leads 5 and 6 by the inverse VHE in the second crossbar. The corresponding nonlocal resistance R NL = V NL /I 3 has been observed previously also near the DP in multiterminal graphene due to an external magnetic field inducing edge states in the quantum Hall regime or the Zeeman spin Hall effect at higher temperatures [6][7][8][9][10], as well as due to the spin Hall effect [10][11][12]. However, none of these mechanisms is operational in the experiment of Ref. [1].
Instead, the physics of graphene on hBN with their crystallographic axes aligned is expected to be governed by the broken spatial inversion symmetry due to different potentials on two triangular sublattices of carbon atoms induced by the hBN substrate. This opens a gap E g at the DP of two valleys K and K in the band structure of an infinite two-dimensional sheet of graphene, where ab initio calculations have estimated E g 58 meV [13]. In addition, the finite Berry curvature [14] of opposite sign at the two valleys was predicted to generate valley-dependent transverse conductivities, σ K xy = e 2 /h and σ K xy = −e 2 /h [2][3][4]. The VH current is characterized by the transverse VH conductivity, σ v xy = σ K xy − σ K xy = 2e 2 /h, and zero transverse charge conductivity, σ xy = σ K xy + σ K xy ≡ 0, within the arXiv:1706.09361v2 [cond-mat.mes-hall] 15 Jan 2018 gap [3,4,15]. Such charge-neutral valley currents, denoted also as "topological" [1][2][3][4] due to the involvement of the Berry curvature hot spots in reciprocal space, are not conserved but are expected to be long-ranged when the intervalley scattering is weak [15]. However, σ v xy is not directly observable, and semiclassical transport theories [16] attempting to connect σ v xy to observable nonlocal resistance [1,4,5] are not applicable [17] to electronic transport near the DP. In the case of G/hBN, this is further emphasized [18] by the presence of the gap forcing electrons to tunnel through the system, which is a phenomenon with no classical analog. The Landauer-Büttiker (LB) formalism [19][20][21], which offers quantum transport framework to compute V NL and R NL directly and has been used for decades to model nonlocal transport measurements [22,23], yields R NL ≡ 0 near the DP in Fig. 2(c) in multiterminal geometries whose channel length is larger than its width (L/W 4 in the experiments of Ref. [1]). The numerically exact result in Fig. 2(c), for a G/hBN system described by the same simplistic Hamiltonian employed [1-100 120 10 20 30 y (Å) 40 60 80 x (Å) 10 20 30 x (Å) 10 20 30 x (Å) 10

4, 15
] to obtain σ v xy = 0, is clearly incompatible with the interpretation of R NL = 0 based on the picture of topological valley currents carried by the Fermi sea states just beneath the gap [3]. Such currents are conjectured to be persistent and circulating in equilibrium [3], but they become mediative VH currents connecting the two crossbars in Fig. 1 under the application of bias voltage, thereby circumventing the absence of electronic states around the DP while demanding a major overhaul of the LB theory [20] in which the absence of states within the gap is a fundamental reason for R NL ≡ 0 in Fig. 2(c). In geometries with L < W , we (as well as Ref. [18]) do get R NL = 0 in Fig. 2(c), but this is trivially explained by transport through evanescent wave-functions able to propagate from first to second crossbar through the gap in such geometries [12,18].
Another long-standing puzzle is the metallic-like ordinary longitudinal resistivity, ρ xx ∼ 10 kΩ, measured experimentally [1] for the G/hBN channel between the two crossbars despite the expected finite gap in its bulk. This suggests the presence of additional conduction pathways, such as edge currents observed very recently [24], which shunt the highly insulating state at low temperatures. However, previous studies have concluded that edge states are either absent [3] or, when they are present due to edges like zigzag [25,26] or chiral [27,28], they become gapped near the DP and dispersionless away from it so they do not carry any current [29]. The latter conclusion is reproduced in Fig. S1(f) in the Supplemental Material (SM) [30] where we compute band structure of G/hBN wire with zigzag edges using the same simplistic Hamiltonian employed in prior studies [1-4, 15, 18, 29].
In this Letter, we resolve both puzzles-R NL ≡ 0 [ Fig. 2(c)] obtained for multiterminal G/hBN whose bulk counterpart exhibits σ v xy = 0 [ Fig. 4(a)]; and metallic-like ρ xx despite presumed gap opening around the DP-by performing density functional theory (DFT) calculation combined with quantum transport simulations, based on both the multiterminal LB formula (Figs. 2 and 3) and the Kubo formula (Fig. 4). We demonstrate that these puzzles are an artifact of the simplistic Hamiltonian [1-4, 15, 18, 29] which inadequately describes G/hBN wires. For example, in contrast to previously obtained [29] gapped band structure for all types of G/hBN wires, the ab initio band structure [ Fig. 3(a)] of G/hBN wires with zigzag edges has no gap. The ab initio Hamiltonian combined with the LB formula yields R NL [ Fig. 2(a)] and ρ xx [ Fig. 2(b)] whose features closely match those measured in Ref. [1]: (i) ρ xx peak is wider than R NL peak whose height decays exponentially with increasing L (Fig. S6(b) in the SM [30]); (ii) aligned hBN is necessary to obtain R NL = 0 in Fig. 2(a)-isolated zigzag graphene wire also has edge states ( Fig. S3(a)-(f) in the SM [30]), but its band structure leads to R NL ≡ 0 (Fig. S4 in the SM [30]).
Since the resolution of the two puzzles relies on an accurate Hamiltonian for G/hBN wires, as well as its combination with a proper quantum theory of observable transport quantities, we first summarize inconsistencies arising in previous theoretical analyses. The seminal arguments [2] for the VHE are based on semiclassical transport theory describing the motion of a narrow wavepacket constructed by superposing [2,44] eigenstates of gapped Dirac HamiltonianĤ D with dispersion k . This assumes that each valley of G/hBN can be described bŷ where v F is the Fermi velocity, σ = (σ x ,σ y ) is the vector of the Pauli matrices corresponding to the sublattice degree of freedom and k is the momentum operator. The wave-packet velocity [14], v k = 1 ∂ k /∂k + dk/dt × Ω k , acquires an anomalous term due to the Berry curvature Ω k hot spot near the apex of the valley described byĤ D . Since Ω k in valley K has opposite direction to that in valley K , electrons belonging to two valleys will be separated [44] in the opposite transverse directions in the presence of an applied electric field E which is required to accelerate electrons according to dk/dt = eE. This gives rise [2] to σ K,K xy = e 2 πh d 2 k f (ε k )Ω k , where f (ε k ) is the Fermi function forcing the integration over the whole Fermi sea, i.e., from the bottom of the band to the Fermi level E F . However, it has already been pointed out in Ref. [18] that nonzero E cannot appear in the linear-response limit of the multiterminal LB formula [19][20][21]. The experiments measuring R NL are carefully kept [1] in the linear-response regime in order to avoid heating of the device and the ensuing thermoelectric effects that can add large spurious contributions to R NL [7]. The multiterminal LB formula, I p = (2e 2 /h) q G pq (V p − V q ), relates the total charge current I p in lead p to voltages V q in all other leads via the conductance coefficients G pq = dE (−∂f /∂E)T pq (E), where the derivative of the Fermi function confines the integration to a shell of states of width ∼ k B T around E F [20]. The transmission functions T pq (E) do not include any effect of E [19][20][21]. We use the multiterminal LB formula implemented in KWANT package [21,45] to compute V NL and R NL in response to an injected current I 3 = −I 4 [12] while keeping I p ≡ 0 in the other four leads. The same procedure allows us to compute ρ xx = R 4T W/L from the four-terminal resistance R 4T = (V 3 − V 4 )/I 1 obtained by injecting current I 1 = −I 2 into the device in Fig. 1 and by imposing voltage probe condition, I p ≡ 0, in leads p = 3-6.
The semiclassical arguments for the origin of σ v xy can be replaced by quantum transport theory based on the Kubo formula [15,46], which requires to first obtain [47,48] the velocity operatorv. The physical consequences of the equation of motion forv = v Fσ defined byĤ D [15], dv/dt are extracted by finding its expectation value in a suitably prepared wave-packet [44,49] injected with initial velocity into G/hBN where it can propagate even in the absence of any external electric field [49]. The first term of dv/dt causes Zitterbewegung motion of the center of wave-packet, while the second one acts on it like a Lorentz force due to an effective magnetic field in the direction of the unit vector e z perpendicular to the graphene plane. For electrons from the K valley, v y = −v Fσy , leading to opposite direction of such force. The gapped Dirac HamiltonianĤ D is a long wavelength limit of a more general tight-binding Hamiltonian (TBH) which accounts for both valleys and can, therefore, be used to capture intervalley scattering effects. The TBH, which is preferred for numerical calculations using the LB [18] or the Kubo formula (as well as for simulations of wave-packet dynamics [44]), is defined on a honeycomb lattice hosting a single p z orbital per sitê Hereĉ † i (ĉ i ) creates (annihilates) an electron at site i; hopping t 1 = 2.7 eV is nonzero between nearest-neighbor (NN) carbon atoms; and the on-site energy ε i = ±∆, responsible for the gap E g = 2∆, is positive on one atom of the graphene unit cell and negative on the other one to take into account the staggered potential induced by the hBN substrate (while neglecting any reorganization of chemical bonding or change in the atomic order of graphene and hBN layers). We use ∆ = 29 meV [13] to compute R NL in Fig. 2(c) for G/hBN channel with zigzag edges, as well as to compute its band structure (Fig. S1(f) in the SM [30]) exhibiting gapped flat bands [29]. Figure 4(a) shows σ v xy , computed using the Kubo formula [46,48] combined with a valley-projection scheme [10,30,50], for a rhomboid supercell of G/hBN with periodic boundary conditions described byĤ TB in Eq. (2). The supercell is either clean or it contains longor short-range disorder (Sec. IV in the SM [30]) as additional terms in the on-site energy ε i . In the clean limit, we confirm [15] that σ v xy = 2e 2 /h is quantized inside the gap [ Fig. 4(a)], as well as that the Fermi sea states just beneath the gap provide the main contribution to it [3]. For long-range disorder that does not mix valleys, σ v xy remains close to the clean limit within a smaller energy range than E g [ Fig. 4(a)] due to disorder-induced broadening of the states [15]. High concentration of valley mixing short-range disorder slightly reduces σ v xy [ Fig. 4(a)]. In order to replaceĤ TB with a more accurate Hamiltonian, we proceed by computing the ab initio band structure of G/hBN wires using local-orbital pseudopotential DFT, as implemented in ATK [51] and OpenMX [52,53] packages, whose Kohn-Sham Hamiltonian can be also easily combined with the LB formula [54]. We assume stacking where one C atom is over a B atom and the other C atom in the unit cell is centered above the hBN hexagon, as energetically most stable configuration found in DFT calculations [13]. The DFT Hamiltonian for a wire-composed of C, B and N atoms, as well as H atoms passivating dangling bonds along the zigzag edges-produces the gapless band structure in Fig. 3(a) and the corresponding zero-temperature twoterminal conductance G 21 = 5 × 2e 2 /h [ Fig. 3(b)] near the DP (at E − E F = 0). This value is insensitive to the bulk nanopore [ Fig. 3(d)], signifying edge transport, but it is reduced to G 21 2e 2 /h in the presence of edge vacancies [ Fig. 3(e)] or quantum interferences generated by edge vacancies in series [ Fig. 3(f)] due to lack of topological protection against backscattering [29,56]. The valley-polarized states [55] above and below the DP are bulk states since G 21 is reduced at those energies by the presence of the bulk nanopore [ Fig. 3(d)]. The edge state are visualized by plotting the local density of states (LDOS) in Figs. 3(c)-(f), which is peaked near the edges but it remains nonzero in the bulk [28]. This can be contrasted with topologically protected edge states in quantum (ordinary [10], spin [57] and anomalous [58]) Hall insulators where LDOS in the bulk vanishes. The corresponding local current distributions (Fig. S5 in the SM [30]) shows that edge currents can survive edge disorder breaking G/hBN wire into short zigzag-edge segments, as suggested also by experiments [24].
Since the usage of the full DFT Hamiltonian is prohibitively expensive for large number of atoms (∼ 10 6 in our LB or Kubo formula calculations), we derive simpler ab initio TBHs [59]. A widely-used approach for this purpose is to transform the DFT Hamiltonian to a basis of maximally localized Wannier functions in a selected energy window around E F [60]. We combine Wannier TBH (truncated [59] to third NN hoppings without loss of accuracy in the energy window considered in Fig. 4) with the Kubo formula in Fig. 4(b) where we find quantized σ v xy in the gap, as well as its surprising resilience to short-range disorder that was not exhibited in Fig. 4(a) for the simplistic TBH in Eq. (2). However, the Wannier TBH applied to G/hBN wires generates much larger group velocity ∂ε kx / ∂k x of edge state bands near the DP (Fig. S1(d) in the SM [30]) than in Fig. 3(a) due to lack of information about atoms (like H) passivating bonds of edge carbon atoms. Nevertheless, such Wannier TBH could be useful for fully encapsulated G/hBN wires where edges are not exposed to the environment [1,24]. Therefore, we derive an alternative ab initio TBH that precisely fits (Fig. S1(c) in the SM [30]) the bands around the DP in Fig. 3(a). This requires up to the fifth NN hoppings and on-site energies in Eq. (2), with adjusted values along the edges (Fig. S2 in the SM [30]). The ab initio 5NN TBH combined with the multiterminal LB formula leads to R NL = 0 exhibiting sharp peak [ Fig. 2(a)] near the DP in the case of ballistic G/hBN channel. The introduction of disorder-edge ( Fig. S5(b) in the SM [30]) or bulk of different range (in the model of Sec. IV in the SM [30])-reduces the height of R NL peak [ Fig. 2(a)], while concurrently generating peak of ρ xx [ Fig. 2(b)]. Although metallic-like ρ xx could arise due to trivial rea-sons, such as charge inhomogeneity induced by chemical or electrostatic doping, recent imaging [24] of proximityinduced supercurrents flowing near the edges of G/hBN wires is in accord with LDOS [ Fig. 3(c)-(f)] or local current distributions (Fig. S5 in the SM [30]).
In conclusion, R NL peak exists in the clean limit [where ρ xx → 0] and is reduced by the disorder, whereas σ v xy is insensitive to disorder at the DP [ Fig. 4(b)]. Furthermore, plotting R NL (E F ) vs. σ v xy (E F ) in Fig. 2(d), computed for a range of E F ≥ 0 independently while using the same long-range or short-range disorder in both the LB and Kubo formula calculations, reveals that these two quantities are not related in a way conjectured by Eq. (1). This strengthens our principal conclusion-the observed R NL is driven by particular gapless band structure and its edge eigenstates carrying current near the DP, rather than by the Fermi sea dissipationless topological valley currents conjectured for the gapped Dirac spectra.