Bulk valley transport and Berry curvature spreading at the edge of flat bands

2D materials based superlattices have emerged as a promising platform to modulate band structure and its symmetries. In particular, moiré periodicity in twisted graphene systems produces flat Chern bands. The recent observation of anomalous Hall effect (AHE) and orbital magnetism in twisted bilayer graphene has been associated with spontaneous symmetry breaking of such Chern bands. However, the valley Hall state as a precursor of AHE state, when time-reversal symmetry is still protected, has not been observed. Our work probes this precursor state using the valley Hall effect. We show that broken inversion symmetry in twisted double bilayer graphene (TDBG) facilitates the generation of bulk valley current by reporting experimental evidence of nonlocal transport in a nearly flat band system. Despite the spread of Berry curvature hotspots and reduced quasiparticle velocities of the carriers in these flat bands, we observe large nonlocal voltage several micrometers away from the charge current path — this persists when the Fermi energy lies inside a gap with large Berry curvature. The high sensitivity of the nonlocal voltage to gate tunable carrier density and gap modulating perpendicular electric field makes TDBG an attractive platform for valley-twistronics based on flat bands.

Supplementary Figure 5: Estimating angle inhomogeneity. Variation of two-probe and four-probe local resistance using different terminals of the device as a function of V TG at V BG = 0 V. In R i j,kl , the indices i, j denote the voltage probes and k, l denote the current probes. The difference between the positions of two moiré peaks is used to estimate the local twist angle. Inset: An optical microscope image of the device with the terminals numbered. The scale bar corresponds to 5 µm.  (a) Color scale plot of local resistance as a function of charge density and electric displacement field. Inset: Optical microscope image of the device with electrodes labeled by numbers. Scale bar corresponds to 4 µm.
(b) Color scale plot of nonlocal resistance as a function of charge density and electric displacement field.
For measuring the nonlocal resistance, the current was sent through the terminals 9 and 4, and voltage was measured across the terminals 11 and 2.  Figure 9: Decay of nonlocal resistance over length. (a-c) The variation of the nonlocal resistance as a function of length for n = −n S (a), n = 0 (b) and n = n S (c) for Device 1 and Device 2. Here, R NL has been normalized with respect to the nonlocal resistance measured in the nearest probe, R max NL , for both the devices.

Supplementary Note 1: Calculation of band structure and valley Hall conductivity in TDBG
As discussed in the theory paper Ref. 3 , moiré bands theory for the moiré pattern superlattice 4 and the accurate continuum models 5 are used to obtain the electronic structure of the twisted double bilayer graphene (TDBG). The continuum model of Bistritzer-MacDonald for the twisted bilayer graphene (TBG) 4 is extended to the case of twisted double bilayer graphene (TDBG), the Hamiltonian of TDBG at the valley K with the interlayer coupling between the twisted layers through a first-harmonic stacking-dependent interlayer tunneling function, and subject to ∆ i intralayer potentials as where h ± t/b = h t/b (±θ /2) such that the relative twist angle between the bilayers is θ .
The Dirac Hamiltonian given by h(θ ) = υ FR−θ p·σ σ σ includes a phase shift due to a rotationR −θ such that e ±iθ p → e ±i(θ p −θ ) , where σ σ σ = (σ x , σ y ) and σ z are the graphene sublattice pseudospin Pauli matrices, and the momentum is defined in the xy plane p = (p x , p y ), where we assume K valley unless stated otherwise. The Fermi velocity υ F = υ 0 defined from υ i = √ 3|t i |a/2h is related to the intralayer nearest-neighbor hopping term t 0 = -3.1 eV that captures the experimental moiré band features better 6 .
The top and bottom bilayer graphene (BG) are labeled through the positive/negative (+/-) rotation signs, while in turn we have top/bottom (t/b) graphene layers within each BG that are coupled through the matrices t ± s . The interlayer coupling model of a bilayer graphene is given by satisfying t s=+ = t † s=− for AB or BA (s = ±1) stacking-dependent interlayer coupling that consists of a minimal coupling term t s = t 1 (σ x − isσ y )/2 plus remote hopping contributions through the terms t 3 = 0.283, t 4 = 0.138 eV, giving rise to trigonal warping and electron-hole asymmetry. The π ± operators include the phases due to ±θ /2 layer rotation. The Hamiltonian of graphene is given by h ± l (θ ) = h ± (θ ) + δ (1 − lsσ z )/2 where the second term adds a δ =0.015 eV sublattice potential at the higher energy dimer sites at the t/b layers l = ± 7 , that depends on AB or BA stacking s = ±, respectively. The site potentials∆ i are mapped on its sublattices through∆ i = ∆ i 1 where i = 1, 2, 3, 4 are the layer labels from top to bottom, and 1 is a 2 × 2 identity matrix. Here, we have an additional control knob to change the electronic structure through a perpendicular external electric field that modifies the interlayer potential ∆ i values in equation (1).
We can identify the interlayer tunneling with the first harmonic expansion coefficient of the interlayer coupling such that t 1 = 3ω 5,7 , and for simplicity we use the same AB stacking tunneling within each Bernal BG and the twisted interfaces. In the small-angle approximation, the interlayer coupling Hamiltonian is given by where the three Q j vectors Q 0 = Kθ (0, −1) and Q ± = Kθ (± √ 3/2, 1/2) are proportional to twist angle θ and K = 4π/3a is the Brillouin zone corner length of graphene, whose lattice constant is a = 2.46 Å, and here the indices l, l label the sublattices of neighboring twisted surface layers. The interlayer coupling matrices between the two rotated adjacent layers are given by using a form that distinguishes interlayer tunneling matrix elements ω = ω BA and ω = ω AA for different and same sublattice sites between the layers. The convention taken here for the T j matrices 5 assumes an initial AA stacking configuration τ = (0, 0) and differs by a phase factor with respect to the initial AB stacking τ = (0, a/ √ 3) 8,9 . The greater interlayer separation c compared to the carbon-carbon distances a CC lead to slowly varying interlayer tunneling function T (r) and the moiré patterns can often be accurately described within a first-harmonic approximation 4,5 . Furthermore, the effects of atomic relaxation in the moiré patterns lead to corrugations that have non-negligible effects in the details of the electronic structure for both intralayer potentials and interlayer coupling 10 . In the case of ω = ω, we proposed a single parameter relation through ω = C 1 ω 2 +C 2 ω +C 3 , where C 1 = −0.5506, C 2 = 1.036, C 3 = −0.02245 as discussed in Ref. 3

