Imaging quantum oscillations and millitesla pseudomagnetic fields in graphene

The exceptional control of the electronic energy bands in atomically thin quantum materials has led to the discovery of several emergent phenomena1. However, at present there is no versatile method for mapping the local band structure in advanced two-dimensional materials devices in which the active layer is commonly embedded in the insulating layers and metallic gates. Using a scanning superconducting quantum interference device, here we image the de Haas–van Alphen quantum oscillations in a model system, the Bernal-stacked trilayer graphene with dual gates, which shows several highly tunable bands2–4. By resolving thermodynamic quantum oscillations spanning more than 100 Landau levels in low magnetic fields, we reconstruct the band structure and its evolution with the displacement field with excellent precision and nanoscale spatial resolution. Moreover, by developing Landau-level interferometry, we show shear-strain-induced pseudomagnetic fields and map their spatial dependence. In contrast to artificially induced large strain, which leads to pseudomagnetic fields of hundreds of tesla5–7, we detect naturally occurring pseudomagnetic fields as low as 1 mT corresponding to graphene twisting by 1 millidegree, two orders of magnitude lower than the typical angle disorder in twisted bilayer graphene8–11. This ability to resolve the local band structure and strain at the nanoscale level enables the characterization and use of tunable band engineering in practical van der Waals devices.

Determining the band structure (BS) and the Fermi surface is a crucial step in understanding and utilizing the electronic properties of materials.The most sensitive canonic method for mapping the BS of bulk metals and semiconductors is measurement of the de Haas-van Alphen (dHvA) oscillations (12).In this quantum mechanical effect, in the presence of magnetic field , electrons coherently circulate in closed electronic orbits along the Fermi surface, giving rise to quantum oscillations (QOs) in the grand thermodynamic potential Ω and in the associated magnetization  = − ⁄  (12).In 2D systems these oscillations are described by the formation of discrete set of Landau energy levels (LLs) with sharp peaks in the density of states (DOS).Charge carriers orbiting in the metallic LL states give rise to diamagnetic response, while groundstate currents flowing in the gapped edge states contribute to paramagnetic magnetization, giving rise to magnetization oscillations with either magnetic field or carrier density (12).Since the measured total magnetic moment scales with sample volume, observation of dHvA effect in 2D systems has been challenging (13,14).Instead, non-thermodynamic Shubnikov-de Haas (SdH) oscillations in the charge transport coefficients is the benchmark characterization tool for atomic layer materials (15).
The recent advances in fabrication of multilayer stacked and twisted atomic layer devices based on van der Waals (vdW) materials have provided a fertile ground for a plethora of novel strongly correlated electronic phases which are rare in bulk materials, including tunable correlated insulators (16), orbital magnetism (17)(18)(19), integer and fractional Chern insulators (20)(21)(22)(23), unconventional superconductivity (24,25), and excitonic superfluidity (26,27).Through material selection, stacking order, and twist angle, a vast variety of vdW structures with different properties can be created.Their BS is strongly dependent on the stacking geometry and can be further manipulated by varying transverse electric field, magnetic field, strain, or pressure.This has ushered in a new era in materials research and offers numerous opportunities for applications in nextgeneration electronics.
The experimental investigation of the structure of the Fermi surface in micron-sized vdW devices is currently based predominantly on the detection of quantum oscillations using SdH effect (15,28) and capacitance oscillations (4,25).However, various types of disorder, such as charge inhomogeneity, twist angle disorder, and strain, widely exist in these samples (9,29,30) resulting in spatial variations in the BS, while the aforementioned methods lack spatial information.While several scanning probe techniques, including scanning tunneling microscopy (31,32) and single electron transistors (22,33), are powerful probes of the local electronic properties, the former requires the electron layers to be exposed to vacuum as in photoemission studies, and neither of them can be used for devices encapsulated with a metallic top gate required for application of displacement fields.Development of a universal tool for measurement of the local BS in the diverse family of 2D quantum materials is thus highly desirable.
Strain is a particularly interesting but difficult to characterize tunable parameter, and it should be ubiquitous in vdW devices due to their high mechanical flexibility.In addition to changing the BS and breaking of crystal symmetries, it creates pseudomagnetic fields (PMFs) due to the valley degree of freedom in hexagonal vdW materials (34).If sufficiently large, these gauge fields may form effective LLs with sharp peaks in DOS that strongly vary in space and profoundly alter the transport properties of the electronic devices (5).PMFs of tens to hundreds of tesla have been observed in artificially strained nanostructures of monolayer graphene (6,7), bilayer graphene (35), twisted homo- (36) and hetero-structures (37), and folding-formed trilayer graphene (38).Yet, the PMFs due to strain naturally forming during the fabrication process have not been observed so far.
Utilizing scanning SQUID-on-tip (SOT) microscopy (39), we imaged the dHvA effect in hBN encapsulated dualgated Bernal-stacked trilayer graphene (TLG) (Fig. 1a).The high magnetic sensitivity of the SOT allows imaging of the QOs at low fields resolving the multi-band electronic structure with sub-meV energy resolution and unperturbed by elevated magnetic fields.The quantitative information provided by the thermodynamic oscillations allows unique high-precision derivation of the tight-binding hopping parameters and accurate reconstruction of the tunable band hybridization induced by the displacement field.Moreover, the nanoscale spatial resolution allows detailed quantitative study of the spatial variations of the QOs over the entire device, revealing the presence of PMFs of mT magnitude in micron-sized domains.

