Imaging the flat bands of magic-angle graphene reshaped by interactions

Cryogenic QTM

All measurements in this work were performed in a cryogenic QTM system operating at a temperature of 4 K, as described previously45. For conductance measurement, voltage biases are applied using a custom-built digital-to-analogue converter array, capable of supplying both d.c. and a.c. signals with 1 μV resolution. Currents were measured using a FEMTO current amplifier, followed by a National Instruments sampler.

Fabrication and characterization of the QTM tip and van der Waals device

The fabrication of both QTM tips and van der Waals devices on flat substrates were described previously45,46. Briefly, we use the standard dry transfer and polymer membrane transfer techniques to fabricate the TBG sample and the QTM tip, respectively. Extended Data Fig. 1a shows an optical image of the TBG heterostructure, comprising bilayer-WSe2/TBG/hexagonal boron nitride (h-BN)/Pt on an approximately 100 μm width and 80 μm tall Si/SiO2 pillar. Extended Data Fig. 1b shows an optical image of the MLG backed by h-BN and graphite placed on a QTM tip.

In Extended Data Fig. 1c, we characterize the TBG sample by conductive AFM scan of the region of magic twist angle carried out at room temperature (Bruker Icon). To compensate for the thermal drift during this AFM scan, we performed a small window scan at rate 10.9 Hz, such that the time for an entire scan is around 20 s. We apply an approximately 1 V d.c. voltage bias between the tip and sample and measure the tunnelling current. The scan image in Extended Data Fig. 1c shows both the moiré lattice and atomic defects in WSe2 (bright spots). To accurately correct the thermal drift, we trace the position of three defects inside the scan window. We then obtain the drift velocity by following the shift of defects positions between consecutive scans and assuming that, within a single 2D scan, the drift velocity is constant. Extended Data Figure 1c,d shows the image and corresponding fast Fourier transform (FFT) after thermal drift correction. The three moiré lattice vectors obtained from this image are: 12.72 nm, 12.73 nm and 13.14 nm. This gives a twist angle of θTBG = 1.1° ± 0.02° and strain of ϵ = 0.03% ± 0.02%.

In Extended Data Fig. 1e, we characterize the QTM tip by in situ imaging of the tip contact area by an atomic defect in the WSe2 barrier as a localized tunnelling channel. At large bias (Vb ≈ −900 mV, well outside the range used in the actual spectroscopy measurements), the defect level becomes energetically accessible and provides an extra, spatially localized current path. When the QTM tip passes over the defect, we observe a small but measurable increase in current. By scanning the tip across such defects and mapping this current enhancement, we obtain a real-space image of the tip’s contact footprint.

Extended Data Fig. 1e presents such a high-bias scan, showing several images of the tip footprint, each produced by a different defect, revealing a contact size of approximately 200 × 50 nm. These real-space dimensions set the momentum-space resolution through the Heisenberg uncertainty relation. The measured footprint corresponds to angular momentum resolutions (in angular units) of roughly δθQTM ≈ 0.05° and ≈ 0.2° along the two principal directions. In Extended Data Fig. 2, we show that the sharpest momentum features—obtained in the remote bands—exhibit a full width at half maximum (FWHM) of δθQTM ≈ 0.1°, consistent with the resolution expected from the measured tip dimensions.

Energy and momentum resolution

We determine the energy resolution of our QTM measurements by analysing a sharp spectral feature at ν = −0.7 in Fig. 4a (reproduced in Extended Data Fig. 2a). A linecut through this feature yields a FWHM of δE ≈ 10 meV by fitting a Lorentzian function. This is an upper bound to the energy resolution of the measurements and it may well be coming from the intrinsic lifetime of the energy bands of MATBG at this temperature. To determine the momentum resolution, we analyse a sharp spectral feature of the remote energy bands of MATBG measured under the same conditions as in Fig. 2a. Extended Data Fig. 2b plots a linecut of dI/dV along the dashed white line, cutting through the dispersive bands, yielding an angular FWHM of δθ = 0.1° by fitting a Lorentzian function. This translates into a δk ≈ 0.03 nm−1 momentum resolution. For the flat bands, we trace the dI/dV peak width versus θQTM at zero bias around the Γ point and identify the peak width δθ ≈ 0.5°. This peak width is composed of two peaks overlapping each other, such that the single peak broadening is δθ ≈ 0.25°.

Electrostatics of the QTM junction

The QTM junction is modelled as a three-plate capacitor consisting of: (1) the MLG on the QTM tip; (2) the TBG layer on the sample; and (3) the metallic back gate. Below, to stay consistent with our previous publication, we use in all of the formulas the notation ‘top layer’ (‘T’) and bottom layer (‘B’) to refer to the two layers between which the electrons perform the QTM’s momentum-resolved tunnelling, not the top and bottom layers of TBG. Namely, the top layer is the MLG layer on the QTM tip and the bottom layer is the TBG on the sample side.