(in its Supplementary information). Our calculations have used a configuration space with variable
cutoff in momentum space of a radius of up to 6G 1 = 24πθ /( √ 3a) using Hamiltonian matrices with sizes as large as 676 × 676 such that θ ω/(12π|t 0 |) to obtain converged results in the limit of small θ and large ω.
The possibility of band gap opening at charge neutrality point (CNP), primary gap (δ p ), through an electric field together with the presence of moiré gaps (secondary gaps, δ s ) with the higher-energy bands leads to well-defined valley Chern numbers. The valley Chern numbers were calculated through by integrating the moiré Brillouin zone for each valley the Berry curvature for the q-th band through 11 where for every k point we take sums through all the neighboring q bands, the |u q are the moiré superlattice Bloch states, and E q are the eigenvalues.
The Hall conductivity that results due to a nonzero Berry curvature (Ω ) at a particular valley of graphene (K or K ) is given by 12 where q indicates the band index and f (ε q (k), E F ) = 1 e (εq(k)−E F )/K B T +1 denotes the Fermi occupation function. At low temperature (shown in orange curve for T = 5 K in Supplementary Fig. 3), it saturates to 2e 2 h in the band gap between the flat bands ( Supplementary Fig. 3a). Away from the CNP gap, it starts decreasing since the Berry curvature in the low energy conduction band and valence band have opposite signs, as seen from the line plots in Supplementary Fig. 2. This decrease of σ xy away from the CNP gap in TDBG is asymmetric in nature, which is in contrary to that obtained for bilayer graphene 13,14 . The Hall conductivity is also shown for an elevated temperature of T = 50 K (blue curve) where we have observed bulk valley transport at the CNP gap. The decrease in σ xy from its low temperature value of 2e 2 h at the gap is due to the thermal excitation of valence band electrons to the conduction band. The valley Hall conductivity, σ VH xy , is obtained by adding the contribution to the Hall conductivity from the individual valleys at K and K , and is given by σ VH xy = 2σ xy 13 .
Supplementary Fig. 3b shows finite non-zero σ xy at the moiré gaps at low interlayer potential of ∆ = 1 meV. Our observation of nonlocal resistance at the moiré gaps is consistent with the fact that nonlocal resistance was also observed at the moiré gaps of other systems like hBN aligned monolayer graphene 1 .