Band structure of ABA graphene
The Bernal-stacked ABA TLG is the minimal graphene structure that requires the full set of parameters in the Slonczewski-Weiss-McClure (SWMc) tight-binding model (40) proposed for describing the BS of graphite (Fig. 1e) with six hopping parameters   (  = 0 to 5), on-site energy difference  between the A and B sublattices, potential difference ∆ 1 between the top and bottom graphene layers induced by the applied displacement field  , and ∆ 2 potential that describes the non-uniform charge distribution between the middle and outer layers.The influence of these individual parameters on the BS is shown in Extended Data Fig. 5.
Due to the mirror symmetry of the crystal structure, the bands can be decomposed into a monolayer graphene-like (MLG) Dirac band and a bilayer graphene-like (BLG) quadratic band, with a relative energy shift between them (Fig. 1f left).A transverse displacement field breaks the mirror symmetry and induces strong hybridization between the two bands, causing multiple Lifshitz transitions.At high  , the Dirac band is divided into three sections, which we label M1, M2 and M3 (Fig. 1f right).The M1 and M3 bands remain separated from the BLG-like bands B1 and B2, while M2 merges with B1 and evolves into gapped mini-Dirac cones (Dirac gullies) at the bottom of B1.The gullies have a three-fold rotational symmetry leading to various possible quantum Hall ferromagnetic and nematic states (41,42).
Previous studies have explored SdH and capacitance oscillations in TLG at elevated fields (3,43,44) to determine the BS and search for broken-symmetry states (4,45,46), resulting in a wide span of derived SWMc parameters (Table 1 in Methods).However, the details of the ABA graphene BS are still under debate.In particular, the size and the sign of the small gap in the Dirac band at  = 0 (Fig. 1f inset), is controversial (46).In addition, it was predicted that the trigonal warping due to non-vanishing  3 breaks the rotational symmetry of the BS with significant consequences on the LL structure, resulting in level anticrossings, which occur between a given MLG LL and every third BLG LL (42,47).Though some single anticrossings were reported (46), the predicted periodicity has not been observed directly.Moreover, some gully-polarized states have been reported in high magnetic fields (4), but gully coherence at low fields remains an open question.
To address these questions, a technique capable of probing the BS with high energy resolution and at low magnetic fields is required.The unique advantage of the dHvA effect is that the QOs are fully described by thermodynamic quantities which can be readily derived directly from BS calculations, thus allowing detailed quantitative analysis, which is not possible in other methods like SdH oscillations.

Nanoscale magnetic imaging and results
The TLG device (Fig. 1c) was fabricated using the common dry transfer method, with the graphene flake sandwiched between two Pt gates, separated by ≅30 nm of hBN (Methods), as shown schematically in Fig. 1a.Transport measurements of   vs. the carrier density  and the applied magnetic field   show a Landau fan with several LL crossings (Extended Data Fig. 1), similar to previous reports (46).The local measurements of the dHvA oscillations were performed in a home-built cryogen-free scanning SOT system at  ≅ 160 mK (Methods).The Indium SOT (Fig. 1b) had an effective diameter of 150 nm and a magnetic sensitivity of 20 nT/Hz We first describe the QOs at a single point by positioning the SOT statically at a height of ℎ ≅ 150 nm above the graphene.Figure 2c shows the dHvA oscillations acquired at   = 320 mT as a function of  and , focusing on low carrier densities between −1.2×10 12 and 2.3×10 12 cm -2 .A line cut of the data at  = 0 V/nm is shown in Fig. 2a.Remarkably, in this relatively small range of  we observe over 100 LLs, in sharp contrast to transport measurements (Extended Data Fig. 1) where almost no SdH oscillations can be discerned at such low   .Moreover, we can resolve dHvA oscillations at fields as low as 40 mT as shown in Extended Data Fig. 3. To the best of our knowledge, this is the lowest magnetic field at which QOs have been observed in 2D systems.As illustrated in Fig. 1g, mapping the very dense LLs provides a unique quantitative tool for reconstruction of the BS with unprecedented precision, allowing us to derive high accuracy SWMc parameters as described in Methods and summarized in Table 1.,e).Moreover, higher M1 LLs display a pronounced negative (dark) signal due to the diamagnetic response in the compressible states (48).Remarkably, as opposed to all other LLs, the 0 1 + and 0 2 − LLs in Figs.2c-e are invisible showing no diamagnetic signal, and their presence is discerned only as a phase shift in the B1 LLs.This is the result of the Dirac Berry curvature giving rise to 0 th LLs pinned to the band extrema with no kinetic energy and no magnetic field dependence, and hence no diamagnetism in their compressible states.In contrast, the incompressible states in the gaps of MLG LLs do show a paramagnetic response (48) determined by the gap Chern number  as shown in Extended Data Fig. 3c.
M2 band.The hybridization of M2 and B1 bands results in a small M2 hole pocket with total DOS that is too low to accommodate even a single LL at elevated magnetic fields.As a result, the M2 LLs could not be identified in previous studies (4,43,46).Our low   and high sensitivity allow clear resolution of M2 LLs as shown in Figs.2c-e and Extended Data Fig. 3.In contrast to M1 and M3 LLs that coexist with B1 and B2 bands with the same carrier type, the M2 hole LLs coexist with B1 electron LLs, reflecting a Mexican hat potential.
Figure 2e shows a gap   0 between the 0 1 − and 0 2 + LLs at  = 0, that is comparable to the gap between two B1 LLs (≅ 1 meV).Upon increasing , the 0 1 − and 0 2 + LLs spread apart rapidly without crossing, clearly indicating that   grows continuously without intermediate gap closing, in contrast to previous suggestions (4,46).

M3 band.
At elevated hole doping the M3 LLs are four-fold degenerate, showing behavior similar to the M1 LLs.At low doping, in contrast, the strong hybridization between M3 and B2 bands gives rise to an unusual valley polarization.This is demonstrated in Fig. 2f by the four-fold degenerate −1 3 LL, which splits into valley polarized, spin-degenerate −1 3 + and −1 3 − LLs at low hole doping.This valley splitting process is accompanied by multiple M3 and B2 LLs crossings and anticrossings as described in Extended Data Fig. 8.In particular, BS calculations have predicted anticrossings between BLG LLs and low-index M3 LLs arising from the trigonal warping term  3 (42,47), giving rise to avoided crossings with enhanced gaps for every third BLG LL.This tripled period of anticrossings with BLG LLs, which has remained unidentified so far, is clearly resolved in our data marked by the white bars in Fig. 2f in excellent agreement with the calculated results in Extended Data Fig. 8.In addition, we resolve also the 0 3 − LL which has no diamagnetism and gives rise to a  shift in the B2 LLs as shown by the red dotted line in Fig. 2f.
LLs in the gullies.Upon increasing , the enhanced hybridization between the MLG and BLG bands and the trigonal warping give rise to formation of three-fold rotationally symmetric Dirac gullies (42,47)  The intricate structure of the QOs is highly sensitive to the shape of the BS and its evolution with .Our fine energy resolution thus allows high precision determination of the SWMc parameters and of the Dingle broadening of the LLs as described in Methods.

