Generation and control of non-local chiral currents in graphene superlattices by orbital Hall effect

Graphene-based superlattices offer a new materials playground to exploit and control a higher number of electronic degrees of freedom, such as charge, spin, or valley for disruptive technologies. Recently, orbital effects, emerging in multivalley band structure lacking inversion symmetry, have been discussed as possible mechanisms for developing orbitronics. Here, we report non-local transport measurements in small gap hBN/graphene/hBN moir\'e superlattices which reveal very strong magnetic field-induced chiral response which is stable up to room temperature. The measured sign dependence of the non-local signal with respect to the magnetic field orientation clearly indicates the manifestation of emerging orbital magnetic moments. The interpretation of experimental data is well supported by numerical simulations, and the reported phenomenon stands as a formidable way of in-situ manipulation of the transverse flow of orbital information, that could enable the design of orbitronic devices.


I. INTRODUCTION
The electronic properties of graphene and other twodimensional (2D) materials with a honeycomb lattice are dictated by the low-energy physics at two inequivalent K and K valleys of the reciprocal space [1].The large separation between these valleys allows to distinguish valley quantum numbers that, likewise the spin degree of freedom, can be used to store and process information [2].Moreover, the valleys in graphene possess opposite orbital magnetic moments, that at K and K are proportional to the inverse of the band gap.Systems with small band gaps can have extremely large magnetic moments.This results in giant Zeeman splittings upon interaction with weak external magnetic fields, lifting valley degeneracy [3,4].
In gapped graphene, the application of an electric field has been predicted to induce a flow of electrons, moving in opposite directions for different valleys, and thus giving rise to a valley Hall effect (VHE) [5,6], that could be detected by non-local transport measurements [5][6][7][8][9][10].However, this is also accompanied by a transverse flow of valley magnetic moments.Consequently, the VHE can also be depicted as an orbital Hall effect [11].The orbital magnetic moments, which are physical quantities defined in the entire momentum space, replace the valley quantum numbers, that depend on the existence of welldefined pockets [11].Different from the VHE, this interpretation leads to a transverse current of a physical observable and can be used to analyze the interaction arXiv:2206.04565v1[cond-mat.mes-hall]9 Jun 2022 with magnetic and electric fields in the same framework.
Furthermore, electronic structure of graphene-based van der Waals (vdW2Ds) heterostructures can be tailored in a remarkable way by varying the twist angle between weakly interacting atomic layers and generating graphene moiré superlattices [5-7, 9, 10, 12-14].In particular, single or doubly aligned graphene/ hexagonal boron nitride (hBN) heterojunctions are very interesting for the study of inversion symmetry breaking in graphene.Indeed, such systems present considerable non-locality [7,9,[15][16][17], whose origin, although frequently associated to the VHE, is currently strongly debated [18][19][20][21].It has been shown that doubly aligned hBN/graphene stacks might generate a super-moiré pattern [22], leading to the presence of small and non-uniform band gaps [20,21].vdW2Ds are, henceforth, a perfect platform for the study of inversion symmetry breaking in graphene.Because of the small gaps, the valley orbital magnetic moments are large and can be manipulated with magnetic and electric fields.
In this work, we report unambiguous formation of chiral non-local currents in doubly aligned hBN/graphene/hBN heterostructure, presenting direct evidence of their orbital magnetic origin.Our alignment design between layers minimizes the band gap, resulting in giant orbital magnetic moments.The interaction with weak magnetic fields lifts their degeneracy, generating chiral non-local currents.Quantum transport simulations in graphene nanoribbons with dispersive edge states and linear response theory calculations support our interpretation of the experimental findings regarding the origin of the chiral non-local resistance in graphene.

Figure 1a
displays an optical image of the heterostructure consisting of a graphite back gate (15 nm thick) and monolayer graphene encapsulated between a top and bottom layer of hBN with a thickness of 10 and 50 nm respectively.The crystals were aligned following their exfoliated straight edges using a micromechanical rotator and, for that reason, the relative twisting angles along the vertical heterostructure are expected to obey m × 30 • , being m = 0, ±1, ±2 . . . .From our electrical data displayed in Fig. 10e of the Supplementary Information we can reasonably discard twisting angles between graphene and hBN that are 0 • or a multiple of 60 • .On the one side, charge neutrality point (CNP) appears as a standalone main Dirac peak with no traces of secondary satellite peaks in the measured carrier density range |n| < 3 × 10 12 cm −2 that arise from the existence of electron-hole pockets at both sides of the main peak at such twisting angles [9,10,12,14,15,17,23,24].On the other side, CNP resistivity of heterostructures where graphene is aligned to 0 • or 60 • with the hBN exhibits a strong thermally activated behaviour with values exceeding hundreds of kΩ at low temperature, indicating a moiré coupling-induced band gap of the order of 30 meV [10,14,15,23,25].In our case, we observe a thermally activated behaviour at low temperature but a CNP resistivity of only ∼ 7 kΩ at room temperature.These characteristics are consistent with band gaps smaller than 10 meV.Charge mobility extracted from magnetotransport Hall measurements rises to 200.000 cm 2 /Vs, as shown in Fig. 10 in the Supplementary Information.
From a careful Raman analysis, we underpin the relative orientation of the flakes from the evolution of full-width at half-maximum of the 2D peak (FWHM2D) as a function of the twisting angle.Finney et al. [24] kept the bottom hBN aligned with graphene at 0 • and varied the relative angle of the top layer from 0 • to 60 • .They showed that the vertical structure exhibits a noticeably high value of the FWHM2D, exceeding the standard one found on isolated graphene by 20 (40) cm −1 if one (both) hBN layer(s) is (are) aligned to the graphene at the commensurate angles of 0 • or 60 • [24].Such broadening results from the moiré-scale relaxation of the graphene lattice, which strongly modifies the band structure [26].We found a value of ∼ 20 cm −1 for the FWHM2D (see Fig. 9 of the Supplementary Information), in good agreement with the value for a standalone graphene flake that would correspond to a twisting angle that is neither 0 • nor 60 • in Ref. [24].The ratio between the 2D and G peaks I(2D)/I(G) > 9 also strengthens the assumption of twisting angles of ±30 • .
Figure 1c shows the non-local resistances R nl for two different configurations (1 and 2) of the injection i and collection c terminals at 1.5 K.This allows us to infer the decay of the non-local signal as a function of the distance ∆x i between the injection (local) and collection (non-local) terminals.For ∆x 1 = 2.5 µm, R nl ∼ 1600 Ω which is consistent with other non-local measurements in graphene/hBN heterostructures [7,10,15].Moreover, the non-local signal gets weaker for increasing distances, reaching R nl ∼ 470 Ω for ∆x 2 = 5 µm.In the absence of external magnetic field, the position of the non-local peaks are centered around the CNP and are symmetric in respect to electron/hole regions.The relation R nl (∆x) = V nl /I 0 πρ xx e −π|∆x|/W , where W = 1.5 µm is the bar width, displays an exponential decay of the Ohmic contribution to R nl as a function of distance between the driving current and the non-local pair of contacts.This relation, already seen in graphene-based devices [15][16][17], serve us to rule out the Ohmic contribution to R nl as its prevailing mechanism.We extracted the ratio between the measured non-local resistances at different distances, obtaining R nl (∆x 2 )/R nl (∆x 1 ) Measured = 0.29 while the purely-Ohmic expression gives us R nl (∆x 2 )/R nl (∆x 1 ) Ohmic ∼ 0.005.From this analysis, one can see that the non-local signal is orders of magnitude higher than the expected Ohmic contribution (see Fig. 12 in the Supplementary Information).

III. CHIRAL NON-LOCAL SIGNAL
To explore the relation between the non-local currents and orbital magnetic moments, we apply positive and negative perpendicular magnetic fields and analyze the non-local signals for different injectioncollection configurations.Unless stated otherwise, all measurements were performed at a fixed electronic temperature of 1.5 K and with an excitation current set to 10 − 20 nA.This low current amplitude was chosen to minimize thermal contributions to the nonlocal transport due to Joule heating and Ettingshausen effects [19] whilst simultaneously maximizing the signalto-noise ratio of the measured voltages.A summary of our results is shown in Fig. 2, which contains 9 different panels divided in three columns.Each column presents a different configuration of the external magnetic field, i.e., −0.5 T, 0 T and 0.5 T. Each of the three rows presents a different injection-collection setups, sketched in the diagram on the left.Each pair of contacts on the diagram has a specific color, while the arrows show the direction of the current for the pair of contacts of a particular injection-collection configuration.Each panel displays the local and non-local resistances as a function of the voltage applied to the graphite backgate.Their color palette matches the ones of the corresponding contacts, irrespective of the specific injection-collection setup.
Let us first comment on the effect of the magnetic field on the local resistances.The CNP is located at V g ∼ 0.6 V, as extracted from the Lorentzian fit of the local resistance, and does not present a noticeable shift when the magnetic field is applied.Still, there is a sharp increase of the local resistance with the magnetic field, which is very pronounced at the CNP.Panel b1 in Fig. 2 displays the evolution of the non-local signal as a function of the distance without external magnetic field.Panel b2 represents a configuration of the nonlocal pairs of contacts that are placed symmetrically to the current flow but at opposite directions.In this case, the homogeneity of the sample is demonstrated because the non-local signals at both sides of the current flow have an expected matching value, as the magnitude of the non-local signal decays with the absolute value of the distance to the injection current.This characteristic, discussed previously, can also be seen in panels a2 and c2 of Fig. 2.
A striking behavior of the non-local signal arises in presence of an external magnetic.We first focus on the case where the current is injected between two non-local contacts.Figure 2.a2 shows the non-local resistances for the B = −0.5 T, where one can see a clear separation between the peaks of opposite contacts.Moreover, they are mostly located either in the electron or hole sectors.Surprisingly, the position of the two peaks is swapped upon magnetic field reversal, which is a clear indication of a chiral behavior of the electronic response.If the two collectors are located at the same side of the injector, as shown in the first and last row of Fig. 2, the situation is different.In these two cases, the two non-local resistance peaks are centered and located either at the electron or hole sector and switch positions with the sign of the magnetic field and the relative orientation with respect to the collector.While both peaks appear at the hole sector in Fig. 2.a1, when the sign of B is reversed, they appear at the electron sector (see Fig. 2.c1).If instead, we switch the position for positive B, as in Fig. 2.a3, the peaks also appear at the electron sector, changing to the hole sector if the field is reversed (see Fig. 2.c3).Furthermore, it is important to mention, that in resemblance to similar experiments, the non-locality is strongly enhanced with the magnetic field in all configurations [12,15,17].
The nine panels in Fig. 2 demonstrate a fully chiral behaviour of the non-local signal at low magnetic fields, which has not been observed in similar heterostructures in monolayer graphene with twisting angle set to different commensurate angles [10,15,17].
To clarify the underlying mechanism, we use the modern theory of magnetism and numerical simulations based on the linear response theory and the Landauer-Büttiker formalism implemented in the KWANT toolkit [27].Addressing the orbital magnetic moment in solids is a nontrivial topic due to the ill-defined behaviour of the r operator in the Bloch basis [28][29][30].Nonetheless, modern theory of magnetism provides an appropriate description of this phenomenon by treating the Bloch electron as a self-rotating wave packet whose magnetic moment is expressed purely in bulk quantities as where n is the band index and U is the periodic part of the Bloch eigenstate.Following Ref. [31], one can show that applying this approach to gapped Dirac materials leads to an expression for the orbital magnetic moment that reads F is the effective mass at the Dirac point (DP), ∆ is the system gap and v F is the Fermi velocity.In this case, the intensity of the orbital magnetic moment is inversely proportional to the system gap while it can generally be seen as inversely proportional to the effective mass as well.This orbital moment in the presence of weak magnetic fields behaves similarly to the spins, i.e. it can couple directly to them.This gives rise to a kdependent Zeeman effect that, in first order perturbation theory, renormalizes the energy spectrum close to the Dirac points [31,32], as shown in Fig. 3a.It can also be interpreted in the usual way as a Zeeman effect with giant g-factor, which is proportional to the inverse of the effective mass, a well known effect in semiconductors.The valleys have opposite magnetic moments, producing a relative shift between the valleys.This effect can be very strong for small gap systems, leading, for example, to situations where the Fermi energy lies inside the gap for one valley while it is located in the electron/hole sector for the other.
Given the proportionality with the inverse of the gap size and coupling with external fields, one can argue that the chiral behavior observed in the non-local resistance measurements appears as a manifestation of the orbital Zeeman effect allowed by the formation of small a gap in the doubly-encapsulated hBN/graphene/hBN heterostructure.Considering the OHE formalism, one should expect a chiral behavior as observed in Fig. 2 (see the Supplementary Information for details).The magnetic moments of each valley flow in opposite directions.Because of the sign of the Berry curvature, the flows invert direction when switching from hole to electron sectors.At the same time, the Zeeman shift between the valleys leads to a relative shift between the non-local peaks that should invert with the change in the magnetic field orientation and collector's locations.Using this reasoning, we can estimate the band gap from non-local peaks of Fig. 2.a2, obtaining ∆ ∼ 5 − 8 meV.
To inquire about this, we calculated the non-local resistance in an eight terminal device containing a graphene nanoribbon with the geometry shown in the inset of panel c of Fig. 3 in the presence of external magnetic fields.The advantage of this approach is that it does not rely on any assumption regarding the orbital Zeeman effect, as the magnetic field is taken into account by the well established Peierls' substitution.Moreover, it can directly connect with experiments, as it allows the calculation of the non-local resistances.Finally, it considers both bulk and edge states.Still, here we are not particularly interested in the controversy related to the specific channels conveying the non-local current because, as it will become clear below, the orbital Zeeman effect is present in bulk and edge states and either of them can, in principle, produce chiral non-local currents.
Previous numerical simulations showed the need of dispersive edge states near the Dirac point for nonlocal transport in gapped graphene [18].They are absent in theories based on the simplistic Hamiltonian considering a single p z orbital.Instead, we used a 6bands tight-binding model that takes into account also the d orbitals [33] and a staggered sublattice potential to break the inversion symmetry of the system.It is important to note that other effects such as non-uniform potential and coupling to hBN layers can also lead to dispersive edge states and similar non-local currents (see Sec. 7 of the Supplementary Information for more details).
Panels b and c of Fig. 3 show the energy bands for a zig-zag graphene nanoribbon with 13 nm of width and a sublattice staggered potential of 5 meV.From panel c, it is noticeable that the energy bands from this multi-orbital model exhibit a dispersive behavior similar to the ones observed from ab-initio calculations of Ref. [18] without a gap in the whole spectrum, but have valence-conduction band separation due to the inversion symmetry breaking and very well-defined electron pockets at opposite sides in the nanoribbon Brillouin zone.Panel b depicts the comparison of the energy bands of the nanoribbon in the cases without considering magnetic fields and with out-of-plane magnetic fields via Peierls' substitution.The most relevant feature displayed by this panel is the stark valley-contrasting coupling with the magnetic field, which is similar to the one observed for the bulk bands.Although the modern theory of magnetism is welldeveloped only for bulk systems, it is clear, from the results in this panel, that the behavior displayed by the energy bands of the nanoribbon is in qualitative agreement with this theory.
Aiming to reproduce the results from the middle columns in Fig. 2, in panels d and e we used the same injection-collection scheme for our non-local resistance simulations.Panel d portrays the case in which the injection occurs at one side of the device.As in the measurements, the simulations results show considerable decay of the non-local signal with the channel length.On the other hand, panel e shows the case in which the current injection occurs at the middle terminals of the device and the non-local resistance is computed in the leads of its sides.Comparing the non-local signal in the lateral terminals, it is clear that the non-local signal at both sides of the sample is overall symmetric.However, the most striking behavior appears when opposite magnetic fields are taken into account, as in the cases displayed in panels f and g.Comparing these figures with the experimental measurements shown in Fig. 2, we find convincing agreement.Although, due to the small sizes of the simulation as compared to the experiment, we need to apply much higher magnetic fields to obtain comparable results.Nevertheless, the most outstanding feature of these results becomes evident when comparing them with the renormalized energy bands in panel b, where the chirality and energy selectivity is directly related to the coupling between the orbital magnetic moment and the external magnetic field, suggesting that the mechanism at play in the generation of these non-local signals observed in the experiment is the orbital valley Hall effect [11].Still, our analysis suggests that Fermi surface edge currents carry the nonlocal signal, once the absence of dispersive edge states destroys the non-local signal.
Figures 4a and 4 b display the contour plot of R nl as a function of both B and V g − V DP for two symmetrical configurations for the local and non-local contacts as sketched below each panel.The electronic temperature at which the curves were recorded was T = 250 mK and the voltage has been centered at the Dirac peak for B = 0.These two panels show a clear chiral behavior and an apparent valley-carrier locking in the non-local signal for low magnetic fields ranging from −0.5 to 0.5 T .In Fig. 4a we can observe a distinct transition from a hole-mediated non-local transport for negative magnetic fields towards an electron-like one when the magnetic field is reversed.Moreover, a strong asymmetry in the non-local curves is clearly visible, with a sudden decay while approaching the DP from the dominant carrier species towards the prohibited one.Figure 4b displays the mirrored configuration for the pairs of contacts, compared to Fig. 4a as sketched in panel Fig. 4c.Consequently, we can argue that the valley-carrier locking, visible in the non-local signal, is a quite robust phenomenon and implies a flow of carriers characterized by different orbital magnetic moments.As the non-local signal should originate from the flow of orbital magnetic moments from a single valley, we compare R nl with the valley filtered OHE, calculated according to Ref. [11], shown in the insets of Figs.4a and 4b.In conclusion, the interpretation of the experimental results are very consistent with the changes in the OHE resulting from the orbital Zeeman effect.The theoretical calculations indicate that the observed Zeeman shift is consistent with a graphene gap of ∆ ∼ 5 − 10 meV.For details on the OHE calculations, see Sec. 6 of the Supplementary Information.
Figure 5 presents the local and non-local signals for a configuration where the excitation current lies symmetrically between two different pairs of contacts for R nl at a fixed perpendicular magnetic field of B ⊥ = 0.5 T and variable in-plane component (B ) measured at 1.5 K. B ranges from 0 to 12 T and its evolution has been marked with an arrow as a guide to the eye.When no parallel magnetic field is applied (B = 0), the charge carrier type is effectively coupled to one of the valleys and the asymmetric chiral behavior is found.Distinctively, while enhancing the in-plane magnetic field component, both non-local and local signals diminish and become symmetric.From these data, it becomes clear that the chiral non-local currents cannot be attributed to a spindependent effect.In-plane magnetic fields can be used to probe spin-polarized currents, as they lead to a nonlocal resistance that oscillates in function of the field.At the same time, they can also be used as a tuning parameter to manipulate the orbital characteristics of 2D multilayers [34,35].They affect the quasi-momentum of each layer differently [36][37][38][39], modifying the effective coupling between the layers and thus the resulting bandstructure.Therefore, similar to strain, they can modify substantially the physics of twisted bilayers and moiré superlattices [34,35].
Still, the application of an in-plane magnetic field B introduces a layer-dependent gauge field A l = B × z l that modifies the electron momentum p → p + (e/c)A l , where l indexes the layer.For graphene encapsulated by two hBN layers and located at z = 0, the magnetic field shifts the momenta of electrons in each hBN layer along opposite directions.Time-reversal symmetry is broken, shifting the momenta of the electrons of different valleys in the same direction.This can have an important impact in the resulting band-structure, modifying the orbital magnetic moments and henceforth modifying local and non-local signals.
In conclusion, we have presented non-local transport measurements on hBN/graphene/hBN narrow gap heterostructures at low magnetic fields, which clearly indicate the presence of chiral effects.Such chiral response is inferred from the non-local resistance when reversing both the magnetic field and the injectioncollection configurations.The interaction between large orbital magnetic moments arising in small gap graphene based superlattices and external magnetic field produces a relative Zeeman shift between the two valleys in both bulk and edge electron states.Furthermore, based on our experimental and theoretical analysis, it is clear that, regardless of the details about the location of current flow, the manifestation of strong chiral effect originates from the interplay between the Zeeman shift and the transverse flow of orbital magnetic moments.Importantly, the analysis of the non-local transport as a function of the magnetic field direction rules out spin effects, whereas its dependence with the distance between contacts clarifies that Ohmic and thermal contributions bring marginal contributions to the main signal.Finally, our computational results show that the valley-orbital Hall effect displays fingerprints in both bulk and edge transport, being of topological nature or not.
All these facts support the interpretation that the origin of the giant non-locality in the studied graphene superlattices is linked to the orbital Hall effect resulting from the valley magnetic moments.Beyond shining a new light on a fierce debate concerning the formation of topological versus non-topological valleydriven phenomena to explain previously reported nonlocal signals [40,41], our findings pave the way towards future developments in room-temperature graphene orbitronics.Methods: The device fabrication of the superlattices follows the standard dry transfer technique with a polycarbonate film fabricated and deposited onto a polydimethylsiloxane stamp.
The relative rotation between the different layers, following their natural edges, was controlled using a heated stage with a micromenchanical rotator with an accuracy better than 0.5 • .The heterostructure rested atop a commercial Si/SiO 2 substrate.The fabricated stack was patterned using electron beam lithography followed by a dryetching process in an ICP-RIE to define the sample geometry.The sample was patterned into a form of a multiterminal Hall bar (optical image in Fig. 1) with the central horizontal bar of width W = 1.5 µm, total length of ∼ 11 µm and a distance between the centres of the contacts of 2.5 µm.Electrical contact to all devices was made by Cr/Au (10nm/50nm) deposited by electron beam evaporation.Supplementary Sections 1 Figure 6 and 7 provide full details of the device fabrication.We extracted the all Raman spectra and their associated FWHM2D from a Lorentzian fit (see Supplementary Section 2 Figures 8 and 9).Transport measurements in the multiterminal device were conducted in two-and four-terminal geometries with a.c.current excitation of 10-20 nA using the standard lock-in technique at 17.7Hz.The graphite layer was gated by applying a direct bias to it in the range of ±12 V.Each flake and the final heterostructure were characterized prior to their stacking by Raman spectroscopy using a micro-Raman spectrometer LabRAM HR Evolution at a wavelength 532 nm and an incident power of ∼10 mW. Figure 8 presents a detailed analysis of the Raman spectra found in the final heterostructure.Panel a) displays an optical image of the final stack prior to the electron beam lithography with a complete 2D Raman map in panel b) obtained at XXX cm −1 ).Panels c) to e) show individual Raman spectra taken at different positions of the stack where only encapsulated graphene was present c, graphite back-gate is found d and the whole final structure where the Hall bar was designed e.
From a careful Raman analysis we can study the relative orientation of the flakes conforming the stack as shown in Ref. [24] where the evolution of the value of the full-width at half-maximum of the 2D peak (FWHM2D) as a function of the twisting angle is presented.In that work, Finney et al. maintained the bottom hBN aligned to the graphene at 0 • and varied the relative angle of the top one ranging from 0 • and 60 • showing that the vertical structure exhibits a noticeably broadened value of the FWHM2D exceeding the standard one found on isolated graphene by 20 (40) cm −1 if one (both) hBN layer(s) is (are) aligned to the graphene at the commensurate angles of 0 • or 60 • .Such broadening result from the moiré-scale relaxation of the graphene lattice which will strongly modify the graphene band structure [26].In our case, a value of ∼ 20 cm −1 for the FWHM2D is found (Figure .S5 Figure.9), in good agreement with the value obtained for a standalone graphene flake that would correspond to a twisting angle that is neither 0 • nor 60 • in Ref. [24].FIG. 9. SUPPORTING INFO 4. a, Raman spectrum of the encapsulated single layer graphene flake.b, 2D peak extracted data from the highlighted region in (a,) and its Lorentzian fit.A value for the full-width at half-maximum of the 2D peak of FWHM2D∼ 20 cm −1 is found which will rule out twisting angles between the layers in the vicinity of θTOP, θBOTTOM ∼ 0 • , ∼ ±60 • as shown in Ref. [24] 3. Magnetotransport at low temperatures.Quantum Hall effect and mobility.
Longitudinal (R XX ) and Hall (R XY ) resistances as a function of the external magnetic field are shown in figure 10.Every spectrum was measured at 1.5 K at different values of the back gate voltage, namely 0 V, 0.3 V and 1.5 V from panel a to c. Values for the mobility have been extracted from the well-known formula from magnetotransport measurements where µ = 1/(R • q • n) where n is the carrier density, q the charge of the carrier and R = R XX(0) W L at a every given value of the backgate voltage.Results for the mobility for both holes and electrons consistently lie in the range of µ ∼ 200, 000cm 2 /Vs, comparable to values obtained in similar vertical heterostructures [10,14,15,17].Figure 11 displays the longitudinal resistance and the filling factor as a function of V g measured at 1.5 K and 12 T. On top of the standard evolution of the Landau levels for single layer graphene ν = ±2, ±6, . . .we can clearly observed unusual plateaus of conductance where both spin and valley degree of freedom have been lifted such as ν = 0, ν = ±1, and eventually fractional quantum Hall plateaus.No traces of the interference with Landau levels arising from secondary Dirac cones have been found thorough the whole study.The Ohmic contribution to the non-local signal as a function of the distance (∆x) to the driving current can be extracted from the relation R nl = V nl /I 0 ∼ (4π)ρ xx e −|∆x|/λ as shown in Refs.[14,15,17] being λ = W/π, W the bar width and ρ xx the longitudinal conductivity.This relation serves us to extract the expected ratio between the non-local resistances at different distances from a pure Ohmic-bulk contribution.The enormous discrepancy between measured and expected ratios in the non-local signals   where L and W are the channel length and width.In zero magnetic field and for a maximal value of R xx = 18 kΩ and being L/W= 2.5/1.5 we would obtain R nl,Ω = 73 Ω which is two orders of magnitude smaller than the measured R nl .Moreover, in absence of external magnetic field only Joule heating effect contributes to the second harmonic of non-local signal R 2f .We were unable to measure significant R 2f nl at zero magnetic field field while for B = 500 mT we obtained a signal smaller than 2% as shown of R nl as shown in Supplementary 15 where R 1f nl and the envelope signal of R 2f nl multiplied for a prefactor ×50 are displayed fro two different pair of non-local contacts.Here, we present a general discussion on perturbation theory to include the effect of magnetic field up to the first-order in the orbital Hall effect calculation.We assume that the intensity of the magnetic field is weak enough and the system is far from the Landau level regime.In this situation, we can treat the effect of the magnetic field in the framework of perturbation theory following Ref.[44][45][46].In graphene systems, this situation is observed in experiments for the magnetic field with intensities | B| 1.0T .We also present the linear response formulas for Hall conductivity (HC) and orbital Hall conductivity (OHC) used to obtain the theoretical contour plots presented in the main text.
Despite the complexity of graphene/h-BN heterostructure used in experiments, near the Dirac point, it is possible to analyze the physics of the system using simply the Hamiltonian of graphene monolayer (ML) with mass [7,[47][48][49].This simplification should not be applied to secondary peaks of non-local resistance measurements, but may be used to understand some features of the central peak.Writing the Hamiltonian of graphene monolayer (ML) on the tight-binding basis, β M L tb = {A, B}, and expanding it near the valleys K = (4π) / 3 √ 3a x and K = − K, we obtain: where, γ ± = v(τ q x ± iq y ) with, q = k − τ K is the wave vector relative to valleys and, τ = ±1 for Dirac cone at valleys K and K respectively.The velocity v = 3at/2 , where t = 2.8eV is the nearest-neighbor hopping amplitude of graphene and a = 1.42 Å is the carbon-carbon distance.∆ is the mass gap term induced by the h-BN substrate.
The mass breaks the spatial-inversion symmetry of the graphene monolayer and opens a bandgap with amplitude E g = ∆ in the electronic spectrum.E We use these unperturbed energies and states to compute the orbital magnetic moment.Considering the application of a weak magnetic field, the energy of electronic bands is corrected by where the second term on the right-hand side of the equation is the orbital magnetic moment m z n,n ( k) of n-th electronic states.The correction in electronic states is In Eq. (A4), N ... means that, it is necessary to normalize the state inside the square brackets after including the perturbative contribution.The equations above give the first-order perturbation theory in the linear order of the magnetic field.In the case of graphene ML, the correction in Eq. (A4) vanishes due to the dimensionality of Hilbert space, i.e., one-dimensional Hilbert space on c(v)-bands.
With the corrected energies and states of Eqs.(A3, A4), we compute the HC and OHC.In the low-temperature limit, the HC is The OHC is given by In Eqs.(A6, A8), the velocity operators are defined by, vx(y) ( k) = −1 ∂H( k)/∂k x(y) .In Eq. (A8), the current with OAM polarized in z-direction that flows in the y-direction is ĵLz y e ) is the atomic Bohr magneton defined using the electron rest mass m e .