The QTM tip and TBG are separated by two layers of WSe2 (1.2 nm) and the TBG to the back gate by a 37 nm h-BN layer. The corresponding interlayer geometric capacitances are denoted as Cg (tip–TBG) and Cbg (TBG–back gate). A voltage bias Vb is applied between the tip and TBG sample layer. The gate voltage Vbg is applied between the TBG sample and the metal back gate. The electrostatic configuration of the measurements can be described by the following46:

$$-e{V}_{{\rm{b}}}={\mu }_{{\rm{B}}}-{\mu }_{{\rm{T}}}-\frac{{e}^{2}{n}_{{\rm{T}}}}{{C}_{{\rm{g}}}},$$

(1)

$$e{V}_{{\rm{bg}}}={\mu }_{{\rm{B}}}+\frac{{e}^{2}({n}_{{\rm{B}}}+{n}_{{\rm{T}}})}{{C}_{{\rm{bg}}}}$$

(2)

Note again that the notation ‘top’ (‘T’) refers to the MLG on the QTM tip and the notation ‘bottom’ (‘B’) refers to the TBG. μB and μT are the chemical potentials of the bottom TBG and top graphene probe, respectively, and nB and nT are their corresponding carrier densities. The back-gate capacitance was experimentally calibrated to be Cbg = 62.3 nF cm−2 based on the gate voltage required to reach filling factor ν = 4 in the magic-angle region (as verified by the conductive AFM measurements; Methods section ‘Fabrication and characterization of the QTM tip and van der Waals device’).

It is important to note that, under a finite voltage bias, TBG becomes doped as a result of the gating from the tip. This means that, at a fixed Vbg, the filling factor ν changes with the tip bias, Vb. To account for this and present our data as a function of the TBG’s filling factor, we determine the gating effect directly from the data. Extended Data Fig. 3 shows the raw measurement of QTM tunnelling conductance, dI/dVb (for simplicity called dI/dV throughout the paper), and the transconductance with respect to the back gate, dI/dVbg, as functions of Vb and Vbg. At integer fillings, and in particular at ν = ±4, dI/dVb shows a suppression and dI/dVbg shows a dipole-like feature. By tracing the ν = 4 line in the dI/dVbg measurements (through a polynomial fitting, overlaid on the data), we extract the trajectory of constant carrier density in the Vb–Vbg plane. We then ‘skew’ the raw data by plotting it as a function of ν that is determined by both Vbg and Vb. The skew-corrected data are shown in Figs. 3 and 4. For Figs. 1 and 2, the scans were taken at a fixed gate voltage without skew correction, so at high bias, the filling factor may deviate from the nominal one.

Along this skewed line, the TBG carrier density nB and chemical potential μB remain fixed. By using the graphene relation \({n}_{{\rm{T}}}={\mu }_{{\rm{T}}}^{2}/\pi {\nu }_{{\rm{f}}}^{2}\), the electrostatic equation can be further simplified to:

$$-e{V}_{{\rm{b}}}={\mu }_{{\rm{B}}}-e{V}_{{\rm{b}}{\rm{g}}}^{\ast }\frac{{C}_{{\rm{b}}{\rm{g}}}}{{C}_{{\rm{g}}}}-\sqrt{e{V}_{{\rm{b}}{\rm{g}}}^{\ast }\frac{{C}_{{\rm{b}}{\rm{g}}}\pi {\nu }_{{\rm{f}}}^{2}}{{e}^{2}}},$$

(3)

in which \({{eV}}_{{\rm{bg}}}^{\ast }={{eV}}_{{\rm{bg}}}-{\mu }_{{\rm{B}}}-{e}^{2}{n}_{{\rm{B}}}/{C}_{{\rm{bg}}}\), giving for every Vb the corresponding shift in Vbg to maintain fixed nB and μB. By fitting the skew lines with this relation, we derive the interlayer capacitance to be Cg = 2.3 μF cm−2, consistent with our previous tunnelling experiment using two layers of WSe2 (ref. 46).

The tunnelling matrix elements

The tunnelling current is modelled by Fermi’s golden rule assuming momentum-resolved tunnelling between the TBG sample and the graphene probe. As in the previous section, to stay consistent with our previous publication, in all of the formulas below, we use the notation ‘top layer’ (‘T’) and ‘bottom layer’ (‘B’) to refer to the two layers between which the electrons perform the QTM’s momentum-resolved tunnelling, not the top and bottom layers of TBG. Namely, the top layer is the MLG layer on the QTM tip and the bottom layer is the TBG on the sample side.

The calculation is described in detail in refs. 46,47 and here we briefly review the formula and tunnelling matrix elements:

$$I=\frac{2\pi e{g}_{{\rm{s}}}{g}_{{\rm{v}}}}{\hbar }\int {\rm{d}}\omega (\,{f}_{{\rm{B}}}(\omega )-{f}_{{\rm{T}}}(\omega +e\phi ))\sum _{n,n{\rm{{\prime} }}}\sum _{{{\bf{k}}}_{{\rm{T}}},{{\bf{k}}}_{{\rm{B}}}}{|\langle {u}_{{\rm{gr}},{{\bf{k}}}_{{\rm{T}}},n}|{T}^{{{\bf{G}}}_{{\rm{T}}},{{\bf{G}}}_{{\rm{B}}}}|{u}_{{\rm{tbg}},{{\bf{k}}}_{{\rm{B}}},n{\rm{{\prime} }}}\rangle |}^{2}{\delta }_{{{\bf{k}}}_{{\rm{T}}}+{{\bf{G}}}_{{\rm{T}}},{{\bf{k}}}_{{\rm{B}}}+{{\bf{G}}}_{{\rm{B}}}}{A}_{{{\bf{k}}}_{{\rm{B}}},\omega }{A}_{{{\bf{k}}}_{{\rm{T}}},\omega +e\phi }$$