LL interferometry and pseudomagnetic fields
Next, we analyze the QOs over the full range of accessible carrier densities || ≲ 9×10 12 cm -2 , which allows resolution of much finer details of the BS and its spatial dependence.Since in this range of  at   = 320 mT there are about 500 BLG LLs in addition to over 100 MLG LLs, we focus on the sparser MLG LLs by applying a larger    excitation, which averages out the BLG QOs and enhances the visibility of the MLG LLs (Methods).Figures 3a-d show the spatial dependence of    (, ) at several carrier densities , while Fig. 3f presents    () vs.  along the dotted line in Fig. 3a revealing QOs due to the MLG LLs.For || ≲ 3×10 12 cm -2 , discussed in Fig. 2, the QOs display relatively high spatial uniformity as shown in Fig. 3d and at the bottom of Fig. 3f, demonstrating high sample quality.At higher carrier densities, however, distinctly different behavior is observed depending on the location as demonstrated in Figs.3g,h showing the QOs at sites A and B indicated in Fig. 3b.Large parts of the sample, exemplified by site A, show continuous evolution of the amplitude of the QOs with  (Fig. 3g), which matches very well the calculated magnetization oscillations considering the finite Dingle broadening and the finite    (Methods).In other parts of the sample, however, striking low-frequency beating of the MLG LLs is found as shown in Fig. 3h at site B. Upon decreasing the field to 170 mT, more beating nodes are observed and their position is shifted to lower LL indices (Fig. 3i).
The low-frequency beating of MLG QOs indicates interference of two close oscillation frequencies, usually originating from two bands with very similar dispersion relations.Since the MLG and BLG bands have very different dispersions, the beating cannot arise due to interference between these two bands.It must therefore originate from small symmetry breaking between the four flavors of the MLG band.In Methods we consider various possible mechanisms of such degeneracy lifting including staggered substrate potential, Kekulé distortions, band shifting, Zeeman effects, and spin-orbit coupling, as well as non-symmetry-breaking disorder, and show that they are inconsistent with the observed behavior.In the following, we demonstrate that the interference of the QOs is well described by the natural presence of strain-induced pseudomagnetic fields   .
Long-wavelength mechanical strain induces an effective gauge field in graphene, which has an opposite sign for the two valleys.Isotropic or uniaxial strains result in zero PMF, whereas non-uniform shear strain gives rise to finite   (34).In the presence of magnetic field, the carriers in the  + and  − valleys experience effective fields   of   +   and   −   which causes a relative shift in the LLs leading to interference as illustrated in Fig. 4a  (black curve) show excellent agreement with the data and in addition well reproduces the secondary nodes at −19 and 57.Moreover, the experimental evolution of the interference of the QOs with  (Fig. 4e) is well reproduced by the simulations (Fig. 4f) both in the position of all the three beating nodes (red arrows) and in their evolution with  using a single fitting parameter   .Since the BS changes profoundly with  while the strain-induced PMF should be independent of the BS, the fact that the observed beating is well described by a -independent   provides an additional strong support of the model.Furthermore, the PMF should affect all the bands.An important self-consistency check is therefore observation of beating also in the BLG LLs.Extended Data Fig. 9 shows a closer examination of the BLG QOs acquired at the same location B. Beating is clearly observed and well reproduced by simulations using the same set of parameters as in Figs. 3 and 4. Finally, a crucial test of the PMF origin of the interference is the predicted linear dependence of   1 on   , which sets it apart from other possible mechanisms (Methods).Figure 3i (red curve) presents the LL interference measured at a lower field of 170 mT showing multiple beating nodes.The calculated QOs (black curve) show an excellent agreement with the data with the same   as derived at 320 mT, confirming the linear   1 dependence on   as shown in Fig. 4c.
By analyzing the LL interference over the entire sample, we derive a map of   as shown in Fig. 3e.In large parts of the sample,   is below our resolution of about 1 mT determined by the largest accessible carrier density.We find three regions with characteristic length  ≅ 1 µm with smoothly varying   reaching up to 6 mT.Two types of lattice distortions-finely tuned triaxial strain or arc-like in-plane bending-have been shown theoretically to produce relatively homogenous   (5,49).In our sample it is very likely that the source of the strain is a small arc-like in-plane twist introduced during the fabrication processes (Fig. 4b inset).
An arc segment of length  with a small twist angle  in a graphene strip generates   =   0   (49), where  = 0.14 nm is the graphene interatomic distance,  ≅ 1 is a numerical constant,  ≅ 2 describes the hopping parameter dependence on ,  0 = ℎ/ is the flux quantum, ℎ is the Planck constant, and  is the elementary charge.The measured 1 <   < 6 mT thus corresponds to twisting by 1 <  < 6 millidegree, or equivalently bending radius of 1 <  < 6 cm and corresponding strain of 8•10 -6 <  � < 5×10 -5 , where  = / and  � = /2 (49).Such minute bending angles and strains should be abundant in exfoliated atomic layer devices.In fact, orders of magnitude larger strains and angle disorder have been reported to occur naturally in twisted and stacked graphene structures (8-11, 32, 50).