Technical details of the non-local Resistance Simulations
In what follows, we give a brief overview of the details of the non-local resistance simulations.As mentioned in the main text, Marmolejo-Tejada et al. demonstrated that the simplistic model that considers only p z electrons in graphene does not reproduce the electronic structure of zigzag Gr/hBN nanoribbons [18].These results motivated us to construct a more involved tight-binding model to provide an appropriate description of the energy states of these nanoribbons.Based on previous first-principles works [50], we took advantage of the D 3h symmetry of the π-bands of graphene near the Brillouin zone corners and constructed a model composed of the p z and the d xz , d yz orbitals.Using the Slater-Koster [51] parameters of reference [33], we built a 6-bands tight-binding model that considers only nearest-neighbour hopping integrals.Figure 17 portrays the energy states of nanoribbons of 50 nm breadth for the p z model and the 6-bands tight-binding model.Upon brief inspection, it is clear that the simplistic p z model, though capable of capturing most of the bulk properties of monolayer graphene and other heterostructures, misrepresents the behaviour of the edge states leading to the appearance of an energy gap when the staggered sublattice potential breaks the inversion symmetry.In contrast, the 6-bands tight-binding model allows the hybridization between p z and d xz , d yz orbitals, promoting the appearance of massive gapless edge-states that gives rise to the non-local signals in our simulations in agreement with the results from Marmolejo-Tejada et al. [18] and supports the formation of orbital magnetic moments that can couple with the external magnetic field.Additionally, the 6-bands model reproduces the width-independent insulating behaviour of armchair graphene nanoribbons reported in Ref. [52].The model up to 5 nearest-neighbours exposed in Ref. [18] has been tested as well in an equivalent device, finding non-local resistances of the same order of magnitude of Fig. 3.However, the profile of the non-local resistance does not resemble the experimental anisotropic peak and the interpretation up to long distance neighbours is not as clear as the multi-orbital model.
In quantum transport simulations, we used an 8-terminal geometry depicted in the inset of panel c of figure 3 in which we fixed the width of the channel and the leads to 13 nm and the separation between each pair of leads perpendicular to the channel direction to 157 nm.Besides this, we also doped the leads perpendicular to the ribbon to counter the insulating nature of the armchair Gr nanoribbons and avoid contact resistance in our simulations.We computed the non-local resistance by obtaining the transmission probabilities for all the contacts and constructed the conductance matrix.Then, we determined the voltages in the leads using the equation I = GV , and to ensure the charge conservation within the system fixed the voltage of the lead from 2 to zero.To include the effects of magnetic fields in our simulation, for the scattering region we have redefined the hopping integrals like t ij µν → t ij µν × exp −i e Ri− Rj A • d l , where A = −Byx is the vector potential associated to the field chosen to not depend on the periodic direction of the ribbon.
in the framework of the project Sym-Break.J. M. C. acknowledges support from the MICINN Ramón y Cajal program (Project No. RYC2019-028443-I).M. V. was supported as part of the Center for Novel Pathways to Quantum Coherence in Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences.K. W. and T. T. acknowledge support from JSPS KAKENHI (Grant Numbers 19H05790, 20H00354 and 21H05233).J. A. D.-N.acknowledges funding from by the CENTERA Laboratories in the framework of the International Research Agendas program for the Foundation for Polish Sciences, co-financed by the European Union under the European Regional Development Fund (No. MAB/2018/9).L. M. C. acknowledges funding from the project "Conversiópn de energía en heteroestructuras de van der Waals" funded by the Agencia Estatal de Investigación (AEI) (PID2019-106684GB-I00/ AEI/10.13039/501100011033)Author contributions: : M. A., E. D. and F. D.-A.developed the concept of the experiment.T. T. and K. W. provided hBN crystals.J. S.-S., A. P.-R., V. C., D. V., and J. A. D.-N.performed device fabrication and carried out Raman spectroscopy.J. S.-S., A. P.-R.E. D., and M. A. performed transport measurements.J. S.-S. and A. P.-R.performed experimental analysis.J. S.-S., A. P.-R., E. D. and M. A. interpreted results with help from L. M. C., T. G. R., S. R., J. M. C., Y. B., R. A. M and F. D.-A.Theoretical calculations were performed by L. M. C., T. G. R., T. P. C., M. V. and S. R. The manuscript and the Supplementary Materials were written by T. G. R., M. A., S. R., L. M. C., J. S.-S., A. P.-R.and F. D.-A. with additional contributions from all authors.

