Device fabrication
The hBN-encapsulated ABA graphene heterostructure was fabricated using the dry-transfer method. The graphene flakes were first exfoliated onto a Si/SiO2 (285 nm) substrate. The number of layers in the graphene flakes was determined using Raman microscopy49. Then, the hBN (about 30 nm thick) and the graphene flakes were picked up using a polycarbonate on a polydimethylsiloxane dome stamp. The stacks were then released onto a pre-annealed Ti (2 nm)/Pt (10 nm) bottom gate, patterned on the Si/SiO2 wafer. The finalized stacks were annealed in a vacuum at 500 °C for strain release50. A Ti (2 nm)/Pt (10 nm) top gate was then deposited on top of the stack. The one-dimensional contacts were formed by SF6 and O2 plasma etching followed by evaporating Cr (4 nm)/Au (70 nm). Then, the device was etched into a Hall bar geometry. Finally, the device was re-annealed at 350 °C in a vacuum. The capacitances per unit area of the bottom and top gates are Cbg = 0.649 × 1012 e cm−2 V−1, Ctg = 0.668 × 1012 e cm−2 V−1. The top and bottom gates are used to control the carrier density n = (CbgVbg + CtgVtg)/e and the effective transverse displacement field D = (CtgVtg − CbgVbg)/2ε0, where ε0 is vacuum permittivity. From fitting the experimental QOs to simulations, we find that D = 1 V nm −1 corresponds to the energy difference between the adjacent graphene layers of Δ1 = 92 meV.
Transport measurements
Transport characterization of ABA graphene devices was carried out using standard lock-in techniques. The Rxx shows a peak along the diagonal charge neutrality line that increases with D, suggesting a gap opening (Extended Data Fig. 1a). The Landau fan shows LL crossings (Extended Data Fig. 1b), consistent with the previous reports3,37,38,40,51,52,53,54,55,56. The QOs from MLG band LLs are visible at low fields, but the BLG LLs can be only resolved above 0.75 T on the electron side and at notably higher fields on the hole doping side (Extended Data Fig. 1c).
SOT measurements and magnetization reconstruction
The local magnetic measurements were conducted in a custom-built scanning SOT microscope in a cryogen-free dilution refrigerator (Leiden CF1200) at a temperature of 160–350 mK (ref. 57). Indium SOT with an effective diameter of about 150 nm and magnetic sensitivity of 20 nT Hz−1/2 was fabricated as described previously33,58,59. The SOT readout circuit is based on SQUID series array amplifier60,61. The SOT is attached to a quartz tuning fork vibrating at about 32.8 kHz (Model TB38, HMI Frequency Technology), which is used as a force sensor for tip height control62. The scanning height was about 150 nm above the ABA graphene. An a.c. voltage \({V}_{{\rm{b}}{\rm{g}}}^{{\rm{a}}{\rm{c}}}\) at a frequency of about 1.8 kHz was applied to the bottom gate to modulate the carrier density by \({n}^{{\rm{a}}{\rm{c}}}={C}_{{\rm{b}}{\rm{g}}}{V}_{{\rm{b}}{\rm{g}}}^{{\rm{a}}{\rm{c}}}/e\). A lock-in amplifier was used to measure the corresponding local \({B}_{z}^{{\rm{a}}{\rm{c}}}\) by the scanning SOT. The \({B}_{z}^{{\rm{a}}{\rm{c}}}\) data were symmetrized with respect to the displacement field D where applicable. In contrast to other scanning techniques, the magnetic signal is transparent to the metallic top gate, enabling the investigation of a wide range of heterostructures and encapsulated devices.
The 2D \({B}_{z}^{{\rm{a}}{\rm{c}}}(x,y)\) images were used to reconstruct the magnetization mz(x, y) using the numerical inversion procedure described in ref. 63 (Extended Data Fig. 2). As the reconstruction of mz requires 2D \({B}_{z}^{{\rm{a}}{\rm{c}}}(x,y)\) information, the QOs at a single location or along the one-dimensional line scans are presented in the main text as the raw data of \({B}_{z}^{{\rm{a}}{\rm{c}}}\).
Magnetic field and modulation amplitude dependence of QOs
The measured signal \({B}_{z}^{{\rm{a}}{\rm{c}}}={n}^{{\rm{a}}{\rm{c}}}({\rm{d}}{B}_{z}/{\rm{d}}n)\) is proportional to the modulation amplitude of the carrier density nac induced by \({V}_{{\rm{b}}{\rm{g}}}^{{\rm{a}}{\rm{c}}}\). It is therefore desirable to use large nac to improve the signal-to-noise ratio. To resolve QOs, however, nac has to be substantially smaller than the period of the oscillations Δn. Extended Data Fig. 2c–e shows the QOs acquired at Ba = 320 mT using \({V}_{{\rm{b}}{\rm{g}}}^{{\rm{a}}{\rm{c}}}=8\,{\rm{m}}{\rm{V}}\), 35 mV and 100 mV rms corresponding to nac of 5.19 × 109 cm−2, 2.27 × 1010 cm−2 and 6.49 × 1011 cm−2 rms, respectively. The four-fold degenerate BLG LLs have a period of Δn = 4Ba/ϕ0 = 3.1 × 1010 cm−2. The lowest \({V}_{{\rm{b}}{\rm{g}}}^{{\rm{a}}{\rm{c}}}=8\,{\rm{m}}{\rm{V}}\) rms was chosen to result in a peak-to-peak value of nac of 1.47 × 1010 cm−2, approximately equal to Δn/2 = 1.55 × 1010 cm−2, which results in an optimal signal-to-noise ratio for detecting the BLG LLs, albeit suppresses the measured \({B}_{z}^{{\rm{a}}{\rm{c}}}/{n}^{{\rm{a}}{\rm{c}}}\) ratio by a factor of π/2. A larger nac washes out the QOs from the BLG LLs, leaving the MLG LLs resolvable as demonstrated in Extended Data Fig. 2d,e. The largest nac also enables observation of the paramagnetic response ∂M/∂μ = C/ϕ0 in the gap between the zeroth and the first MLG LLs dictated by the Chern number C = 2 on the electron side and C = −2 on the hole side (Extended Data Fig. 2e).
Extended Data Fig. 2f–h shows the QOs at Ba = 40 mT, 80 mT and 170 mT. At these low fields, the Dingle broadening greatly suppresses the QOs due to BLG LLs (Extended Data Fig. 4) and reduces the visibility of the MLG LLs at large displacement fields because of the reduction in the gap energies. At 170 mT and \({V}_{{\rm{b}}{\rm{g}}}^{{\rm{a}}{\rm{c}}}=8\,{\rm{m}}{\rm{V}}\) rms, the M2 LLs and the 12-fold degenerate LLs in the gullies are resolved as seen in Extended Data Fig. 2h.
BS calculations
The BS of ABA graphene was calculated in the tight-binding model following refs. 2,36 based on SWMc parameterization34. On the basis of {A1, B1, A2, B2, A3, B3}, where Ai and Bi are the two sublattice sites in the ith layer, the low-energy effective Hamiltonian can be written as
$${H}_{0}=\left(\begin{array}{cccccc}{\varDelta }_{1}+{\varDelta }_{2} & {v}_{0}{\pi }^{\dagger } & {v}_{4}{\pi }^{\dagger } & {v}_{3}\pi & {\gamma }_{2}/2 & 0\\ {v}_{0}\pi & \delta +{\varDelta }_{1}+{\varDelta }_{2} & {\gamma }_{1} & {v}_{4}{\pi }^{\dagger } & 0 & {\gamma }_{5}/2\\ {v}_{4}\pi & {\gamma }_{1} & \delta -2{\varDelta }_{2} & {v}_{0}{\pi }^{\dagger } & {v}_{4}\pi & {\gamma }_{1}\\ {v}_{3}{\pi }^{\dagger } & {v}_{4}\pi & {v}_{0}\pi & -2{\varDelta }_{2} & {v}_{3}{\pi }^{\dagger } & {v}_{4}\pi \\ {\gamma }_{2}/2 & 0 & {v}_{4}{\pi }^{\dagger } & {v}_{3}\pi & -{\varDelta }_{1}+{\varDelta }_{2} & {v}_{0}{\pi }^{\dagger }\\ 0 & {\gamma }_{5}/2 & {\gamma }_{1} & {v}_{4}{\pi }^{\dagger } & {v}_{0}\pi & \delta -{\varDelta }_{1}+{\varDelta }_{2}\end{array}\right)$$
where Δ1 = −e(U1 − U3)/2 and Δ2 = −e(U1 − 2U2 + U3)/6, with Ui the potential of layer i. Δ1 is determined by the displacement field, whereas Δ2 describes the asymmetry of the electric field between the layers. The band velocities vi (i = 0, 3, 4) are related to the tight-binding parameters γi by \({v}_{i}\hbar =\frac{\sqrt{3}}{2}{a}_{{\rm{c}}}{\gamma }_{i}\), where ac = 0.246 nm is the crystal constant of graphene, π = ξkx + iky, and ξ is the valley index (ξ = ±1 for valley K+ and K−, respectively).
On a rotated basis (A1 − A3)/\(\sqrt{2}\), (B1 − B3)/\(\sqrt{2}\), (A1 + A3)/\(\sqrt{2}\), B2, A2, (B1 + B3)/\(\sqrt{2}\), the Hamiltonian can be rewritten as
$${H}_{{\rm{TLG}}}=\left(\begin{array}{cccccc}{\varDelta }_{2}-\frac{{\gamma }_{2}}{2} & {v}_{0}{\pi }^{\dagger } & {\varDelta }_{1} & 0 & 0 & 0\\ {v}_{0}\pi & {\varDelta }_{2}+\delta -\frac{{\gamma }_{5}}{2} & 0 & 0 & 0 & {\varDelta }_{1}\\ {\varDelta }_{1} & 0 & {\varDelta }_{2}+\frac{{\gamma }_{2}}{2} & \sqrt{2}{v}_{3}\pi & -\sqrt{2}{v}_{4}{\pi }^{\dagger } & {v}_{0}{\pi }^{\dagger }\\ 0 & 0 & \sqrt{2}{v}_{3}{\pi }^{\dagger } & -2{\varDelta }_{2} & {v}_{0}\pi & -\sqrt{2}{v}_{4}\pi \\ 0 & 0 & -\sqrt{2}{v}_{4}\pi & {v}_{0}{\pi }^{\dagger } & \delta -2{\varDelta }_{2} & \sqrt{2}{\gamma }_{1}\\ 0 & {\varDelta }_{1} & {v}_{0}\pi & -\sqrt{2}{v}_{4}{\pi }^{\dagger } & \sqrt{2}{\gamma }_{1} & {\varDelta }_{2}+\delta +\frac{{\gamma }_{5}}{2}\end{array}\right)$$
For Δ1 = 0, the Hamiltonian can be block-diagonalized into MLG-like and BLG-like blocks, that is, HTLG = HMLG ⊕ HBLG. A finite displacement field hybridizes the two blocks.
In an external magnetic field, in the Landau gauge, the canonical momentum π can be replaced by π − eA, where A is the vector potential. π obeys the commutation relation [πx, πy] = −i/lB, where \({l}_{B}=\sqrt{\left({\hbar }/{eB}\right)}\) is the magnetic length. As in the usual one-dimensional harmonic oscillator, on the basis of LL orbital |n⟩, the matrix elements of π, π† are given by
$$\begin{array}{c}{{\rm{K}}}^{+}:\pi | n\rangle =\frac{i\hbar }{{l}_{B}}\sqrt{2(n+1)}| n+1\rangle \\ {\pi }^{\dagger }| n\rangle =-\frac{i\hbar }{{l}_{B}}\sqrt{2n}| n-1\rangle \\ {{\rm{K}}}^{-}:\pi | n\rangle =\frac{i\hbar }{{l}_{B}}\sqrt{2n}| n-1\rangle \\ {\pi }^{\dagger }| n\rangle =-\frac{i\hbar }{{l}_{B}}\sqrt{2(n+1)}| n+1\rangle \end{array}$$
Therefore, the new Hamiltonian can be written on the basis of LL orbitals. Using matrix elements of π and π† operators, the momentum operators are replaced by raising and lowering the diagonal matrix of dimensions Λ × Λ, where Λ is the cutoff number for the infinite matrix, restricting the Hilbert space with indices n ≤ Λ. All the other nonzero elements γi are substituted by γiIΛ, where IΛ is the identity matrix with dimensions Λ × Λ. As our measurements were performed in low magnetic fields and high-index LLs are often involved, a large cutoff was used so that it spans the energy range significantly larger than in the experiment. We also removed false LLs caused by imposing the cutoff, which usually have very large indices. In the simulations, Λ was set to 400 for small carrier-density ranges (Fig. 2) and to 800 for calculations over larger ranges (Figs. 3 and 4).
Evolution of the BS and LLs with displacement field
Extended Data Fig. 3 shows the calculated BS of ABA graphene using the derived SWMc parameters and the evolution of the LLs with D and Ba. At D = 0 (Δ1 = 0), there is essentially no hybridization between the MLG and BLG bands. All the LLs are valley (and spin) degenerate except for the zeroth LLs of the MLG and BLG bands that are valley polarized because of the Berry curvature (Extended Data Fig. 3a). With increasing Δ1, the gaps of the MLG and BLG bands increase and the hybridization between the bands grows resulting in the formation of mini-Dirac cones (gullies) and in LL anticrossings (Extended Data Fig. 3b–d). At our highest accessible Δ1 ≈ 50 meV, the lowest LLs in the gullies are well isolated from the rest of the LLs as shown in Extended Data Fig. 3e,f. As the BLG bandgap \({\varDelta }_{{\rm{G}}}^{0}\) is characterized by C = 0, it has no magnetization. The six-fold degenerated compressible zeroth LLs \({0}_{{\rm{G}}}^{+}\) and \({0}_{{\rm{G}}}^{-}\) in the gullies also have no magnetization at low fields, M = −∂ε/∂B = 0, because of their zero kinetic energy. As a result, zero magnetization is observed around the CNP over a width of δn = 12Ba/ϕ0 in carrier density as indicated in Fig. 2c,f. The first paramagnetic signal appears when the Fermi level reaches the C = ±6 gaps \({\varDelta }_{{\rm{G}}}^{1}\) and \({\varDelta }_{{\rm{G}}}^{-1}\) between the zeroth and the first gully LLs as shown in Fig. 2f. At elevated magnetic fields, the six-fold gully degeneracy of the zeroth LLs is partially lifted4,53.
Reconstruction of BS parameters
Several experimental studies3,4,37,38,40,51,52,54,64 have investigated the tight-binding parameters of ABA graphene as shown in Extended Data Table 1. The high resolution of our data and the fine features attained at low magnetic fields allow high-precision reconstruction of SWMc parameters as follows. We set γ0 to the standard literature value of 3,100 meV, which corresponds to Fermi velocity of graphene \({v}_{{\rm{F}}}=\frac{\sqrt{3}}{2\hbar }{a}_{{\rm{c}}}{\gamma }_{0}=1{0}^{6}\,{\rm{m}}\,{{\rm{s}}}^{-1}\). The γ0 sets the overall energy scale, whereas the value of the remaining seven parameters, relative to γ0, determine the BS. The fitting of the parameters was performed manually. We first determined the effect of the individual parameters on particular features of the BS as shown in Extended Data Fig. 4, which then guided us in the iterative fitting process. In particular, in the absence of displacement field, Δ1 = 0, the MLG band is affected by only γ0, γ2, γ5 and δ, with the gap at the Dirac point given by \({E}_{{\rm{g}}}^{0}=\delta +\frac{{\gamma }_{2}-{\gamma }_{5}}{2}\). The BLG band is strongly dependent on γ0, γ1 and γ3, weakly dependent on γ4 and essentially independent of γ5 and δ. The BLG gap size is mainly governed by γ2 and Δ2. The relative energy shift between the MLG and BLG bands is mainly governed by γ2.
The dependence of the measured QOs on n and D at low Ba provides a very sensitive tool for determining the SWMc parameters. After developing an understanding of the influence of the individual parameters on the relative position of the LLs in specific regions in the (n, D) plane, an initial set of parameters was chosen to attain an approximate fit to the data. Then fine-tuning of the parameters is achieved by calculating the QOs for each set of parameters and comparing with the data at D = 0 V nm−1. This process is repeated manually adjusting the different parameters in an iterative manner. After attaining a good fit at D = 0, additional fine-tuning was performed to fit the entire range of D. As the different parameters have a distinctive effect on the relative positions of the LLs, this manual procedure is readily manageable. The error bars were determined by the values of the individual parameters for which a visible deviation from the data was observed.
The following attributes were particularly informative for the fitting processes:
-
1.
The number of BLG LLs between the adjacent MLG LLs
-
2.
The relative energy shift between MLG and BLG bands
-
3.
LL anticrossings in the gullies
-
4.
The gap size of MLG band
Attribute 1 is determined by the DOS ratio of the two bands, which is predominantly governed by γ1. By adjusting γ1 to fit the relative number of BLG and MLG LLs along with optimization of other parameters we obtain γ1 = 370 ± 10 meV.
Attribute 2 is then used to determine γ2. The energies of the band extrema and hence the relative position of the zeroth LLs can be calculated analytically. In particular, for Δ1 = 0, the \({0}_{{\rm{M1}}}^{-}\) LL at the bottom of M1 band is positioned at energy Δ2 − γ2/2, whereas the top of BLG valence band is at Δ2 + γ2/2. Thus, the relative position between MLG and BLG bands is determined by γ2 and Δ2. As the LL spectrum is quite sensitive to Δ2, γ2 is determined first. We use the relative position between −1M3 and the nearby BLG LLs to fit γ2, and we get γ2 = −19 ± 0.5 meV.
Attribute 3 is governed by γ3, which induces trigonal warping of the BLG bands. As shown in Extended Data Fig. 6, this results in the anticrossings between the BLG LLs and MLG \({0}_{{\rm{M3}}}^{-}\) and \({-1}_{{\rm{M3}}}^{+}\) LLs. From fitting to the experimental data, we obtain γ3 = 315 ± 10 meV.
Attributes 2 and 4 are used to derive δ and γ5. The MLG band gap at D = 0 V nm−1 is \({E}_{{\rm{g}}}^{0}=\delta +\left({\gamma }_{2}-{\gamma }_{5}\right)/2\), whereas the gap centre is located at 2Δ2 + δ − (γ2 + γ5)/2. In our experimental data, one BLG LL fits within the MLG gap and 20 BLG LLs reside between \({0}_{{\rm{M1}}}^{-}\) and \({-1}_{{\rm{M3}}}\), from which we attain δ = 18.5 ± 0.5 meV and γ5 = 20 ± 0.5 meV. Note that \({E}_{{\rm{g}}}^{0}\) can be either positive or negative. We find that \({E}_{{\rm{g}}}^{0}\) is negative, which means that the zeroth K− LL (\({0}_{{\rm{M1}}}^{-}\)) resides at the bottom of the M1 band and the zeroth K+ LL (\({0}_{{\rm{M2}}}^{+}\)) is at the top of M2. In this case, the Dirac gap Eg increases with Δ1 and the \({0}_{{\rm{M1}}}^{-}\) and the \({0}_{{\rm{M2}}}^{+}\) LLs spread apart with the displacement field as shown in Extended Data Fig. 3f, consistent with experimental data in Fig. 2c,e and calculations in Fig. 2d. If \({E}_{{\rm{g}}}^{0}\) is positive, the zeroth K− LL will reside at the top of M2, whereas the zeroth K+ LL will be at the bottom of M1. In this case, on increasing D, the Dirac gap closes and then reopens with the crossing of the two zeroth LLs, such that Eg is always negative at high D with zeroth K− LL at the bottom of M1. Extended Data Table 1 shows that the value and the sign of \({E}_{{\rm{g}}}^{0}\) varies notably in the literature. However, only in ref. 40 and in the present work the Dirac gap is reported directly. For the rest of the references, the \({E}_{{\rm{g}}}^{0}\) values presented in the table are calculated from the reported values of δ, γ2 and γ5.
Δ2 mainly affects the gap of the BLG bands and as −1M3 resides closely to the BLG band gap, we use the number of BLG LLs between −1M3 and −2M3 to fit Δ2 and get Δ2 = 3.8 ± 0.05 meV. γ4 plays the most negligible role, slightly adjusting the shape of the BLG bands. The fitting procedure is to choose these parameters such that the inaccuracy of the number of BLG LLs between any pair of MLG LLs is no more than one. By optimizing all parameters for best fit to the experimental data, we derive γ4 = 140 ± 15 meV, as shown in Extended Data Table 1.
Orbital magnetization calculations
Oscillations in orbital magnetization M from the LLs can be calculated analytically for either parabolic or Dirac bands as shown previously65. However, there is no analytical expression for the LL spectrum in ABA graphene; therefore, the magnetization oscillations have to be calculated numerically. We follow the method described in ref. 14 to derive the magnetization M(n) and then calculate its derivative ∂M/∂n.
We first consider the case with zero LL broadening. For an arbitrary LL spectrum Ei with degeneracy Di (i is the Landau-level index), the DOS N0(ε) of the system is
$${N}_{0}\left(\varepsilon \right)=\sum _{i}{D}_{i}\delta \left(\varepsilon -{E}_{i}\right).$$
Ei describes spin-degenerate LLs from both valleys with degeneracy \({D}_{i}=2\frac{eB}{h}\). The grand thermodynamic potential Ω0(μ, B) is then given by
$${\varOmega }_{0}=-kT{\int }_{-\infty }^{\infty }{N}_{0}\left(\varepsilon \right){\rm{ln}}\left(1+{{\rm{e}}}^{\left[\left(\mu -\varepsilon \right)/kT\right]}\right){\rm{d}}\varepsilon ,$$
where k is the Boltzmann constant, T is the temperature and μ is the chemical potential.
Now we consider LL broadening of width Γ (Dingle parameter) with a Lorentzian form
$$L\left(\varepsilon \right)=\frac{1}{\pi }\left(\frac{\varGamma }{{\varepsilon }^{2}+{\varGamma }^{2}}\right).$$
The DOS and the grand potential are then described by
$$N\left(\varepsilon \right)=\sum _{i}{D}_{i}L\left(\varepsilon -{E}_{i}\right),$$
$$\Omega =-kT{\int }_{-\infty }^{\infty }N(\varepsilon ){\rm{ln}}(1+{{\rm{e}}}^{[(\mu -\varepsilon )/kT]}){\rm{d}}\varepsilon =-kT{\int }_{-\infty }^{\infty }\sum _{i}{D}_{i}L(\varepsilon -{E}_{i}){\rm{ln}}(1+{{\rm{e}}}^{[(\mu -\varepsilon )/kT]}){\rm{d}}\varepsilon .$$
Then \(M\) is given by
$$M=-\frac{\partial \Omega }{\partial B}=kT{\int }_{-\infty }^{\infty }\sum _{i}(\frac{\partial {D}_{i}}{\partial B}L(\varepsilon -{E}_{i})-{D}_{i}{L}^{{\prime} }(\varepsilon -{E}_{i})\frac{\partial {E}_{i}}{\partial B}){\rm{ln}}(1+{{\rm{e}}}^{[(\mu -\varepsilon )/kT]}){\rm{d}}\varepsilon ,$$
where \({L}^{{\prime} }\left(\varepsilon \right)=\partial L\left(\varepsilon \right)/\partial \varepsilon \). In the zero-temperature limit (T → 0), M can be simplified:
$$M={\int }_{-\infty }^{\mu }\sum _{i}(\frac{\partial {D}_{i}}{\partial B}L(\varepsilon -{E}_{i})-{D}_{i}{L}^{{\prime} }(\varepsilon -{E}_{i})\frac{\partial {E}_{i}}{\partial B})(\mu -\varepsilon ){\rm{d}}\varepsilon .$$
To compare with our experiment, we need to calculate
$$\frac{\partial M}{\partial n}\left(n\right)=\frac{\partial M}{\partial \mu }\left(n\right)\frac{\partial \mu }{\partial n}\left(n\right),$$
where \(\frac{\partial \mu }{\partial n}\left(n\right)\) is the inverse of the DOS as a function of the carrier density and \(n(\mu )={\int }_{-\infty }^{\mu }N(\varepsilon ){\rm{d}}\varepsilon \).
Extended Data Fig. 5 shows the calculated n(μ), \(\frac{\partial n}{\partial \mu }\left(\mu \right)\), \(\frac{\partial n}{\partial \mu }\left(n\right)\), \(\frac{\partial M}{\partial \mu }\left(\mu \right)\), \(\frac{\partial M}{\partial \mu }\left(n\right)\) and \(\frac{\partial M}{\partial n}\left(n\right)\) versus Δ1 at Ba = 320 mT using the derived SWMc parameters and Dingle broadening Γ = 0.3 meV. The modulation in DOS, ∂n/∂μ, is well resolved in Extended Data Fig. 5b,c, but it is relatively small because of the LL broadening, except near CNP, in which large gaps with vanishing DOS open between the lowest LLs in the gullies at elevated Δ1.
The calculated ∂M/∂μ versus μ in Extended Data Fig. 5d shows that the crossing of the MLG and BLG LLs does not cause any phase shift. By contrast, in ∂M/∂μ versus n in Extended Data Fig. 5e, the BLG LLs show a 2π shift on crossing the four-fold degenerate MLG LLs and a π shift on crossing the two-fold degenerate zeroth LLs. This arises from the fact that the filling an MLG LL delays filling the next BLG LL versus total n, but not versus μ. As the DOS modulation \(\frac{\partial n}{\partial \mu }\left(n\right)\) is quite small, \(\frac{\partial M}{\partial n}\left(n\right)\) in Extended Data Fig. 5f looks very similar to \(\frac{\partial M}{\partial \mu }\left(n\right)\) except near CNP.
Derivation of the Dingle parameter
The energy bands are broadened by the intrinsic broadening Γ, given by the quantum scattering time τq = ħ/2Γ. Hence Γ sets the finest meaningful energy resolution with which the band energy can be described. To attain this energy resolution experimentally, we need to use the lowest Ba for which the LL energy gaps are comparable to Γ. In this limit, the amplitude of the QOs is rapidly suppressed with increasing Γ. Extended Data Fig. 4c–g shows the calculated QOs for various Dingle parameters Γ = 0.2–0.8 meV. As the energy spacing of the BLG LLs in the conduction band is about 1 meV, the amplitude of their QOs is suppressed by about two orders of magnitude over this range of Γ, whereas in the valence band, in which the LL gaps are about 0.6 meV, the QOs are completely quenched with the higher Γ. By contrast, the amplitude of the QOs of the MLG LLs, which have an order of magnitude larger gaps at low carrier densities, is much less affected by these Γ values. As a result, the relative amplitude of the MLG and BLG QOs is strongly dependent on Γ, enabling its accurate determination. By fitting to the experimental data in Fig. 2, we obtain Γ = 0.3 ± 0.05 meV, which also provides a very good agreement in quantitative comparison between the amplitudes of the measured \({B}_{z}^{{\rm{a}}{\rm{c}}}\) and the calculated mz taking into account the 2D magnetization reconstruction.
The finite nac modulation by \({V}_{{\rm{b}}{\rm{g}}}^{{\rm{a}}{\rm{c}}}\) also causes a suppression of the apparent amplitude of the QOs. It can be shown that if the peak-to-peak value of the carrier-density modulation is less than half of the LL degeneracy, nac < 2Ba/ϕ0, which is the case in our high-resolution measurements, the suppression is less than a factor of π/2. For larger nac, the visibility is suppressed rapidly as shown in Extended Data Fig. 2d,e. In particular, in Figs. 3 and 4 we have intentionally used larger nac to suppress the QOs due to BLG LLs and to improve the signal-to-noise ratio for detections of the MLG LLs. As this type of suppression of the apparent amplitude of QOs is harder to simulate in our BS calculations, we have used Γ = 0.3 meV for the calculations presented in all the figures except in Figs. 3 and 4, where Γ = 0.6 meV was used instead for suppression of the BLG QOs artificially. This larger Γ does not affect the shape of the calculated MLG QOs appreciably but reduces their amplitude.
Our derived Γ = 0.3 meV with corresponding local quantum scattering time τq = ħ/2Γ ≈ 1 ps, is about four times lower than the value reported based on global SdH oscillations40. This is consistent with the observation that the lowest magnetic field for detection of QOs in our local dHvA measurements is substantially lower than what is required for detection of the SdH oscillations (Extended Data Fig. 1). The large Γ reported based on SdH oscillations is probably because of sample inhomogeneity, such as charge disorder and the PMFs (BS). Hence, the measurement of the local dHvA QOs enables the determination of the local BS with energy resolution set by the intrinsic broadening Γ of the energy bands. This is of key importance for the study of BS of twisted vdW materials that are particularly prone to strain and spatial inhomogeneities.
LL anticrossings
The hybridization between the BLG and MLG bands on increasing Δ1 with the displacement field gives rise to partial lifting of valley degeneracy of the LLs. This effect is particularly pronounced near the top of the BLG valence band at intermediate values of Δ1 as shown in Extended Data Fig. 6c,d. Here, when MLG and BLG LLs in the same valley intersect, the strong band hybridization and non-vanishing γ3 leads to avoided crossing between the LLs as marked by the open symbols. Interestingly, the anticrossing occurs between the MLG LLs and every third BLG LL. Our derived SWMc parameters provide an excellent fit to the experimentally observed anticrossings as demonstrated in Extended Data Fig. 6a,b. Moreover, the strong hybridization lifts the valley degeneracy of the first MLG LL in the M3 sector as shown by the pronounced splitting between \({-1}_{{\rm{M3}}}^{-}\) and \({-1}_{{\rm{M3}}}^{+}\) in Extended Data Fig. 6b–d. This splitting is resolved experimentally in Extended Data Fig. 6a.
Interference of BLG LLs
The interference of the LLs can be observed also in the BLG bands at the same locations at which it is present in the MLG bands. Extended Data Fig. 7 shows the QOs acquired at site B as in Figs. 3h and 4e, but using lower \({V}_{{\rm{b}}{\rm{g}}}^{{\rm{a}}{\rm{c}}}=8\,{\rm{m}}{\rm{V}}\) rms that enables resolving the BLG LLs. The beating nodes at around 0.5 × 1012 cm−2and 1.8 × 1012 cm−2 are seen (Extended Data Fig. 7b), which can be well reproduced by the simulations using BS = 4.2 mT (Extended Data Fig. 7c).
Resolution of the PMF by LL interference
The minimal PMF that can be measured using the interference method is determined by the highest accessible LL index of the beating node \({N}_{{\rm{b}}}^{1}\). At Ba = 320 mT in the accessible range of n, the highest MLG LL index in ABA graphene is ±70, and hence the minimal \({B}_{{\rm{S}}}={B}_{{\rm{a}}}/\left(4{N}_{{\rm{b}}}^{1}\right)=1.14\,{\rm{mT}}\). For comparison, the lowest PMF that has been recently resolved by scanning tunnelling microscope is BS ≈ 0.5 T (ref. 66).
PMFs on different length scales
In moiré 2D materials, notable lattice relaxation occurs, giving rise to periodic strain and PMFs up to tens of tesla within moiré unit cell67,68,69. This short-range periodic PMF is part of the periodic potential that determines the BS70,71, but does not affect the usual LLs. By contrast, the strain that we probe varies gradually on a much larger length scale (about 1 µm). This strain gives rise to smooth PMFs, which shift the LLs in the presence of Ba and form strain-induced LLs at zero magnetic field6,72,73,74,75.
Towards characterization and use of PMFs
Strain engineering has been proposed to realize programmable PMFs leading to topological phases and various electronic devices5,48. Although large, short-range PMFs have been widely observed6,7,67,68,69,72,73,74,75, long-range homogeneous and controllable PMFs required for the development of new functionalities and valleytronics have not been realized5,43,44. Several methods have been proposed to induce variable mesoscale strain, including bending, MEMS, piezoelectric devices and polyimide deformation47,76,77,78,79, but the generated PMFs could not be detected. Our method enables the integration of such in situ controllable strain engineering, transport measurements and high-resolution local PMF imaging, laying the groundwork for investigation and use of PMFs.
Discussion of possible alternative mechanisms of interference of QOs
We consider below several other possible mechanisms that can alter the BS and induce degeneracy lifting, which may lead to interference of the LLs, and show that they are incompatible with the experimental data.
Band shifting
Spin–orbit coupling as well as the Zeeman effect at elevated fields can lift flavour degeneracy producing an energy shift between the bands of opposite spin or valley. Both the intrinsic spin–orbit coupling in graphene and the Zeeman contributions at our low magnetic fields result in a negligible energy shift of the order of µeV (refs. 80,81), which cannot account for the experimental data. Nevertheless, we explore whether a generic rigid shift between bands can reproduce the revealed LL interference pattern. In Fig. 3h, the first node in the interference of the MLG LLs occurs at an index N ≈ 19. The corresponding LL energy gap is \(\triangle {E}_{N}={E}_{N+1}-{E}_{N}=\sqrt{2e\hbar {v}_{{\rm{F}}}^{2}{B}_{{\rm{a}}}}\left(\sqrt{N+1}-\sqrt{N}\right)\approx 2.5\,{\rm{meV}}\). For the destructive interference, the LLs of the two bands have to be out of phase, namely, shifted by δEN ≈ 1.25 meV. Extended Data Fig. 8a shows the BS with a rigid shift of 1.25 meV between the K+ and K− bands with the corresponding calculated QOs presented in Extended Data Fig. 8b. The main resulting feature is that the MLG LLs are split into two, which is markedly different from the experimental QOs. This points out that to reproduce the observed QOs, the energy shift δEN between the interfering LLs has to grow with the LL index rather than being constant or decreasing with N. This is the behaviour in the case of PMF, where \(\delta {E}_{N}=\sqrt{2e\hbar {v}_{{\rm{F}}}^{2}N}\left(\sqrt{{B}_{{\rm{a}}}+{B}_{{\rm{S}}}}-\sqrt{{B}_{{\rm{a}}}-{B}_{{\rm{S}}}}\right)\approx {B}_{{\rm{S}}}\sqrt{2e\hbar {v}_{{\rm{F}}}^{2}N/{B}_{{\rm{a}}}}\) grows as \(\sqrt{N}\).
Staggered substrate potential
The possible alignment between the hBN and ABA graphene can cause an on-site potential difference between the A and B sublattices. Here we consider the simplest situation in which one of the graphene layers (bottom) is aligned with the hBN giving rise to a staggered substrate potential. In this case, the Hamiltonian can be written on the basis of {A1, B1, A2, B2, A3, B3} as
$${H}_{0}=\left(\begin{array}{cccccc}{\varDelta }_{1}+{\varDelta }_{2} & {v}_{0}{\pi }^{\dagger } & {v}_{4}{\pi }^{\dagger } & {v}_{3}\pi & {\gamma }_{2}/2 & 0\\ {v}_{0}\pi & \delta +{\varDelta }_{1}+{\varDelta }_{2} & {\gamma }_{1} & {v}_{4}{\pi }^{\dagger } & 0 & {\gamma }_{5}/2\\ {v}_{4}\pi & {\gamma }_{1} & \delta -2{\varDelta }_{2} & {v}_{0}{\pi }^{\dagger } & {v}_{4}\pi & {\gamma }_{1}\\ {v}_{3}{\pi }^{\dagger } & {v}_{4}\pi & {v}_{0}\pi & -2{\varDelta }_{2} & {v}_{3}{\pi }^{\dagger } & {v}_{4}\pi \\ {\gamma }_{2}/2 & 0 & {v}_{4}{\pi }^{\dagger } & {v}_{3}\pi & -{\varDelta }_{1}+{\varDelta }_{2}+{\delta }_{A3} & {v}_{0}{\pi }^{\dagger }\\ 0 & {\gamma }_{5}/2 & {\gamma }_{1} & {v}_{4}{\pi }^{\dagger } & {v}_{0}\pi & \delta -{\varDelta }_{1}+{\varDelta }_{2}+{\delta }_{B3}\end{array}\right).$$
For concreteness, we choose δA3 = 2 meV and δB3 = −2 meV. The resulting BS is shown in Extended Data Fig. 8c (red) in comparison with the original BS (black). The staggered substrate potential increases the gaps of the MLG and BLG bands but does not lift the valley degeneracy and therefore does not lead to beating. Extended Data Fig. 8d presents the calculated QOs showing no LL beating.
Kekulé distortion
Kekulé distortions are the bond density waves that have been observed in graphene epitaxially grown on copper82 or in the presence of strain83. In contrast to the O-type Kekulé distortion that opens a gap at the Dirac point, we find that the Y-type84 distortion can result in LL interference. The Y-shaped modulation of the bond strength, parametrized by the hopping parameters γ0 and \({\gamma }_{0}^{{\prime} }\) (Extended Data Fig. 8e), gives rise to valley-momentum locking and to inequivalent Fermi velocities for both the MLG and BLG bands. Hence, it lifts the valley degeneracy of the LLs resulting in chiral symmetry breaking. In the SWMc model, γ0 is the sole parameter that controls the Fermi velocity vF of the MLG band (\({v}_{{\rm{F}}}=\frac{\sqrt{3}}{2}\frac{a{\gamma }_{0}}{\hbar }\)). The energy difference between the LLs from the two valleys with the same index N is \(\delta {E}_{N}=\sqrt{2e\hbar N{B}_{{\rm{a}}}}{\Delta v}_{{\rm{F}}}\), where \({\Delta v}_{{\rm{F}}}=\frac{\sqrt{3}}{2}\frac{a}{\hbar }({\gamma }_{0}-{\gamma }_{0}^{{\prime} })\). The first beating node appears when δEn is equal to half of the gap size: \(\sqrt{2e\hbar N{B}_{{\rm{a}}}}{\Delta v}_{{\rm{F}}}\,=\)\(\sqrt{e\hbar {v}_{{\rm{F}}}^{2}{B}_{{\rm{a}}}/2N}/2\), which yields \({N}_{{\rm{b}}}^{1}={v}_{{\rm{F}}}/\left(4{\Delta v}_{{\rm{F}}}\right)\) as shown in Extended Data Fig. 8f. In Fig. 3h, \({N}_{{\rm{b}}}^{1}=19\), which corresponds to a very weak Kekulé distortion with ΔvF/vF = 1.4 × 10−2. However, the Kekulé distortion results in \({N}_{{\rm{b}}}^{1}\) that is independent of Ba as corroborated by the calculated QOs for Ba = 320 mT and 170 mT in Extended Data Fig. 8h,i. This is because the LLs shift in the same proportion in the two valleys with Ba. This is in sharp contrast to beating due to PMF for which \({N}_{{\rm{b}}}^{1}={B}_{{\rm{a}}}/\left(4{B}_{{\rm{S}}}\right)\) is proportional to Ba. The experimental data points in Extended Data Fig. 8g (circles) are consistent with PMF and incompatible with the Kekulé distortion.
Disorder in BS parameters
The BS can vary in space because of various types of disorder. Focusing on the Dirac bands, for example, the energy of the Dirac point or vF could be position dependent without breaking the valley symmetry. If the parameters change gradually in space on lengths scale larger than our spatial resolution of about 150 nm, the LLs will shift gradually in space following the variations in the BS without showing interference at any location. Let us now consider the opposite case of sharp boundaries between domains with different BS. In this situation, at the boundaries, the finite size of our SOT may result in the simultaneous detection of LLs originating from the two neighbouring domains giving rise to apparent interference. In such a case, we expect to observe interference along a network of grain boundaries with width comparable to our SOT size. Instead, Fig. 3e shows well-defined domains of typical width of 1 µm and length of up to 2 µm, much larger than the SOT size, over which the interference is rather uniform. Furthermore, most of the domains showing beating are located at the ends or corners of the device, so they do not have two neighbouring domains that can cause the apparent interference. Finally, if there is a relative shift in the Dirac point between the neighbouring domains, the apparent interference patterns at the boundary would evolve similar to that calculated in Extended Data Fig. 8b, whereas if vF changes between the domains the beating node \({N}_{{\rm{b}}}^{1}\) of the apparent interference would be independent of Ba as calculated in Extended Data Figs. 8f–i. Both these possibilities are inconsistent with the experimental data. More generally, the Ba dependence of the LL interference due to variations in BS is distinctly different from the one caused by BS. We therefore conclude that disorder that causes spatial variations in BS without creating PMFs cannot explain the observed LL interference.