Does filling-dependent band renormalization aid pairing in twisted bilayer graphene?

Magic-angle twisted bilayer graphene (MATBG) exhibits a panoply of many-body phenomena that are intimately tied to the appearance of narrow and well-isolated electronic bands. The microscopic ingredients that are responsible for the complex experimental phenomenology include electron–electron (phonon) interactions and nontrivial Bloch wavefunctions associated with the narrow bands. Inspired by recent experiments, we focus on two independent quantities that are considerably modified by Coulomb interaction-driven band renormalization, namely the density of states and the minimal spatial extent associated with the Wannier functions. First, we show that a filling-dependent enhancement of the density of states, caused by band flattening, in combination with phonon-mediated attraction due to electron-phonon umklapp processes, increases the tendency towards superconducting pairing in a range of angles around magic-angle. Second, we demonstrate that the minimal spatial extent associated with the Wannier functions, which contributes towards increasing the superconducting phase stiffness, also develops a nontrivial enhancement due to the interaction-induced renormalization of the Bloch wavefunctions. While our modeling of superconductivity (SC) assumes a weak electron-phonon coupling and does not consider many of the likely relevant correlation effects, it explains simply the experimentally observed robustness of SC in the wide range of angles that occurs in the relevant range of fillings.


INTRODUCTION
The origin of electron pairing in magic angle twisted bilayer graphene (MATBG) has been at the heart of the discussion since the discovery of superconductivity (SC) [1][2][3] in MATBG, while the phase diagram and the associated experimental phenomenology has continued to evolve dramatically 4 . At the time of writing of this letter, the following observations related to SC are universally accepted in MATBG: (i) multiple pockets of SC are present over an extended range of fillings, −4 < ν < 4 (ν ≡ electron filling in the moiré "flat" bands). The location of these SC regions are not simply tied to either the near vicinity of the correlation-induced insulators at commensurate fillings 2,3,5 , or to the van Hove singularities (vHs) associated with the noninteracting bandstructure. (Note that although different models for the noninteracting bandstructures may predict vHs to occur at different fillings, no single model predicts multiple vHs to be present at fillings corresponding to locations of the SC pockets.) (ii) The SC regions are more resilient to external screening and deviations away from magic-angle [6][7][8][9] , i.e., even when the sharp insulating gaps in the limit of low temperatures are no longer observable at the various commensurate fillings, SC continues to remain robust with only minor changes to the transition temperatures, T c .
Inspired by these experimental facts, we focus here on the following interesting theoretical scenario, where the sole effect of the electron-electron (Coulomb) interaction is to renormalize the bare noninteracting bandstructure in a filling-dependent fashion (see Fig. 1a, b), while the attraction required for pairing stems from electron-phonon interactions. We capture the effects of these renormalizations on both the bandstructure and Bloch wavefunctions at the level of a Hartree-Fock (HF) approximation. As explained below, we match various qualitative aspects of the filling-dependent, renormalized bandstructure to recent experimental observations 10 . Within this setup, we will demonstrate that the superconducting phase diagram as a function of density and twist-angle is markedly different from the one derived from a model of noninteracting bandstructure. With this plan in place, we are led to a number of important questions, that we address in this letter: (i) What controls the propensity towards pairing at angles away from magic-angle? (ii) To what extent does the bare electronic bandstructure influence the SC phase diagram and its pairing tendencies? (iii) What are the modified properties of the Bloch functions associated with the renormalized Hamiltonian and their possible effect on superconducting properties?
A number of recent theoretical works have focused on the role of bandstructure renormalizations in MATBG on the possible symmetry-broken insulating phases at commensurate fillings at the level of an HF approximation. However, the role of these renormalizations and especially the band "flattening" behavior (to be made precise below; see Fig. 1a-h) on the pairing tendencies has not been analyzed explicitly. In MATBG, as a result of the uneven real-space charge distribution within the unit cell that reflects the effective triangular symmetry of the TBG lattice (Fig. 1i, j), the Hartree corrections become prominent [11][12][13][14][15] . The exchange effects, which can lead to gap openings and change the topological properties, are directly accessible in transport experiments (e.g., Landau fan) [16][17][18][19][20] . On the other hand, the Hartree corrections only alter qualitative aspects of the bandstructure and the underlying Bloch wavefunctions, leaving the topological properties unaltered and thereby making them harder to detect in transport. Interestingly, these changes can be imaged directly in local-probe experiments 10 . Here we build on this recently seen mechanism relying on the Hartree-correction, that modifies the bandstructure and Bloch wavefunctions, and analyze their role on enhancing the tendency towards pairing in MATBG.
We note at the outset that we intentionally do not include the effects associated with the "cascade transitions" at integer fillings near magic-angle 21,22 , which though relevant for concrete experimental observations will inevitably complicate further the discussion of SC. As already indicated above, it is at present unclear to what extent the cascade is tied to the origin of SC; nevertheless, in what follows the doping dependence of SC near magic-angle will be modified in the vicinity of integer fillings where a cascade would be expected to occur. On the other hand, away from magic-angle where the effects of the cascade become less pronounced, the results for SC are less likely to change qualitatively (even though cascade-related signatures were observed down to low twist angles 6 ). For simplicity, we leave a careful analysis of the SC phase diagram, including the effects of the cascade transitions, to a future study.