Discussion
Our measurements of dHvA effect in low fields have revealed multiple fine features of the BS in a trilayer graphene device, which allow high-precision derivation of the hopping parameters in the SWMc model, while the nanoscale spatial resolution enables mapping the local BS over the entire device.The revealed PMFs in a honeycomb lattice, cause beating of the LLs, which is usually ascribed to sub-band splitting or spin-orbit coupling.The presence of PMFs has important implications for our comprehension of various types of disorder and their impact on strongly correlated states in a wide range of vdW structures.In particular, magic angle twisted bilayer and multilayer graphene are known to be prone to twist angle disorder which strongly affects their properties.Spatial variations in twist angle and strain are understood to cause fluctuations in the bandwidth of the flat bands, in electron interactions, and in formation of symmetry broken states (51).Yet, the effects of the accompanying spatially varying PMF have not been investigated systematically.Our findings imply that the typical reported twist angle disorder of 0.1° (8)(9)(10)(11) will generate   ≅ 0.1 T, which should greatly affect the magnetotransport behavior.Indeed, resolving LLs in transport at low fields has proven to be challenging in twisted devices, which could be attributed to such highly spatially varying PMFs.Moreover, in the presence of magnetic field, the PMF breaks the valley symmetry resulting in different DOS in the two valleys.  fluctuations may thus affect the local Stoner instabilities and symmetry breaking mechanisms that lead to the quantum anomalous Hall effect, Chern insulators, and inhomogeneities in the spontaneous orbital magnetization (18)(19)(20).Finally, the recent development of controllable in-plane bending of graphene ribbons (52) opens the door to microscale engineering and exploitation of PMFs.The derived method of high precision determination of the local band structure and imaging of the pseudomagnetic fields provides a powerful tool for characterization and optimization of tunable electronic bands and calls for further investigation of the role of strain-induced gauge fields in the formation of symmetry broken strongly correlated states of matter.
32.8 kHz (Model TB38, HMI Frequency Technology), which is used as a force sensor for tip height control (66).
The scanning height was about 150 nm above the ABA graphene.An ac voltage    at frequency of about

Magnetic field and modulation amplitude dependence of QOs
The measured signal    =   (  /) is proportional to the modulation amplitude of the carrier density   induced by    .It is therefore desirable to use large   to improve the signal-to-noise ratio.In order to resolve QOs, however,   has to be substantially smaller than the period of the oscillations ∆ .Extended Data Figs.3a-c show the QOs acquired at   = 320 mT using    = 8, 35, and 100 mV rms corresponding to   of 5.19×10 9 cm -2 , 2.27×10 10 cm -2 , and 6.49×10 11 cm -2 rms.The four-fold degenerate BLG LLs have a period of ∆ = 4   0 ⁄ = 3.1×10 10 cm -2 .The lowest    = 8 mV rms was chosen to result in peak-to-peak value of   of 1.47×10 10 cm -2 , approximately equal to ∆ 2 ⁄ = 1.55×10 10 cm -2 , which results in optimal signal-to-noise ratio for detecting the BLG LLs, albeit suppresses the measured    /  ratio by a factor of /2 .A larger   washes out the QOs due to BLG LLs, leaving the MLG LLs clearly resolvable as demonstrated in Extended Data Figs.3b,c.The largest   also allows clear observation of the paramagnetic response   ⁄ = / 0 in the gap between the 0 th and 1 st MLG LLs dictated by the Chern number  = 2 on the electron side and  = −2 on the hole side (Extended Data Fig. 3c).Hamiltonian can be rewritten as:
In an external magnetic field, in the Landau gauge, the canonical momentum  can be replaced by  − , where  is the vector potential. obeys the commutation relation �  ,   � = −/  , where   = �(ℏ/) is the magnetic length.As in the usual one-dimensional harmonic oscillator, in the basis of LL orbital |⟩, the matrix elements of ,  † are given by: (3) Therefore, the new Hamiltonian can be written in the basis of LL orbitals.Using matrix elements of  and  † operators, the momentum operators are replaced by an up (down) diagonal matrix of dimensions  × , where  is the cutoff number for the infinite matrix, restricting the Hilbert space with indices  ≤ .All the other nonzero elements  are substituted by   , where   is the identity matrix with dimensions  × .
As our measurements were performed in low magnetic fields and high-index LLs are often involved, large cutoff was used so that it spans 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 range (Fig. 2) and to 800 for calculations over larger range (Fig. 3,4).

Evolution of the band structure and Landau levels with displacement field
Extended Data Fig. 4 shows the calculated BS of ABA graphene using the derived SWMc parameters and the evolution of the LLs with  and   .At  = 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 0 th LLs of the MLG and BLG bands which are valley polarized due to the Berry curvature (Extended Data Fig. 4a).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 Figs.4b-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 Figs.4e,f.Since the BLG bandgap ∆  0 is characterized by  = 0, it has no magnetization.The six-fold degenerated compressible 0 th LLs 0  + and 0  − in the gullies also have no magnetization at low fields,  = −   ⁄ = 0, due to their zero kinetic energy.As a result, zero magnetization is observed around the CNP over a width of  = 12   0 ⁄ in carrier density as indicated in Figs.2c,f.The first paramagnetic signal appears when the Fermi level reaches the  = ±6 gaps ∆  1 and ∆  −1 between the 0 th and the 1 st gully LLs as marked in Fig. 2f.At elevated magnetic fields the six-fold gully degeneracy of the 0 th LLs is partially lifted (4,57).

Reconstruction of band structure parameters
A number of experimental studies (3,4,43,44,46,55,56,58,68) have investigated the tight-binding parameters of ABA graphene as summarized in 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 3100 meV and fit the remaining seven parameters.Extended Data Fig. 5 shows the effect of the individual parameters on the BS.
In the absence of displacement field,  1 = 0, the MLG band is affected only by  0 ,  2 ,  5 and , with  Attribute (i) is determined by the DOS ratio of the two bands.For a given  0 , the relative DOS 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 attain  1 = 370 ± 10 meV.
Attribute (ii) is then used to determine  2 .The energies of the band extrema and hence the relative position of the 0 th LLs can be calculated analytically.In particular, from Eq. ( 2) for  1 = 0, the 0 1 − LL at the bottom of M1 band is positioned at energy  2 −  2 /2, while 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 .Since the LL spectrum is quite sensitive to  2 ,  2 is determined first.We use the relative position between −1 3 and the nearby BLG LLs to fit  2 , and we get  2 = −19 ± 0.5 meV.
Attribute (iii) is governed by  3 , which induces trigonal warping of the BLG bands.As shown in Extended Data Fig. 8, this results in the anticrossings between the BLG LLs and MLG 0 3 − and −1 3 + LLs.From fitting to the experimental data we attain  3 = 315 ± 10 meV.
Attributes (ii) and (iv) are used to derive  and  5 .From Eq. ( 2), the MLG band gap at  = 0 V/nm is   0 =  + ( 2 −  5 )/2, while the gap center is located at 2 2 +  − ( 2 +  5 )/2.In our experiment data, one BLG LL fits within the MLG gap and 20 BLG LLs reside between 0 1 − and −1 3 , from which we attain  = 18.5 ± 0.5 meV and γ 5 = 20 ± 0.5 meV.Note that   0 can be either positive or negative.We find that   0 is negative, which means that the 0 th  − LL (0 1 − ) resides at the bottom of M1 band and the 0 th  + LL (0 2 + ) is at the top of M2.In this case the Dirac gap   increases with  1 and the 0 1 − and the 0 2 + LLs spread apart with the displacement field as shown in Extended Data Fig. 4f, consistent with experimental data in Figs.
2c,e and calculations in Figs.2d.If   0 is positive, the 0 th  − LL will reside at the top of M2 while the 0 th  + LL will be at the bottom of M1.In this case, upon increasing , the Dirac gap closes and then reopens with crossing of the two 0 th LLs, such that   is always negative at high  with 0 th  − LL at the bottom of M1.Table 1 shows that the value and the sign of   0 varies significantly in the literature.Note, however, that only in Datta(46) and in the present work the Dirac gap is reported directly.For the rest of the references the   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 −1 3 resides closely to the BLG band gap, we use the number of BLG LLs between −1 3 and −2 3 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 1.By optimizing all parameters for best fit to the experimental data we derive  4 = 140 ± 15 meV, as summarized in Table 1.This work 3100 370 ± 10 -19 ± 0.5 315 ± 10 140 ± 15 20 ± 0.5 18.5 ± 0.5 3.8 ± 0.05 -1.0