(4)

gs/v is the spin/valley degeneracy, kB/T is the electron momentum of the bottom/top layer, n and n′ are the band indexes, fB/T(E) is the Fermi–Dirac distribution function and ϕ is the electrostatic potential that shifts the relative position of energy bands between TBG and graphene, derived from the electrostatic model in Methods section ‘Electrostatics of the QTM junction’. The chemical potential μB/T is defined relative to their charge neutrality (Dirac point of TBG and graphene probe). \({A}_{{{\bf{k}}}_{{\rm{B}}/{\rm{T}}},E}\) is the spectral function of the TBG and the graphene probe. \({T}^{{{\bf{G}}}_{{\rm{T}}},{{\bf{G}}}_{{\rm{B}}}}\) is the tunnelling matrix in sublattice space between the graphene probe wavefunction \(|{u}_{{\rm{gr}},{{\bf{k}}}_{{\rm{T}}}}\rangle \) and the TBG wavefunction \(|{u}_{{\rm{tbg}},{{\bf{k}}}_{{\rm{B}}}}\rangle \) projected on the top layer. GT and GB are the reciprocal lattice vectors of the graphene probe and the top layer graphene of TBG, at which it satisfied the relation kT + GT = kB + GB. We assume that the tunnelling happens within the first BZ of graphene probe and it limits the Bragg scattering process of {GT, GB} choices to be three. For example, when GT = GB = 0, the matrix is given by: \({T}^{{{\bf{G}}}_{{\rm{T}}}=0,{{\bf{G}}}_{{\rm{B}}}=0}=I+{\sigma }_{x}\), in which I is the unity matrix. The other two Bragg scattering processes are related by C3. Here we focus on \({T}^{{{\bf{G}}}_{{\rm{T}}}=0,{{\bf{G}}}_{{\rm{B}}}=0}\) and the matrix elements are given by:

$$\begin{array}{c}{|\langle {u}_{{\rm{g}}{\rm{r}},{{\bf{k}}}_{{\rm{T}}}}|{T}^{{{\bf{G}}}_{{\rm{T}}}=0,{{\bf{G}}}_{{\rm{B}}}=0}|{u}_{{\rm{t}}{\rm{b}}{\rm{g}},{{\bf{k}}}_{{\rm{B}}}}\rangle |}^{2}\\ \,=\,{|{\langle ({u}_{{\rm{g}}{\rm{r}},{{\bf{k}}}_{{\rm{T}}},A}^{\ast }+{u}_{{\rm{g}}{\rm{r}},{{\bf{k}}}_{{\rm{T}}},B}^{\ast })\rangle }_{FS}({u}_{{\rm{t}}{\rm{b}}{\rm{g}},{{\bf{k}}}_{{\rm{B}}},A}+{u}_{{\rm{t}}{\rm{b}}{\rm{g}},{{\bf{k}}}_{{\rm{B}}},B})|}^{2}\\ \,\propto \,{|{u}_{{\rm{t}}{\rm{b}}{\rm{g}},{{\bf{k}}}_{{\rm{B}}},A}+{u}_{{\rm{t}}{\rm{b}}{\rm{g}},{{\bf{k}}}_{{\rm{B}}},B}|}^{2}=|{\langle I+{{\sigma }}_{x}\rangle }_{{u}_{{\rm{t}}{\rm{b}}{\rm{g}},{{\bf{k}}}_{{\rm{B}}}}}|,\end{array}$$

(5)

in which \({u}_{{{\rm{gr}},{\bf{k}}}_{{\rm{T}}},A/B}^{\ast }\) is the graphene probe wavefunction with the sublattice component A/B. The graphene probe wavefunction is averaged over the Fermi surface and we assume that, within the graphene Fermi surface range, the TBG wavefunctions \({u}_{{\rm{tbg}},{{\bf{k}}}_{{\rm{B}}},A/B}\) remain the same. The above equation shows that the tunnelling current is proportional to the \(|{\langle I+{\sigma }_{x}\rangle }_{{u}_{{\rm{tbg}},{{\bf{k}}}_{{\rm{B}}}}}|\), in which the \(|{u}_{{\rm{tbg}},{{\bf{k}}}_{{\rm{B}}}}\rangle \) contains both the layer polarization and sublattice information.

In Fig. 2d, we show the matching between the matrix element \(|{\langle I+{\sigma }_{x}\rangle }_{{u}_{{\rm{tbg}},{{\bf{k}}}_{{\rm{B}}}}}|\) and the dI/dV intensity in the flat-band region. It shows an excellent agreement at ratio w0/w1 = 0.6. Here we further show that the naive layer polarization quantity (Extended Data Fig. 4b,d) is inconsistent with the data, regardless of the w0/w1 ratios.

