Coulomb interactions and migrating Dirac cones imaged by local quantum oscillations in twisted graphene

Main

Moiré materials offer a combination of strongly correlated phases together with high experimental tunability. They provide a controllable window into the mechanisms and origins of interacting quantum materials. The example of magic-angle twisted bilayer graphene (MATBG)^{1,2,3,4,5,6,7,8,9,10} led the way. It has flat bands (FBs) whose electronic interactions produce prominent correlated insulating states and unconventional superconductivity. Subsequently, the same basic phenomena were discovered in alternating-twist trilayer graphene (tTLG)^{11,12,13,14,15,16,17,18}. In that material, the top and bottom layers are twisted with the same angle (theta) relative to the middle layer (Fig. 1a). This similarity is a consequence of mirror symmetry, which splits the tTLG band structure (BS) into two decoupled sectors^19,20, MATBG-like FBs and a bystander monolayer-graphene-like Dirac cone. The extra Dirac cone sector can be used to probe the correlated physics in the FB sector¹⁷. Additionally, the application of a transverse displacement field D breaks the mirror symmetry and hybridizes the two sectors, thus providing a different degree of in situ tunability. Indeed, possibly due to the larger magic angle, the correlated insulators and superconductors in tTLG are more stable than in MATBG²¹, making it a prime platform for exploring the nature of these phases.

Coulomb interactions and migrating Dirac cones imaged by local quantum oscillations in twisted graphene — **Fig. 1: Transport measurements in tTLG.**

However, despite extensive theoretical and experimental efforts, the identity of many phases—including the phase at the charge neutrality point (CNP)—remains unknown. Theoretical studies^{22,23,24,25,26,27,28,29,30,31,32,33}, including self-consistent Hartree–Fock (HF) and strong coupling analytics, have considered symmetry-broken states of magic-angle FB systems at integer fillings. Numerous ground states have been proposed, including the valley polarized state, the valley Hall state and Kramer’s intervalley coherent (KIVC) state²³, all of which are gapped, and a gapless nematic semimetal (NSM) that spontaneously breaks the ({C}_{3}) symmetry²⁴. However, there is little direct experimental evidence distinguishing between these orders, especially at the CNP. Transport studies have revealed Chern insulators at some integer fillings^3,8,10, but many of the underlying symmetries of the various many-body states remain unclear. More broadly, transport measurements are less well suited for identifying ground states as they do not directly probe the density of states (DOS), and merely a few techniques directly probe the BS. One example is scanning tunnelling microscopy, which has recently revealed an incommensurate Kekulé spiral order at (nu =pm 2) in twisted graphene systems^18,34,35. From a theoretical perspective, the ground state at the CNP plays a central role in determining the entire phase diagram. Yet, barring observations of (mathrm{C}_{3})-symmetry breaking by scanning tunnelling microscopy in some systems^{18,36,37,38,39}, there has been little experimental exploration of CNP ground states. This provides strong motivation for developing experimental techniques that can probe charge neutrality.

Here we study the local thermodynamic de Haas–van Alphen quantum oscillations (QOs)^40,41 in tTLG slightly away from the magic angle. We detected QOs in the Dirac sector down to a magnetic field of 56 mT. This was facilitated by the low DOS of the Dirac sector and the lower sensitivity of the local measurements to the effects of disorder. In addition to providing high-energy resolution, these oscillations allowed us to probe the ground-state physics and to directly visualize the symmetry breaking and BS renormalization due to the direct and exchange Coulomb interactions in the FBs.

By employing self-consistent mean-field Hartree calculation, we found that the FB dispersion was substantially renormalized by Coulomb repulsion and that the Hartree interaction term describes well the evolution of the bandwidth with filling^15,42,43,44. At half-filling of the conduction FB, we found clear evidence of spontaneous flavour-symmetry breaking^6,45, which suggests that the Stoner-polarized symmetry-broken phase is the parent state of the pervasive (nu =2) correlated insulator in MATBG^46,47. Most importantly, near the CNP we found strong evidence that the CNP ground state is an exchange-interaction-driven NSM that spontaneously breaks (mathrm{C}_{3}) symmetry. This phase is not generally favoured theoretically²³. However, even a small amount of heterostrain can make it energetically favourable⁴⁸. Our findings of a semimetallic ground state provide invaluable insights into the long-standing puzzle regarding the absence of an observable gap at the CNP in MATBG.