Renormalized bands
We begin with the single-particle continuum Hamiltonian 23 : ( where the explicit form ofĤ ðξ;σÞ appears in the Methods section. The spinor, ψ γ , is written in the basis of (A 1 , B 1 , A 2 , B 2 ) sites of the original two layers (l = 1, 2) and we use the shorthand notation, γ ≡ {ξ( = ±1), σ( = ±1)}, for the valley/spin degrees of freedom. The real space integration is over a moiré unit cell Ω. In what follows, any reference to the "noninteracting model" corresponds to a calculation that is based solely on the eigenvalues and eigenstates of this Hamiltonian in Eq. 2.
The Coulomb interaction is given by, δρðrÞ ¼ X γ ψ y γ ðrÞψ γ ðrÞ À ρ CN ðrÞ: Here δρ(r) tracks the density relative to that at charge neutrality, ρ CN ðrÞ, and V c ðr À r 0 Þ is the Coulomb potential with a Fourier transform, V c ðqÞ ¼ 2πe 2 =εq. For reasons to be made clear below, the dielectric screening by the substrate (denoted ε) is treated as a free parameter. We approximate the above interaction term using a selfconsistent HF approximation as, where the many-body renormalization, Σ HF (ν), will lead to a modified electronic bandstructure due to either the Hartree (H H ) or Fock (H F ) terms, respectively. The Hartree correction is given by, where 〈. . .〉 H denotes a summation over occupied states measured from CNP (ν = 0) 11 . As a function of increasing doping relative to charge neutrality, there is a preferential buildup of charge at AA sites in real space (Fig. 1j), corresponding to electronic states near K points of the mini-Brillouin zone. The nonuniform spatial charge distribution generates an electrostatic potential that prefers an even redistribution of the electron density. In contrast, the real space charge distribution corresponding to electronic states near Γ point is more uniform in the unit cell. The effect of the electrostatic Hartree potential and the associated charge redistribution thus leads to a lowering of the energy of the electronic states near the Γ point compared to the energy of states near the K points ( Fig. 1a, b). The effect of the Hartree potential becomes increasingly pronounced as a function of decreasing twist-angle, especially near the magic-angle where the noninteracting bandwidth is Fig. 1 Hartree-Fock renormalized bandstructures. TBG bandstructures as a function of filling for a θ = 1.2°, and b θ = 1.06°after including the HF corrections. One of the important features is related to band flattening and eventual inversion at the Γ point of the MBZ. The energy dependence of the density of states (c, d, g, h), demonstrating maximal enhancement when the band gets flattened, is shown next to each bandstructure. TBG bandstructures including HF corrections as a function of twist angle for e ν = 3, and f ν = −3, respectively. In (a, b, e, f) we plot bandstructures for the ξ = −1 valley. i The noninteracting energy landscape in the extended zone scheme of an electron band for ξ = −1 at θ = 1.06 ∘ . One MBZ is shown as a white hexagon. With red energy contour we denote Fermi surface for ν = 3. j Charge density for a noninteracting bandstructure corresponding to E F = 4.92 meV corresponding to the red energy contour in (i) and a filling of ν ≈ 3. Note that even at this large filling majority of charge density is located on the AA sites of the effective triangular lattice giving rise to the large Hartree potentials. Here charge density is normalized by the highest charge density at full filling.
minimal. There is an increasing tendency towards band-inversion near the Γ point 13,14 , a feature that has not been observed in experiments till date 10 . However, it is important to note that the Fock term, H F , inherently acts against this tendency towards band-inversion via two key mechanisms 24 : (i) by increasing the overall bandwidth and (ii) by contributing an opposing correction to the self-energy as compared to the Hartree term, Eq. 7.
The Fock term, H F , is given by, Note that the Fock potential, unlike Hartree, does not contain a summation over valley/spin degrees of freedom and as a result of the block-diagonal nature of the noninteracting Hamiltonian (Eq. 2) that does not contain any inter-flavor terms. We explicitly forbid any inter-flavor terms to be generated spontaneously 24,25 , since our goal here is to focus on the qualitative changes to the band structure and not on determining the precise nature of the correlated insulators 26 . Finally, the notation h i F corresponds to a summation over occupied states; see Methods for a more detailed discussion of the subtleties and the various conventions adopted in earlier works regarding the Fock term. Our modeling of the bandstructure is motivated by recent experiments 10 ; see Methods for details. In particular, we determine the microscopic parameters for the model by matching our theoretical bandstructures to the experimental results sufficiently far away from the magic angle. These parameters are kept fixed for all twist angles and as a result, we do not capture the subtle lattice-relaxation effects near magic-angle 27,28 . Note also that as a result of our adopted procedure, the location of the magic angle is θ ≈ 1 ∘ , which is different from the value encountered most often in literature, θ ≈ 1.1 ∘ . At the same time, it is worth noting that in the continuum model, varying the ratio of the interlayer parameters does not drastically alter the location of vHs. It is thus possible, in principle, to disentangle the effects due to twist-dependent interlayer hopping ratio from those due to band-flattening physics (see Supplementary Materials for further discussion). For general agreement with the experimental results, we found it necessary to use a dielectric constant ε larger than that set by the substrate, in accordance with similar observations made in earlier studies 11,13,24 . In spite of these simplifying approximations, our modeling captures the qualitative behavior exhibited by the measured MATBG bandstructure as the twist angle is brought closer to the magic-angle condition (See ref. 10 and discussion in Supplementary Materials).The final renormalized bandstructures at fixed angles of θ = 1.20 ∘ , 1.06 ∘ are shown as a function of filling in Fig. 1a, b. Similarly, for a fixed filling, the bandstructures for increasing twist angles are shown in Fig. 1e, f. Most notably, we see that as the twist angle approaches the magic angle, the bandstructure becomes locally flat near the Γ point beyond a certain filling. At these fillings, the location of the vHs changes from that set by the noninteracting model parameters as a result of the HF renormalization, c.f. the density of states panels in Fig. 1. We note that ref. 29 has also studied the onset of band-flattening in MATBG and commented on its possible relevance for enhancing the effects of interactions. For our ensuing discussion of phonon-mediated attraction, the nontrivial band flattening and associated interaction-induced shift of the vHs will play a crucial role in determining the shape of the superconducting dome.

Phonon-mediated attraction
With the role of electron-electron interactions limited only to the corrections discussed above, we now take these bandstructures and accompanying wavefunctions to investigate phonon-mediated pairing. Earlier works have highlighted the importance of a purely phonon-driven mechanism, within various approximations, when the electronic bandstructure is limited to the noninteracting model [30][31][32][33][34][35][36][37][38][39] . We employ here the framework and notation of an earlier work by us, focusing exclusively on an intervalley, spinsinglet gap function with zero center of mass momentum for simplicity 37 . The effective electron-electron interaction, after integrating out the phonons, has the form: where the projected density operators, ρ ξ (q, iω), include the formfactors, Λ γγ 0 ðp; kÞ ¼ δ ξξ 0 δ σσ 0 p; fγg e iðpÀkÞÁr k; fγ 0 g . Here iω is the fermionic Matsubara frequency and k; fγ 0 g j i denotes a Bloch wavefunction of the mean-field Hamiltonian, which includes the HF renormalizations due to the Coulomb interactions. The phonon mediated interaction vertex is given by, where ω ph (q) = c s q is the acoustic phonon dispersion for graphene. The electron-phonon coupling constant, g ¼ D 2 =ρ m c 2 s , is related to the deformation potential, D = 25 eV, the speed of sound in graphene, c s = 12,000 m/s, and atomic mass density, ρ m = 7.6 × 10 −8 g/cm 2 40,41 . Here we have redefinedg ¼ g=N with a large-N prefactor to obtain a controlled theoretical limit in which the Eliashberg equations for the pairing gap function are asymptotically exact; see ref. 37 for details of our earlier large-N framework. In the large-N formulation, the HF corrections of interest to us here are subleading in 1/N 37 . On the other hand, partly inspired by the emerging phenomenological considerations pointing towards the importance of these corrections, we employ a two-step calculation procedure where the HF modifications to the bandstructure are included first at leading order in the interaction strength, followed by the second stage where the pairing computation is carried out using the Eliashberg approach ignoring vertex corrections.
The Eliashberg equation for the gap function, Δ(iω, k), is given by with a kernel K(. . .), where ∫dΩ p denotes integration over the angle between vector k and p for a fixed direction of k. Importantly, E ξ;p is the electronic dispersion including the density-dependent HF renormalization. Note that it is important to include a summation over the moiré umklapp processes in Eq. 13, that were shown to lead to an increase in the pairing scale as in ref. 37 . We summarize the role of the form-factors Λð Þ and the origin of this effect in the Supplementary Discussion and demonstrate its effect on pairing tendencies in Supplementary Fig. 3.

Robustness of SC away from magic-angle
Within the framework of Eqs. 12, 13, T c is determined by the temperature at which the linearized Eliashberg equation (Eq.12) has an eigenvalue of 1. One of the central results of this paper appears in Fig. 2, which shows T c for (i) a noninteracting model without HF corrections (Fig. 2a) and (ii) a model that includes the density-dependent HF corrections (Fig. 2b), for a range of fillings and twist-angles. The results are qualitatively distinct; in particular, the intricate structure for T c as a function of θ and ν in Fig. 2b is seemingly unrelated to properties of the noninteracting bandstructure. For the latter, T c peaks at the magic angle and then rapidly falls off with a varying twist angle. The peak of the SC dome as a function of ν is also pinned to be at the same filling. On the other hand, T c for the HF corrected bandstructure does not abruptly fall off with changing twist angle and its maximal value is not pinned at a fixed filling. The qualitative behavior for ν < 0, both with and without HF interactions, is similar to the corresponding results for ν > 0 (Fig. 2); the quantitative differences arise from the particle-hole asymmetry that is present between the electron and hole bands (see Supplementary Discussion for further details). Within our Eliashberg analysis of the "weak" electron-phonon coupling (controlled by large-N), T c is ultimately controlled by the density of states at the Fermi surface, N(E F ), as shown in panels Fig. 2 c, d for the same range of fillings and twist-angles. The interesting structure associated with N(E F ) in Fig. 2d arises from the Hartree contribution to the bandstructure. The phenomenon of band-flattening shifts the van Hove singularity in the density of states away from the location dictated by the parameters of the noninteracting band structure (bright feature indicated with arrows in Fig. 2c) and enhances the density of states at the Fermi level as shown in Fig. 1a, b. As a function of decreasing twist angle (especially approaching magic-angle), the filling beyond which we see an onset of band flattening near Γ point (see Fig. 1a, b) also decreases. This is evident from the "three-prong" like feature in Fig. 2d. Remarkably, the inclusion of the many-body renormalization to the bandstructure offers a plausible explanation for the robustness of SC over a broader range of twist angles, in line with experimental results 1,2,4,6-8 and unlike the prediction of the bare unrenormalized bandstructures.
Before proceeding any further, we remind the reader that the apparent appearance of a two-peak-like structure at magic angle ( Fig. 2b) will be masked by the onset of "cascade-like" transitions occurring near the integer fillings 21,22 ; these latter effects have not been included in this study. In particular, this will lead to a suppression of band-flattening and, more crucially, induce a sequence of flavor selective cascade transitions that will alter the profile of density of states (and thereby T c ). Importantly, the Hartree-induced band-inversion near Γ becomes suppressed for ν > 1 after a cascade transition. This in turn will likely increase the density of states near ν ≈ 2 for twist angles near the magic-angle compared to the results of our current analysis (Fig. 2d), possibly resolving the apparent contradiction with the experimental results 4 . We also anticipate that Cooper pairing will be unlikely for ν ≳ 3 due to the same underlying reasons.
In our discussion thus far, we have identified T c with the onset of a pairing amplitude. In two-dimensions, the actual T c is determined by the superconducting phase stiffness. Thus even though our simplified analysis predicts an enhanced tendency towards SC near ν~0 in the vicinity of magic-angle (Fig. 2b) and a peak that is shifting towards higher fillings at increasing twist angles (θ ≳ 1.15 ∘ ), the SC phase stiffness will drop precipitously as the bands become nearly filled (ν ≲ 4) or approach charge neutrality (ν~0). Nevertheless, we already see that bandflattening aids in the onset of the simplest phonon-induced attraction as a result of a renormalized density of states at the Fermi energy at larger twist angles.

Fubini-Study metric
As noted above, T c is ultimately determined by the superconducting phase stiffness. For large bandwidths, the higher electronic kinetic energy contributes to a larger phase stiffness. On the other hand, for narrower bands, the problem becomes inherently non-perturbative and the exact mechanism that leads to a finite phase stiffness is presently unclear. However, if the Wannier functions are non-localizable 42 (e.g., as in topological bands), or have a finite geometric extent, the local Cooper pairs can contribute to the phase stiffness even in a perfectly flat-band limit 43 . Within a BCS mean-field description of the projected interaction, the stiffness in the flat-band limit is proportional to the minimal spread of the Wannier functions [44][45][46][47] , which would also be the case for our weak-coupling computation. To investigate how the geometrical properties of the Bloch functions of the renormalized HF Hamiltonian evolve with density and twist-angle, we study the Fubini-Study (FS) metric, evaluated for the same flavor γ. The trace of the FS metric, Trðg αβ Þ, controls the minimal spread of the associated Wannier functions via M ¼ R k Trðg αβ Þ. The results are as shown in Fig. 3. The Bloch wavefunctions and the associated FS metric (Eq. 14) for the renormalized Hamiltonian, H 0 þ H c , undergo a qualitative change from the corresponding quantities for the single-particle Hamiltonian, H 0 , as twist angle approaches magic-angle. Sufficiently far away from magic-angle (e.g., θ = 1.20°), there is a negligible difference between the metric for the noninteracting (Fig. 3a) vs. HF-modified bands (Fig. 3b). The subtle bandflattening features present in the bandstructures, Fig. 1a, do not qualitatively alter the momentum dependence of the metric. Closer to the magic-angle (e.g., θ = 1.06°) wherein the HFmodified bands develop band inversions (Fig. 1b), the modified Bloch wavefunctions lead to an appearance of distinct new features in the momentum dependence of the metric (see orange arrows and dashed circles in Fig. 3c, d). The locations of these new features correspond to the local maxima of the HF-corrected bandstructure. As a result, the integrated metric, M, is also strongly modified in the vicinity of magic-angle (Fig. 3e) and might contribute towards an independent source of enhancement of the superconducting phase stiffness. Note that within the general weak-coupling setup, the contribution to the stiffness due to the enhanced spatial extent of Wannier-functions is independent of the underlying details of the pairing symmetry and phonon-mediated attraction.
It is important to note here that irrespective of whether the bands are topological or not, the minimal spatial extent of the Wannier functions is related to the integral of the trace of the FS metric. Regardless of the actual form, shape, and other properties of the Wannier functions, the strong non-monotonic enhancement of this quantity as a function of twist angle near the magic-angle is an important consequence of the band-flattening mechanism (Fig. 3c).
We note that while the presence of cascade can partly obscure the large enhancement of the integrated metric near magic-angle, c.f. Fig. 3e, we anticipate M to increase even when effects of the cascade are taken into consideration. The enhancement of M is tied to incipient band-inversions near the Γ point, hints of which have been observed in recent experiments 10 . Interestingly, the distribution of the Berry curvature in the MBZ is also affected due to the HF corrections even though the valley projected Chern invariant remains unchanged; see Supplementary Discussion for further discussion. A recent experimental work 48 highlights the role played by redistribution of Berry curvature in the MBZ as a key-element driving appearance of new insulating states.

DISCUSSION
In this work, we have focused on HF-driven corrections to the narrow electronic bands and their effect on phonon-induced pairing in the s-wave channel. However, phonons can also induce (weaker) attraction in other channels (e.g., d-wave) that will nevertheless exhibit qualitatively similar behavior as a function of filling and twist-angles (albeit with lower T c ), especially if the short distance part of the Coulomb repulsion suppresses s-wave pairing 33 . Beyond the weak-coupling limit considered here, it is difficult to speculate on the robustness of our conclusions and on the possibility of other competing instabilities in the absence of a non-perturbative analysis. Nonetheless, given the experimental evidence for filling-dependent band-flattening and our explicit demonstration of the dramatic many-body effects on the electronic bands and Bloch wavefunctions, it is clearly necessary to include them for a complete understanding of SC in TBG, and moiré materials in general. (See refs. 49,50 for a more recent discussion of the role of band-flattening in other settings).
Let us now comment on the subtle effects of an external gateinduced screening on the interplay between electron-electron interaction induced modifications to the bands and phononmediated SC for the same device without simultaneously varying twist-angle-a question that was partly addressed in a recent experiment 9 . Although the overall scale of both Hartree and Fock terms is controlled by the strength of the dielectric screening (ε), the absolute magnitude of Hartree and Fock contributions can differ as they are dominated by scattering processes on different momentum scales. While the Hartree correction is controlled only by momenta on scales comparable to and larger than the moiré reciprocal lattice vector, 11,13 (L M ≡ moiré period), the Fock term receives a contribution from small (less than MBZ size) momentum processes as well. As such, when the screening gate is at a distance, d ≲ 1/2G M , the Fock term will be more strongly suppressed relative to the Hartree term. Since the former is responsible for preventing band-flattening (and incipient bandinversion), we expect screening to enhance pairing tendencies. On the other hand, even if future experiments can resolve such changes in T c , it might be difficult to disentangle these effects It is worthwhile to comment on the regime of validity of some of the approximations we have made, specifically with regard to the temperature-dependent corrections in the HF analysis and the possible effects of including cascade-like transitions in the analysis here. The HF corrections are included at zero temperature, which is justified and self-consistent within our setup. Specifically, the leading temperature-dependent corrections to the Hartree selfenergy are small in T/E F , where E F is defined as the chemical potential difference of the minority carriers measured with respect to the bottom/top of the band. At weak-coupling, the T c obtained by us is already much smaller than E F , and including the additional corrections to the Hartree self-energy will not alter the results significantly. The origin of the Hartree corrections (and the associated band-flattening) is expected on more general grounds even at finite, but small, temperatures due to the inhomogenous charge distribution in real space.
The effects of the cascade transitions, which have been neglected, can be included, in principle, in a future extension of the framework discussed here. The actual implementation would inevitably have to rely on a detailed input of the precise ground state near the integer fillings (see e.g., refs. 25,26,51,52 )and the nature of the superconducting state in the vicinity of such a state. Our results for the FS metric are likely less sensitive to the cascades, since the resulting real space inhomogeneous charge distribution and the associated bandflattening should continue to be relevant.
Finally, given the extent to which interactions can modify the Bloch wavefunctions and the associated quantum geometric tensor, it is desirable to design experiments suited towards moiré materials that can probe these quantities directly (regardless of whether, e.g., M determines the actual phase stiffness). In principle, transport or optical measurements can be used to infer the distribution of Berry curvature [53][54][55][56] . Similarly, it would be interesting to extract the FS metric by analyzing corrections to the predictions of semiclassical equations of motion in the presence of nonuniform electric fields away from magic-angle 57 . Given that many of these techniques are not restricted to operate only at low temperatures, hopefully, future experiments can extract these quantities and probe the physics of band-flattening well above the temperatures associated with cascade transitions and SC, thereby providing further direct insight into the complexities of TBG.

TBG continuum model
Let us review the continuum model used to capture the noninteracting bandstructure of TBG. The parametrization of the noninteracting Hamiltonian and the mean-field treatment of the Coulomb interactions follows the procedure outlined in ref. 10 , which we reproduce here for completeness.
For the noninteracting Hamiltonian, H 0 , we employ the continuum model introduced in ref. 23 . As described in Eq. 2 of the main text, the noninteracting Hamiltonian is given by where the explicit form of the elements of the operatorĤ ðξ;σÞ are detailed below. The spinor, ψ γ , is written in the basis of (A 1 , B 1 , A 2 , B 2 ) sites of the original two layers (l = 1, 2) and we use the shorthand notation, γ ≡ {ξ( = ±1), σ( = ±1)}, for the valley/spin degrees of freedom. The real space integration is over a moiré unit cell Ω. The intralayer elements of theĤ ðξ;σÞ are where k is a momentum in the BZ of the original graphene layers, R φ ð Þ is the 2 × 2 two-dimensional counterclockwise rotation matrix that accounts for rotation of the BZs for the original graphene layers, and the signs ± in Eq. 16 correspond to the layers l = 1 and 2, respectively. For all twist angles studied in the paper, we use ℏv/a = 2.1354 eV as the energy scale for the Hamiltonians H ξl . In Eq. 16 above the k ⋅ p expansion is taken around the two vectors K ðlÞ 1 , K ðlÞ À1 , which denote the Dirac points K and K 0 of the two layers, respectively: Within the convention of ref. 23 the moiré superlattice BZ is defined using two reciprocal lattice vectors with L M ¼ a=2sinðθ=2Þ being the moiré effective lattice period and a = 0.246 nm corresponding to the lattice constant of graphene. We denote the reciprocal lattice vector length as The operator U ξ (r) in Eq. 15 encodes the interlayer hopping and is given by: We treat the interlayer couplings u and u 0 as fitting parameters for the band structure according to the procedure introduced in ref. 10 and summarized below. To determine the energy bands and the eigenstates of both the noninteracting model and the one that includes mean-field corrections, we take the Bloch wavefunction ansatz for valley ξ as Ψ j ξ;n;k ðrÞ ¼ where j = A 1 , B 1 , A 2 , B 2 label the spinor components, n is a band index, and k is the Bloch wave vector in the BZ of the original graphene layers.
In the above ansatz, the G sum runs over all possible integer combinations of the reciprocal lattice vectors G ¼ m 1 G M 1 þ m 2 G M 2 with integer m 1 and m 2 . In practice we constrain −15 ≤ m 1 , m 2 ≤ 15 thus ensuring that the mean-field band-flattening and pairing calculations do not suffer from any noninteracting model cut-off effects.

Mean-field treatment of interactions
In this section, we provide explicit forms for the mean-field Hartree and Fock potentials defined in the main text. The Hartree self-energy, H H of Eq. 7, when expressed in the basis used for the Bloch ansatz from Eq. 20, takes the form γ 00 ;k 0 ;G 00 ;i 00 C i 00 γ 00 k 0 G 00 ð Þ Ã C i 00 γ 00 k 0 ðG 00 À G þ G 0 Þ: In the above expression P 0 denotes summation over occupied states for a given filling measured with respect to the CNP as explained in the main text (denoted h i H ). As explained in refs. 11,13,25 , since the Hartree potential is controlled primarily by the contribution of the first-star of reciprocal lattice vectors, it is sufficient to consider combinations of G, G 0 satisfying G À G 0 ¼ mG M 1 þ nG M 1 with (m, n) = {(±1, 0), (0, ±1), (±1, ±1)}. As argued in refs. 13,25 , going beyond the first-star of reciprocal vectors does alter results drastically. We note however that this is one explicit cut-off used in the calculation. As usual, the G ¼ G 0 ¼ 0 contribution is excluded since it is canceled by the positive (jellium) ionic background.
Similarly, the Fock self-energy, H F , in Eq. (9) becomes where ∑″ denotes summation over occupied states, i.e., h i F of Eq. 9 from the main text, in a manner explained in Section 3 below. Importantly, the self-energy explicitly depends on the crystal momentum k due to the nonlocal nature of the Fock potential. This dependence imposes significant numerical complexity for self-consistent calculations, unlike the Hartree form of Eq. 21, since the self-energy due to the Fock term has to be determined separately for each momentum k.
C. Lewandowski et al.

Procedure for fitting TBG continuum model to scanning tunneling microscopy (STM) experiments
Here we elaborate further on the modeling and analysis introduced in Section 1, 2, specifically focusing on the relation to experiments. The specific fitting procedure follows closely the one adopted in ref. 10 and is meant to be an accompanying reference. Our key objective here is to provide experimentally inspired modeling of the qualitative behavior of the bandstructure and corresponding Bloch wavefunctions as a function of twist angle and filling, rather than provide a complete and full treatment of the HF problem in TBG. We also note the large body of existing work in the literature that analyzes HF contributions (e.g., refs. 13,[24][25][26]58,59 ). Finally, as already emphasized in the main text, we do not address here the detailed nature of cascade states, and the presence or absence of Fock-induced gaps at integer fillings. The noninteracting continuum model as introduced in Section 1 for a given twist angle has two free parameters-the interlayer couplings u; u 0 . Although their dependence on twist angle has been studied through ab initio methods 27 , here we choose a simpler approach intended to highlight the important interaction-driven qualitative changes to the band structure. We assume that the two parameters u and u 0 , corresponding to samesublattice and opposite-sublattice interlayer tunneling, have fixed values for all twist angles. This approximation misses the subtle role relaxation physics plays on increasing the ratio of these parameters η ¼ u=u 0 as the twist angle is brought closer to the magic angle 27 . To fix u and u 0 , we focus on the measurements at the largest available angle of θ = 1.32 ∘ in ref. 10 , where the role of interactions is least important. By matching the measured Landau-level (LL) spectrum to that obtained numerically from this continuum model, we fix u 0 ¼ 90 meV and u ¼ 0:4u 0 . As noted in the main text, as a result of this approximation scheme the magic-angle of the noninteracting model occurs at θ ≈ 0.99 ∘ , which differs from the typical values θ ≈ 1.1 ∘ quoted in the literature. At the same time, however, the twist-dependent ratio of η ¼ u=u 0 will not qualitatively alter our results. Specifically, the location of vHs in each band does not change with η for values of η ≲ 0.8 (See also discussion in ref. 60 ) Combining this observation with the property that the peak of the pairing dome (within weak-coupling limit) is dictated by the location of the van Hove singularity, we argue that band-flattening corrections to pairing can be disentangled from any of those coming from varying ðu; u 0 Þ.
We now proceed to parametrize the strength of the dielectric screening ε in the Coulomb potential V c ðqÞ ¼ 2πe 2 =εq. Nominally the value set by the substrate, i.e., typical of an hBN encapsulated graphene is ε ≈ 5. This value however massively overestimates the role of Hartree and Fock processes, leading to band structures with large Γ point inversions that are not observed experimentally. To overcome this unwarranted behavior, earlier works 11,13,24 use a range of values for ε ranging from 5 to 66. In the same spirit, we choose ε = 15 to quantitatively capture the following three experimental characteristics seen in the LL spectra of ref. 10 : (i) the energy spacing between LLs arising from the Γ point band structure at θ = 1.32 ∘ (Fig. 1 in ref. 10 ), (ii) the energy spacing between the highest energy LL from the flat band and the lowest energy LL from the dispersive bands at θ = 1.32 ∘ (Fig. 1 in ref. 10 ), and (iii) the critical angle at which largest energy LL from the flat-band joins the vHs (Fig. 2 in ref. 10 ). These criteria are met with a choice of ε = 15, which we then keep constant for all values of θ. With this parametrization, it was found in ref. 10 that including Hartree-only correction adequately captures experimental observations at θ = 1.32 ∘ , suggesting that the Fock term plays a subdominant role at least at this large twist angle.
As mentioned in the previous Section 2, the Fock potential is nonlocal and thus carries a high computational cost, making a parameter sweep like that in Fig. 2 prohibitive. Motivated by the physical intuition that the role of Fock is to oppose the Hartree potential, further substantiated by recent works 24 , we approximate Eq. 22 as In particular, we replace the nonlocal dependence with a constant, twistangle dependent prefactor g(θ) with the characteristic energy scale of the Fock interaction being set by the Hartree potential V c ðG À G 0 Þ. We stress that despite the similarity to a local approximation, e.g., V c ðrÞ / δðrÞ, the form of Eq. 23 carries the additional dependence on reciprocal momenta G; G 0 . In analogy to the Hartree potential, we limit the reciprocal momenta included in the self-energy to a difference G À G 0 residing in the first star of reciprocal lattice vectors. Unsurprisingly, the Fock potential so constructed acts in opposition to the Hartree potential due to the overall minus sign above. We reiterate here that this is a phenomenological approximation, which nonetheless captures qualitative behavior of the Fock term reported by other authors, e.g., refs. 13,24,25 . Namely, it demonstrates how the Fock term can counteract band inversion stemming from the Hartree term; at charge neutrality point it also predicts an increase in the van Hove to van Hove separation and a broadening (splitting) of each flat-band's van Hove peak-both features seen in STM experiments 10,61,62 . For a further discussion of this approach together with other possible mechanisms (e.g., strain) that can similarly prevent Hartree-driven band inversions we refer to ref. 10 . In ref. 10 several forms for the prefactor of the Fock interaction g(θ) were considered. It was observed that a twist-angle-independent prefactor does not manage to capture the qualitative band structure behavior seen in experiments: A constant value that is too large prevents unphysical band inversion near the magic angle, yet it also prevents band flattening for ν = 3, 4 seen in the experiment at intermediate angles that could be recovered with a Hartree-only analysis; A constant value that is unusually small allows for band flattening seen in experiments, but does not prevent unphysical band inversions near the magic angle. To remedy these issues in ref. 10 , the following g(θ) dependence was chosen which we use in our analysis as well: we take g(θ) = 0 for θ ≥ 1.14 ∘ and θ ≤ 0.84 ∘ , whereas for 0.84 ∘ ≤ θ ≤ 1.14 ∘ range we assume a triangular profile with a maximum of g(θ = 0.99 ∘ ) = 1.875.
A final element of the Fock term description is the specification of what bands are included in the summation ∑″ in Eqs. 22 and 23 (or as schematically indicated with h i F in the main text). Several summation schemes are present in the literature with the key difference being whether only flat bands 13,25 or also dispersive bands 26,59 are included. Qualitatively, at ν = 4 the Fock term arises predominantly from the flat bands; on the other hand, at ν = −4 the dominant contribution arises from the dispersive valence bands. The contribution from the dispersive valence bands is also opposite in sign and effect to that of the flat bands as seen in ref. 24 . We verified that this behavior holds within our approximation of Eq. 23, but due to the local form of the Fock potential contribution from dispersive bands at ν = −4, one finds gap opening near the K; K 0 points beyond what is experimentally plausible. To qualitatively capture the above sign trend in the Fock term while mitigating the preceding gap issue, we include in the summation ∑″ all flat band states both at ν = 4 and ν = −4, but change the sign of the contribution at ν = −4. For other fillings, motivated by ref. 24 , we interpolate our solution as explained in the next paragraph.
As the contribution of the Hartree and Fock terms to the self-energy depends linearly (in the absence of cascade and gap opening at the CNP) on filling 13,14,24 , for ease of computation we linearly interpolate between the solutions at ν = 4 and ν = −4. We verified that the results do not qualitatively change if a self-consistent approach is used at every filling. To that end, the mean-field interacting HF Hamiltonians for filling ν are taken to be where H H ðνÞ, H F ðνÞ denotes Hartree or Fock correction for a filling ν. For the bandstructures and calculations presented in this paper, we evaluate the solutions at ν = 4 and ν = −4 self-consistently until convergence is reached. We set the convergence threshold for all self-consistent parameters as 0.1% total relative error (difference between successive self-consistent steps). A grid of 441 k-points was used for the analysis of the self-consistent potentials (see next section for parametrization of the SC calculation), where the convergence was reached after a few (typically in less than 5) iterations.

Numerical solution of the linear gap equations
Here we provide additional technical details of the Eliashberg procedure we adopt in this paper. We note again that we have used the same framework as detailed in our earlier work ref. 37 . For a given twist angle and filling we start the calculation by pre-computing a 2D MBZ mesh of points, with their associated Bloch wavefunctions and energies. We note that because of the HF-corrections that enter at a mean-field level into this calculation, a different 2D mesh of points has to be precomputed for C. Lewandowski et al. each filling. In the calculations, we use a mesh of 10201 MBZ points for each ν and θ. To carry out an angle-average of the kernel from Eq. 13, we first fix a particular direction of vector k (upon verifying that the conclusions are not dependent on the specific direction) and then identify all p points that are of magnitude p (within a resolution admitted by the mesh). We then estimate the angular average Eq. 13 by averaging over these p points-in practice,~50−100 points are used for each pair of k and p momentum values.
To determine the critical temperature, we seek the temperature T for which Eq. 12 has an eigenvalue of unity. In practice, we carry out a bisection method search for a T giving an eigenvalue within ±0.01 of unity.
In the calculations, we use a linearly spaced grid of 30 k points ranging from the center of the MBZ Γ to the K point. For the Matsubara grid, in analogy with the procedure of ref. 37 we employed an approximate scheme that consists of ten first Matsubara frequencies followed by 20 linearly spaced frequencies starting from the 11th Matsubara frequency to the UV cutoff, where the cutoff is chosen to be 8W (W is the flat-band bandwidth).

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.