For the tunnelling current calculation in Fig. 1, we input the single-particle band structure of TBG at θQTM = 1.2°, the electrostatics from Methods section ‘Electrostatics of the QTM junction’ and calculate the tunnelling current based on equation (4).

Spectroscopy at different momenta: its Gaussians decomposition and spectral function peak and lifetime analysis versus quasiparticle energy

In Extended Data Fig. 5, we show the tunnelling spectroscopy at different momenta along the flat-band region. The momenta at which the spectroscopy images are measured are marked by dashed black lines on a zoom-in of the bands from Fig. 2b, shown in the top panel. The range of momenta covers angles from θQTM = −0.5° to 0.8°. Apart from an overall magnitude, the spectroscopy in all momenta looks very similar. In Fig. 3c, we use this invariance on momentum in the flat bands to make Gaussian fit to the averaged data over the flat-band region.

The Gaussian fitting procedure is performed as follows. We first identify the peak centres directly from the spectroscopy data by tracing the local maxima in the contour representation of Fig. 4a; these trajectories are shown as dashed black lines in Extended Data Fig. 6a. Specifically, we identify two classes of features: (1) the ‘Hubbard’ bands that exhibit the characteristic cascading behaviour and (2) further excitations whose energy (about ±15 meV) is independent of filling. For each filling, the number of peaks identified in this manner determines the number of Gaussians used in the corresponding spectral fit. The peak positions extracted from the dashed-line trajectories serve as initial guess for the fit, although we allow these positions to vary during the fitting. The free parameters in the fitting routine are the Gaussian centres, widths, amplitudes and the background offset.

In Figs. 3b and 4a, we observed many spectral features originating from the Hubbard bands of heavy electrons. The spectral features typically appear at high energy as smeared ‘plumes’ and become stronger the closer they approach the Fermi level.

In a simple lifetime-broadening picture, the spectral function peak corresponding to a quasiparticle at an energy E acquires a width inversely proportional to the lifetime and the peak amplitude of the spectral function decreases correspondingly. In the present experiment, several bands coexist within a narrow energy range, making it challenging to extract the linewidths of the bands quantitatively. Instead, we can reliably extract the peak height, Apeak, which is a robust experimental observable and reflects the combined effects of quasiparticle residue and lifetime. Inspection of Fig. 4 shows a clear and systematic trend: when a band lies far from the Fermi energy, its Apeak is small and this maximal spectral weight increases continuously as the band approaches the Fermi level.

To address this more quantitatively, in Extended Data Fig. 6b, we track the maximum of the spectral peak (dashed line) from the highest energy that we can reliably identify a peak to the point at which it reaches the Fermi energy, and the peak becomes less well defined. In panel c, we plot the one over the extracted peak height 1/Apeak versus peak energy E and compare it with two functional forms: \(\frac{1}{{A}_{{\rm{peak}}}}=a(E-{E}_{{\rm{F}}})+b\) and \(\frac{1}{{A}_{{\rm{peak}}}}=a{(E-{E}_{{\rm{F}}})}^{2}+b\). Overall, we see that a linear dependence of 1/Apeak on E − EF better fits the observations.

The height of a spectral-function peak reflects both the quasiparticle residue Z(E − EF) and the lifetime τ(E − EF). In principle, both quantities may evolve as a Hubbard band moves to higher energies. If, however, we interpret the dominant trend as arising mainly from lifetime effects, the observed scaling would imply that it is more likely that τ(E −E F) ∝ 1/(E − EF) and not the Fermi-liquid expectation τ(E − EF) ∝ 1/(E − EF)2. Such dependence is often seen in strongly correlated/Hubbard systems, but because we do not measure this quantity in separation, we cannot make a decisive claim.

Fitting the transition between light and heavy electrons

Figure 5 shows dI/dV at Vb = 0 mV as a function of θQTM and filling factor ν. To determine the boundary between heavy and light electrons, we plot horizontal linecuts of dI/dV versus θQTM for each filling and perform a fit of a sigmoid-type function:

$$\frac{{\rm{d}}I}{{\rm{d}}V}={c}_{1}+\frac{{c}_{2}-{c}_{1}}{1+{e}^{-\frac{{\theta }_{{\rm{QTM}}}-{\theta }_{0}}{\sigma }}}$$

(6)

On the basis of the fitting, when the heavy-electron spectral intensity decays by 80%, we mark the position θhl(ν) as a proxy for the transition between heavy and light electrons. The parameter σ characterizes the width of this transition. We then fit θhl(ν) to a third-order polynomial and plot this as the dashed white lines in Fig. 5.

Determine the larger twist angle of θ
TBG = 1.2°

We determine the angle of larger-angle TBG from the measurement in Fig. 1d. By tracing the dI/dV intensity change at the conduction-band edge of flat bands around the Γ point, we performed a step function fitting of the form:

$$y=C+\frac{A}{2}(1+\tanh (B\ast (x-{x}_{0}))),$$