Orbital magnetization calculations
Oscillations in orbital magnetization  due to LLs can be calculated analytically for either parabolic or Dirac bands as shown previously (69).However, there is no analytical expression for the LL spectrum in ABA graphene, and therefore the magnetization oscillations have to be calculated numerically.We follow the method described in (14) to derive the magnetization () and then calculate its derivative /.
We first consider the case with zero LL broadening.For an arbitrary LL spectrum   with degeneracy   ( is the Landau level index), the density of states  0 () of the system is: describes spin-degenerate LLs from both valleys with degeneracy   = 2  ℎ . The grand thermodynamic potential  0 (, ) is then given by: Here  is the Boltzmann constant,  is the temperature and  is the chemical potential.Now we consider LL broadening of width  (Dingle parameter) with a Lorentzian form The DOS and the grand potential are then described by Then  is given by where '() = ()/.In the zero-temperature limit (0),  can be simplified: To compare with our experiment, we need to calculate

Derivation of the Dingle parameter
At our low   the intrinsic LL broadening strongly suppresses the amplitude of the measured QOs.Extended Data Fig. 7 shows the calculated QOs for various Dingle parameters  = 0.2 to 0.8 meV.Since 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 , while in the valence band, where the LL gaps are about 0.6 meV, the QOs are completely quenched with the higher .In 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 , allowing its accurate determination.By fitting to the experimental data in Fig. 2, we attain  = 0.3 ± 0.05 meV, which also provides a very good agreement in quantitative comparison between the amplitudes of the measured    and the calculated   taking into account the 2D magnetization reconstruction.
Note that the finite   modulation by    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,   < 2  / 0 , which is the case in our high-resolution measurements, the suppression is less than a factor of /2.For larger   the visibility is suppressed rapidly as shown in Extended Data Fig. 3.In particular, in Figs. 3 and 4 we have intentionally used larger   in order to suppress the QOs due to BLG LLs and to improve the signal to noise ratio for detections of the MLG LLs.Since 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 for Figs. 3 and 4, where  = 0.6 meV was used instead for suppression of the BLG QOs artificially.This larger  does not affect appreciably the shape of the calculated MLG QOs but reduces their amplitude.
Our derived  = 0.3 meV with corresponding local quantum scattering time   = ℏ/2 ≅ 1 ps, is about four times lower than the value reported based on global SdH oscillations (46).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 likely due to sample inhomogeneity, such as charge disorder and the pseudomagnetic fields (  ).

LL anticrossings
The hybridization between the BLG and MLG bands upon 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   = �2ℏ  2   �√ + 1 − √� ≅ 2.5 meV.For the destructive interference, the LLs of the two bands have to be out of phase, namely shifted by   ≅ 1.25 meV.Extended Data Fig. 10a shows the BS with a rigid shift of 1.25 meV between the  + and  − bands with the corresponding calculated QOs presented in Extended Data Fig. 10b.The main resulting feature is that the MLG LLs are split into two which is markedly different from the experimental QOs.This clearly points out that in order to reproduce the observed QOs, the energy shift   between the interfering LLs has to grow with the LL index rather than being constant or decreasing with .This is indeed the behavior in the case of PMF, where   = �2ℏ Disorder in band structure parameters.The BS can vary in space due to various types of disorder.Focusing on the Dirac bands, for example, the energy of the Dirac point or   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 simultaneous detection of LLs originating from the two neighboring 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, Figure 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.In addition, most of the domains showing beating are located at the ends or corners of the device, so they do not have two neighboring domains that can cause the apparent interference.Finally, if there is a relative shift in the Dirac point between the neighboring domains, the apparent interference patterns at the boundary would evolve like calculated in Extended Data Fig. 10b, while if   changes between the domains the beating node   1 of the apparent interference would be independent of   as calculated in Extended Data Figs.10f-i.
Both these possibilities are inconsistent with the experimental data.More generally, the   dependence of the LL interference due to variations in BS is distinctly different from the one caused by   .We therefore conclude that disorder that causes spatial variations in BS without creating PMFs cannot explain the observed LL interference.

Fig. 1 .
Fig. 1.Experimental setup and ABA graphene band structure.a, Schematic sample structure with trilayer graphene (TLG) encapsulated by hBN and top (TG) and back (BG) Pt gates, with scanning SQUID-on-tip (SOT).b, Scanning electron microscope image of the Indium SOT.c, Optical image of the TLG device.d, The stacking geometry of the device with indicated voltages applied to the gates.e, Atomic structure of ABA graphene with indicated SWMc parameters.f, The BS of ABA TLG for zero displacement field (Δ 1 = 0 meV) and for Δ 1 = 50 meV.Inset: a small Dirac gap   0 is present in the MLG band at Δ 1 = 0 meV which grows rapidly with Δ 1 .g, 3D rendering of the BS reconstructed from the dHvA oscillations with overlaid contours of the calculated LLs.The LLs are shown for   = 1 T for clarity.At our   = 320 mT the LLs are three times denser.The color map represents the wavefunction projection onto the MLG (red) and BLG-like (blue) bands.