Supplementary Note 2: Nonlocal transport measurement
The transport measurement reported in the main manuscript and all other sections in the supplementary was measured using a low frequency (∼ 17 Hz) lock-in technique by sending a current ∼ 10 nA and measuring the voltage after amplifying using SR560 preamplifier or preamplifier model 1021 by DL instruments, Ithaca. Since the nonlocal resistance is appreciable when the system is gapped, spurious signal can be measured for these high resistivity states 13,15 . To verify the nonlocal resistance we measured is not a measurement artifact, we employ Keithley 2182 nanovoltmeter to measure dc voltage while sending current using Keithley 6221 current source. The nanovoltmeter has input impedance > 10 GΩ while the SR560 or DL 1021 preamplifier has an input impedance of 100 MΩ. In Supplementary Fig. 4a, we have plotted the nonlocal resistance as a function of top gate voltage V TG around the CNP at V BG = 28 V measured both by the lock-in method aided by the preamplifier and the dc measurement using nanovoltmeter. For measurement using the nanovoltmeter, we measured resistance both in forward and reverse directions and took the average resistance to nullify any dc voltage drop due to thermo-electric effect at various junctions of the current path inside the cryostat. In Supplementary Fig. 4b, we plotted the nonlocal resistance of the moiré peak at the hole side as a function of the perpendicular electric field using the two schemes. As seen from both Figs. 4(a) and 4(b), the resistance values are independent of the measurement schemes. All the measurements subsequently were done using lock-in with the preamplifier.

Supplementary Note 3: Angle homogeneity for nonlocal detection of valley current
The twist angle homogeneity is an important prerequisite for measuring nonlocal valley transport in twisted graphene devices. This is because the generation of valley current in the injection probes happens when the Fermi energy lies in a Berry curvature hotspot, i.e., the Fermi energy lies near a gap. This requires the charge density to be tuned by the gates to specific values: n = 0 (CNP gap) or n = ±n S (moiré gap), where To estimate the local twist angle near various probes we measure two-probe and four-probe resistance using different combinations of current and voltage probes. In Supplementary Fig. 5, we present such different plots of local resistance as a function of V TG for V BG = 0 V. All the curves have two moiré peaks corresponding to n = ±n S . The difference in the positions of the two peaks on the n-axis corresponds to 2n S , which is used to estimate the local twist angle. We find that the twist angle to vary between 1.14°to 1.19°, establishing that the device has good angle homogeneity.

Supplementary Note 4: Reciprocity of the nonlocal resistance measurement
For the Onsager reciprocal relations to be valid, the nonlocal resistance we measure should be the same if one swaps the injection and the detection terminals 13 . We verify this in Supplementary Fig. 6 where we present nonlocal resistance as a function of V TG with V BG = 0 V for two reciprocal combinations of injection and detection probes. We find that the nonlocal resistance does not change if we swap the terminals.

Supplementary Note 5: Additional data on scaling
A tell-tale signature of bulk valley current is cubic scaling between the nonlocal and local resistances. In  Fig. 7i shows cubic scaling in the high electric field end, while the scaling deviates from cubic and saturates towards the low electric field end where the band gap at n = n S is higher 2 . This saturation at high band gap regime is attributed to large valley Hall angle physics 16 and is also seen in earlier studies at the CNP gap of bilayer graphene 13 . In Supplementary Fig. 7f we find that the CNP shows cubic scaling at high temperatures. As the temperature is lowered (brown and black curve in Supplementary Fig. 7f), the scaling deviates from cubic to higher exponents in the high electric field end. This is similar to that in Fig. 3f of the main manuscript, where high D and low T shows the same departure.

Supplementary Note 6: Data from the second device
In Supplementary Fig. 8, we present local and nonlocal resistance from Device 2 which has a twist angle of 1.24°.

Supplementary Note 7: Decay of Nonlocal resistance with length
The dependence of the nonlocal resistance on the length of the current path is governed by the valley diffusion length l v , as seen from equation (8), Here, σ VH xy is the valley Hall conductivity. L and W represent the length and width of the Hall bar channel, respectively. To extract l v we plot the decay of the nonlocal signal along the sample length at the moire gaps and the CNP in Supplementary Fig. 9. We fit the exponential decay using equation (8) and extract out l v for both the devices. The data for Device 1 is plotted at T = 35 K to negate low temperature effects. For moiré gaps, we find l v to be 0.9 µm and 1.8 µm for Device 1 and Device 2 respectively. At the CNP, l v is 1.15 µm for Device 1. These values are similar to those obtained in bilayer graphene 13 .