in which the x0 marks the Γ point position, which is half the value of the tunnelling intensity decay. By measuring the angle distance between the Γ point to the KT point, at which the Dirac point crosses, we determine the twist angle to be θTBG = 1.2°.

Spectroscopy of flat bands at different locations

Extended Data Fig. 7a–f presents tunnelling spectroscopy at a momentum point within the flat part of the bands (KT) measured at several locations across the sample, separated by several microns. These separations are large enough such that each region effectively behaves as independent devices. In all cases, we observe excited states whose energy is independent of filling, similar to the measurement in Fig. 4a.

For each location, we use room-temperature conductive AFM to image the real-space moiré structure (Methods section ‘Fabrication and characterization of the QTM tip and van der Waals device’) and plot its corresponding FFT. From the latter, we extract the local twist angle, θTBG, and heterostrain, ϵTBG, of the moiré. From the spectroscopy, we determine the excitation energy, ΔE. These three locally measured quantities are indicated above each panel. The centre panel shows a real-space map of the measurement locations, with the colour of each point representing the corresponding ΔE (see colour bar).

Although we cannot definitively rule out heterostrain as the origin of the observed approximately 15 meV excitation, the data strongly indicate that strain alone cannot account for this feature, for several reasons:

1. 1.

   Universality across widely separated regions of the same device. We observe nearly identical excitation energies (13–16 meV) at locations separated by several microns. It is highly unlikely that the same level of strain would persist uniformly over such large distances.

2. 2.

   Universality across different samples. Similar excitation energies are measured in two distinct MATBG devices, measured in separate cooldowns, indicating that the feature is neither device-specific nor related to a particular thermal cycle.

3. 3.

   Lack of correlation with independently measured strain. The heterostrain extracted from the moiré FFT varies by a factor of 3 across the measured positions (from ε = 0.02% to ε = 0.06%). Despite this substantial variation, the excitation energy remains essentially unchanged.

4. 4.

   Quantitative inconsistency with theoretical expectations. For the least-strained region (ε = 0.02%), the observed excitation is about 15 meV. Existing theoretical estimates predict a strain-induced level splitting of about 10 meV per 0.1% strain51, which would correspond to only about 2 meV splitting at ε = 0.02%, an order of magnitude smaller than what we observe.

QTM probing of the K and K′ valleys

In Extended Data Fig. 8, we explain why the time-reversal symmetry of the QTM probe dictates that it measures the combined contribution of the K and K′ valleys.

Panel a shows the inequivalent Fermi surfaces of the K and K′ valleys (red and blue) folded into the same mini-BZ. For illustration, we depict triangularly warped Fermi surfaces to emphasize that the two valleys can, in principle, be different.

To understand what the QTM actually probes, it is useful to unfold the picture back to the full BZ (panel b). There the K and K′ Fermi surfaces sit at opposite corners. In the experiment, the Dirac points of the probe (purple) move along a circular trajectory (dashed line) and current is detected whenever the Dirac point of the probe intersects a Fermi surface. Crucially, if both the probe and the sample obey time-reversal symmetry, the K and K′ valleys produce identical spectra as a function of the QTM angle θQTM (panel c). If we take the illustration in panel b as an example, although the two triangular Fermi surfaces are mirror reflections, the rotating Dirac points intersect them in exactly the same sequence—first the flat edge and then the triangle tip—and therefore yield identical spectroscopic measurement.

We might expect that spontaneous valley symmetry breaking (for example, a valley-polarized Chern insulator) would manifest as clear energy splitting between the K and K′ bands and that, consequently, it should be easy to identify flavour symmetry breaking from this band splitting. However, this expectation is incorrect—recent dynamical mean-field theory calculations performed in both the symmetric and flavour-symmetry-broken34,35 states show that very similar band splittings appear in both situations. These splittings originate primarily from the formation of Hubbard bands, which occur regardless of whether flavour symmetry is broken. As a result, the spectroscopic signatures of flavour symmetry breaking are subtle and cannot be identified simply by looking for energy-split bands.

Hartree-driven band stretching and its connection to wavefunctions at different momenta within the flat band

In this section, we provide a systematic explanation of the momentum-dependent band stretching observed near the Γ point. To highlight the effect, we plot in Extended Data Fig. 9h,i the measured bands from Fig. 3a alongside schematic guides to the eye marking the stretched areas indicated by arrows. Specifically, the stretching appears as a peak at Γ at ν = −4, which continuously evolves to a dip at Γ at ν = 4.

This stretching originates from the momentum-dependent response of the electronic wavefunctions to the Hartree potential that builds up with increasing carrier density. Early theoretical works (for example, refs. 21,22) already noted that, although moiré-scale electrostatic potentials are negligible at charge neutrality, a spatially varying Hartree potential VH(r) develops as electrons or holes are added. The energy of a state at momentum k then shifts according to its real-space overlap with this potential,

$$\delta {E}_{{\rm{H}}}(k)=\int {{\rm{d}}}^{2}r{V}_{{\rm{H}}}(r){|{\psi }_{k}(r)|}^{2},$$