Fig. 2 .
Fig. 2. Measurement of the dHvA effect in ABA graphene.a, The measured local magnetic QOs signal    as a function of  at 160 mK at a fixed location in the interior of the sample at   = 320 mT and  = 0 V/nm using    = 8 mV rms.Some indices of the LLs in M1, B1, and M3 bands are indicated.b, The calculated dHvA differential magnetization   at  = 0 V/nm using the derived BS. c, The measured    vs.  and .Crossings between the four-fold degenerate BLG LLs (horizontal yellow/blue lines) and the 0 th  − valley MLG LL, 0 1 − results in a  shift (red dotted line), while crossing with the higher MLG LLs introduces a 2 shift (white dotted line).The white vertical bars indicate the 12-fold degeneracy of the LLs in the gullies.d, Calculated   (, ) using the fitted SWMc tight-binding parameters with Dingle broadening of 0.3 meV.The MLG LLs in M1, M2, and M3 bands and the gully LLs are labeled.e, Zoom-in of the measured    (, ) in the vicinity of MLG Dirac gap   with marked MLG LLs (red).The  shifts at crossings between BLG and 0 2 + and 0 1 − valley-polarized MLG LLs and 2 shift at crossing of valley degenerate −1 2 LL are indicated.f, Zoom-in of    (, ) near the top of the BLG valence band B2 showing anticrossings between the MLG and BLG LLs.The white horizontal bars indicate enlarged anticrossing gaps for every third BLG LL (see Extended Data Fig. 8).
and to highly intriguing evolution of the LLs.The low energy BLG LLs, which are mostly valley degenerate at  = 0 , undergo valley polarization and regrouping with , forming valley-polarized six-fold degenerate LLs in the gully pockets as shown in Extended Data Fig. 4f.The 0 th electron and hole gully LLs, 0  − and 0  + , display no diamagnetism due to the Berry phase, and the gap between them, ∆  0 has zero Chern number  = 0.As a result, a zero magnetism strip of width corresponding to 12-fold degeneracy (incompressible gap and 0  − and 0  + compressible states) is observed around the charge neutrality point (CNP) in Figs.2c,d,f at elevated .The positive and negative (yellow and blue) signals above and below the zero-magnetization strip are the paramagnetic responses in the LL gaps ∆  1 and ∆  −1 , as marked in Fig. 2f and Extended Data Figs.4e,f.

Fig. 3 .
Fig. 3. Imaging the beating of QOs and mapping the pseudomagnetic fields.a-d, Spatial imaging of    (, ) at   = 320 mT and  = 7.26×10 12 (a), 6.41×10 12 (b), 4.08×10 12 (c) and 2.81×10 12 (d) cm -2 corresponding to dotted lines in (f).The black rectangle indicates the boundaries of the TLG.e, Map of the derived pseudomagnetic field   across the sample.Regions with   below our resolution of 1 mT are shaded in white.f, Line scans of    () vs.  measured along the dotted line marked in (a) showing QOs due to MLG LLs in M1 band.g, The measured QOs due to MLG LLs (upper panel) at location A indicated in (b) and the calculated QOs (bottom panel).h, The measured QOs (upper panel) at location B and the calculated QOs (bottom panel) with   = 4.2 mT.The LL indices at the beating nodes are indicated.i, Same as (h) at   = 170 mT.The applied larger    = 10 mV rms in (a-d) and (i), and 20 mV rms in (f-h) averages out the QOs due to BLG LLs intensifying the visibility of MLG LLs.

Figure
Figure 4b shows the theoretical   1 dependence on   for   = 320 mT, and the calculated beating patterns of the QOs vs.   are presented in Fig. 4d with overlaid traces (red) of the beating nodes   1 to   3 .At the nodes, the  + and  − LLs are out of phase, giving rise to amplitude suppression and to a barely visible frequency doubling.In Fig. 3h (red curve) we find that the first node occurs at   1 = 19, from which we derive   =   /�4  1 � = 4.2 mT.The theoretically calculated QOs in the presence of   = 4.2 mT

Fig. 4 .
Fig. 4. LL interferometry of strain-induced pseudomagnetic field.a, Schematic of LL beating in the presence of PMF (  ) in graphene.Only the MLG Dirac bands with LLs are shown for clarity.b, The LL index of the first beating node   1 in the Dirac band vs.   at   = 320 mT.Inset: schematic of a graphene strip with arclike bent section (red) of length  and bend angle  generating strain-induced   .The illustrated  is greatly exaggerated compared to the maximal derived  ≅ 6×10 -3 degree.c, Calculated dependence of   1 on   for   = 4.2 mT (red).The open circles show the measured   1 at   = 170 and 320 mT.d, Calculated QOs in the MLG band vs.  and   at   = 320 mT and  = 0 V/nm.The locations of the beating nodes are highlighted in red.e, The measured QOs in the MLG bands vs.  and  using    = 20 mV rms showing beating nodes (red arrows).f, Calculated QOs vs.  and  at   = 320 mT and   = 4.2 mT.

1. 8 Extended Data Fig. 2 |
kHz was applied to the backgate to modulate the carrier density by   =      /.A lock-in amplifier was used to measure the corresponding local    by the scanning SOT.The    data were symmetrized with respect to the displacement field  where applicable.The 2D    (, ) images were used to reconstruct the magnetization   (, ) using numerical inversion procedure described in Ref. (67) as shown in Extended Data Fig. 2. Since the reconstruction of   requires 2D    (, ) information, the QOs at a single location or along 1D line scans are presented in the main text as the raw data of    .Reconstruction of the local magnetization from    (, ) .a, Example of the measured    (, ) at  = 5.16×10 12 cm -2 ,   = 320 mT, and    = 10 mV rms.b, Reconstructed differential magnetization   = ∂/ ∂.

