Nodal band-off-diagonal superconductivity in twisted graphene superlattices

The superconducting state and mechanism are among the least understood phenomena in twisted graphene systems. Recent tunneling experiments indicate a transition between nodal and gapped pairing with electron filling, which is not naturally understood within current theory. We demonstrate that the coexistence of superconductivity and flavor polarization leads to pairing channels that are guaranteed by symmetry to be entirely band-off-diagonal, with a variety of consequences: most notably, the pairing invariant under all symmetries can have Bogoliubov Fermi surfaces in the superconducting state with protected nodal lines, or may be fully gapped, depending on parameters, and the band-off-diagonal chiral p-wave state exhibits transitions between gapped and nodal regions upon varying the doping. We demonstrate that band-off-diagonal pairing can be the leading state when only phonons are considered, and is also uniquely favored by fluctuations of a time-reversal-symmetric intervalley coherent order motivated by recent experiments. Consequently, band-off-diagonal superconductivity allows for the reconciliation of several key experimental observations in graphene moiré systems.

Tunneling conductance measurements taken within the superconducting state reveal V-shaped density of states (DOS) [48,49] which can become U-shaped at other electron concentrations [49].Setting aside the possibility of thermal fluctuations as origin [50], this is most naturally interpreted as a transition from nodal to fully gapped superconductivity.For a consistent microscopic theoretical understanding, this provides the following challenges: (i) electron-phonon coupling-a widely discussed [33,[35][36][37][38][39][40] pairing mechanism in TBG and TTG-will typically mediate an entirely attractive interaction in the Cooper channel, with leading pairing state that transforms trivially under all symmetries and is thus fully gapped [51,52].(ii) Even when the low-energy interactions favor an irreducible presentation (IR), e.g., E of C 3 , with nodal basis functions (p-or d-wave), the generically fully gapped chiral configuration wins over the nodal nematic one within mean-field.(iii) Even if we assume that the nodal state is energetically favored, e.g., due to sig-nificant corrections beyond mean-field [27,[53][54][55], one is still left to explain why there is a transition to another, fully gapped superconductor upon changing the filling.
In this work, we show that the combination of flavor polarization and the representations of the symmetries in the flat bands of TBG and TTG allow for pairing channels that are completely off-diagonal in the flat bands and that such band-off-diagonal states can naturally reconcile all three key challenges (i-iii).More specifically, we find two distinct band-off-diagonal states: one of them transforms under the trivial representation A of the system's point group C 6 (or one of A 1,2 of D 6 if we set the displacement field to zero) but can nonetheless have symmetryprotected nodal lines, akin to Bogoliubov Fermi surfaces discussed in [56,57], see Fig. 1(a-c) for an intuitive visual explanation.The surprising possibility of the existence of such Bogliubov Fermi surfaces without external magnetic field is unique to twisted graphene systems in that it follows as a direct consequence of both the symmetry and relative flatness of their normal-state bands.The second off-diagonal state transforms under a two-dimensional IR (E 2 of C 6 ).Its associated chiral state, E 2 (1, i), which is favored in mean-field over the nematic one, has the unique property of exhibiting nodal lines or being fully gapped depending on the filling fraction, even when the order parameter is kept fixed.We supplement our general symmetry arguments and phenomenological models with Hartree-Fock (HF) calculations on the continuum model, studying a variety of different pairing mechanisms.We find that nodal band-off-diagonal pairing is favored by the optical A 1 and B 1 phonon modes and by fluctuations of a time-reversal symmetric intervalley coherent (T-IVC) state (the T-IVC state has Kekulé order on the graphene scale [58][59][60]).Evidence for the former has been provided by a recent photoemission study [61] while evi-dence for the latter comes from recent STM experiments [7].Furthermore, also fluctuations of a time-reversalsymmetric sublattice polarized state (SLP+) are attractive in the band-off-diagonal channel (see Table II for a formal definition of the order parameters).We also show that fluctuations of both T-IVC and of a nematic, time-reversal symmetric IVC order [62] favor either the band-off-diagonal A or an E 1 state with band-diagonal components, which may also be nodal; the winner is determined by the relative amount of nematic IVC and T-IVC fluctuations.

A. Possible pairing states
Let us begin by classifying the superconducting instabilities in graphene moiré systems in the limit where the low-energy bands are spin polarized but allowing for multiple bands.We denote the spinless low-energy fermionic creation operators by c † k,α,η with momentum k in valley η = ±, and of band index α labeling the upper (α = +) and lower (α = −) quasi-flat bands.As a result of twofold rotational symmetry, C 2z , along the out-of-plane (z) direction or effective spinless time-reversal symmetry, Θ, the non-interacting band structure ξ k,α,η obeys ξ k,α,η = ξ −k,α,−η ≡ ξ η•k,α and intervalley pairing is expected to dominate.A general pairing order parameter in the inter-valley channel couples as where the order parameter ∆ k,η = −∆ T −k,−η is a matrix in band space.The physical spin texture of the superconductor is entirely determined by the form of the underlying normal state's polarization: if the spins are aligned in the two valleys, the superconductor is a nonunitary triplet, while anti-alignment [24,28] leads to a singlet-triplet admixed state [13,27,28].In both cases, all of the following states are well defined, with the aforementioned spin structures and symmetries given by appropriate combinations of spinless operations and spin rotations (see Appendix A1).
We will classify the pairing states according to the irreducible representations (IRs) of the system's point group D 6 , which is generated by six-fold rotations (C 6z ) along the z axis and two-fold rotation symmetry (C 2x ) along the in-plane x axis.Note a displacement field (D 0 ̸ = 0) breaks the in-plane rotations leading to the point group C 6 .Importantly, all IRs of D 6 and C 6 are either even or odd under C 2z .Choosing the phases of the Bloch states such that C 2z acts as c k,α,η → c −k,α,−η , it holds The last column shows which states merge when D0 ̸ = 0, reducing the point group from D6 to C6.
This immediately implies that the pairing states in all IRs even under C 2z (A 1 , A 2 , E 2 of D 6 ) must be antisymmetric in band space and, thus, entirely band offdiagonal, whereas the order parameters of the other IRs (B 1 , B 2 , E 1 ) are symmetric and can contain both band-diagonal and band-off-diagonal components.While superconducting order parameters with finite band-offdiagonal components are rather common in multi-band systems, the existence of pairing states that are constrained to be entirely band-off-diagonal is rather unique and follows from the combination of C 2z symmetry and the spin polarization in the normal state.Importantly, this is unaffected by strain or nematic order breaking C 3z as long as C 2z remains, which guarantees that there are IRs with entirely band-off-diagonal order parameters.
Choosing the phase conventions of the Bloch states such that C 2x and C 3z act as c k,α,η → (σ z ) αα c (kx,−ky),α,η and c k,α,η → c C3zk,α,η , respectively, the resulting candidate order parameters are summarized in Table I.Note that a momentum-independent representation of C 2x must be σ z due to the bands' eigenvalues at the Γ-M line, which in turn are connected to the topological obstruction of the flat bands [63].The reality (Hermiticity) constraint in Table I on χ, X, and Y ( χ, X, and Ŷ ) comes from the residual spinless time-reversal symmetry Θ of the normal state [64,65].The two two-dimensional IRs E 1,2 are each associated with three pairing statestwo nematic phases E 1,2 (1, 0), E 1,2 (0, 1) and one chiral state E 1,2 (1, i).

B. Spectral properties
We here have the rather unique situation that there are pairing channels, associated with the IRs A 1,2 and E 2 , where the pairing is constrained by C 2z to be entirely band off-diagonal.One immediate very unusual consequence is that the superconducting order parameter transforming under the trivial representation (A 1 ) has a symmetry-imposed line of zeros along the Γ-M line, and hence a nodal point in the spectrum.This is related to the topology-induced non-trivial representation of C 2x in band space.We refer to Ref. 39 for the discussion of other topological nodal points for pairing in obstructed TBG bands.As we will show next, band-off-diagonal pairing leads to additional unusual spectral properties with far reaching consequences for graphene moiré systems.To this end, consider the following effective Hamiltonian, , where the scalar function ∆ k describes the form of pairing.We will here study two cases which are conventionally considered to be fully gapped, (i) a momentumindependent "s-wave state" (A 2 or A pairing in Table I) where ∆ k = ∆ 0 and (ii) a "chiral p-wave" state, or more precisely an E 2 (1, i) state, where ∆ k = ∆ 0 (X k + iY k ) with (X k , Y k ) being smooth, MBZ-periodic functions transforming as (x, y) under C 3z .Furthermore, we parameterize the dispersion, ξ η•k,α , of the two flat bands (α = ±) in valley η = ± as ξ k,α = ϵ k − µ + α δ k , where ϵ k and δ k are C 3z (and, for D 0 = 0, C 2x ) symmetric functions.
The Bogoliubov spectrum of H σy has four bands, given by ±δ k ± (ϵ k − µ) 2 + |∆ k | 2 .Consequently, the excitation gap at momentum k reads as which is shown in Fig. 1(d), and therefore exhibits nodes where As long as the band structure has Dirac points, there are points k D in the Brillouin zone with δ k D = 0, associated with the blue cross in Fig. 1(d).Furthermore, for a metallic normal state, µ must be within the bandwidth and, hence, there must be a region R in momentum space where |δ k | > |ϵ k −µ|.For the momentum-independent A 2 state, ∆ k = ∆ 0 , this implies that there exists ∆ c 0 > 0 such that there is k * ∈ R with parameters (such as the blue circle) above the red solid line in Fig. 1(d) as long as |∆ 0 | < ∆ c 0 .By continuity, this means that there must be a nodal point on any line connecting k D and k * .Consequently, for µ within the bandwidth and δ k D = 0 for some k D , the A 2 will always have a nodal line if |∆ 0 | is sufficiently small, consistent with the intuitive picture based on the Bogoliubov spectrum in Fig. 1(a-c).
We illustrate this further in Fig. 1(e) using a toy model with T .This leads to the second unexpected conclusion that, for any pairing mechanism, including conventional electron-phonon coupling, the leading instability either has nodal lines in a finite region below T c or transforms non-trivially under the symmetries of the normal state.For electron-phonon pairing (or pairing mediated by the fluctuations of any time-reversal-symmetric order parameter [52], such as the T-IVC state) this is particularly unexpected since it is generally believed to always lead to a fully gapped state that transforms trivially under all symmetries.In fact, this can be proven in general terms [51,52], even for spinorbit-split Fermi surfaces and beyond mean-field theory [52].The crucial difference to these works, however, is that spinfull time-reversal is broken in our case such that the Fermi-Dirac constraint is inconsistent with a nonsign-changing, band-diagonal pairing state.This leads to the unique situation that although electron-phonon coupling will lead to entirely attractive interactions in the Cooper channel, the superconducting energetics is frustrated: the dominant pairing state is determined by whether the energetic loss due to non-resonant band-offdiagonal Cooper pairs (A 2 pairing) or the costs from sign changes of the order parameter (such as B 1 ) are less harmful.We will demonstrate this explicitly by a model calculation in Sec.II E below, where either A 2 or B 1 is dominant, depending on the form of the electron phonon coupling.
Let us first, however, discuss the general spectral properties of the "chiral p-wave" state which is canonically expected to be fully gapped as long as the Fermi surfaces do not cross the zeros of X k + iY k .Three of these zeros have to be at the Γ, K, and K ′ points as a consequence of C 3z symmetry.In the absence of fine-tuning, X k + iY k will have vortices at these points with vorticity v = +1.As can be seen in Fig. 1(f), where we show the phase of X k + iY k using an admixture of the two lowest-order terms, the net vorticity of +3 at these highsymmetry points has to be compensated by anti-vortices at generic momenta.The lowest possible number is three C 3z -related vortices, which appear near the M points in Fig. 1(f).If it holds |δ k | > |ϵ k − µ| at any of these zeros k = k j , we obtain a point above the red line in Fig. 1(d) and, thus, a nodal point along any contour between that k j and k D ; as opposed to the A 2 state, this holds irrespective of the value of ∆ 0 and therefore all the way to zero temperature.In summary, we find that also the E 2 (1, i) "chiral p-wave" state is not generically fully gapped but instead will exhibit a nodal line encircling any zero k j of X k + iY k with |δ kj | > |ϵ kj − µ|.This leads to an interesting filling dependence of the superconducting gap, as we illustrate in our toy model in Fig. 1(g at the vortices at Γ, K/K', and near M. Depending on µ, D k is positive only near the Γ point or only in a region surrounding the vortices close to the M points, leading to nodal lines encircling Γ and near the M points, respectively, as shown in the inset of Fig. 1(g).These regimes are separated by a fully gapped region where D k < 0 for all k, which could explain the fully gapped to nodal transition seen in tunneling experiments [49] when the filling fraction is changed.Note that D kj = −|ϵ kj − µ| ≤ 0 for k j at the K and K' points.In Fig. 1(g), D K = D K ′ vanishes close to the top of the band, which simply means that the Fermi surfaces cross the K, K' points and the superconductor has nodal points for this fine-tuned value of the chemical potential.

C. Fluctuation-induced pairing
Having discussed the unique energetics of pairing and spectral properties of the resulting superconductors in spin-polarized quasi-flat-bands with Dirac cones on a general level, we next study these aspects more explicitly by solving the superconducting self-consistency equations in the flat bands common to alternating-twist graphene systems.We will start with pairing induced by fluctuations of a nearby symmetry-broken phase.To this end, we will couple the low-energy electrons introduced in Eq. ( 1) to a collective bosonic field ϕ j (q) = ϕ † j (−q) via where the Hermitian matrices λ j capture the nature of the correlated insulating phase; we here choose and normalize λ j such that (λ j ) 2 = 1.Both for twisted bi- [9] and trilayer graphene [14,15,29], the stable phases emerging out of the U (4)×U ( 4) [9] manifold in the chiralflat (decoupled) limit are natural candidates.Integrating out the bosonic modes, we obtain an effective electronic interaction which in the for superconductivity relevant intervalley Cooper channel reads as with vertex t ϕ = ±1 encoding whether the order parameter is even or odd under time-reversal, Θϕ j (q)Θ † = t ϕ ϕ j (q), and χ q > 0 denoting the (static) susceptibility of ϕ j .Before discussing numerical results for the full model, we first focus on perfectly flat bands.In this limit, the leading superconducting instability within mean-field theory is given by the largest eigenvalue of V in Eq. ( 6) viewed as a matrix in the multi-index (η, α, β).Furthermore, if there is an anti-symmetric, valley-off-diagonal matrix D obeying (see Methods) the associated leading superconducting order parameter in Eq. ( 1) is given by (∆ k,η ) α,α ′ = δ k (Dη x ) α,η;α ′ η with δ k > 0; here η j denote Pauli matrices in valley space and the precise form of δ k is determined by χ(q).

D. T-IVC fluctuations
Motivated by recent experiments [7] providing direct evidence for T-IVC order, we start with T-IVC fluctuations as a pairing glue.In the U (4) × U (4) symmetric limit, the T-IVC state is associated with λ j = σ 0 η j , j = x, y, within our conventions.Since t ϕ = +1, we are looking for Dη x that commutes with λ j .Interestingly, there is a unique anti-symmetric, valley-off-diagonal matrix D ∝ σ y η x with that property, implying that the leading pairing state has the form ∆ k,η = σ y δ k , δ k > 0. This is exactly the A 2 state in Table I, which, as discussed above, will have nodal lines at least in the vicinity of T c when a finite band dispersion is taken into account.Intuitively, the fact that A 2 pairing is favored can be understood by noticing that the valleyoff-diagonal form of λ j leads to an attractive interaction across the valleys, which penalizes the B 1 state with its sign change between the two valleys.In fact, it holds To go beyond the flat-band limit, we solve the superconducting mean-field equations numerically.We take the flat TBG bands from the continuum model [66] as the starting point.To capture the spin polarized normal state, we supplement it with Coulomb repulsion and a perform HF calculation (see Appendix A for details).As can be seen in the resulting band structure shown in Fig. 2(a) with interaction renormalization assuming filling fraction ν = 2, this not only pushes one spin flavor below the Fermi level but also induces significant band renormalizations.For our subsequent study of superconducticity, we project onto the two bands at the Fermi level and associate them with the creation operators c k,α in the interactions in Eqs. ( 4) and (5).In our numerical computations, we choose where A m is the real space area of a moiré unit cell, and take α = 0.05 for concreteness, although we checked our main conclusion do not crucially depend on this form.In all of our numerics, we work at doping ν = 2.5.
As expected, we indeed find that the A 2 state dominates, both right at the critical temperature T c , obtained from the linearized gap equation, and at T = 0 as we show by iteratively solving the full self-consistency equation (see Appendix C).One crucial effect of the finite dispersion and splitting between the bands is that a finite interaction strength, V > V c,1 , is required to stabilize the superconducting phase, as can be seen in the plot of T c in Fig. 2(b).Superconductivity ceases to be a weak-coupling instability as the Bloch states (k, α, η) and (−k, α ′ , −η) are not degenerate for α ̸ = α ′ , cutting off the logarithmic divergence known from BCS theory.The quasi-particle spectrum and order parameter of superconductivity from T = 0 numerics are shown in Fig. 2(c,d).In accordance with our general discussion above, we observe that the order parameter only has finite components proportional to σ y , which do not mix with the band-even contributions ∝ σ 0,x,z as a result of C 2z symmetry.Furthermore, it does not change sign as a function of k and, for sufficiently small V but still with V > V c,1 , the nodal lines in the superconducting spectrum persist all the way to T = 0, while the nodal line is gapped out at low The interaction-strength-dependence of the superconducting gap can be more clearly seen in Fig. 2(e), where we show the DOS for the self-consistent solution at T = 0.For large V , the superconductor becomes fully gapped at T = 0, leading to a U-shaped DOS.With smaller V , the magnitude of the order parameter decreases and the superconductor eventually exhibits nodal lines, as explained above.In the regime just before these nodal lines appear, there is an increase in the DOS near the Fermi level, roughly when the order parameter and the maximal band splitting are comparable, leading to a V-shaped DOS (green line).The lifetime parameter used to compute the DOS is 0.3 meV; this choice was based on our k-grid spacing.While it is not necessarily small with respect to the tunneling gap (which vanishes at V c,2 ), it is small with respect to ∆(k) which is of order 5 meV just as the state is becoming fully gapped for our choice of normal state.This behavior of the DOS with interaction strength may offer a natural explanation for the Ushaped tunneling conductance measurements near ν = 2 and V-shaped tunneling conductance measurements near ν = 3 observed in TTG [49]; if we are considering T-IVC fluctuations of the insulator at ν = 2, then it may be reasonable to expect the coupling to these fluctuations could grow weaker as we dope towards ν = 3, in line with the experimentally observed ν dependence.
Note that the regime we call V-shaped here is strictly speaking fully gapped.However, the crucial difference to the BCS state is that the gap is much smaller than the order parameter magnitude as a result of the different Bogoliubov spectrum in Eq. ( 3).This is why, depending not only on the magnitude of the pairing but also on the precise form of the normal state, the resulting tunneling spectra can resemble those observed experimentally [48,49], such as the green curve in Fig. 2(e), making the A 2 state an attractive candidate.The regime of small V where stable superconductivity with true Bogoliubov Fermi surfaces is observed can further exhibit a peak at ω = 0 which is due to a Van Hove singularity crossing the Fermi level, see blue curve in Fig. 2(e); while this peak has not been observed experimentally, its presence crucially depends on details of the normal state band structure and is only found to be energetically favored in a very small regime of V in our model.

E. Electron-phonon coupling
To illustrate that the off-diagonal A 2 state is more generally favored beyond just T-IVC fluctuations, we next discuss electron-phonon coupling, which is frequently considered as a plausible pairing mechanism for twisted moiré systems [33,[35][36][37][38][39].Similar to Ref. 35, we use that the optical A 1 , B 1 , and E 2 phonon modes are known [68] to dominate the electron-phonon coupling in single-layer graphene.As these are optical phonons, we further assume that the impact of the interlayer coupling on the phonons can be neglected and arrive at for the electron-phonon coupling, where v µ encode the layer structure of the modes (see Methods).Symmetry dictates that the vertices Λ g are given by Λ A1 = η x ρ x , Λ B1 = η y ρ x , and Λ E2 = (η z ρ y , −ρ x ) where ρ acts on the microscopic sublattice basis.Integrating out the phonons and projecting to the flat bands, we obtain an effective electron-electron interaction (see Methods) where the coupling constants V g of the three different phonon modes g = A 1 , B 1 , E 2 are estimated to obey V A1 = V B1 ≃ 1.33V E2 for parallel spins in the two valleys, while V A1 = V B1 = 0 for anti-parallel spins.From Eq. ( 9), it is clear that the induced interaction would be always completely attractive if we focused on intra-band pairing, α = α ′ = β = β ′ , which in spinful systems generically favors the trivial pairing channel [51,52].In our case, the combination of two energetically close bands and the trivial pairing being purely band-off-diagonal leads to the competition between different superconductors, even with electron-phonon coupling alone.
To demonstrate this, we study intra-valley pairing within the mean-field approximation and parametrize the relative strength of the different phonon modes with an angle variable θ ph according to The results of the mean-field calculation are summarized in Fig. 3.We see that the A 2 pairing state is favored by the intervalley phonons (θ ph = 0) inspite of its band-off-diagonal nature leading to a suppressed gap [see Fig. 3(a)].This is natural as these phonons mediate an attractive interaction between the two valleys which disfavors the B 1 state, similar to T-IVC fluctuations.In fact, focusing on the leading, momentum independent term, In accordance with symmetry, the A2 (B1) state only has order-parameter components ∝ σy (∝ σ0,x,z).We took ν = 2.5 and V0 = 250 meV • (nm) 2 with a continuum model bandwidth ≃ 2 meV.We point out that if A1 phonons are dominant, as suggested by recent experimental work [61] and past theoretical study in mono-layer graphene [67], we would expect our A2 pairing to dominate assuming the pairing potential is sufficiently large.We also emphasize that although the pairing functions for A2 pairing (b) when θ = 0 and B1 pairing when θ = π/2 (e,f,g) are roughly equal, the excitation spectra shows the B1 state with a band gap on the order of the pairing strength (d) while the A1 state's band gap is nearly zero, see (a).
the chiral limit (see Appendix D3).This maps the problem exactly to that of T-IVC fluctuations, immediately explaining why the order parameter has a fixed sign in Fig. 3(b).As θ ph is increased, the B 1 state is favored (roughly for θ ph > π/4) as can be seen in Fig. 3(c).This is expected since the intravalley E 2 phonon mediates an attractive interaction within each valley such that the energy gain due to the enhanced gap [Fig.3(d)], associated with the band-diagonal matrix elements of the B 1 state, will overcompensate the energetic loss due to the sign change of B 1 's order parameter between the two valleys.This picture is consistent with the dominance and nonsign-changing nature of the band-diagonal components of the B 1 state, see Fig. 3(e-g).Finally, this behavior can also be understood by applying the commutator criterion in Eq. ( 7) in the microscopic sublattice basis, see Appendix D1.This shows that, as opposed to the conventional scenario [51,52], there are two possible leading superconducting states and the superconducting pairing state does not transform trivially under the symmetries of the system even when phonons alone provide the pairing glue.We have checked in our T = 0 numerics that a 60-70 meV•(nm) 2 coupling to A 1 and B 1 phonons (based on Ref. 68) is roughly of the order needed to stabilize the A 2 pairing, assuming the normal state is the flat bands of the un-renormalized continuum model, which in our case has a bandwidth of 2 meV.However, we note that if the interaction-renormalized band splitting is much larger than the continuum model band width, or if the normal state has anti-parallel spins in either valley, addi-TABLE II: Leading superconducting states in the flat-band limit, following from Eq. ( 7), for pairing mediated by fluctuations of the indicated orders, defined by using λ j in Eq. ( 4).Here δ k > 0 and states separated by commas are degenerate.The couplings in the microscopic basis, used in Fig. 4 for the respective orders, are listed under λj .Except for SLP−, the leading superconducting states for λ j and λj are the same (cf.Fig. 4 and Appendix D1).

Fluctuating Order
Leading Superconductor type tional particle-hole fluctuations, such as those of T-IVC order, will also be required for pairing.An interesting scenario arises for anti-parallel spins in the two valley as a magnetic field will cant the spins and, hence, increase the projection of the intervalley phonon matrix elements to the flat bands.At least in TTG, with the suppressed orbital coupling, this could give rise to re-entrant superconductivity at high fields [25].

F. Other particle-hole fluctuations
Finally, we discuss pairing induced by fluctuations of other particle-hole instabilities.In Table II, we list the resulting leading superconductors taking λ j in Eq. ( 4) to be any of the different strong-coupling candidate order parameters [9,[13][14][15]29].In particular, in addition to the T-IVC, we will consider the time-reversal-odd Kramers intervalley coherent state (K-IVC), and time reversal-odd and -even sublattice polarized states (SLP− and SLP+).To analyze how sensitive our conclusions are to the precise form of the coupling of the strong-coupling fluctuating orders to the electrons, we also perform numerics by projecting momentum-independent coupling vertices in the microscopic basis with the correct symmetries (see, e.g., Table II in [13]), listed as λj in Table II, to the flat bands.In the band basis, this leads to momentumdependent coupling vertices, cf.Eq. ( 9).Motivated by recent experiments [7], we will also consider fluctuations of an additional nematic, time-reversal symmetric, layerodd, intervalley coherent state (N-IVC) [62] which is not a candidate ground state in the strong coupling limit; unlike the other strong-coupling ground states, the N-IVC has no momentum independent representation in the flat band basis but does have a momentum-independent matrix order parameter in the sublattice basis which takes the form λ (j,j ′ ) = (η x , η y ) j (ρ 0 , ρ z ) j ′ .The results for fluctuations of the projected strong-coupling orders λj in Table II are shown in Fig. 4, where we use the angle θ fluc.to tune the relative strength between T-IVC and any of the other type of fluctuation-induced interactions by multiplying the T-IVC interaction potential with cos(θ fluc. ) and the other fluctuation potential with sin(θ fluc.).In our microscopic numerics, we have taken a potential form χ(q) = 1 Am V α 2 +|q| 2 /k 2 θ again with α = 0.2 and with V = 4200 meV•(nm) 2 .We chose the value of V such that the transitions between the different pairing states are clearly visible in Fig. 4 when varying θ fluc. .In accordance with the prediction for λj in Table II, SLP+ fluctuations further stabilize the A 2 superconductor, see Fig. 4(a).As such, the band-diagonal B 1 superconducting channel, where SLP+ fluctuations are also attractive, can become the leading channel (favored over A 2 as a result of the finite bandwidth) only very close to θ fluc.= π/2.K-IVC fluctuations, however, are repulsive for A 2 pairing and favor the B 1 state more strongly.
So far, the strong-coupling (λ j ) and sublattice ( λj ) form of the couplings in Table II lead to the same conclusions.This is different for SLP− fluctuations [Fig.4(d)], where the projection-induced momentum-dependence in the band basis can stabilize the E 1 superconductor.This can be understood by applying Eq. ( 7) in the sublattice basis (see Appendix D1).We also find the E 1 state when fluctuations of the N-IVC state of Ref. 62 dominate.Examples of the E 1 nematic and B 2 order parameters which emerge for SLP− fluctuations or N-IVC fluctuations are shown in Appendix F. We point out the nematic E 1 pairing is also an interesting candidate given that despite having nonzero pairing in the σ 0 , σ x , σ z channels, it will be nodal as long as the σ x components do not gap out the nodes in the band-diagonal parts.

III. DISCUSSION
Taken together, we see that the proposed band-offdiagonal A 2 superconductor is an especially attractive candidate for TBG and TTG: first, it can lead to both V-shaped or U-shaped DOS, depending on lifetime parameters, the normal state, and the coupling strength V , see Fig. 2(e).As these parameters might vary from sample to sample and within a sample (e.g., V is expected to decrease upon doping further away from the insulator), this can naturally explain the tunneling data of [48,49].We emphasize however that at least at the level of our mean-field numerics, we only expect a Vshape in the regime where the superconducting pairing is of the order of the bandwidth; this is the regime, where although the pairing is finite and can be quite large, the gap in the superconducting spectrum is either just closing or very small relative to the pairing.Increasing the pairing further will lead to an evolution from V to U shaped while decreasing the pairing will eventually lead to a nodal Fermi surface and presumably a peak at zero energy in the DOS.Second, despite its interband nature, A 2 is the unique pairing state that is favored by fluctuations of two out of the four strong-coupling-candidates we consider for the correlated insulator, see Fig. 4(a-c).What is more, this includes the T-IVC state, signatures of which are observed in recent experiments [7].Finally, it is also favored by the likely dominant [61,67] optical intervalley phonon modes.We emphasize that, both in the case of fluctuating correlated insulators and phonons, the minimum attractive coupling needed to stabilize a purely band off diagonal state depends on the energy splitting between the two flat bands in the normal state; if the bands of our normal state are closer to degenerate, irrespective of the total bandwidth, the needed coupling to stabilize the A 2 pairing in mean-field will decrease.
The other band-off-diagonal superconductor we identify transforms under the IR E 2 , i.e., can be thought of as a p-wave state.Its spectral properties also agree well with experiment as the chiral configurations, E 2 (1, i), which is favored within mean-field theory over a nematic E 2 state, can also have nodal regions, depending on filling.As can be seen in Fig. 1(g), this can lead to a transition from gapped to nodal when increasing the electron filling starting at ν ≃ 2. However, as opposed to the A 2 state, E 2 does not naturally appear as leading instability when considering optical phonons or fluctuations of any of the strong-coupling order parameters of the correlated insulator.While this makes it energetically less natural than A 2 , we cannot exclude it since its phenomenology agrees well with experiment and since the precise form of the coupling of the dominant low-energy collective excitations are not known-significant momentum dependencies beyond λ j and λj in Table II could stabilize E 2 pairing as well.We also find in our numerics a nematic E 1 state which may be preferred over its chiral version in the presence of sufficient strain or due to fluctuation corrections [27,[53][54][55].We find the E 1 state is the leading instability of nematic IVC fluctuations and SLP− fluctuations, and is a subleading instability of T-IVC fluctuations.The E 1 state is interesting in its own right, as it can also be nodal.
As superconductivity might further coexist with T-IVC order [7], we have checked (see Appendix E) that this does not alter our main observation: the preserved C 2z symmetry still allows for entirely band-off-diagonal states, with transitions from nodal to full gapped, which are stabilized (among other fluctuations) by intervalley phonons.
For the future, it will be interesting to go beyond mean-field and analyze the competition of our band-offdiagonal states with odd-frequency pairing, which we study in a follow-up work [69].It also seems promising to study Andreev reflection [48,49] for our interband pairing scenario.On a more general level, our work shows that the observation of nodal pairing in twisted graphene systems does not immediately exclude a chiral superconducting state nor an entirely electron-phonon-based pairing mechanism.It illustrates that a microscopic understanding of the superconducting states in graphene moiré systems requires taking into account their intrinsically multi-band nature.
Note added.Just before posting our work, Ref. 70 appeared online, which discusses pairing induced by A 1 phonons in spinful TBG bands.
As the Frobenius inner product ⟨A, B⟩ F = tr[A † B] reaches its maximum (minimum) at fixed ⟨A, A⟩ and ⟨B, B⟩, if A = cB with c > 0 For the ansatz ∆k = δ k Dη x (and assuming for now that δ k has a fixed sign for all k), this is obeyed if We state Eq. ( 12) as the (anti)commutator condition (7) in the main text [equivalent if (λ j ) 2 = 1], not only because it highlights the simple algebraic and basis independent nature of the condition but also since it emphasizes the similarities to the generalized Anderson theorem of [71,72].
If we can find a solution to Eq. ( 12), we know that the maximum (or at least one of the possibly degenerate maxima) of F[ ∆k ] is of the form of ∆k = δ k Dη x where δ k is obtained as the maximum of the reduced functional or equivalently as the largest eigenvector of χ k−k ′ viewed as a matrix in k and k ′ .As χ k−k ′ > 0 (due to stability), the Perron-Frobenium theorem then immediately implies δ k > 0, in line with out assumption above and as stated in the main text.
Electron-phonon coupling.To present more details on the electron-phonon coupling, the associated displacement operators in Eq. ( 8) can be expressed in terms of canonical bosons, b g,α,µ,q , (u g,µ (r)) j = q b g,j,µ,q e iq•r + H.c.
where j refers to the two components for the E 2 phonon (is idle for A 1 , B 1 ), M is the carbon mass, and ω g (q) is the phonon dispersion, characterizing the phononic part of the Hamiltonian, H P = q ω g (q)b † g,j,µ,q b g,j,µ,q .As for (v µ ) ℓ in Eq. ( 8), ℓ = 1, 2 refers to the physical graphene layer in the case of TBG.One can, in principle, choose any orthonormal basis; we will find it convenient to use the layer-exchange even and odd states, For TTG, the situation is more involved (see Appendix D2), but our arguments about which phonons are attractive in which pairing channels will hold for both systems.
We project H EP in Eq. ( 8) onto the two flat bands (α = ±) in each valley η of the spin polarized continuummodel, leading to a coupling term similar to Eq. ( 4) with momentum-dependent coupling matrices, λ j → λ g,j,µ k,k ′ .Investigating the matrix elements λ g,j,µ k,k ′ , we notice that they almost vanish for the layer-odd intervalley (A 1 , B 1 ) phonons, which can be understood as a consequence of chiral and particle-hole symmetry (see Appendix D3).The situation is the reverse for the intravalley (E 2 ) phonons, where the layer-even matrix elements are numerically small and the layer-odd matrix elements dominate.We therefore focus on layer-even (odd) intervalley (intravalley) phonon couplings.
Neglecting the momentum dependence in the phonon frequencies and retardation effects, the resulting electronelectron interaction in the inter-valley Cooper channel obtained by integrating out the phonons is given by Eq. ( 9).Here, V g = g 2 g /(2N ω 2 g ) > 0 and V A1 = V B1 ≃ 1.33V E2 results from g A1 = g B1 ≃ g E2 and the phonon frequencies estimated in Ref. 68.Importantly, this only holds for parallel spins in the two valleys.For antiparallel spins, the projection of the coupling matrices to the flat bands vanishes for the intervalley phonon modes A 1 and B 1 such that V A1 = V B1 = 0.
The situation is more non-trivial in case (ii).Let us assume, for notational simplicity, that the spin polarization of the active flat bands in valley η = + is s =↑ and in valley η = − is ↓.Accordingly, we define which are symmetries of the system and obey the same algebraic relations as the symmetries in the main text, In fact, their representation on the fermions defined in Eq. ( A2) is exactly the same as that of Θ and C 2z in the main text, As such, for case (ii), the time-reversal symmetry Θ and two-foldrotational symmetry C 2z in the main text can be identified with Θ and C 2z in Eq. (A3).To illustrate this further and also explicitly discuss the spin structure of the order parameter, we transform the superconducting order parameter back to the d-fermions via Eq.(A2), which shows that we obtain an admixture of singlet and (unitary) triplet pairing.To demonstrate the action of Θ and C 2z more explicitly and provide a consistency check, let us focus on ∆ k = 2∆σ y , where Eq. (A5) becomes From Eq. (A3), we find the representations C 2z : d k → η x s x d −k and Θ : d k → η x s x d −k ; applying this in Eq. (A6), we find that exactly as in the main text.
For case (i), Eq. (A5) instead becomes i.e., a non-unitary triplet state-as expected [27] since this is the "Hund's partner" of the singlet-triplet admixed state in Eq. (A5), obtained by an independent spin-rotation in the two valleys [SU(2) − × SU(2) + ].For ∆ k = 2∆σ y this yields Again in accordance with the spinless formulation of the main text, we get C 2z : ∆ → ∆ and Θ : ∆ → −∆ * .We finally note that the normal-state polarization also determines the spin-structure of the fluctuating orders in Table II and Table III: switching between the two scenarios (i) and (ii) requires replacing an order parameter for the correlated insulator by its "Hund's partner" (see, e.g., Table II in [13] for a complete list).As the system is believed to be close to the SU(2) − × SU(2) + symmetric limit (the intervalley Hund's coupling was estimated to be smaller than 0.1 meV in [24]), the strength of fluctuations of Hund's partners is expected to be roughly the same.As such, both scenarios (i) and (ii) are consistent with a mechanism based on fluctuations of an order parameter of a correlated insulator.As mentioned in the main text, this is different for phonons, where only scenario (i) allows for intervalley phonons providing the pairing glue.In this appendix, we will describe how we compute solutions to the linearized gap equation at T c .As in App.B, we will assume a spin polarized normal state and only consider superconducting instabilities within a single spin flavor.We recall that for the case of fluctuation-mediated superconductivity, we couple electrons to bosonic modes (j = 1, 2, . . . ) as, e.g., in Eq. ( 4), with λ j α,η;α ′ ,η ′ capturing the symmetries broken by the corresponding order parameter.In order to compactly write down the linearized gap equation, it is convenient to express λ j as Here we also include the momentum dependence of the matrix elements, which arises when we study phonons and order parameter fluctuations projected from the sublattice basis to the band basis.In Eq. (C1), A k,k ′ are the valley off diagonal pieces of the form factor λ j α,η;α ′ ,η ′ and B k,k ′ is the valley diagonal pieces.With this notation in hand, the linearized gap equation we solve is where the Greens function Here ∆ αβ k denote the pairing in band space where α, β = 0, 1 label the upper and lower flat band.Finding a solution to the above equation then amounts to computing the right-hand side of Eq. (C2), diagonalizing it in the space of momenta, Nambu index, and band index, and looking at the eigenvectors which attain eigenvalue 1 for some value of T .To enforce Fermi-Dirac statistics, we solve the above equation on half of the moiré Brillouin zone.We also exclude the edge points in our linearized gap equation computations for phonons and projected order fluctuations.We expect including these points would reduce the needed coupling to obtain a finite T c (or reduce T c for fixed coupling) but not change the leading instabilities.

T-IVC & SP order
Given the current insights from experiment, the most natural alternative scenario is that the normal state exhibits both T-IVC [7] and spin polarization [24,28] simultaneously.The projection to the remaining two active flavor degrees of freedom is given by Increasing ν beyond ν = 2 will lead to a metallic state with two non-degenerate bands α = ± coming from the original flat-band manifold.Let us denote the associated creation operators by c † k,α , which have one index less than the associated operators discussed in the main text since valley is not a good quantum number anymore.The superconducting order parameter is a 2 × 2 matrix, coupling to the electrons as ., and thus has to obey ∆ k = −∆ T −k .As the projector in Eq. (E1) commutes with C 2z (in fact, also with C 2x and C 3z ), all pairing states must still be either even or odd under C 2z (transform under one of the IRs of D 6 or C 6 ).Since Eq. (E1) projects onto the subspace where η which is the analogue of Eq. ( 2) of the main text.As before, all C 2z -even states must be entirely band-off-diagonal, ∆ k = δ k σ y .However, since the number of active degrees of freedom is reduced, there are more restrictions: all C 2z -odd superconductors must have zeros in the Brillouin zone due to For completeness and to conveniently address energetics, we extend the discussion to the microscopic sublattice basis.Let ∆k be the corresponding superconducting order parameter-an 8 × 8 matrix in sublattice, valley, and spin space.Then pairings are constrained to obey The order parameters which are compatible with Eq. (E3) will all be spin triplets.The C 2z -even states, i.e., order parameters transforming under A 2 , E 2 , or A 1 , will have the form (suppressing k-dependencies) ∆k ∼ P ν=2 s x η z ρ z ; in line with our symmetry arguments above, one can check that they will go as σ y in band space and thus be purely band off diagonal in the subspace defined by P ν=2 .The pairings which are odd under C 2z include the B 1 and B 2 pairings with ∆k ∼ P ν=2 s x ρ 0 , and E 1 pairings with ∆k ∼ P ν=2 s x (ρ x , ρ y η z ) previously discussed in our main text; however, as pointed out above and unlike in the main text, the C 2z -odd pairings in both the band basis and sublattice basis are no longer allowed to have a component without a sign change since only the momentum odd components of the B 1 , B 2 , and E 2 pairings survive projection P ν=2 .Since only the band-off-diagonal A 2 state can have a non-sign-changing order parameter, a superconducting state satisfying the criterion around Eq. ( 7) of the main text can only be this state (or none).We have studied which of the pairing mechanisms survive the projection and whether they favor or disfavor A 2 pairing, see Table IV.We find that A 1 phonons, T-IVC fluctuations, and spin fluctuations all provide an attractive pairing potential, and if any of these have large enough couplings to overcome the normal state band splitting, the A 2 triplet pairing is the leading instability, as in the main text.Furthermore, due to the fact that the remaining bands after reconstruction, as described by the projector P ν=2 , are not degenerate (there is no remaining spin symmetry to guarantee degeneracy), a Bogoliubov Fermi surface or a fully gapped state and, thus, a transition from nodal to gapped as a function of filling are possible depending on parameters (similar to our discussion in the main text).

TABLE IV:
We list the possible pairing glues which are compatible with a T-IVC+SP normal state (i.e., the interactions survive projection to the space of the upper T-IVC bands of a single spin flavor).We denote interactions which will generate an attractive interaction for the A2 pairing with a ✓ and interactions which will generate a repulsive interaction with a ✗.

T-IVC normal state
We will now consider a simpler normal state which leaves twice the number of degrees of freedom as the previous normal state we considered.In particular, we can consider a strong coupling T-IVC normal state with projector of the form: In contrast to the case for a normal state with coexisting T-IVC and spin-polarized order, there are now more possible pairing options and singlet pairing is once again possible.We can classify the possibilities as pairings which are triplet, singlet, and by IRs of the point group.We find the possible pairings include triplet A 1 and A 2 pairings and singlet B 1 and B 2 pairings with: triplet B 1 and B 2 pairings and singlet versions of our A 1 and A 2 states with: ∆ k ∼ P ν=2 s x η 0 ρ 0 ∆ k ∼ P ν=2 s 0 η 0 ρ 0 (E6) and triplet E 1 pairing and singlet E 2 pairing with: ∆ k ∼ P ν=2 s x (η 0 ρ x , η z ρ y ) ∆ k ∼ P ν=2 s 0 (η 0 ρ x , η z ρ y ) (E7) Of the above, the only options which are not enforced to have a sign change are our purely inter-band A 2 triplet pairing, the A 1 singlet pairing, and the E 2 singlet pairing.Since these pairings do not have a sign change, they are the only possible candidates for the criterion around Eq. ( 7) of the main text and we have enumerated the possible pairing glues for these s-wave states in Table V.

TABLE V:
We list the possible pairing glues which are compatible with a T-IVC normal state (ie the interactions survive projection to the space of the upper spin-degenerate T-IVC bands).We denote interactions which will generate an attractive interaction with a ✓and interactions which will generate a repulsive interaction with a ✗.
We find in this case that all of the pairing glues which are attractive for our A 2 triplet pairing are also attractive for one of the singlet pairings, except for spin fluctuations.Therefore, we can say that if the pairing is triplet for a spin-degenerate T-IVC normal state, the leading instability is likely to be our A 2 pairing provided the pairing glue interaction is sufficiently strong and spin fluctuations may play an important role in energetically favoring this state.In this case, we expect that the phenomenology of Bogoliubov Fermi surfaces and a nodal to gapped transition as a function of interaction strength will again apply.

Appendix F: More Superconducting Instabilities
In this appendix, we will discuss the superconducting instabilities we find beyond the A 2 and B 1 states shown in Figs. 2 and 3 of the main text and focus on the other leading instabilities we find in the presence of fluctuations of different particle hole orders.For SLP− fluctuations, we find the B 2 state can be favored over the B 1 when the strength of T-IVC fluctuations are on the same order as SLP− fluctuations, as shown in Fig. 4. We show the B 2 state for parameter value θ f luc.≃ π 4 in Fig. 5.For N-IVC fluctuations as well as for SLP− fluctuations, we find the E 2 is the leading instability, as shown in Fig. 6.We show the two components of the E 2 state for parameter value θ f luc.≃ π 2 in Fig. 4.

FIG. 1 :
FIG. 1: Spectral properties of interband pairing.While for band-diagonal pairing a small superconducting order parameter can immediately open up a gap as time-reversal symmetry guarantees that the associated avoided crossings [gray regions in (a)] in the Bogoliubov spectrum are at the Fermi level, this is not the case for band-off-diagonal pairing (b).Here, a sufficiently strong order-parameter value is required to establish a full gap, see (c).Its k dependence according to Eq. (3) is shown in (d), where the red line indicates nodal points.If the band structure has Dirac points, there will be a point on the horizontal axis (blue cross).Consequently, if there is another momentum point located above the red line (blue circle), continuity of the Hamiltonian implies a nodal point on any path connecting the two momenta.(e) Gap of the isotropic A2 state and δ k , ϵ k (zeros indicated in red) for the normal-state toy model defined in the text.BEC/BI refers to the Bose-Einstein condensate/band insulator limit.(f) Complex phase φ k = arg(X k + iY k ) for leading basis function with small subleading corrections.(g) Shows the gap of the chiral p-wave E2(1, i) state with ∆0 = 1.5t and the value of D k j := |δ k j | − |ϵ k j − µ| for kj at the three symmetry-in-equivalent vortices in (f) as a function of µ.We took t ′ = −2.2t,t > 0, in (b,d).

FIG. 2 :
FIG. 2:Pairing mediated by T-IVC fluctuations.We show (a) the band structure of the normal state with spin polarization (K, K', and Γ label the high-symmetry points of the moiré scale Brillouin zone) and (b) the critical temperature Tc (in units of the maximum band splitting W0 ≃ 9.4 meV) as a function of coupling strength V measured in units of the critical coupling Vc,1 = 105 meV • nm 2 obtained from the linearized gap equation.The band structure (with color indicating the band-projected value of the anomalous correlator) of the A2 state and its order parameter are shown in (c) and (d).The DOS of the T = 0 superconductor for several different values of coupling strength V is plotted in (e).The DOS was computed as k δ (E k − ω), replacing the δ function with Lorentzians with half width at half max 0.3 meV (much smaller than the typical superconducting order parameter).The critical coupling Vc,2 where the nodal lines disappear is Vc,2 ≃ 1.4Vc,1.

FIG. 3 :
FIG.3: Pairing from electron-phonon coupling.We show (a) the band structure and (b) the self consistent order parameter of the A2 pairing for θ ph = 0 and T = 0.The eigenvalues corresponding to the A2 and B1 pairings in the linearized gap equation at T = 5 K, which is close to their Tc, are shown in (c) as a function of θ ph .We show an example of the band structure (d) of the B1 pairing and its order parameter (e,f,g).In accordance with symmetry, the A2 (B1) state only has order-parameter components ∝ σy (∝ σ0,x,z).We took ν = 2.5 and V0 = 250 meV • (nm) 2 with a continuum model bandwidth ≃ 2 meV.We point out that if A1 phonons are dominant, as suggested by recent experimental work[61] and past theoretical study in mono-layer graphene[67], we would expect our A2 pairing to dominate assuming the pairing potential is sufficiently large.We also emphasize that although the pairing functions for A2 pairing (b) when θ = 0 and B1 pairing when θ = π/2 (e,f,g) are roughly equal, the excitation spectra shows the B1 state with a band gap on the order of the pairing strength (d) while the A1 state's band gap is nearly zero, see (a).

TABLE I :
Summary of pairing states in spin-polarized flat bands.Here χ k ( χk ) is a real-valued (real and symmetric 2 × 2 matrix-valued) MBZ-periodic function invariant under C3z.Furthermore, X k and Y k ( Xk and Ŷk ) transform as x and y under D3, generated by C3z and C2x, while also being real (and symmetric).The third column indicates the type of nodes-line (ln), point (pt), or none (n)-on a generic Fermi surface for sufficiently small/large order-parameter magnitudes; options separated by "or" indicates that this depends on the normal-state band splitting, see main text.

TABLE III :
Generalization of TableIIof the main text, where we also indicate the dominant superconducting orders ( ∆) in the microscopic basis, obtained by applying Eq. (7) in the sublattice basis.The phonon modes refer to the sublattice-basis form λj of the coupling, cf.Eq.(8), and "g-nematic" stands for the (intravalley) graphene nematic state of Ref. 62, which has the same coupling as the E2 phonon.