and because different k-states have distinct real-space charge distributions, they experience different Hartree shifts. This naturally produces a stretching of the bands, even in the absence of strong correlations. Although indirect signatures of this effect have been reported previously (for example, through Landau-level spectroscopy in ref. 16), Fig. 3 presents the first, to our knowledge, direct observation of this Hartree-driven band stretching.

Notably, our measurements do not only show the overall magnitude of the band stretching but they directly provide the momentum-dependent energy shifts δEk. From the expression above, it is evident that these shifts encode detailed information about the real-space structure of the wavefunctions at different k. This is especially valuable in MATBG, in which the flat bands are known to have strongly momentum-dependent wavefunctions arising from their nontrivial topology.

In fact, one of the most notable theoretical predictions for TBG, already present in the BM single-particle model, is that the flat-band wavefunctions change qualitatively across the mini-BZ. This effect does not need electronic correlation and happens also when the bands are not extremely flat—for example, for TBG with twist angles larger than the magic angle. In the BM continuum model, more than 90% of the momentum states in the band have a localized wavefunction in the AA regions of the moiré cell. However, near the Γ point, this structure is reversed: the states become more extended and feature a suppression of charge density at the AA sites.

Although this momentum-dependent real-space structure was used as the basis for formulating the c-electron and f-electron decomposition in the THF model, it is not specific to that framework. Rather, it is a robust and very general feature of the non-interacting flat bands themselves. Our ability to measure δEH(k) with momentum resolution therefore provides direct access to these underlying wavefunction characteristics, offering a powerful new experimental probe of the topology-driven structure of the MATBG bands.

Extended Data Fig. 9a–g provides a step-by-step illustration of how the observed Hartree stretching arises and how it is directly connected to the wavefunction structure of the flat bands. Panel a sketches the real-space Hartree potential VH(r) that develops on electron filling. Because more than 90% of the flat-band wavefunctions in the BM model reside on the AA sites and have an f-electron Wannier function, the added carriers generate a Hartree potential that is strongly peaked at these sites and therefore highly repulsive for f-electrons. We emphasize that here we are discussing single-particle wavefunctions, so the c/f decomposition is fully valid in this context. As shown in panel b, the shape of VH(r) remains fixed as a function of filling, whereas its overall amplitude grows approximately linearly with ν (ref. 21).

Panel c plots the BM charge densities |Ψf(r)|2 (red) and |Ψc(r)|2 (blue). The f-electron density has far greater overlap with VH(r) than the c-electron density and thus experiences a substantially larger Hartree shift. This is captured schematically in panel d, in which both c-electron and f-electron energies increase linearly with ν but with a much steeper slope for the f-electrons.

Panel e translates these energy shifts into momentum space. States near the Γ point—dominated by c-character—shift only weakly with filling, whereas most of the states in the rest of the band—dominated by f-character—shift strongly. Because most of the flat-band states are f-like, the Fermi energy EF (dashed line) remains effectively pinned to the flatter regions of the band as ν varies.

Our experiment, however, measures energies relative to EF. Referencing to EF has the effect shown schematically in panel f: the strongly shifting f-electron states remain near zero energy, whereas the weakly shifting c-electron states acquire an apparent downward slope with filling. This behaviour matches precisely the trend observed in Fig. 4 (here we focus only on the Hartree component; further interaction terms responsible for Hubbard-like band splitting and cascades appear on top of this effect in the experiment).

Finally, panel g shows how the ν-dependent bands of panel e appear once referenced to EF: the bands exhibit a clear stretching near Γ—exactly as observed experimentally in Fig. 3a. This stretching provides compelling evidence for the distinct wavefunction character of states near Γ, a key prediction of the topological structure of the BM bands and a foundational ingredient of the THF model.

Toy model to understand the essence of the Dirac revivals

In Extended Data Fig. 10, we present a minimal toy model that separates the low-energy spectrum into two sectors: (1) light carriers with a Dirac-like dispersion (‘c-like’ electrons) and (2) heavy, nearly flat bands (‘f-like’ electrons). The aim is to capture, in a transparent way, the key experimental features, namely the cascade phenomenology and the Dirac revivals—without committing to microscopic details.

The components/assumption of the model are as follows:

1. 1.

   We assume two decoupled electronic components, which we call ‘c-like’ and ‘f-like’. They are not necessarily the original c-electrons and f-electrons but they do carry their two important aspects: light/heavy, different Hartree couplings.

2. 2.

   Specifically, we model the f-like electrons with a simplified quantum dot model, which assumes completely flat bands but with a phenomenological ‘lifetime’ broadening that gives a finite width to the energy bands. We compute the f-electron spectral function to capture the Hubbard bands filling evolution.

3. 3.

   We assume that there is a finite Coulomb coupling between c-like and f-like electrons, W, and that W < U.