Transport measurements

The tTLG sample was fabricated using the dry-transfer method and encapsulated in hexagonal boron nitride (hBN; Fig. 1a,b and Methods). The d.c. voltages ({V}_{rm{tg}}^{;rm{dc}}) and ({V}_{rm{bg}}^{;rm{dc}}) applied to the top and bottom gates allowed us to control the carrier density (n) and (D) simultaneously. Transport measurements of ({R}_{{xx}}) and ({R}_{{yx}}) at constant D = 0, were performed at a temperature T = 300 mK as a function of the applied out-of-plane magnetic field ({B}_mathrm{a}) and (n), respectively (Fig. 1c,d). At D = 0, the mirror symmetry preserved the decoupling of the FB sector and the Dirac cone (Fig. 2h)^19,20. As such, two Landau fans emanate from the CNP: a dense fan emerging at large ({B}_mathrm{a}) reflecting Landau levels (LLs) in the FBs and a highly dispersing fan at low ({B}_mathrm{a}) reflecting the LLs arising from the Dirac sector (Fig. 1c and Extended Data Fig. 1a). As the Dirac band is populated in parallel with the FBs, it shunts any gapped FB state, including the fourfold filling of the FBs (which is highly resistive in MATBG). Therefore, extracting information from transport demands a careful analysis of the two sectors (Supplementary Information Section I)¹⁷.

**Fig. 2: Imaging Dirac LLs and Hartree interactions.**

The Chern number sequence in the Dirac band is ({C}_mathrm{D}=4{N}_mathrm{D}+2), and the energies of the corresponding LLs are

$${varepsilon }_mathrm{N}^mathrm{D}=operatorname{sgn}left({N}_mathrm{D}right){v}_mathrm{F}sqrt{2e{hbar }leftvert{N}_mathrm{D}rightvert{B}_mathrm{a}}.$$

(1)

Here, e is the elementary charge, ℏ the reduced Planck constant, ({v}_mathrm{F}) is the Dirac Fermi velocity and ({N}_mathrm{D}=0,pm 1,pm 2,ldots) is the LL index in the Dirac sector. The ({N}_mathrm{D}=0) LL resides at the CNP with ({C}_mathrm{D}=2), and the half-filled compressible ({N}_mathrm{D}=1) LL (dotted line in Fig. 1c) crosses the edge of the FB at ({n}_{0}=4.03times {10}^{12}) cm⁻² and ({B}_mathrm{a}=2.25) T (blue circle in Fig. 1c). At this point in phase space, the ({N}_mathrm{D}=0) LL is filled along with half of the ({N}_mathrm{D}=1) LL. Therefore, the density in the Dirac sector is ({n}_mathrm{D}=4{B}_mathrm{a}/{phi }_{0}=0.22times {10}^{12}) cm⁻², where ({phi }_{0}=h/e) is the flux quantum. The density in the FBs is then ({n}_mathrm{F}={n}_{0}-{n}_mathrm{D}=3.81times {10}^{12}) cm⁻², which we assigned to be the filling factor (nu =4) in the FB. This allowed us to determine the twist angle of the tTLG device, (theta =1.3^{circ}) (Supplementary Information Section I), which is slightly below the magic angle.

Inserting ({B}_mathrm{a}=2.25) T and ({N}_mathrm{D}=1) (Fig. 1c, blue circle) into equation (1), we found ({varepsilon }_mathrm{F}=W=49.4) meV, where (W) is the bandwidth of the conduction FB. This is substantially larger than (W) calculated from the non-interacting continuum model of (Wapprox24) meV (Fig. 2h), indicating that Coulomb interactions are occurring as discussed below. Note that the hole spectrum has no pronounced features within the FB, whereas the electron spectrum has a distinct dip in ({R}_{{xx}}) above (nu =2), along with non-trivial ({R}_{{yx}}) (Fig. 1d). Additionally, the slope of the Dirac LLs changes as they cross the (nu =2) peak in ({R}_{{xx}}) (black dotted lines in Fig. 1c and also Extended Data Fig. 1a), indicating strongly correlated behaviour in the FBs^11,12,17, as elaborated below.

Imaging QOs and the Dirac LLs

