Nematic superconductivity in magic-angle twisted bilayer graphene from atomistic modeling

Twisted bilayer graphene (TBG) develops large moir\'e patterns at small twist angles with flat energy bands hosting domes of superconductivity. The large system size and intricate band structure have however hampered investigations into the superconducting state. Here, using full-scale atomistic modelling with local electronic interactions, we find at and above experimentally relevant temperatures a highly inhomogeneous superconducting state with nematic ordering on both atomic and moir\'e length scales. The nematic state has a locally anisotropic real-valued d-wave pairing, with a nematic vector winding throughout the moir\'e pattern, and is three-fold degenerate. Although d-wave symmetric, the superconducting state has a full energy gap, which we tie to a {\pi}-phase interlayer coupling. The superconducting nematicity is further directly detectable in the local density of states. Our results show that atomistic modeling is essential and also that very similar local interactions produce very different superconducting states in TBG and the high-temperature cuprate superconductors.

While the superconducting transition temperature T c ≈ 3 K is low in TBG [2,5], the ratio to the Fermi temperature is large T c /T F ≈ 0.1, exceeding the weak coupling regime [12]. In fact, with superconductivity appearing in close proximity to correlated insulators and a pseudogap state with reduced DOS above the superconducting domes [7][8][9] with accompanied strange metal behavior [14,15], the phase diagram of TBG shares striking similarities to the high-temperature cuprate and pnictide superconductors [12,16,17]. This points to the possibility of strong local electronic interactions being responsible for superconductivity, although electron-phonon pairing is also a plausible candidate [14,[18][19][20].
The large scale of the moiré pattern (∼ 10 4 carbon atoms at the magic angle) and the non-trivial band topology have severely hampered studies of superconductivity in TBG. Continuum models reproduce the normal-state band structure [21,22], but are harder to reconcile with the possibility of strong local electronic interactions. Effective lattice models, using e.g. the moiré pattern or otherwise rescaled lattice schemes, have been used in Hubbard-like tight-binding studies, but must strike a difficult balance between accuracy and orbital proliferation [23]. Many of these studies have proposed a time-reversal symmetry (TRS) breaking chiral d-wave state [23][24][25][26][27][28][29][30], as also found in heavily doped graphene [31,32].
Finding chiral d-wave superconductivity is at first sight not surprising, considering that the d x 2 −y 2 -wave state of the cuprates naturally transforms into the two d x 2 −y 2and d xy -wave states with a symmetry enforced degeneracy at T c , on lattices with three-or six-fold symmetry, such as graphene's honeycomb lattice and the TBG moiré lattice [32][33][34][35]. On the honeycomb lattice, the chiral d x 2 −y 2 + id xy -wave combination then becomes fully gapped, while any real-valued nematic d-wave state always has a nodal spectrum, thus making the chiral state energetically favored and the ground state at zero temperature [31][32][33]. A similar favoring could naively also be expected to hold in TBG. However, in apparent contrast to these basic expectations, different nematic superconducting states have very recently been proposed for TBG. Some of these require higher order coupling to additional coexisting orders or pairing channels to energetically favor the nematic state over the chiral state [35][36][37][38]. More intriguing are direct findings, i.e. without other orders or pairing channels, of a nematic phase in a phase diagram region using either rescaled or atomistic lattice models, although with seemingly varying properties [23,24,39]. Experimentally, recent magnetotransport measurements have also identified intrinsic nematicity in the superconducting state of TBG [40]. Taken together, these results all highlight the question of if and how nematic superconductivity appears in TBG, as well as its measurable consequences.
In this work, to accurately capture superconductivity in TBG, we use full-scale atomistic modeling, including all carbon atoms, a dense k-point sampling, and electron interactions in agreement with both current TBG experiments and known interactions in graphene(-like) systems, as well as consistent with the cuprates. By solving both the mean-field self-consistent and linear gap equations using realistic transition temperatures, we find a highly inhomogeneous and real-valued nematic d-wave superconducting state, with both atomic-and moiré-scale nematicity, dominating the phase diagram and at the experimentally observed critical temperatures. We further show directly measurable signatures of the nematicity in the local density of states. Unexpectedly, the nematic d-wave state also has a full energy gap, which we tie to a strong π-phase Josephson locking between the graphene layers.

II. RESULTS
A. Normal state properties TBG consists of two graphene layers rotated by a twist angle θ, which produces a periodic moiré interference pattern of length L M = a/(2 sin (θ/2)) [21,41], far exceeding the graphene lattice constant a = 2.46Å, see Fig. 1(a). We model commensurate angle TBG with a full-scale tight-binding model of all carbon atoms in the large periodic moiré unit cell. The out-of-plane p z carbon orbitals hybridize in-plane through standard nearest neighbor hopping t = 2.7 eV, while the interlayer hybridization is captured by an exponentially decaying hopping together with a Koster-Slater angular dependence [22,42], see Methods.
In reciprocal space, the Dirac cones from the graphene layers are displaced by the twist and sit at the corners of the smaller moiré Brillouin zone (MBZ), see Fig. 1(b). As the twist angle decreases, the Fermi velocity of the Dirac cones is reduced. The layer and valley degrees of freedom then form four spin-degenerate narrow bandwidth moiré bands that separate from the remaining band structure [21], see Fig. 1(c). As a consequence, TBG has a large DOS around zero energy, peaking in two van-Hove singularities (VHSs) that correspond to saddle points in the band structure [3,21]. At the magic angle, the moiré bands become essentially completely flat and the two VHSs merge, see Figs. 1(d,e).
The states of the moiré bands are primarily localized to the AA region of the moiré cell, where the carbon atoms of the two layers are aligned, as shown by the inhomogeneous and three-fold symmetric charge density N m (x) in Fig. 1(f). Resolving the charge density with respect to energy, Fig. 1(g) further shows the local density of states (LDOS) along a periodic path passing the AB, AA, and BA regions (black dotted line in Fig. 1(f)).

B. Modeling superconductivity
Many properties of the superconducting state in magicangle TBG are still unknown. We therefore employ a general model for the superconducting pairing, guided by only a few constraints: the weak interlayer van der Waals coupling motivates intralayer pairing, while the observed suppression of superconductivity in in-plane magnetic fields and no evidence for spin-polarized Cooper pairs, restrict us to consider spin-singlet pairing [2,43]. The strong on-site repulsion in graphene, graphite [44], and TBG [6] also points to spin-singlet pairing, since in the strong coupling limit of the Hubbard model, the resulting t-J model gives spin-singlet nearest-neighbor bond interactions [31,45]. In fact, preference for bond spinsinglet configurations, over double and single occupancies, was already central in early treatments of pπ-bonded planar organic molecules (of which graphene is the infinite extension) [46], and subsequently also used for cuprate superconductors [47]. Spin-singlet bond pairing has also recently been derived from spin-fluctuations in TBG [24] and also used in smaller rescaled lattice models [23,39]. Consequently, we model superconducting pairing in TBG by spin-singlet order parameters on every in-plane nearest-neighbor carbon bond, using a uniform coupling strength J, see Methods. Together these bond order parameters forms an order parameter field∆ throughout the moiré cell.
We solve for the superconducting order parameter field with each bond order parameter treated fully independently using both the full non-linear, self-consistent and the linearized gap equations, see Methods. The linear gap equation is valid at T c and has the same symmetry as the normal state, leading to its solutions always belonging to one of the irreducible representations of the symmetry point group D 3 of TBG: the one-dimensional A 1 and A 2 , or the two-dimensional E representation. Below T c , non-linear contributions enter, possibly further breaking the symmetry, which we capture by iteratively solving the full non-linear and self-consistent gap equation at zero temperature. This combined approach enables us to completely characterize the superconducting state.

C.
Moiré-scale nematicity Starting with the results from the linear gap equation, we show in Fig. 2(a) the four highest critical temperatures T c as a function of coupling strength J at the magic angle and with the Fermi level aligned with the upper VHS, corresponding to ν ≈ 1. Across the wide range of all investigated coupling strengths, the solution with highest T c (green highlight) belongs to the E irreducible representation and is therefore two-fold degenerate and spanned by two real-valued order parameter fields:∆ x and∆ y . In Fig. 2 (left) and lower (right) graphene layer in the moiré cell, while the amplitude of∆ y in the upper layer is plotted in Fig. 2(c). Both∆ x and∆ y clearly break the C 3 rotation symmetry of the normal state and are instead enhanced in the AA regions and along the C 2 nematic axes aligned with the x-and y-axis, respectively.
In contrast to pristine graphene where superconductivity only develops for J/t above the quantum critical point (QCP) J c = 1.76t [31,48,49], TBG shows substantially enhanced ordering in Fig. 2(a), with finite a T c even at very weak coupling. For weak coupling, T c is approximately proportional to J/t, as expected for an isolated flat band [50][51][52][53]. For slightly stronger coupling, the growth in T c however accelerates, especially when approaching the graphene QCP. This both indicates contributions from dispersive bands [39] and suggests that the low energy moiré flat bands have a catalytic effect that at least partially triggers the QCP cascade boosting T c at couplings also below the QCP.
Alongside T c extracted from the linear gap equation, we also plot in Fig. 2(a) (blue line, right y-axis) the spatially averaged amplitude of the self-consistent order parameter field, obtained from the full non-linear gap equation at zero temperature (in Kelvin). As seen, T c and the averaged order amplitude are directly comparable, despite the strong amplitude inhomogeneity over the moiré unit cell. After fixing J to match the experimentally observed T c ≈ 4 K [2,5], we tune the chemical potential across all moiré band fillings. The resulting zero temperature self-consistent order amplitude is shown in the inset of Fig. 2(a). The superconducting amplitude predictably drops sharply at the two moiré band edges. In between, the amplitude is rather stable but has two local maxima . Plot marker color labels the symmetry of each solution, gray lines guide the eye. The leading two-fold degenerate solution is highlighted in green. Highlighted in blue are the spatial averaged amplitude order parameter, |∆(xi)| k −1 B , (right y-axis) obtained from the non-linear self-consistency (SC) equation at zero temperature (T = 10 −7 t). Inset: Self-consistent order parameter as a function of chemical potential µ at weak coupling J = 0.42t. Vertical lines mark upper VHS (dashed red), lower VHS (dashed black), and half-filling (solid black). (b) Plot of the order parameter field∆x throughout the moiré cell in the top (left) and bottom graphene layer (right). Underlying color intensity plot shows the amplitude, clearly depicting a global, moiré-scale, nematicity. The streamlines of the local vectors field χ(x) illustrate the intricate atomic-scale nematicity in the local d-wave bond order parameters, as well as the superconducting π-phase-shift between the two layers. (c) Same as (b) but for∆y in the top graphene layer. (d) Color intensity plot of the order parameter field amplitude for the chiral∆x + i∆y solution in the top graphene layer. (e) Relative free energy landscape at T = 0 for all independent linear combinations∆(Θ, ϕ) of the∆x,y solutions. (f ) Cut along the real-valued linear combinations of (e), which due to gauge invariance are located at the lines ϕ = 0, 2π (blue) and ϕ = π (green) in (e). The line is a least square fit, showing a three-fold degenerate nematic ground state. (g) Energy difference between the most and least favorable nematic combinations in (f) as a function of tuning the averaged order parameter amplitude, with dot marking value for J = 0.43t. Inset: Energy difference between most favorable nematic state and the time reversal symmetry (TRS) breaking chiral∆x + i∆y combination. (h), (i), (j) Low-energy density of states (DOS) for most favorable nematic, least favorable nematic, and chiral combination, respectively, extracted at the points of the same color in (e). Plots show that the ground state is fully gapped, while the other combinations show a (partially) filled energy gap reflecting lower superconducting condensation energies (given in figures). For all but panels (a) and (g): J = 0.43t and chemical potential aligning the Fermi level with the upper VHS (band filling ν ≈ 1), but conclusions are unchanged for other J.
near the upper and lower VHS (dashed lines) and a local minima at half filling (solid black line), largely replicating the features of the normal-state DOS in Fig. 1(c). Further analyzing the zero temperature self-consistent solutions at different experimentally relevant couplings and band fillings, we always find that they are exclusively real-valued, even when self-consistency is achieved iteratively from an initial complex number field configuration.
In fact, we find that the self-consistent solutions are always real-valued linear combinations of the two leading and degenerate solutions∆ x and∆ y of the linear gap equation, see Methods. This makes TBG a nematic superconductor [54], with nematic symmetry breaking on the moiré length scale. As both the normal-state and electronic interactions are fully isotropic, this nematicity is entirely due to spontaneous symmetry breaking in the superconducting state.
D. Degenerate ground state.
At T c all the states∆(Θ, ϕ) are degenerate solutions, but due to non-linear contributions of a finite order field only specific combinations are expected to be energetically favorable at zero temperature. Demonstrating this symmetry lifting, we show in Fig. 2(e) the relative free energy of all possible∆(Θ, ϕ) at zero temperature, confirming the self-consistent results, the free energy minima are achieved for the real-valued linear combinations. The free energy maxima (i.e. smallest condensation energy) are instead attained around the rotationally symmetric (see Fig. 2(e)) TRS breaking complex combinationŝ ∆(π/4, ±π/2) =∆ x ± i∆ y . The free energy maxima (i.e. smallest condensation energy) are instead attained around the rotationally symmetric (see Fig. 2(e)) TRS breaking complex combinations∆ ± =∆(π/4, ±π/2) = ∆ x ± i∆ y .
On closer inspection, the energy splitting among the real-valued combinations is six-fold symmetric as a function of the nematic angle Θ, see Fig. 2(f). As a result, ∆ x is one of the three gauge-inequivalent nematic ground states, all with the nematic C 2 -axis directed towards one of the next-nearest-neighbor AA regions. In contrast,∆ y is one of the three least stable gauge-inequivalent nematic states that instead has the C 2 -axis directed towards one of the nearest-neighbor AA regions, as seen in Fig. 2 The amplitude of the energy splitting among the realvalued nematic states depends non-monotonically on the field norm (or equivalently on J or T c ), as seen in Fig. 2(g), with notably large values, up to 10% of k B T c , around the experimentally observed T c . In contrast, the energy of the TRS breaking chiral state∆ ± relative to the nematic ground state is large around the experimentally observed critical temperatures and also increases monotonically with the field norm (or T c or J/t). We thus find that the nematic state is heavily favored in this regime, and even more so at stronger coupling. Based on these results it is also not surprising that nematic states have previously also been reported in the strong coupling regime [23,24]. While any degeneracy lifting of the E manifold must necessarily vanish as the field norm goes to zero, we do find that the chiral state,∆ ± , just barely becomes the ground state (beyond the resolution used in Fig. 2(g)) for very small field amplitudes, before transitioning to the nematic state for more realistic temperatures.
The energy split among the different superconducting states is further reflected in the different gap structures among the∆(Θ, ϕ) states, as revealed by the low energy DOS in Figs. 2(h)-(j) at an experimentally relevant temperature. The ground state∆ x =∆(0, 0) is fully gapped with a resulting large condensation energy achieved by pushing occupied VHS states down in energy in to sharp coherence peaks. In contrast, the least favored nematic states, including∆ y , are not fully gapped, in turn limiting the condensation energy. For the least favored of all the states in the E manifold, the TRS breaking chiral state∆ ± , we find a substantially reduced gap with coherence peaks that are closer together, both features substantially reducing the resulting condensation energy.
Based on the results above, we have demonstrated that TBG is a nematic superconductor at experimentally relevant temperatures with a strong dependence on the orientation of the C 2 nematic axis affecting both the condensation energy and gap structure. In particular, the realvalued nematic ground state is energetically favored due to a fully gapped spectrum, which is in sharp contrast to superconductivity in doped graphene where the nematic state is always nodal while the chiral d-wave state is energetically favored and has a full energy gap [31,32,55]. In fact, our findings of a d-wave nematic state with a full energy gap is overall unexpected as d-wave symmetry traditionally is associated with a nodal energy spectrum. However, it has recently been shown that the situation can be more complicated in multiorbital or multiband systems, where completely nodeless and gapped d-wave superconducting states has recently been found [56][57][58][59][60]. This points towards the large scale moiré pattern and the multiple nearly degenerate moiré bands as having an important role in determining the properties of the superconducting state. In particular, we find that the gap structure sensitively depends on the fine distinctions that arise with different C 2 nematic axis orientations, and we also note that the gap structure of previously reported nematic d-wave states have been found to be nodal in rescaled lattice and continuum model while unreported in recent atomistic model [23,24,27]. Moreover, we note that our results show that the real-valued nematic ground state is favored over the complex chiral combination directly in the E pairing channel of a simple pairing model.
Our nematic ground state is therefore produced without any higher order coupling to additional coexisting orders or nearly degenerate pairing channels that have additionally been shown to otherwise favor nematicity [35][36][37].

E. Atomic-scale d-wave nematicity
Having established a nematic superconducting state in TBG on the moiré length scale, we next characterize the superconducting state on the atomic scale, where the symmetry is well described by the point group D 6h of the graphene honeycomb lattice. We accomplish this by introducing a three-dimensional local order vector ∆(x) for the three nearest-neighbor bond order parameters of each carbon and projecting this vector onto the complete set of form factors f for the D 6h irreducible representations: an extended s-wave and two d-waves, d xy and d x 2 −y 2 , see Methods. Because the nematic state has a pure local d-wave character with a negligible s-wave component (< 0.025%) and is also completely real, each local order vector is uniquely expressed as a real-valued linear combination: captures the spatially varying d-wave orientation. As shown by the streamlines of the vector field χ(x i ) in Figs. 2(b,c), there exists a strong atomic-scale variation in the nematicity, which is aligned with the moiré-scale nematic axis in the central AA-region, but then also forms two vortex-antivortex pairs of opposite circulation outside the AA region. The two antivortices are always pinned to the center of the AB and BA regions, independently of the C 2 nematic axis orientation, while the two vortices stay close to the AA-region on opposite sides of the nematic axis, following its orientation. This atomicscale nematic pairing vortex pattern consistently appears for all investigated filling factors, coupling strengths J/t, and in both the zero temperature self-consistent linear and the non-linear gap equation solutions, making it a very robust feature of the nematic state.
We also note that this nematic pairing vortex pattern is distinct from superficially similar patterns of spontaneous supercurrents found in TBG with a chiral ground state [23,24,27], since the time-reversal invariant nematic ground state is necessarily free from supercurrents. Still, a vortex pattern has also been found in the normalstate Dirac mass term in a TBG continuum model [61], which together might suggest an overall tendency for vortex formation within large moiré patterns.

F. Nematic signatures in LDOS
The rotational symmetry breaking of the nematic superconducting state is clearly visible in the superconducting order parameter: the moiré-scale nematic C 2 axis, the intricate atomic-scale d-wave nematicity pattern, and the three-fold degenerate ground state. Nematic superconductivity is also known to exhibit magnetic field directional dependencies in various properties tied to superconductivity, such as the upper critical field, as also recently found in superconducting TBG [40].
In Fig. 3, we show that the nematic superconducting state in TBG additionally gives rise to measurable signals in the electronic local density of states (LDOS). In the main panels, (a) to (d), we plot the LDOS on sublattice A in the top layer along the corresponding four different line cuts shown in (e) to (h). In the normal state, the three-fold rotational symmetry leads to equivalent LDOS on the cuts of (a),(c) and (b),(d), respectively. In the superconducting state, however, the four cuts in (e) to (h) form an progressively increasing angle against the nematic C 2 axis (dashed line), and the corresponding LDOS in (a) to (d) show an increasingly pronounced asymmetry around the moiré cell center, unambiguously demonstrating the nematic superconducting order. In fact, the nematic superconducting state induces shifts, primarily in the coherence peaks, that are up to almost half the maximal moiré state LDOS in size. The sublattice LDOS polarization, i.e. the difference in LDOS between the A and B sublattices, is shown in (i) to (l) along the cuts, demonstrating that all shifts occur in opposite directions between the two sublattices. This means the observed anisotropy washes out in a joint LDOS, which further highlights the importance of atomistic modeling and might also explain why no nematic anisotropy in the LDOS has been found in the superconducting state in rescaled lattice models [23]. To summarize, our results show that the nematic superconducting state in TBG is clearly observable in sublattice-resolved scanning tunneling spectroscopy/microscopy (STS/STM) measurements.

G. Interlayer π-locked Josephson coupling
The superconducting order parameters at both T c and zero temperature always show a rigid interlayer π-shift in the superconducting phase and the atomic-scale nematic vectors in Fig. 2(b) consequently reverse direction between the two layers. This result is consistent with earlier continuum model results predicting a π-shift due to a large layer counterflow velocity [21,27], and was also recently reported numerically using atomistic modeling [24], albeit seemingly missing for a rescaled model [23]. This rigid π-shift suggests a strong and unconventional interlayer Josephson coupling that we here explore by manually tuning the interlayer phase of the selfconsistent order parameter, as is standard in Josephson setups.
Adding a phase factor ∆ → e iφ ∆ in the bottom layer, we plot in Fig. 4(b) the resulting superconducting condensation energy as a function of the interlayer phase difference φ. For roughly half of the available phase space, the superconducting state is unstable with a total energy higher than in the normal state. This extraordinarily stiff energy-phase relationship proves that the π-shift is central for the existence of superconductivity.
The underlying reason for the strong interlayer Josephson coupling is that the π-shift is responsible for the entire superconducting gap of the nematic d-wave superconducting state of TBG, as seen in the DOS color intensity plot in Fig. 4(a), with specific spectra extracted in Figs. 4(c-e). At the self-consistent solution, φ = π, a full superconducting gap is present, whereas it decreases and eventually disappears for φ approaching zero (or equivalently 2π). This result is in notable contrast to any realvalued d-wave state in graphene or AA stacked bilayer graphene, which are always nodal, regardless of interlayer phase difference. We thus conclude that the moiré band structure generates the interlayer π-shift, which in turn results in the nematic d-wave state of TBG being fully gapped and the ground state solution.

III. CONCLUSION
Using full-scale atomistic modeling of magic-angle TBG, we find a highly inhomogeneous and nematic dwave superconducting state with a full energy gap at experimentally observed critical temperatures. The large real-space inhomogeneity demonstrates that atomic-scale modeling is crucial for TBG. Besides atomic spatial resolution, we also find that a dense k-point sampling of the MBZ and the moiré flat band structure is required to correctly capture the superconducting order. In fact, if sampling only the MBZ center, we instead find a chiral d-wave ground state [24] at the experimentally relevant temperatures.
The only possible uncertainty left in our work is the exact nature of the electronic interactions, where we assume intralayer spin-singlet bond interactions. This choice is well motivated both by experimental evidence [2,6,43] and theoretical work on TBG [24], graphene [31,44,55], and pπ-bonded organic molecules [46]. Similar interactions are also present in the cuprate superconductors [45], with which TBG shares many phase diagram features [12]. Moreover, even if additional longer-range interactions were considered, results from the honeycomb lattice indicate that our results will remain qualitatively correct [32,55].
The only unknown parameter is then the coupling constant J/t. Our results are however remarkably stable with the same nematic ground state for all experimentally relevant coupling strengths J/t, further supporting their validity. Quantitatively, we also find experimentally observed T c :s already at very weak coupling J/t ∼ 0.4, which clearly illustrates the strong tendency of the moiré flat bands to electronically order [2,40]. Thus, the required coupling strength J/t is easily attained and also compares favorably with order of magnitude estimates of the super-exchange J/t = 4t/U ∼ 1 from graphene based on ab-initio results and a strong coupling scenario [6,44]. We also note that we calculate the superconducting mean-field pair amplitude ∆, while superconductivity in two dimensions is reached only at the Berezinskii-Kosterlitz-Thouless (BKT) transition, marking the onset of phase coherence and requiring a finite superfluid weight. Here our use of the full TBG band structure with its non-trivial topology guarantees a finite geometric superfluid weight [39,[62][63][64]. While the BKT temperature is always substantially lower than the mean-field temperature, they have been shown to correlate directly in several TBG models [39]. Thus, the mean-field pair amplitude calculated in this work should be a good measure of superconductivity, albeit it requires us to consider appreciably higher mean-field critical temperatures than those measured experimentally.
With the insulating states surviving to higher temperatures, superconductivity then primarily appears as domes flanking the insulating regions when tuning the filling. In fact, this is a universal feature of flat band systems: superconductivity always thrives further away from high DOS peaks compared to any particle-hole (insulating) order, thus naturally forming domes [53]. This universal picture is consistent with current TBG experiments and notably independent of the mechanisms behind the insulating behavior and superconductivity being the same or competing in nature. As such, our results are not sensitive to the exact nature of the insulating states. In conjunction with this, it is interesting to point out that, despite similar phase diagram features and electronic interactions for the cuprate superconductors and TBG [12], we find a remarkably different (fully gapped, highly inhomogeneous, and nematic) d-wave superconducting state in TBG. We might expect any insulating state to also display similar atomic-scale inhomogeneity, opening up for a wide range of different behaviors in TBG and other moiré systems.

IV. METHODS
Moiré lattice structure. In the first layer of TBG, the two graphene unit vectors are a 11 = ax, a 12 = (a/2)x + (3a/2)ŷ, while in the second layer, the lattice is rotated by the relative twist angle θ and spanned by the rotated unit vectors a 2i = R(θ)a 1i , where R(θ) is the rotation matrix in two dimensions. With two carbon atoms per graphene unit cell, there are in both layers two sublattices, A and B, of carbon atoms placed at {na l1 +ma l2 +δ SB η l +(d 0 δ l /2)ẑ} for n, m ∈ Z, l = 1, 2, and S ∈ A, B, where η 1 = acŷ, η 2 = R(θ)η 1 , ac = a/ √ 3,ẑ is the unit vector perpendicular to the layers, d 0 = 3.35Å, and δ l = (−1) l . The two layers produces a periodic moiré interference pattern of period length L M = a/(2 sin θ/2).
We perform calculations using a large and periodically repeated moiré unit cell with the unrelaxed lattice positions [65][66][67]. This requires using commensurate twist angles, for which the lattices of the two layers periodically match up: n 1 a 11 +n 2 a 12 = m 1 a 21 +m 2 a 22 , for integers m i , n i . Twist angles satisfying this condition are given by cos (θ) = (3q 2 − p 2 )/(3q 2 + p 2 ) and parametrized by a relative prime integer pair (gcd(p, q) = 1, greatest common divisor), with p, q ∈ N, p ≥ q > 0 [42]. Specifically, the number of carbon atoms in the moiré unit cell is Nc(p, q) = 4 gcd(3, p) p 2 + 3q 2 / gcd(p − 3q, p+3q) 2 . We choose p = 1 with q odd, as this results in the fewest number of carbon atoms within the moiré unit cell for a given twist angle. In particular, we focus in the main text on the commensurate unit cell (p, q) = (1, 55), which gives the twist angle θ ≈ 1.2 • closest to the magic angle, defined as where the Fermi velocity vanishes. We here also note that the commensurate condition is independent of the twist center. Since we use a gcd(3, p) = 1 commensurate moiré lattice constructed from an AA stacked bilayer and twisted around an axis passing through two carbon atoms, the resulting lattice has a D 3 point symmetry group. Due to the long wave length moiré pattern at small twist angles, this lattice also has an an approximate D 6 symmetry that even becomes exact if the twist axis is instead taken through the center of a honeycomb lattice hexagon [68]. Normal state Hamiltonian. We employ a standard tightbinding model including all carbon atoms: The intralayer Hamiltonian is given by graphene: with in-plane nearest neighbor hopping t between the carbon pz electrons, created by the operator c † i,σ,l on the site i, layer l, and spin σ. The overall occupancy is regulated by the chemical potential µ. Further, the interlayer hybridization has been shown to be well-captured by hopping amplitudes decaying exponentially with distance and with Koster-Slater factors for the bond angle dependence [22]: where r ij is the displacement between the carbon sites i and j and cos β z ij =r ij ·ẑ. The input parameters are fixed by matching to the electronic band structure single and AB-stacked bilayer graphene [22]: t = 2.7 eV, tσ = 0.48 eV, and λ = 0.184ac. To keep the Hamiltonian matrix sparse, we cutoff t ij for distances d > 6a, resulting in ∼ 250 interlayer bonds per atom. This long-ranged cutoff is needed to preserve all symmetries of the TBG lattice, while a larger cutoff does not affect our results. Superconducting pairing. We model superconductivity by: with a homogenous coupling strength J. We solve for the superconducting spin-singlet bond order parameters ∆ ij using both the full non-linear self-consistent and the linearized gap equations: where The non-linear self-consistency gap equation, Eq. (5), is equivalent to ∂F/∂∆ ij = 0, and thus finds the order parameters ∆ ij which minimize the Free energy F . This equaqtion is solved iteratively until convergence is reached. The linear gap equation, Eq. (6), only has solutions for the infinitesimal order field δ∆ (defined up to a complex constant) when the stability matrix S rs ij (T ) has at least one eigenvalue equal to 1, which implicitly defines the critical temperature Tc. This equation is equivalent to the Free energy Hessian ∂ 2 F/∂∆ ij ∂∆rs = ∂ s ij /∂∆rs + δ ir δ js /J changing signature with a zero eigenvalue in the normal state (∆ = 0) and the emerging order therefore lowers the Free energy below Tc. Rational pole expansion. Solving the gap equations, Eqs. (5)-(6), using the standard approach of matrix diagonalization is prohibitively expensive for TBG due to large degrees of freedom. We instead calculate the Fermi operator F β (H) = (e βH + 1) −1 using a rational pole expansion [69,70]. Specifically, we use the minimax rational approximation of J. E. Moussa that minimizes the uniform norm = max z∈[−βE min ,∞] |F β (z) − Np n=1 Rn/(βz − Pn)| for Np poles at Pn with residues Rn, because of its rapid convergence at low temperatures [71].
To illustrate the method, we first show that all single particle (anomalous) expectation values are elements of F β (H BdG (k)), within the Bogoliubov de-Gennes (BdG) formalism. In the 2Nc dimensional block Nambu-spinor basis X k = ({c k,↑ }, {c † −k,↓ }) T , the complete Hamiltonian H = H N + H SC for the Nc carbon atoms and spins has the BdG bilinear form with the constant energy shift C = −µNc + i,j |∆ ij | 2 /J. The BdG bilinear form is diagonalized by a unitary transformation with the block structure The accompanying canonical transformation defines the fermionic Bolgoliubov quasiparticles of opposite energy ±E, In this diagonal basis, the Fermi operator is where the blocks are: From the transformations of Eq. (10), we find that all single particle expectation values are indeed given by the elements of F β (H BdG (k)): where we introduce 2Nc dimensional basis vectors for the electron [e i ] j = δ ij and hole [h i ] j = δ (i+Nc)j blocks of the BdG Matrix, respectively. Next, we apply the rational pole expansion. Specifically, the (anomalous) expectation values for where we define the "propagated" vectors [βH BdG (k) − Pn] −1 q i = Π n k q i = p n i . Thus, by solving the set of linear equations for the propagated vectors we can calculate the anomalous expectation values in Eq. (5) within the rational pole expansion of the Fermi operator using Np terms. Finally, we also need access to derivatives, or generally a static response, in Eq. (6). Given a perturbation λ to H BdG , the static response of any expectation value, can also be computed using the rational pole expansion. Namely, using that ∂ λ A −1 In particular, for derivatives with respect to the superconducting order parameter ∆rs, only the off-diagonal blocks∆(k) are nonzero. In addition, H BdG (k) is block diagonal and contains only the normal state Hamiltonian when evaluated at∆ = 0, as is the case for the linear gap equation. As a consequence we find where x n i and y n j are Nc dimensional vectors satisfying the linear equations: where [i]s = δ is and [j]s = δ js . Together, Eq. (14) and Eq. (17), show how both (anomalous) expectation values and their static response are computed using the rational pole expansion. Thus solving the gap equations Eqs. (5)- (6) is simply reduced to solving the respective sets of linear equations of Eqs. (15)- (18). The main advantage of this approach is that all linear equations can be solved in parallel. Additionally, we solve these equations with the minimal residual method [72], which, based on Lanczos iterations, takes advantages of the sparseness of H BdG , and thus only requires sparse matrix-vector multiplications. The total computational time has therefore a close to linear scaling in the problem dimensions. For many Lanczos iterations, a known risk is that roundoff errors may give numerical instabilities from loss of orthogonality [73]. We have verified by comparing with direct diagonalization of H BdG for rotation angles down to θ ≈ 1.5 • , that all calculated expectation values are within the expected accuracy with no evidence of numerical instability. Solving for and analyzing the superconducting order. We solve the gap equations, Eqs. (5) and (6), using the rational pole expansion of the Fermi operator, Eqs. (14) and (17), with a maximal error less than 10 −8 %. More specifically, to solve the linear gap equation we first compute the static response matrix Γ rs ij (T ) = −∂ s ij /∂∆rs at the fixed temperatures Tc = 6.25 × 2 n 10 −5 t ≈ 2 n K for n = 0, 1, ...6, with a uniform ks × ks grid sampling of reciprocal space (ks = 4). Then the largest eigenvalues λn(Tc) of Γ were found using the eigenvalue solver PRIMME [74]. Finally, for each eigenvalue, the coupling strength of Jn(Tc) = 1/λn(Tc) ensures that Eq. (6) is satisfied. That is, at Jn(Tc) the n th (subordinate) order has the critical temperature of Tc. In Fig.2(a) we present the (Jn(Tc), Tc) pairs for the fourth largest eigenvalues at the fixed Tc.
The eigenvectors corresponding to each solution (Jn(Tc), Tc) are symmetry-classified on the moiré length scale. For this, we construct the irreducible representation (Irrep) projection operatorŝ P I = (d I /|G|) g∈G χ * I (g)Γ(g) of the point group G = D 3 , where I index the Irreps of dimension d I with group characters χ I and whereΓ(g) is the action of the symmetry operations g on all the superconducting order parameters in the moiré unit cell. Because an eigenvector of a symmetric static response matrix necessarily belongs to only a single Irrep, only the corresponding projection operator leaves the eigenvector invariant, while all other Irrep projection operators annihilates it, thereby allowing the symmetry of the eigenvector to be easily identified. Nonetheless, the same classification is also obtained with more insight by considering only the linear subspace consisting of the six order parameters on the nearest-neighbor bonds of two layer-aligned A-sites, since together these six order parameters form a closed set under the symmetry operationsΓ(g) of the D 3 symmetry group and necessarily with the same symmetry as the whole eigenvector. In particular, this linear subspace is the direct sum of the D 3 Irreps: A1, A2 and two E, with Irrep basis vectors (up to normalization) obtained from the non-zero eigenvectors of the projection operators and listed in Table I. We then easily classify the moiré-scale symmetry of each eigenvector by projecting each subset of bonds on these Irrep basis vectors.
As previously noted, the moiré lattice model also has an approximate D 6 symmetry, besides the exact D 3 symmetry [68]. We have verified that also the leading superconducting order in Fig. 2(a) has an approximate D 6 symmetry by calculating that it has a 99% weight in the E 2 Irrep of D 6 . Here the weight P I∆ 2 / ∆ 2 in the D 6 Irrep I is calculated from the projection operators with the D 6 geometric transformationsΓ(g) that interchanges all bonds of TBG with a zeroth order interpolation around the nearest hexagon center to the AA region.
Complimentary to the linear equations Eq. (6), we also solve the fully non-linear self-consistency equations Eq. (5) using Eq. (15) at an effective zero temperature (T = 10 −7 t) and with the same k-point sampling (ks = 4) as the linear equation Eq. (6). We find that the zero temperature self-consistent solutions∆ SC are always completely contained in the eigenspace Λ 1 of the leading solutions of the linear equation. Formally, this is established by the projection P Λ 1∆ SC 2 / ∆ SC 2 being as high as 0.99993 for J = 0.43t and 0.99992 for J = 1.0t, where ∆ = i,j |∆ ij | 2 is the order parameter field norm. Additional physical properties. We additionally compute a number of different physical properties using the following methods.
Band plot. For the band plot in Fig. 1(c), the 16 lowest eigenvalues of the the normal state were computed evenly along the high symmetry cuts shown in Fig. 1(a) with PRIMME [74].
Density of states. From the lowest eigenvalues computed with PRIMME, the density of states, DOS(E) = n δ(E − En), was computed as the kernel density estimation with Gaussian width Kσ, using on a regular ks × ks grid in reciprocal space for both the normal and the superconducting states, see Table II for parameters of each figure. Furthermore, in Figs. 1(d),(e), the DOS was calculated for the commensurate angles parametrized as (p, q) with p = 1 and q = 21, 23, ..., 71 and then normalized by the area of the moiré unit cell. For q < 31 we used Kσ = 1 meV and ks = 32, while for 31 ≤ q, Kσ = 1/2 meV and ks = 16. The energy of the VHSs shown in Fig.1(d) were defined as the position of the two local maximum in the resulting density estimation. Similarly, the value at the two local maxima define the VHS DOS in Fig. 1(e). For Fig. 4(a), we used a two-dimensional Gaussian Kernel density estimation with a kernel widths of 1/10 meV for the energy and π/10 for the interlayer phase.
Local density of states. From the lowest eigenvalues together with eigenvectors computed with PRIMME, the normal state local density of states (LDOS) of Fig. 1(g) was computed using where U λi (k) is the i-site amplitude of the λ eigenvector. For the superconducting state in Fig. 3 we instead used the electronic part of the LDOS: Figure(s) Spectrum ks Kσ (meV) Fig. 1(c) Normal 60 1/20 Fig. 1(d),(e) Normal 16 or 32 1/2 or 1 Fig. 2 where u iλ (k) and v iλ (k) are the particle and hole amplitudes of the Nambu-spinor eigenvector. For both cases, we computed the LDOS as the weighted Kernel density estimation with the same Gaussian kernel width Kσ = 1/8 meV and ks = 24. Charge density. The charge density of the moiré flat band Nm(x) in Fig. 1(f) was obtained by first computing the particle occupation N i (µ) = σ c † iσ c iσ using Eq. (14) at each site i of the top graphene layer for two different uniform chemical potentials: µ 1 = 7.5 meV and µ 2 = 21.4 meV. The lower (upper) µ is in the band gaps below (above) the moiré bands. The change in particle occupation between the two µs, Nm(x i ) = N i (µ 2 ) − N i (µ 1 ), is thus the charge density of the moiré bands and is also equivalent to the integrated LDOS: Nm(x) = µ 2 µ 1 n(x, E)dE. To obtain the joint charge density of both A and B sublattices with in each unit cell, the density on each sublattice was first linearly interpolated as a function of position and then summed on an hexagonal grid in Fig. 1(f).
Relative energy differences. The relative energy differences in Figs. 2(e),(f), were calculated by reinserting the linear combinations∆(Θ, ϕ) of the T = 1.25×10 −4 t eigenvectors∆x,y in to the BdG Hamiltonian with same field amplitude ∆ as the zero temperature self-consistent solution of the corresponding coupling strength J(Tc) ≈ 0.43t. The zero temperature energy difference between the∆(Θ, ϕ) state and the ground state at∆(0, 0) were then computed as ∆E(Θ, ϕ) = n,k θ(−E n k,Θ,ϕ )E n k,Θ,ϕ − θ(−E n k,0,0 )E n k,0,0 for the eigenenergy bands E n k,θ,ϕ indexed by n. Note that all other terms in the free energy cancel between the two solutions. For the k-vector sum, the same k-point sampling (ks = 4) was used as when solving the gap equations. For the results in Fig. 2(g) we calculated in the same way the relative energy differences for∆(Θ, ϕ) but over a range of field norms ∆ , i.e. not constrained to the self-consistent values at the appropriate coupling strength. Since the used basis vectors∆x,y(Tc) might change with Tc and thus field norm, this is technically an approximation. However, we find that basis vectors∆x,y(Tc) of the leading solution are remarkably stable, and that the resulting energy differences are insensitive to the∆x,y(Tc) used for the variation, a result we checked by using both eigenvectors from T = 10 −3 t and T = 1.25 × 10 −4 t, as well as a higher (ks = 6) k-point sampling. To accurately capture the energy difference also in the very small ∆ regime, we computed the energy difference using a very high k-point sampling (ks = 20).
Condensation energies. The zero temperature condensation energies in Fig. 2(h-j), as well as the the maximum condensation energy of the self-consistent φ = π state in Fig. 4(b), were computed using all the eigenvalues of the normal state and of the BdG Hamiltonian attained with Armadillo [75,76]. In terms of the eigenvalues E n of Hamiltonian, the normal state energy is E N = N −1 k n,k θ(−E n k )E n k , while the energy in the superconducting state E(∆) = N −1 k n,k θ(−E n k (∆))E n k (∆) − 2µNc + i,j ∆ ij | 2 J, also contains the constant energy shift C of the BdG form, see Eq. (7). The zero temperature condensation energy is then just the difference E cond = E N − E(∆).
Streamlines. The streamlines in Figs. 2(b),(c) were constructed by numerically integrating the flow equation ∂tr(t) = χ(r) both forwards and backwards starting from initial points on a hexagonal grid throughout the moiré cell. Here the vector field χ(x i ) = cos τ (x i )x + sin τ (x i )ŷ was the linearly interpolated from the angle τ (x i ) = arctan ∆(x i ) · f d x 2 −y 2 /∆(x i ) · f dxy , with separate accounting for the appropriate quadrant for the angle.

V. DATA AVAILABILITY
Data are available from the corresponding authors upon request.

VI. CODE AVAILABILITY
Computer program source codes used for all calculations and analysis are available for review from the corresponding authors upon request.

VII. ACKNOWLEDGMENT
We acknowledge financial support from the Swedish Research Council (Vetenskapsrådet Grant No. 2018-03488) and the Knut and Alice Wallenberg Foundation through the Wallenberg Academy Fellows program. The computations were enabled by resources in project SNIC 2020/5-18 provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX, partially funded by the Swedish Research Council through grant agreement No. 2018-05973.

VIII. AUTHOR CONTRIBUTIONS
A.B.S. and T.L. initiated the project. T.L. developed the computational framework and carried out the simulations with help from J.S. and F.P.. All authors analyzed and interpreted the results. A.B.S. and T.L. wrote the manuscript with input from all authors.