The f-electron quantum dot model Hamiltonian is: H = Unf(nf − 1)/2, in which U is the effective charging energy and nf is the f-electron occupation. This quantum dot has eight flavours such that 0 ≤ nf ≤ 8. We compute the spectral function given by this Hamiltonian: \(A(\omega )={Z}^{-1}{\sum }_{n,m}\delta (\omega -{\varepsilon }_{n}+{\varepsilon }_{m})({e}^{-\beta {\varepsilon }_{m}}+{e}^{-\beta {\varepsilon }_{n}}) < l|{c}_{\alpha ,{\rm{f}}}^{\dagger }{|m > |}^{2}\), in which Z = Tr(e−β(H − μn)) is the partition function, β = 1/kBT and εm is the energy of the |m> eigenstates. We compute the spectral function A(ω, nf) and then add smearing in energy to mimic the lifetime effect seeing in the measurement. From the smeared spectral function, we compute the relation μf(nf).

As well as the quantum dot, we add Dirac (c-like) electrons. Let nf and nc denote the f and c fillings (n = nf + nc is fixed externally). The interaction term can be written as Wnfnc = Wnf(n − nf). The total energy is:

$${E}_{{\rm{tot}}}({n}_{{\rm{f}}};\,n)={E}_{{\rm{f}}}({n}_{{\rm{f}}})+{E}_{{\rm{c}}}({n}_{{\rm{c}}})+W{n}_{{\rm{f}}}(n-{n}_{{\rm{f}}}),$$

and the stationarity condition

$$\frac{{\rm{d}}{E}_{{\rm{tot}}}}{{\rm{d}}{n}_{{\rm{f}}}}=0$$

gives a simple self-consistency equation:

$${\mu }_{{\rm{f}}}({n}_{{\rm{f}}})-{\mu }_{{\rm{c}}}({n}_{{\rm{c}}})+W[n-2{n}_{{\rm{f}}}]=0,$$

with μf/c ≡ dEf/c/dnf/c. We solve this equation graphically to obtain (nf, nc) at each total filling n and then compute the corresponding spectral functions in momentum space for both sectors.

Extended Data Fig. 10a,b presents the calculated spectral functions for doping range from n = 0 to n = 1.5. Panels a and b show the f and c spectral functions in momentum space versus filling, reproducing: (1) the Hubbard-band evolution of heavy f-electrons and (2) the Dirac revival of light c-electrons. Specifically, in panel b, we track the Dirac point evolution: it first shifts downward with increasing ν owing to the Hartree potential and then when f states start to populate, it shifts upward (solid white guideline); this is precisely the ‘Dirac revival’ phenomenology and, within this model, we can see that it is driven by the W Coulomb term. The schematic in Fig. 4 is based on this calculation.

The above model does not include c–f hybridization. To show that the Dirac revivals physics is robust and appears also in the case of strong hybridization (γ/U > 1), we solve a more realistic THF model that also includes hybridization (details in the next section). Extended Data Fig. 10c plots the calculated filling factor dependence of the spectral function. The black line traces the position of the c-electrons Dirac point, clearly demonstrating that the Dirac revivals phenomenology, observed in the experiment, appears also in the presence of hybridization.

THF model reproducing the Dirac revivals

We consider a THF model38 that treats U, the on-site f–f interaction, non-perturbatively within the Hubbard-I approximation, while treating all other interactions within the self-consistent Hartree approximation. We use the following parameters for the calculations: v⋆ = −1 eV nm; v′⋆ = 366 meV nm; γ = −70 meV; M = 3.7 meV; U = 60 meV; W = 45 meV; Y = 45 meV.

To calculate the Green’s function of f-electrons and c-electrons within Hubbard-I, we first consider the non-hybridized f-electrons propagator:

$${G}_{{\rm{f}},\eta }^{0}({\bf{k}},\omega )=\left(\frac{\frac{1}{2}+\frac{{\nu }_{{\rm{f}}}}{N}}{\omega +\delta \mu +\frac{U}{2}+i{\tau }^{-1}}+\frac{\frac{1}{2}-\frac{{\nu }_{{\rm{f}}}}{N}}{\omega +\delta \mu -\frac{U}{2}+i{\tau }^{-1}}\right){I}_{2\times 2}.$$

νf = sign(ν)⌊|ν|⌋ is the integral filling of each f site relative to charge neutrality, within the non-hybridized theory, taken to be the total filling rounded towards zero. δμ = μf − ν · U sets the asymmetry between the Mott-band energies. N = 8 is the number of electronic states (flavours × orbitals) per f site and μf = μ − W · ⟨δnc⟩ is the effective electro-chemical potential felt by the f site owing to repulsion from c-electrons, with δnc measured relative to charge neutrality. τ is a finite quasiparticle lifetime introduced by hand for numerical stability. We take τ−1 = 1 meV.

The non-hybridized c-electrons propagator is given by

$${G}_{{\rm{c}},\eta }^{0}({\bf{k}},\omega )={[\omega +{\mu }_{{\rm{c}}}-{H}^{({\rm{c}},\eta )}({\bf{k}})+i{\tau }^{-1}]}^{-1},$$

in which H(c,η)(k) is the single-particle Hamiltonian of c-electrons in valley η as given by Song and Bernevig38 and μc = μ − W · ⟨δnf⟩ − V · ⟨δnc⟩ is the effective electro-chemical potential felt by the c-electron, with δnf measured relative to charge neutrality. Finally, the Hubbard-I propagator is given by

