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

Magic-angle twisted bilayer graphene 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 non-trivial 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 non-trivial enhancement due to the interaction induced renormalization of the Bloch wavefunctions. While our modeling of superconductivity 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 superconductivity in the wide range of angles that occurs in the relevant range of fillings.

Introduction.-The origin of electron pairing in twisted bilayer graphene (TBG) has been at the heart of the discussion since the discovery of superconductivity [1][2][3] in magic-angle twisted bilayer graphene (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 superconductivity (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 non-interacting bandstructure. (Note that although different models for the non-interacting 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 magicangle [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 * cyprian@caltech.edu † debanjanchowdhury@cornell.edu the sole effect of the electron-electron (Coulomb) interaction is to renormalize the bare non-interacting bandstructure in a filling-dependent fashion (see Fig. 1a, b), while the attraction required for pairing stems from electronphonon 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 non-interacting 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 a HF approximation. However, the role of these renormalizations and especially the band "flattening" behavior (to be made precise below; see i) The non-interacting 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 non-interacting bandstructure corresponding to EF = 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.
the unit cell that reflects the effective triangular symmetry of the TBG lattice (Fig. 1e, f), 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 modify 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 superconductivity. 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 H (ξ,σ) 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 "non-interacting 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, Here δρ(r) tracks the density relative to that at charge neutrality, ρ CN (r), and V c (r − r ) 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 self-consistent 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. 1f), corresponding to electronic states nearK points of the mini-Brillouin zone. The non-uniform 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 theK 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 non-interacting bandwidth is 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 bandinversion 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)) does not contain any inter-flavour terms. We explicitly forbid any inter-flavour 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 . . . 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 magicangle. These parameters are kept fixed for all twist angles and as a result we do not capture the subtle latticerelaxation effects near magic-angle [27,28]. Note also that as a result of our adopted procedure, the location of 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 van Hove singularities. 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 Supp. Mat. Sec. 3 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 Supp. Mat. Sec. 3). 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. 1c, d. Most notably, we see that as the twist angle approaches 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 non-interacting model parameters as a result of the HF renormalization, c.f. 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 effects of interactions. For our ensuing discussion of phonon-mediated attraction, the non-trivial 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 phononmediated pairing. Earlier works have highlighted the importance of a purely phonon driven mechanism, within various approximations, when the electronic bandstructure is limited to the non-interacting 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 inter-action, after integrating out the phonons, has the form: where the projected density operators, Here iω is the fermionic Matsubara frequency and |k, {γ } 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 (= 12000 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 Hartree-Fock 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 Hartree-Fock 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 materials and demonstrate its effect on pairing tendencies in the Supplementary Figure 3.
Robustness of SC away from magic-angle.-Within the framework of Eq. 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 materials for further details).
Within our Eliashberg analysis of the "weak" electronphonon 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 non-interacting 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, 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][7][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 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 band-flattening 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 nonperturbative 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 = 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 negligible difference between the metric for the non-interacting (Fig. 3a) vs. HF-modified bands (Fig. 3b). The subtle band-flattening 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 HF-modified bands develop band inversions (Fig.1b), the modified Bloch wavefuctions 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 HFcorrected bandstructure. As a result, the integrated metric, M, is also strongly modified in the vicinity of magicangle (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 imporatant 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 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 materials for further discussion. A recent experimental work [48] highlights the role played by redistribution of Berry curvature in the MBZ as a keyelement driving appearance of new insulating states.
Discussion.-In this work, we have focused on HFdriven 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 weakcoupling 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 superconductivity 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 gate-induced screening on the interplay between electronelectron interaction induced modifications to the bands and phonon-mediated 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, G M = 4π/ √ 3L M , [11,13] (L M ≡ moiré period), the Fock term receives 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 band-inversion), 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 from a more conventional source of enhancement arising from suppression of the Coulomb repulsion.
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 Hartree-Fock analysis and the possible effects of including cascade-like transitions in the analysis here. The Hartree-Fock corrections are included at zero temperature, which is justified and self-consistent within our setup. Specifically, the leading temperature dependent corrections to the Hartree self-energy 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 selfenergy will not alter the results significantly. The origin of the Hartree corrections (and the associated bandflattening) 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. [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 inhomogenous charge distribution and the associated band-flattening 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 non-uniform 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 superconductivity, thereby providing further direct insight into the complexities of twisted bilayer graphene.
Supplementary materials for "Does filling-dependent band renormalization aid pairing in twisted bilayer graphene?" Let us review the continuum model used to capture the non-interacting bandstructure of TBG. The parametrisation of the non-interacting 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 non-interacting Hamiltonian, H 0 , we employ the continuum model introduced in Ref. 23. As described in Eq. 2 of the main text, the non-interacting Hamiltonian is given by where the explicit form of the elements of the operator H (ξ,σ) 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 H (ξ,σ) 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 , which denote the Dirac points K and K 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/2 sin(θ/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 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) = G C j ξ,n,k (G)e i(k+G)·r , where j = A 1 , B 1 , A 2 , B 2 labels 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 non-interacting 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 In the above expression denotes summation over occupied states for a given filling measured with respect to the CNP as explained in the main text (denoted . . . 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 As argued in Ref. [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 contribution is excluded since it is cancelled 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. . . . F of Eq.(9) from the main text, in a manner explained in Sec. 3 below. Importantly, the self-energy explicitly depends on the crystal momentum k due to the non-local 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.

Procedure for fitting TBG continuum model to scanning tunneling microscopy (STM) experiments
Here we elaborate further on the modeling and analysis introduced in Sec. 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 an experimentally inspired modelling of the qualitative behavior of the bandstructure and corresponding Bloch wavefunctions as a function of twist angle and filling, rather then provide a complete and full treatment of Hartree-Fock problem in TBG. We also note the large body of existing work in the literature that analyzes Hartree-Fock 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 non-interacting continuum model as introduced in the Sec. 1 for a given twist angle has two free parametersthe interlayer couplings u, u . 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 , corresponding to same-sublattice and oppositesublattice 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 as the twist angle is brought closer to the magic angle [27]. To fix u and u , 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 = 90 meV and u = 0.4u . As noted in the main text, as a result of this approximation scheme the magic-angle of the non-interacting 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 will not qualitatively alter our results. Specifically, the location of van Hove singularities 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 ).
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 parametrisation 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 Sec. 2, the Fock potential is non-local 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 non-local dependence with a constant, twist-angle dependent prefactor g(θ) with the characteristic energy scale of the Fock interaction being set by the Hartree potential V c (G − G ). We stress that despite 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 . In analogy to the Hartree potential, we limit the reciprocal momenta included in the self-energy to a difference G − G 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 qualative 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 twistangle-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 . . . 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 theK,K 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 Hartree-Fock 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 parametrisation of the SC calculation), where the convergence 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 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 theK point. For the Matsubara grid, in analogy with the procedure of Ref. [37] we employed an approximate scheme that consists of 10 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).

Comparison of Hartree vs. Hartree-Fock and role of electron-phonon umklapp
Here we present additional results for the calculation of T c as a function of filling and twist angle; we focus specifically on (i) a comparison of a non-interacting, Hartree-only and Hartree-Fock modified bandstructure on pairing tendencies, (ii) the role played by Umklapp processes, and (iii) the shape of the T c dome for negative fillings ν < 0.
In Supplemental Figure 1 we show the pairing dome profile as a function of twist angle and filling for three different "extents" of incorporating Coulomb interactions. Supplemental Figure 1a shows a pairing profile for a non-interacting bandstructure and wavefunctions which, as discussed in the main text, is sharply peaked near the magic angle and falls precipitously away from it. As expected of a weak-coupling calculation, the dome is peaked at the location of the vHs of the non-interacting bandstructure that remains fixed for all twist angles. When Hartree-only (Supplemental Figure 1b) and Hartree-Fock (Supplemental Figure 1c) corrections are introduced into the SC calculation, we find the appearance of a qualitatively similar-looking three-prong feature stemming from band-flattening physics analysed in the main text. The Hartree-only correction however presents an overall lower maximum T c near the magic angle as compared to the Hartree-Fock results. This simply is a consequence of the Hartree-only bandstructure featuring a largeΓ point inversion that increases overall bandwidth at the magic angle beyond that of non-interacting or Hartree-Fock (See also Supplemental Figure 2 for further comparison of non-interacting and Hartree-Fock results.). In Supplemental Figure 1, 2, as in the main text, we normalize the pairing temperature T c by its maximum value for the non-interacting bandstructure of T c,max = 0.0294 meV.
In Supplemental Figure 3 we study the effect of including umklapp phonon processes on the pairing temperature. As we argued in Ref. [37], in TBG due to the small size of the effective MBZ, scattering processes that involve an exchange of multiples of reciprocal lattice vectors (umklapp processes) are favourable. We denote a computation which incorporates an exchange of up to m reciprocal lattice vectors G M 1 , G M 2 as an "mG" calculation. Due to the nature of the moiré potential that couples both graphene layers to each other, to construct a low energy effective theory it is necessary to use several graphene states that are coupled by these reciprocal lattice vectors G. The specific number of such states imposes a critical value of m umklapp processes that can be exchanged above which there is no significant increase in pairing temperature -we refer interested reader to our earlier work Ref. [37]. We stress that again the highest m above which no further increase in pairing temperature (typically mG = 3G) occurs is well below any cut-off imposed by the finite truncation of the continuum Hamiltonian (which would correspond to mG = 15G as specified in the Methods section of the main text. Indeed this physics finds representation in our results of Supplemental Figure 3a-c, where inclusion of higher number of umklapp processes enhances pairing temperature. We also note in passing that this effect becomes less pronounced the larger is the twist angle as the MBZ size increases, thereby suppressing umklapp processes. Finally, note that Supplemental Figure 3a-c is evaluated for negative fillings to provide an explicit demonstration of behavior that is qualitatively similar to ν > 0 (i.e. we do not find discernible qualitative differences between results of Fig.2b and Supplemental Figure 1c beyond some small quantitative differences stemming from particle-hole asymmetry that is present in the TBG Hamiltonian).