Supplementary Note 8: Effect of band flatness on Berry curvature hotspot
To understand the effect of a flat band on the Berry curvature hotspot, we consider a simple band structure of gapped monolayer graphene. We choose the monolayer graphene Hamiltonian since, for a flat band, in any twisted system, there is no reason in general for the linear term in momentum to be zero. The Hamiltonian for the gapped monolayer graphene is given by, H =hv F σ σ σ .k + ∆ g σ 3 . Here 2∆ g is the band gap between two energy bands given by E ± = ± (hv F k) 2 + ∆ 2 g . We incorporate the band flatness by renormalizing the term v F in the Hamiltonian, whereas in the case of monolayer graphene v F = 10 6 m/s = v F0 . The Berry curvature is given by With the charge density n given by n = k 2 F /π and δ = ∆ g /hv F , the Berry curvature at the Fermi wave vector k F is given by Ω (k F ) = δ /[2(nπ + δ 2 ) 3/2 ]. In Supplementary Fig. 10a we plot the Berry curvature as a function of charge density for monolayer graphene with a gap of 2∆ g = 10 meV and for different values of the renormalized velocity, v F . As v F is smaller, i.e., the band is flatter, we find that the Berry curvature is more delocalized away from the gap. In Supplementary Fig. 10b we plot the contribution of the conduction band to the valley Hall conductivity given by, σ xy = (e 2 /2h) Ω (k)kdkdθ = (e 2 /2h)(1 − δ / √ nπ + δ 2 ). We again find that in the case of flat band the Hall conductivity saturates to its asymptotic value much slower. The Berry curvature delocalization and slow saturation of Hall conductivity in flat band systems, as represented in Supplementary   Fig. 10, can be reproduced even by starting with a gapped bilayer graphene Hamiltonian.
Supplementary Note 9: Comparison of extent of nonlocal resistance peak with a system without flat bands Here we compare the extent of the nonlocal signal in the charge density axis from our TDBG device with that of hBN aligned graphene device from Gorbachev et al. 1 . We compare the charge neutrality peak at T = 20 K for both the devices. Upon application of a perpendicular electric field (shown for D/ε 0 = −0.4 V nm −1 in Supplementary Fig. 11), a gap opens up at charge neutrality within the flat bands in the TDBG device 2 .
These bands have Berry curvature hotspots, as shown in Supplementary Fig. 2. For the case of hBN aligned graphene, the superlattice potential results in opening up a gap between the valence band and conduction band at the charge neutrality 17 . This, too, results in Berry curvature hotspots at the band edges near the gap opening. However, unlike the case in small-angle TDBG, the resulting low energy bands are not flat. The broader extent of the nonlocal resistance in our data (the FWHM being ∼5 times larger) can be attributed to the spreading of Berry curvature hotspot due to flat bands in the TDBG system.

Supplementary Note 10: Modulation in FWHM of nonlocal peak by tuning Fermi velocity
We now further demonstrate that the broadening of the nonlocal resistance peak in the charge density axis is indeed due to Berry curvature hotspot spreading as Fermi velocity, v F , decreases. We can tune v F by changing electric field within the same TDBG device having a fixed twist angle. This is because the perpendicular electric field can modulate the bandwidth of the flat bands in TDBG, as shown in the inset of Supplementary   Fig. 12a and also discussed in Chebrolu et al. 3 . As the bandwidth reduces, that is the band becomes flatter, the v F reduces. We have plotted the nonlocal resistance in the 1.18 • TDBG as a function of charge density across the CNP for three different electric fields in Supplementary Fig. 12a. We find that as we increase the electric field |D| (i.e., as the bandwidth of flat bands increases), the width of the nonlocal resistance peak becomes narrower, thus confirming that lower Fermi velocity at low electric field results in wider nonlocal resistance peak.
We systematically extract the FWHM and plot its variation with the electric field for two devices with twist angle 1.18 • (Supplementary Fig. 12b) and 1.24 • (Supplementary Fig. 12c). We find that this feature of increase in FWHM of R NL with decrease in v F of flat bands (by decreasing |D|) is repeatable (left axis corresponding to cyan colored curve in Supplementary Fig. 12b and Supplementary Fig. 12c). To further rule out any effect of local resistance, we also plot the FWHM of nonlocal resistance peak after normalizing with the FWHM of local resistance (right axis corresponding to purple colored curve in Supplementary   Fig. 12b and Supplementary Fig. 12c), which shows the same trend. An increase in the ratio by decreasing electric field |D| signifies that a decrease in Fermi velocity of the flat bands spreads the Berry curvature hotspots in k-space-this translates to expanding the FWHM width of the nonlocal resistance peak at CNP.
This reinforces our understanding of the physics of Berry curvature spreading due to flat bands with low v F .