Under small ({B}_mathrm{a}), the dense LLs in the Dirac band serve as a powerful built-in ruler and spectrometer for probing the FB physics. We used a scanning superconducting quantum interference device fabricated on the apex of a sharp pipette (SQUID-on-tip, SOT)⁴⁹ to detect the de Haas–van Alphen QOs associated with Dirac LLs at ({B}_mathrm{a}lesssim250) mT. This approach, in addition to offering a fine energy scale, provides position-dependent spectroscopy on the nanoscale⁴⁰. An indium SOT of about 160 nm diameter was scanned at a height of (happrox180) nm above the sample surface (Fig. 1a) at (T=300) mK (Methods). In addition to the d.c. voltages applied to the gates, a small a.c. voltage ({V}_{rm{bg}}^{;rm{ac}}) of 85 mV r.m.s. was applied to the bottom gate, which modulated the carrier density by ({n}^{rm{ac}}) corresponding to ({nu }^{rm{ac}}=0.03) r.m.s., and the resulting ({B}_{z}^{rm{ac}}left(x,yright)={n}^{rm{ac}}(mathrm{d}{B}_{z}/mathrm{d}n)) was imaged across the sample.

We explored first the valence FB upon hole doping. A simpler picture was expected, as no symmetry breaking was observed in transport. We measured the evolution of ({B}_{z}^{rm{ac}}left(xright)) with carrier density at (D=0) (Fig. 2a) by repeated scanning along the black dotted line in Fig. 1b, while incrementing (nu) in 0.02 steps at ({B}_{a}=251) mT. The oscillations in ({B}_{z}^{rm{ac}}) reflect the de Haas–van Alphen QOs in the magnetization due to alternating diamagnetic and paramagnetic equilibrium currents in the compressible and incompressible LL states⁵⁰. Owing to the high DOS in the FB sector and the low ({B}_mathrm{a}), only oscillations from Dirac LLs with larger energy gaps are discerned^41,50. As (|nu |) is increased, the frequency of the ({B}_{z}^{rm{ac}}left(xright)) oscillations increases, and at (left|nu right|approx4.2), there is an abrupt jump in the frequency, as the Fermi energy ({varepsilon }_mathrm{F}) enteres the dispersive bands. The spatial variation of the oscillations reflects the inherent disorder in the sample.

Figure 2b shows a cross section of ({B}_{z}^{rm{ac}}) as a function of doping along the dashed white line in Fig. 2a. Each maximum (red triangle) marks the centre of the ({N}_mathrm{D}) compressible LL whereas each minimum (blue triangle) marks the incompressible state above the same LL. In the range (left|nu right|lesssim4.2), where the Dirac LLs and FBs coexist, nine maxima were observed, (left|{N}_mathrm{D}right|=9). Using the Dirac LLs as an energy ruler (equation (1)), we found that the FB is fully occupied when (left|{varepsilon }_mathrm{F}right|) reaches ({varepsilon }_{9}^mathrm{D}=49.5) meV. This large derived (W) of the FB is in good agreement with transport data and is much larger than the calculated single-particle (Wapprox24) meV, emphasizing the key role of interactions. Similar behaviour can be seen in samples 2 and 3 in Extended Data Fig. 2 (Supplementary Information Section II).

Each QO period reflects the filling of another Dirac LL and an increase in the Dirac sector carrier density ({n}_mathrm{D}) by (Delta {n}_mathrm{D}=4{B}_mathrm{a}/{phi }_{0}=2.4times {10}^{10}) cm⁻². Therefore, by counting the periods upon increasing the total carrier density (n), we directly determined ({n}_mathrm{D}) as a function of doping (Fig. 2c). In the FB region ((left|nu right|lesssim4.2)), ({n}_mathrm{D}) grows monotonically with filling, reaching ({n}_mathrm{D}approx0.22times {10}^{12}) cm⁻², in agreement with the value extracted from transport. Upon doping the dispersive bands ((left|nu right|gtrsim4.2)), there was a sharp increase in the slope. This is seen clearly for the derivative (mathrm{d}{n}_mathrm{D}/mathrm{d}n) in Fig. 2e, which can be understood as the fraction of electrons that populate the Dirac sector as carriers are added to the system. Such large values of ({n}_mathrm{D}) and (mathrm{d}{n}_mathrm{D}/mathrm{d}n) cannot be accounted for in a single-particle picture and, therefore, interactions must be considered.

Charging self-energy and doping-dependent BS