FIG. 1 .
FIG.1.a, Optical image and schematics of a complete device, consisting on a heterostructure based on a graphite back gate (15 nm) and monolayer graphene encapsulated between two layers of hBN with thicknesses 10 nm (top) and 50 nm (bottom).b, Art view schematic of the lattice configuration with twist angles (θTOP and θBOTTOM) between graphene and the top/bottom h-BN, assuming θTOP ≈ 30 • .A moiré wavelength of λ ∼ 0.47 nm can be extracted from the relation found in Refs.[10,42], in stark contrast with the usual λ ∼ 14 nm present in fully-aligned samples.c, Non local resistance dependence with the distance to the injection source for two different sets of contacts at T = 1.5 K in absence of external magnetic field.The two pairs of contacts are separated 2.5 and 5.0 µm from the local signal.The solid line serve as the direction of the driving (local) current.Dashed lines represent the position of the maximum of the two non-local resistances as a guide to the eye.

8 FIG. 3 .
FIG. 3. a, Comparison between the energy bands of gapped graphene systems with and without out-of-plane magnetic fields with B = 0.2 T and a gap ∆ = 14 meV.b, Comparison of the energy bands of zigzag nanoribbons in for B = 0 and |B| = 30 T and sublattice staggered potential V ab = 5 meV.c, Energy bands of zigzag gaped graphene nanoribbons with V ab = 5 meV.The inset shows the schematic representation of the device used in the simulations.Non-local resistance for current injection across terminals 0 − 5 (panel d) and 1 − 4 (panel e) in the absence of external magnetic fields.Non-local resistance for current injected across terminals 1 − 4 for B = −30 T (panel f ) and B = 30 T (panel g).