Fig. 3 |
Comparison of    (, ) at different magnetic fields and    .a, The measured QOs at   = 320 mT and    = 8 mV rms.The induced peak-to-peak carrier density modulation,   =1.47×10 10 cm -2 , is about half of the LL degeneracy 4   0 ⁄ allowing clear resolution of the BLG and MLG LLs.b, QOs at   = 320 mT and    = 35 mV rms.The BLG LLs are washed out by the large carrier modulation while MLG LLs and the 12-fold degenerate LLs in the gullies are well resolved.c, At    = 100 mV rms the MLG QOs are smeared out.The paramagnetic response in the MLG LL gaps with Chern numbers  = ±2 are clearly visible as indicated.d, QOs at   = 40 mT and    = 20 mV rms.At this low field most of the LLs are smeared by the intrinsic LL broadening and only the lowest MLG LLs in sections M1 and M2 can be resolved at low displacement fields.e, Same as (d) at   = 80 mT.f, At   = 170 mT and    = 8 mV rms the BLG LLs cannot be resolved, but all the MLG LLs and the 12-fold degenerate LLs in the gullies are very prominent.Extended Data Figs.3d-f show the QOs at   = 40, 80, and 170 mT.At these low fields the Dingle broadening greatly suppresses the QOs due to BLG LLs (see Extended Data Fig. 7), and reduces the visibility of the MLG LLs at large displacement fields due to the reduction in the gap energies.At 170 mT and    = 8 mV rms the M2 LLs and the 12-fold degenerate LLs in the gullies are resolved very clearly as seen in Extended Data Fig. 3f.Band structure calculations The band structure of ABA graphene was calculated within the tight-binding model following Ref.(2, 42) based on SWMc parameterization (40).In the basis { 1 ,  1 ,  2 ,  2 ,  3 ,  3 }, where   ,   are the two sublattice sites in the  th layer, the low energy effective Hamiltonian can be written as (1) where  1 = −( 1 −  3 )/2 and  2 = −( 1 − 2 2 +  3 )/6 , with   the potential of layer  . 1 is determined by the displacement field, while  2 describes the asymmetry of the electric field between the layers.The band velocities   ( = 0,3,4) are related to the tight-binding parameters   by   ℏ = √3 2     , where   = 0.246 nm is the crystal constant of graphene,  =   +   , and  is the valley index ( = ±1 for valley  + and  − respectively).With a rotated basis ( 1 −  3 )/�(2), ( 1 −  3 )/�(2), ( 1 +  3 )/�(2),  2 ,  2 , ( 1 +  3 )/�(2) , the

Fermi velocity 𝑣𝑣 0 Extended Data Fig. 4 |
given by  0 and the gap at the Dirac point by   0 =  +  2 − 52 .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 .Evolution of ABA graphene band structure with displacement field.a-d, Projected 3D band structure of ABA graphene (left panels) and the corresponding evolution of LLs with   at different values of potential difference ∆ 1 = 0, 10, 25, and 50 meV between the top and bottom graphene layers (right panels).Red and blue lines denote the LLs in  + and  − valleys respectively.The tight-binding parameters used for the calculations are given in Table1(bottom row).With increasing ∆ 1 , the MLG and BLG bandgaps grow and band hybridization is enhanced forming mini-Dirac cones (gullies) near CNP.e, Zoomed-in view of the evolution of LLs with   at ∆ 1 = 50 meV.At low fields the 0 th LLs 0  + and 0  − are six-fold degenerate (including spin) and have vanishing magnetization.The six-fold degeneracy of the gully LLs is partially lifted at high   .The Chern numbers  in the large gaps are indicated.f, Evolution of LLs with ∆ 1 at   = 320 mT.Extended Data Fig. 5 | The dependence of the band structure on the SWMc parameters   to   , , and   at   =  .The panels show the effect of individual parameters on the BS, calculated using the parameters of Ref. (46) (see Table 1 ) and multiplied by a factor 0.8, 1.1, 1.4, 1.7 and 2 denoted by the different colors from black to red.The dependence of the measured QOs on  and  at low   provides a very sensitive tool for determining the SWMc parameters.The following attributes are particularly informative: (i) The number of BLG LLs between the adjacent MLG LLs.(ii) The relative energy shift between MLG and BLG bands.(iii) LL anticrossings in the gullies.(iv) The gap size of MLG band.

.Extended Data Fig. 6 |
is the inverse of the density of states as a function of the carrier density and () = ∫ ()  −∞ Calculations of the orbital magnetization.a, Calculated carrier density  as a function of chemical potential  and displacement-field-induced potential difference ∆ 1 .b, / vs.  and ∆ 1 .c, / vs.  and ∆ 1 .d, Calculated differential magnetization / vs.  .e-f, Calculated / (e) and differential magnetization   = / (f) vs. .  = 320 mT and  = 0.3 meV in all the calculations.Extended Data Fig.6shows the calculated () , vs. ∆ 1 at   = 320 mT using the derived SWMc parameters and Dingle broadening  = 0.3 meV.The modulation in DOS,   ⁄ , is well resolved in Extended Data Figs.6b,c, but it is relatively small due to the LL broadening, except near CNP where large gaps with vanishing DOS open between the lowest LLs in the gullies at elevated ∆ 1 .The calculated   ⁄ vs.  in Extended Data Figs.6d shows that the crossing of the MLG and BLG LLs does not cause any phase shift.In contrast, in   ⁄ vs.  in Extended Data Figs.6e, the BLG LLs show a 2 shift upon crossing the four-fold degenerate MLG LLs and a  shift upon crossing the two-fold degenerate 0 th LLs.This arises from the fact the filling a MLG LL delays filling the next BLG LL vs. total , but not vs. .Since the DOS modulation   () is quite small,   () in Extended Data Figs.6c looks very similar to   () except near CNP.