To account for interactions, we began by considering the charging self-energy, which was modelled by a self-consistent calculation of the Hartree term (Supplementary Information Section IV.b). As the wavefunctions of the FB carriers are mostly localized to the AAA regions of the moiré pattern (Fig. 2g), an electron added to the system sees the charge build-up on the AAA sites as a background periodic potential^43,51,52. The strength of the periodic Hartree potential ({V}_mathrm{H}) (equation (11) in Supplementary Information Section IV.b) is proportional to (n), which leads to a doping-dependent BS, as shown in Fig. 2h–j (Supplementary Video 1 and Supplementary Information Sections IV and VI for BS calculations and parameters). Particle–hole symmetry allowed us to refer to the conduction band for convenience. At the CNP, the charging self-energy is zero, and the states at (Gamma) are at the highest energy in the conduction FB (Fig. 2h, purple arrow). Upon increasing the carrier density, the charging self-energy increases, and states away from (Gamma), which are mostly localized to the AAA regions, have their energies increased due to electronic build-up, eventually causing an inversion of the (Gamma) point to reside at the bottom of the conduction FB (Fig. 2j, purple arrow). This effect is accompanied by a large increase of the FB bandwidth relative to the non-interacting single-particle (W) (Fig. 2h)^15,52.

From the Hartree calculation, we extracted the evolution of ({varepsilon }_mathrm{F}) with filling. It crosses the top of the flat conduction band at an energy ~49 meV consistent with (W) extracted from the experiment. As the Dirac cone is static, ({varepsilon }_mathrm{F}) directly translates into the Dirac cone density, ({n}_mathrm{D}={4S}_mathrm{D}/4{uppi }^{2}=uppi {k}_mathrm{F}^{2}/{uppi }^{2}=(1/uppi){left({varepsilon }_mathrm{F}/{hslash }{v}_mathrm{F}right)}^{2}), where ({S}_mathrm{D}) is the (k) space area of the Dirac cone, ({k}_mathrm{F}) is the Fermi momentum and factor 4 comes from degeneracy. ({n}_mathrm{D}) versus doping is plotted in Fig. 2d. The data match the experimental results well (Fig. 2c). Furthermore, Fig. 2f shows the derivative (mathrm{d}{n}_mathrm{D}/mathrm{d}n), which is the fraction of carriers that populate the Dirac band upon further doping. (mathrm{d}{n}_mathrm{D}/mathrm{d}n) increases as a function of (n), reaching values of 0.1 for (left|nu right|lesssim4.2). In a non-interacting system, the BS is static and (mathrm{d}{n}_mathrm{D}/mathrm{d}n) would be given by the ratio of the Dirac DOS to the total DOS and be at most about 0.005 due to the extremely high DOS of the FBs, which causes ({varepsilon }_mathrm{F}) to change slowly with doping. When accounting for the charging self-energy, however, the DOS of the FBs plays only a minor role, as the BS evolves with doping. The charging energy causes ({varepsilon }_mathrm{F}) to increase with doping at a much faster rate, as the FBs are pushed up in energy, transferring more carriers into the Dirac sector (Fig. 2h–j and Supplementary Video 1). For (|nu |gtrsim4.2), the jump in (mathrm{d}{n}_mathrm{D}/mathrm{d}n) is reproduced well by the theory in Fig. 2f, reflecting the drop in the total DOS at ({varepsilon }_mathrm{F}) upon entering the dispersive bands. The agreement between the highly sensitive theoretical and experimental (mathrm{d}{n}_mathrm{D}/mathrm{d}n) for such an extended range in (nu), through the FBs and into the dispersive bands, is striking. This emphasizes the importance of the charging self-energy for providing a complete description of the mirror-symmetric twisted graphene family, which has the strongest superconductivity of all moiré systems²¹.

Symmetry breaking at half-filling

We now explore the electron-doped spectrum, which shows signatures of symmetry breaking in transport. Figure 3a shows ({B}_{z}^{rm{ac}}(x)) versus (n) at a lower ({B}_mathrm{a}=131) mT, which provides higher energy resolution (see Extended Data Fig. 3 for ({B}_mathrm{a}=251) mT data). For (nu lesssim2), the ({B}_{z}^{rm{ac}}) oscillations are like those of their hole counterpart, whereas for (2lesssim nu lesssim 2.6), a bunching of QOs is observed along with increased intensity. This is seen more clearly in the cross section of ({B}_{z}^{rm{ac}}) along the dashed white line in Fig. 3a (Fig. 3b) and in the derived (mathrm{d}{n}_mathrm{D}/mathrm{d}n) (Fig. 3c). For both low ((nu lesssim2)) and high fillings ((nu gtrsim2.6)), (mathrm{d}{n}_mathrm{D}/mathrm{d}n) behaves like a valence band (Fig. 2e). However, at half-filling of the conduction FB, there is a clear jump in (mathrm{d}{n}_mathrm{D}/mathrm{d}n) (shaded region), which is absent in the valence band (Fig. 2e). This signifies that a larger fraction of carriers are populating the Dirac sector and is a clear signature of interaction-driven degeneracy lifting.