0. 8 0 5 FIG. 4 .
FIG. 4. Panels and a and b show the heat maps of the sample non-local resistance for two symmetrical configurations scanned in the space of backgate voltage, corrected by the position of the Dirac peak at B = 0 (VG − VDP ), and magnetic field measured at 250 mK.Insets on each graph show the numerical results for the OHC using Eq.(A7) in the Supplementary Information, showing extraordinary agreement with the experimental results.c, Sketch for the two different injection and collection configurations for the experimental and theoretical results shown in panels a and b above.

FIG. 8 .
FIG. 8. SUPPORTING INFO 3. a, Optical image of the vertical heterostructure.b, Raman spectra of the whole designated area.c,d,e, Full Raman spectra of the positions marked with the coloured symbols in (b) where encapsulated graphene is present (c), the graphite back gate is found (d) and both graphene and graphite coincide in the underlying structure (e).

FIG. 10 .
FIG. 10.SUPPORTING INFO 5. a, b, c, Spectra for the local longitudinal RXX (in grey) and Hall resistances RXY (black) as a function of the external magnetic field measured at constant gate voltage of 0 V (a), 0.3 V (b) and 1.5 V (c) respectively.Each panel includes the extracted mobility value.d, Sketch of a typical Hall-bar-like configuration of the contacts where the solid line represents the longitudinal glocal driving current.e, Local resistance vs. backgate voltage between three different sets of contacts at B=0T.Charge neutrality point (CNP) appears as a Dirac peak at VG ∼0.6V out trivial Ohmic contributions to R nl .Black dashed lines show the expected evolution for R nl (∆x) following the Ohmic-contribution-relation.Chiral non-local resistances are robust as a function of the electronic temperature as shown in Fig.13where R nl and R local as a function of V g are displayed for different temperatures under external low fields of B = 0.1 T. While the maximum value for the different resistances slightly drops with T , chirality is clearly preserved.Finally, Figure14shows the evolution of the normalized R nl as a function of V g measured at different external magnetic fields at a electronic temperature of 250 mK.Every non-local resistance has been normalized by its local counterpart for clarity and they have been recorded following the sketch in the figure.As we can see from both panels the normalized non-local signals preserve a strong chirality and are enhanced by the application of small to moderate external magnetic fields.As the orbital magnetic moment in graphene interacts with the external magnetic field the interplay between the Zeeman shift and the transverse flow of orbital magnetic moments increases giving rise to a strengthened non-local signal.