$$\begin{array}{c}[\begin{array}{cc}{G}_{{\rm{c}},\eta } & {G}_{{\rm{cf}},\eta }\\ {G}_{{\rm{fc}},\eta } & {G}_{{\rm{f}},\eta }\end{array}]({\bf{k}},\omega )={[\begin{array}{cc}{({G}_{{\rm{c}},\eta }^{0}({\bf{k}},\omega ))}^{-1} & {({H}^{({\rm{fc}},\eta )}({\bf{k}}))}^{\dagger }\\ {H}^{({\rm{fc}},\eta )}({\bf{k}}) & {({G}_{{\rm{f}},\eta }^{0}({\bf{k}},\omega ))}^{-1}\end{array}]}^{-1},\end{array}$$

with H(fc,η)(k) the f–c single-particle hybridization term as given by Song and Bernevig38.

The densities ⟨δnf⟩ and ⟨δnc⟩ are calculated by integrating the spectral functions over negative frequencies and subtracting their values at charge neutrality:

$$\langle \delta {n}_{\alpha }\rangle =4{\int }_{-\infty }^{0}{\rm{d}}\omega \int \frac{{{\rm{d}}}^{2}{\bf{k}}}{{A}_{\mathrm{BZ}}}{{\mathcal{A}}}_{\alpha }({\bf{k}},\omega )-{n}_{\alpha }^{0}.$$

The spectral function is given by

$${{\mathcal{A}}}_{\alpha }({\bf{k}},\omega )=-\frac{1}{\pi }{\rm{ImTr}}[{G}_{\alpha ,\eta }({\bf{k}},\omega )],$$

with α = f,c, \({n}_{\alpha }^{0}\) is the density at charge neutrality and the factor 4 corresponds to spin and valley degeneracy, for which we use the fact that the state is flavour symmetric by assumption. For each value of the total density, we solve for ⟨δnf⟩, ⟨δnc⟩ and the chemical potential μ self-consistently.

Comparing the energy spectrum at partially full and completely full flat bands

In Extended Data Fig. 11, we plot the spectroscopy data as in Fig. 4a and compare linecuts A and B, taken when the chemical potential lies inside and outside the flat bands, respectively (dashed lines in panel a). For partial filling of the flat bands, we observe a characteristic feature at ΔE = −15 meV, as indicated by the dashed black lines in linecut A at ν = 1.5. By contrast, when the flat bands are fully filled, we instead find a single broad peak, as shown in linecut B at ν = 4.3.

One possible interpretation of the ΔE feature is a single-particle splitting induced by strain. However, the strain in the regions in which we perform the measurements is generally small (Methods section ‘Spectroscopy of flat bands at different locations’), making it unlikely to account for the observed 15 meV splitting. Moreover, a single-particle splitting should also be present when the flat bands are fully filled. In panel c, we attempt to fit the broad peak in linecut B with two symmetric Gaussians constrained to have the same lifetime as at partial filling (FWHM of about 10 meV, comparable with the energy resolution) and a fixed splitting of 15 meV and obtain a poor fit. This further suggests that simple single-particle strain-induced splitting is unlikely to be the origin of the ΔE feature.

Discussions of possible ground states at integer fillings

At charge neutrality, we observed a semimetallic phase with states at Γ at the Fermi level. There are two physically distinct mechanisms to produce a semimetal: the thermally disordered Mott semimetal27,34,35,39 and the strain-induced or spontaneously C3 symmetry-broken semimetal29,52,53. The two states give very similar spectra, with Hubbard-like nearly flat bands away from the Γ point and a band touching at or near Γ. There should be small differences in the spectra between the two states—for example, the Dirac points of the strain-induced semimetal are generally not precisely at Γ and do not have to occur at the same energy—but we believe that these differences may be too small for us to observe within experimental resolution. The main factor limiting the resolution is probably the tip size, limiting the resolution in momentum space to discern the fine details of the dispersion near Γ.

At non-zero integer fillings, our experiment is not directly sensitive to symmetry breaking, as the tunnelling from the QTM tip is not sensitive to the electron spin and valley and we do not have real-space sensitivity that can identify the spatial modulation of the spectrum that occurs in intervalley coherent or Kekulé spiral states. Experiments do not show any obvious signatures of symmetry breaking in the momentum-resolved electronic spectrum. In particular, there is no clear gap in the electronic spectrum at the Fermi level at any density (Fig. 3a). Our results are incompatible with a gap-opening broken-symmetry state at integer fillings at the temperature of the experiment.

Transport experiments typically see quantum oscillations emanating from ν = 2 and −2. Our measurements, on the other hand, do not show any signatures of Fermi surface at these filling factors. It is important, however, to keep in mind the experimental conditions. Our measurements are performed at T = 4 K and zero magnetic field. By contrast, quantum oscillations are necessarily measured at finite magnetic fields and it is therefore difficult to exclude the possibility that the oscillatory features observed at finite field reflect a field-stabilized ground state that differs from the zero-field state examined here.