Following the data, we took as an ansatz a simple Stoner instability model, where for (2le nu le 2.6), one flavour was increased in energy by ({varDelta}_{rm{HF}}) and the other was decreased by the same amount (Fig. 3e–g and Supplementary Information Section III). Based on HF calculations at even integer fillings in MATBG^20,23,24,32, we took ({varDelta}_{rm{HF}}=20) meV, which we add to the Hartree-calculated bands (Fig. 3f). The resulting occupation of the two flavours in this simple Stoner model is sketched in Fig. 3g. The shaded region marks the symmetry-broken state. Figure 3f shows that in this state, ({varepsilon }_mathrm{F}) is much larger than in the symmetric state (Fig. 3e), and consequently, more carriers are transferred into the Dirac sector. This results in an increase in the calculated (mathrm{d}{n}_mathrm{D}/mathrm{d}n), as shown in Fig. 3d (shaded region), in agreement with the data in Fig. 3c. Additionally, this mechanism is consistent with the increased amplitude of the ({B}_{z}^{rm{ac}}) oscillations, as in the symmetry-broken state there are fewer states in the FB sector to tunnel to, resulting in lower Dirac LL broadening and a larger signal.

Displacement field dependence

To explore the BS evolution with displacement field (D), we measured ({B}_{z}^{rm{ac}}) as a function of (nu) and (D) at a fixed SOT position (red circle, Fig. 1b) at ({B}_mathrm{a}=251) mT (Fig. 4a) and at ({B}_mathrm{a}=56) mT, which has a better energy resolution with denser LLs (Fig. 4b). At ({B}_mathrm{a}=251) mT and (Dapprox 0) (dashed box), the QOs as a function of (nu) are similar to those in Figs. 2a and 3a, as expected. However, the ({B}_mathrm{a}=56) mT data at (Dapprox 0) (dashed box) do not scale as expected. From equation (1), it follows that the ({N}_mathrm{D}=1) LL at 56 mT should appear at a lower energy by a factor of the square root of the ratio of ({B}_mathrm{a}), that is (sqrt{251/56}approx 2.1). For a high and uniform FB DOS, the ({N}_mathrm{D}=1) LL that appears at (nu =0.95) and (Dapprox 0) in Fig. 4a should appear at (nu =0.45) at 56 mT. Strikingly, Fig. 4b shows that the ({N}_mathrm{D}=1) LL appears at (nu =0.27), almost a factor of two below the above estimate. This strongly indicates that considering only the charging self-energy, which was used to elucidate our results thus far, cannot fully describe the correlated physics in this system close to the CNP.

Next, we analyse the dispersion of the Dirac LLs with (left|Dright|), or equivalently, as a function of the outer-layer potential difference (U). Finite (D) breaks the mirror symmetry and hybridizes the Dirac and FB sectors. In models where the graphene and FB Dirac cones overlap at the mini-Brillouin zone corners (K₊ and K₋), such as single-particle and Hartree-only systems, the hybridization occurs at the CNP (Extended Data Fig. 4). With increasing (U), the Dirac node gets pushed to higher energy. As a result, the 0th Dirac LL along with all the higher LLs are pushed away from the CNP (Fig. 4c, grey dashed line), as occurs for ABA trilayer graphene⁴¹. Figure 4a indeed shows that the higher LLs disperse away from the CNP as expected. However, there are two peculiarities. First, the 0th LL stays pinned at the CNP, independent of (D). Second, the ({N}_mathrm{D}=pm 1) LLs behave contrary to expectation, curving towards the CNP.