Extended Data Fig. 7 |
Suppression of QOs by LL broadening.a, Calculated () at  = 0 V/nm,   = 320 mT, and Dingle parameter  = 0.3 meV.The V-shaped dip in  (black arrow) corresponds to the McClure peak at the Dirac point of the MLG band (14).b-f, Calculated differential magnetization   = / for different  = 0.2 to 0.8 meV.The QOs due to BLG LLs with small energy gaps are suppressed much stronger by  than the MLG LLs with large gaps.
∆ 1 as shown in Extended Data Figs.8c,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 Figs.8a,b.Moreover, the strong hybridization lifts the valley degeneracy of the 1 st MLG LL in the M3 sector as shown by the pronounced splitting between −Figs.8b-d.This novel splitting is clearly resolved experimentally in Extended Data Figs.8a.Extended Data Fig. 8 | Landau level anticrossings.a, Measured QOs in   c vs.  and  at   = 320 mT near the top of the BLG valence band reproduced from Fig. 2f.b, Calculated   using the derived SWMc parameters providing an excellent fit to the experimental data.c, A zoom-in of the region marked in (b) with overlaid schematic LLs.d, Calculated LLs in the  − (blue) and  + (red) valleys.The open symbols denote the points of valley splitting and LL anticrossings between the MLG and BLG LLs.The splitting of −1 3 LL into valley polarized −1 3 − and −1 3 + LLs is clearly resolved in calculations and in the experimental data.

(
2 ��  +   − �  −   � ≅   �2ℏ  2 /  grows as √.Staggered substrate potential.The possible alignment between the hBN and ABA graphene can cause an onsite potential difference between the A and B sublattices.Here we consider the simplest situation where 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 in the basis of { 1 ,  1 ,  2 ,  2 ,  3 ,  3 } as: For concreteness, we choose  3 = 2 meV and  3 = −2 meV.The resulted BS is shown in Extended Data Fig. 10c (red) in comparison to 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. 10d 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 copper (73) or in the presence of strain (74).In contrast to the O-type Kekulé distortion which opens a gap at the Dirac point, we find that the Y-type (75) can result in LL interference.The Y-shaped modulation of the bond strength, parametrized by the hopping parameters  0 and  0 ′ (Extended Data Fig. 10e), 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.Within the SWMc model,  0 is the sole parameter that controls the Fermi velocity   of the MLG band (  = The energy difference between the LLs from the two valleys with the same index  is   = �2ℏ  ∆  ,where ∆  =  0 −  0 ′ ) .The first beating node appears when   is equal to half of the gap size:�2ℏ  ∆  = �ℏ  2   /2/2, which yields   1 =   /(4∆  ) as shown in Extended Data Fig. 10f.In Fig. 3h   1 = 19, which would correspond to a very weak Kekulé distortion with ∆    ⁄ = 1.4×10 -2 .Note, however, that the Kekulé distortion results in   1 that is independent of   as corroborated by the calculated QOs for   = 320 and 170 mT in Extended Data Figs.10h,i.The reason being that the LLs shift in the same proportion in the two valleys with   .This is in sharp contrast to beating due to PMF for which   1 =   /(4  ) is proportional to   .The experimental data points in Extended Data Fig. 10g (circles) are clearly consistent with PMF and incompatible with the Kekulé distortion.

Data Fig. 10 |
Alternative mechanisms of LL interference.a, Band structure of ABA graphene with a relative energy shift of  = 1.25 meV between the  + and  − bands.b, The corresponding calculated QOs show splitting of the lowest MLG LLs inconsistent with the experimental data.c, Band structure of ABA graphene with (red) and without (black) staggered sublattice potential  3 = 2 meV,  3 = − 2 meV.d, Corresponding calculated QOs showing no LL interference.e, Kekulé-Y distortion with  0 ′ hopping parameter along the bonds emphasized in black.f, The dependence of the first beating node   1 on the difference of the Fermi velocities ∆    ⁄ of the two valleys.Inset: schematic of the MLG band dispersions of the two valleys.g, The dependence of   1 on   for Kekulé-Y distortion (blue), PMF (red), and the experimental data (circles).h, The calculated QOs with ∆    ⁄ = 1.4% at   = 320 mT.i, Same as (h) at   = 170 mT.
Figures 2b,d present the dHvA oscillations calculated from the BS derived using the fitted SWMc parameters (Methods), showing remarkable qualitative and quantitative agreement with the experimental data.Such accurate comparison between the theory and experiment is possible because the thermodynamic QOs in the magnetization can be calculated quantitatively from the BS.The observed QOs can be classified by five prominent sets of LLs as follows: Due to the low DOS of the MLG bands and approximate linear dispersion, the MLG LLs are much sparser and disperse as ~±�|  |  .Displacement field opens a large gap   between the MLG sections (Figs.1f,g) resulting in parabolic-like upturn of M1 LLs with  in Figs.2c,d.When the B1 LLs are crossed by an M1 LL, there is a shift in the position of B1 LLs vs.  because the M1 LL has to be filled before the following B1 LLs can be occupied.The size of the shift depends on the degeneracy of the M1 LL.All the higher M1 LLs are four-fold degenerate, like the B1 LLs.As a result, at the crossing points, the B1 LLs show a phase shift of 2 as marked by the dotted white line in Figs.1c,e.The 0 th MLG LLs are fundamentally different.Due to the topological nature of the Dirac point and the associated Berry phase, the 0 th LL are valley polarized with the  − 0 th LL, 0 1 − , residing at the bottom of M1 band while the  + 0 th LL, 0 2 + , is pinned to the top of M2 band.Since the 0 th LLs are two-fold spin degenerate, their crossing with the four-fold degenerate B1 LLs gives rise to a phase shift of  rather than 2, as shown by the red dotted lines in Figs.1c,e (see also Methods and Extended Data Figs.6d . In the MLG Dirac band the LL energies are given by   = �2ℏ  2   , and hence for   ≪   the energy shift between the LLs in the two valleys is   = �2ℏ  2 ��  +   − �  −   � ≅   �2ℏ  2 /  .Since   scales with √, the lowest LLs remain essentially fully degenerate.For higher LLs, the relative shift between the valley-polarized LLs grows continuously with energy, giving rise to beating.The first destructive interference occurs when   = ∆  /2 , where ∆  =  +1 −   = �2ℏ  2   �√ + 1 − √� ≅ �ℏ  2 /(2) is the LL energy spacing, resulting in the first beating node at  =   1 =   /(4  ).

Table 1 | Comparison of SWMc parameters in different works.
The parameters are in units of meV.