FIG. 12 .
FIG. 12. SUPPORTING INFO 7. a,b, Non local resistance peaks for two different configurations depicted on the schematic inset for each graph.The Ohmic contribution to the non-local signal has been calculated following the exponential dependence (black dashed line) with the distance R nl = V nl /I0 ∼ (4π)ρxxe −|∆x|/λ , where we have taken the experimental value of λ = W/π being W = 1.5 µm the bar width and |x| the distance between local and non-local pairs of contacts as found in references [14, 15, 17].

FIG. 15 . 6 .
FIG. 15.SUPPORTING INFO 10.Comparison of the non-local resistance and second harmonic signals for two different configurations as indicated in the bottom sketchs.Data were recorded at 1.5 K and with an external applied magnetic field B = 500 mT.The intensity of the second harmonic has been adjusted in order to verify the match between the signals, which indicates a negligible influence of thermal phenomena to which the second harmonic would be sensitive, including Nernst and Ettinghsausen effects.
energies and eigenstates of these unperturbed Hamiltonian.For ML we have n = v(c) to index the state of valence (conduction) band and, for BL we have n = 1v(c), 2v(c) to index the states of the 2-dimensional subspace of the valence (conduction) band.
FIG. 16.SUPPORTING INFO 11.Orbital Hall conductivity a and Hall conductivity b of graphene ML for three values of magnetic fields, B = 0.0 (black), 0.1 (red), and 0.2 (blue) T. Panels c and d show the contributions for conductivities of different valleys K (solid) and K (dashed).Here we used an energy gap parameter ∆ = 10meV.

FIG. 17 .
FIG. 17. SUPPORTING INFO 12.Comparison of the energy bands of zigzag Gr/hBN nanoribbons with 50 nm breadth and staggered sublattice potential ∆ = 5 meV for: a model with pz orbitals, and b model with pz,dxz and dyz orbitals.c Energy bands of Armchair Gr/hBN nanoribbons with the same width and staggered sublattice potential.
FIG. 18. SUPPORTING INFO 13. non-local resistance for current injection across 1 − 4 in 8-terminal devices with 13 nm breadth and 157 nm distance between the injection and collection leads for: a the model with pz orbitals and b the 6-bands model.