At lower ({B}_mathrm{a}=56) mT (Fig. 4b), we confirm that the ({N}_mathrm{D}=0) LL does not disperse, and in addition, all the detectable (left|{N}_mathrm{D}right| > 0) LLs curve towards the CNP. Combining the information from Fig. 4a,b, we conclude that the Dirac LLs residing at (left|nu right|lesssim1.5) curve towards the CNP, whereas LLs appearing at (left|nu right|gtrsim1.5) curve away from the CNP, in sharp contrast to the Hartree prediction in Fig. 4c (grey dashed line). The behaviour of the higher carrier density is, thus, apparently captured well by considering only charging self-energy, as in Figs. 2 and 3. However, as the charging self-energy becomes small near the CNP, it is necessary to consider exchange interactions to gain insights into the underlying physics.

Exchange interactions

The observed (D) dependence reveals a lack of hybridization at the CNP, which means that graphene and FB Dirac cones do not overlap at zero energy. This can arise either due to gapping of the FB Dirac cones, namely an insulating FB state forms, or because they move away from the mini-Brillouin zone corners (K₋ and K₊) to maintain a semimetallic FB phase. Note that the latter automatically breaks (mathrm{C}_{3}) rotational symmetry, independent of the exact nature of the ground state. At the CNP, the charging self-energy is zero but strong exchange interactions in the FBs, modelled by the mean-field Fock interaction, may drive symmetry-broken orders. HF calculations at the CNP in tTLG allow many candidate symmetry-broken ground states, including the valley polarization state, the valley Hall state, KIVC and the NSM²⁰. We first analyse the NSM ground state, in which exchange energy causes the FB Dirac cones to migrate towards the (Gamma) point (Extended Data Fig. 5b and Supplementary Information Section IV.d), spontaneously breaking (mathrm{C}_{3}) and no longer overlapping the graphene Dirac cones. As a result, at (U=0), the Dirac cones and FBs intersect at an energy equal to ({varDelta}_{rm{HF}}approxpm20) meV (Fig. 4d (blue) and Extended Data Fig. 5d). Consequently, at (Une 0), the hybridization occurs at energies around ({varDelta}_{rm{HF}}) (Fig. 4d, red curve), rather than near the graphene Dirac nodes as in the Hartree-only or single-particle cases (Extended Data Fig. 4). As (U) is increased, the hybridization strength increases, effectively lowering ({v}_mathrm{F}) (or equivalently the slope) of the Dirac cone below the hybridization energy, as can be seen in Fig. 4d (red).

This picture successfully explains the low-energy behaviour. First, the graphene Dirac node never hybridizes, so the ({N}_mathrm{D}=0) LL does not disperse. Second, increasing (left|Dright|) decreases ({v}_mathrm{F}) below ({varDelta}_{rm{HF}}), thereby reducing the energy of the Dirac LLs (equation (1)). This can be seen in Fig. 4c (purple), which shows the simulated Dirac LLs at ({B}_mathrm{a}=251) mT, by calculating the renormalized ({v}_mathrm{F}) as a function of (left|Dright|) in the NSM (Extended Data Fig. 5a). Finally, the experimentally observed boundary of (left|nu right|approx1.5) separating the upward- and downward-dispersing Dirac LLs corresponds to a hybridization energy of about 19 meV, in good agreement with ({varDelta}_{rm{HF}}) from the simulation.

Note that other HF ground states similarly result in a renormalization of ({v}_mathrm{F}) with (U). These states, such as KIVC, have gapped FB Dirac cones, which also leads to hybridization at higher energies ({varDelta}_{rm{HF}}) (Extended Data Fig. 5c). Therefore, the measurements of the (D) field dependence provide irrefutable evidence that the FB Dirac cones either break (mathrm{C}_{3}) and migrate away from the mini-Brillouin zone corners, as in the NSM case, or have been gapped, as for KIVC.

Ground state at the CNP

We return focus on the (Dapprox 0) segment of Fig. 4b. Figure 4e shows a cross section of ({B}_{z}^{rm{ac}}) averaged around (Dapprox 0) in the dashed black box. The dense LLs at ({B}_mathrm{a}=56) mT carry a wealth of information about the correct ground state. Note that at this ({B}_mathrm{a}), the ({N}_mathrm{D}=1) LL resides at energy ({varepsilon }_{1}^mathrm{D}=7.8) meV, observed as a diamagnetic peak in the QOs (blue arrow in Fig. 4e). Furthermore, the first paramagnetic peak, due to equilibrium currents flowing in the incompressible gap between the ({N}_mathrm{D}=0) and 1 LLs (red arrow in Fig. 4e), can be approximated to be at energy ({varepsilon }_{1}^mathrm{D}/2=3.9) meV. Figure 4e shows that this paramagnetic peak appears at density (nu approx0.1) corresponding to (n=9times {10}^{10}) cm⁻², whereas only ({n}_mathrm{D}=2{B}_mathrm{a}/{phi }_{0}=2.7times {10}^{9}) cm⁻² carriers reside in the Dirac band. This means that at this energy, most of the carriers reside in the FB, namely the FB has a substantial DOS below 3.9 meV. Hence, if a gap exists, it must be less than a few milli-electronvolts. In fact, a more careful detailed analysis of the QO data allowed us to place an even lower bound of <1 meV on the possible gap at the CNP (Extended Data Fig. 6 and Supplementary Information Section V). This bound seems incompatible with a correlated gapped phase, where the gap is expected to be comparable to the Coulomb repulsion energy scale ~20 meV (refs. ^23,24).

Furthermore, the separation between the upward- and downward-dispersing LLs occurs at (left|nu right|approx1.5), which corresponds to a hybridization energy and, hence, a K-point FB gap of about ±19 meV. There is no generic theoretical reason for an order of magnitude difference between the K- and (Gamma)-point gaps in the FBs. In the vast theoretical literature on mirror-symmetric twisted graphene, there are no ground-state candidates with a K-point gap ~19 meV and a tiny global gap. These theoretical and experimental considerations provide strong evidence for the lack of a gap in the FBs, and they effectively rule out the KIVC or any other gapped ground state at the CNP. We, thus, conclude that the CNP ground state is a (mathrm{C}_{3})-breaking NSM.

To verify this conclusion, we extracted (mathrm{d}{n}_mathrm{D}/mathrm{d}n) (Fig. 4f) from the QOs in Fig. 4e at ({B}_mathrm{a}=56) mT. For the lowest Dirac LLs, (|{N}_mathrm{D}|le 3), we neglected the charging self-energy and assumed a static BS renormalized by the exchange energy. Thus, (mathrm{d}{n}_mathrm{D}/mathrm{d}n) is simply the relative DOS of the Dirac sector in the static BS, which can readily be calculated, as shown in Fig. 4g for the Hartree term only case (grey) and for the HF NSM case (purple). For large (nu) ((left|{N}_mathrm{D}right| > 3)), the charging self-energy becomes dominant and the data follow (mathrm{d}{n}_mathrm{D}/mathrm{d}n) calculated with the Hartree term only. Close to the CNP, however, the HF NSM has the correct FB DOS needed to qualitatively match the experimental values. This provides further corroboration that the ground state at the CNP is, indeed, an NSM.

Discussion

Generally, HF analytics find the KIVC or valley Hall state to be energetically favourable with respect to the NSM state²³. However, even small amounts of heterostrain drive a transition that stabilizes the NSM⁴⁸. Twisted systems in general have been shown to have substantial strain⁵³, including tTLG, for which scanning tunnelling microscopy measurements show that the angle mismatch often relaxes into the mirror-symmetric configuration¹⁴. The apparent strain is probably related to the spatial variations in the QOs observed in Figs. 2a and 3a. Our finding of a NSM ground state at the CNP is consistent with these considerations and consistent with most experiments that find no evidence of a gap at the CNP^2,4,7,8.

The regime of slightly off magic-angle graphene, where correlated-electron effects are weaker, provides opportunities for studying symmetry-breaking instabilities as well as their hierarchy. Our finding that the Stoner transition persists even in the absence of insulating states (Fig. 3) establishes the Stoner polarized states as the parent state for the emergence of correlated insulators. By reducing the interaction strength, we found the transition near half-filling to be most robust, in agreement with the appearance of the strongest correlated insulator state in MATBG at the same filling.

The high sensitivity of the low-magnetic-field thermodynamic QOs to the BS makes our measurement technique a powerful probe of low-energy interaction effects and fragile ground states that have yet remained unsolved in highly correlated systems. The emerging framework isolates the interaction effects at high doping, which are governed by the charging self-energy, from those at low doping where the exchange energy becomes dominant. This technique with thermodynamic QOs can readily be generalized to interacting systems that do not naturally contain a Dirac cone in the BS. One can add a Dirac band to essentially any van der Waals system by adding another monolayer graphene sheet twisted at a large angle^54,55, such that the graphene and the system of interest share the charge density but are effectively isolated from each other at low energy.

Methods

Device fabrication

The hBN-encapsulated tTLG devices were fabricated using the dry-transfer method. The flakes were exfoliated onto a Si/SiO₂ (285 nm) substrate and picked up using a polycarbonate on polydimethylsiloxane dome stamp. The number of graphene layers was determined by Raman microscopy, and the crystallographic orientations of the hBN and tTLG were determined from their straight edges. The WITec alpha300 R Raman Imaging Microscope was used to carry out Raman measurements using a wavelength of 532 nm and to cut monolayer graphene into three pieces using a 1,064 nm laser beam. During the dry-transfer process, the crystal axes of each layer of the tTLG stack were aligned by a mechanical rotation stage. After encapsulation with hBN, the stacks were released onto a pre-annealed graphite local bottom gate (~7–10 nm) patterned on an Si/SiO₂ wafer for devices 1 and 2. For device 3, p-doped Si was used as the bottom gate. The finalized stacks were annealed in vacuum at 350 °C to release the strain. A Ti (2 nm)/Pt (12 nm) top gate was then deposited on top of the stack for device 1. For devices 2 and 3, the top gate was Ti (4 nm)/Au (14 nm). The one-dimensional contacts were etched by SF₆ and O₂ plasma, followed by deposition of contacts of Cr (4 nm)/Au (70 nm) in an angle-rotated e-gun evaporator. Subsequently, the Hall bar geometry was etched using SF₆ and O₂ plasma. Finally, the surface resist and etching residues were swept off by atomic force microscopy in contact mode.

Device summary

Device 1 (presented in the main text): Graphite bottom gate, Ti/Pt top gate; bottom hBN thickness ~55 nm, bottom gate capacitance ({C}_{rm{bg}}approx 3.4times {10}^{11}e) cm⁻² V⁻¹, top hBN thickness ~32 nm, top gate capacitance ({C}_{rm{tg}}approx 5.8times {10}^{11}e) cm⁻² V⁻¹ and twist angle ({theta }_mathrm{M}approx 1.3^circ). The capacitances were calculated by fitting the slopes of Chern insulator lines in the magneto-transport measurements.

Device 2: Graphite bottom gate, Ti/Au top gate; bottom hBN thickness ~36 nm, ({C}_{rm{bg}}approx 5.2times {10}^{11}e) cm⁻² V⁻¹, top hBN thickness ~60 nm, ({C}_{rm{tg}}approx 3.1times {10}^{11}e) cm⁻² V⁻¹ and twist angle ({theta }_mathrm{M}approx 1.5^circ).

Device 3: p-doped Si bottom gate, Ti/Au top gate; bottom hBN thickness ~36 nm, ({C}_{rm{bg}}approx 7.5times {10}^{10}e) cm⁻² V⁻¹, top hBN thickness ~60 nm, ({C}_{rm{tg}}approx 3.4times {10}^{11}e) cm⁻² V⁻¹ and twist angle ({theta }_mathrm{M}approx 1.5^circ).

SOT fabrication and QO measurements

The a.c. magnetic field measurements were done using an indium SOT of 160 nm effective diameter, which was fabricated as described previously^49,56,57. The imaging was performed at (T=300) mK using a cryogenic SQUID series array amplifier⁵⁸. The SOT had a magnetic field sensitivity down to 10 nT Hz⁻^1/2 and was attached to a quartz tuning fork excited at its resonance frequency of ~33 kHz for height control as described in ref. ⁵⁹. The scanning was performed at a height of 180 nm above the sample surface.

The ({B}_{z}^{rm{ac}}) images were obtained with pixel size of 90 nm and acquisition time of 1 s per pixel. The measured signal ({B}_{z}^{rm{ac}}={n}^{rm{ac}}left(mathrm{d}{B}_{z}/mathrm{d}nright)) is proportional to the modulation in the carrier density ({n}^{rm{ac}}) induced by a small a.c. voltage ({V}_{rm{bg}}^{;rm{ac}}) of 85 mV r.m.s. applied to the back gate at a frequency of (fapprox5) kHz. To obtain the best signal-to-noise ratio from the QOs in the Dirac band, the ({V}_{rm{bg}}^{;rm{ac}}) amplitude was optimized for the lowest ({B}_mathrm{a}=56) mT, for which the period of the QOs (Delta n) was the smallest.