Introduction

Precise twist angle and carrier density control have made it possible to map the rich phase diagram of twisted bilayer graphene (TBG)1,2. Around the magic twist angle θ ≈ 1.1° four spin-degenerate topological moiré flat bands are formed around zero energy and the associated density of states (DOS) peaks make TBG extremely prone to electronic ordering1,2,3,4,5,6,7,8,9,10,11,12. Both correlated insulating states at integer flat band fillings (ν = ±1, 2, 3) and superconductivity have been found, with the latter appearing as superconducting domes both flanking the ν = ±2 insulating states2,13 and more broadly across the entire moiré bands5.

While the superconducting transition temperature Tc ≈ 3 K is low in TBG2,5, the ratio to the Fermi temperature is large Tc/TF ≈ 0.1, exceeding the weak coupling regime12. In fact, with superconductivity appearing in close proximity to correlated insulators and a pseudogap state with reduced DOS above the superconducting domes7,8,9 with accompanied strange metal behavior14,15, the phase diagram of TBG shares striking similarities to the high-temperature cuprate and pnictide superconductors12,16,17. This points to the possibility of strong local electronic interactions being responsible for superconductivity, although electron-phonon pairing is also a plausible candidate14,18,19,20.

The large scale of the moiré pattern (~104 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 structure21,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 proliferation23. Many of these studies have proposed a time-reversal symmetry (TRS) breaking chiral d-wave state23,24,25,26,27,28,29,30, as also found in heavily doped graphene31,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}^{2}}\)- and dxy-wave states with a symmetry enforced degeneracy at Tc, on lattices with three- or six-fold symmetry, such as graphene’s honeycomb lattice and the TBG moiré lattice32,33,34,35. On the honeycomb lattice, the chiral \({d}_{{x}^{2}-{y}^{2}}+i{d}_{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 temperature31,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 state35,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 properties23,24,39. Experimentally, recent magnetotransport measurements have also identified intrinsic nematicity in the superconducting state of TBG40. 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.

Results

Normal state properties

TBG consists of two graphene layers rotated by a twist angle θ, which produces a periodic moiré interference pattern of length \({L}_{{{{{{{{\rm{M}}}}}}}}}=a/(2\sin (\theta /2))\)21,41, far exceeding the graphene lattice constant a = 2.46 Å, see Fig. 1a. 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 pz 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 dependence22,42, see Methods section.

Fig. 1: Lattice structure and normal state properties of twisted bilayer graphene (TBG).
figure 1

a Carbon lattice in one moiré cell of TBG with twist angle θ ≈ 3.9 and moiré period LM. b Illustration of TBG in reciprocal space, where θ rotates the two hexagonal Brillouin zones (BZ) of the graphene layers (BZ1 and BZ2), displacing the Dirac cones at the BZ corners (K1 and K2) by \({{\Delta }}K=2K\sin (\theta /2)\), the reciprocal length of the moiré Brillouin zone (MBZ). Depicted are also the high symmetry lines and points (γ, μ, κ) in the MBZ. c Low-energy band structure of TBG at θ ≈ 1.2 along the high symmetry lines in b, consisting of four narrow spin-degenerate moiré bands (solid lines) separated from other bands (dotted lines) by finite energy gaps. The corresponding large density of states (DOS), peaking at the upper and lower van Hove singularities (VHSs), is plotted on same energy axis to the right. Side inset: DOS for a larger energy range, highlighting the massive moiré flat band (red) DOS peak. Inset: DOS as a function of the moiré band filling factor ν in units of electrons per moiré unit cell with VHSs at ν ≈ ± 1 (dashed lines). d Energy of the upper and lower VHSs as a function of twist angle θ. e Maximum peak height of the two VHS peaks as a function twist angle θ, peaking sharply at the magic angle. f Intensity plot of the top layer charge density Nm(x) of the completely filled moiré bands, equivalent to integrating the LDOS in g. The scale bar shows a distance of 20a, where a = 2.46 Å is the graphene lattice constant. g Local density of states (LDOS) n(x, E) as a function of energy E and position along the dashed line in f, centered on the AA region.

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. 1b. 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 structure21, see Fig. 1c. 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 structure3,21. At the magic angle, the moiré bands become essentially completely flat and the two VHSs merge, see Fig. 1d, 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 Nm(x) in Fig. 1f. Resolving the charge density with respect to energy, Fig. 1g further shows the local density of states (LDOS) along a periodic path passing the AB, AA, and BA regions (black dotted line in Fig. 1f).

Modeling superconductivity

Many properties of the superconducting state in magic-angle 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 pairing2,43. The strong on-site repulsion in graphene, graphite44, and TBG6 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 interactions31,45. In fact, preference for bond spin-singlet 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 superconductors47. Spin-singlet bond pairing has also recently been derived from spin-fluctuations in TBG24 and also used in smaller rescaled lattice models23,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 section. Together these bond order parameters forms an order parameter field \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}\) 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 Tc 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 D3 of TBG: the one-dimensional A1 and A2, or the two-dimensional E representation. Below Tc, 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.

Moiré-scale nematicity

Starting with the results from the linear gap equation, we show in Fig. 2a the four highest critical temperatures Tc 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 Tc (green highlight) belongs to the E irreducible representation and is therefore two-fold degenerate and spanned by two real-valued order parameter fields: \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x}\) and \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{y}\). In Fig. 2b we plot in color scale the normalized amplitude of the \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x}\) order parameter field in the upper (left) and lower (right) graphene layer in the moiré cell, while the amplitude of \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{y}\) in the upper layer is plotted in Fig. 2c. Both \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x}\) and \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{y}\) clearly break the C3 rotation symmetry of the normal state and are instead enhanced in the AA regions and along the C2 nematic axes aligned with the x- and y-axis, respectively.

Fig. 2: Properties of the nematic superconducting state in twisted bilayer graphene (TBG).
figure 2

a Critical temperatures Tc (left y-axis) obtained from the linear gap equation as a function of the coupling strength J/t with the Fermi level at the upper van Hove singularity (VHS). 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, \(\langle | {{{{{{{\mathbf{\Delta }}}}}}}}({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})| \rangle {k}_{{{{{{{{\rm{B}}}}}}}}}^{-1}\), (right y-axis) obtained from the non-linear self-consistency (SC) equation at zero temperature (T = 10−7t). 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 \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{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 \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{y}\) in the top graphene layer. d Color intensity plot of the order parameter field amplitude for the chiral \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x}+i{\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{y}\) solution in the top graphene layer. e Relative free energy landscape at T = 0 for all independent linear combinations \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}({{\Theta }},\varphi )\) of the \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{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 \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x}+i{\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{y}\) combination. hj 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.

In contrast to pristine graphene where superconductivity only develops for J/t above the quantum critical point (QCP) Jc = 1.76t31,48,49, TBG shows substantially enhanced ordering in Fig. 2a, with finite a Tc even at very weak coupling. For weak coupling, Tc is approximately proportional to J/t, as expected for an isolated flat band50,51,52,53. For slightly stronger coupling, the growth in Tc however accelerates, especially when approaching the graphene QCP. This both indicates contributions from dispersive bands39 and suggests that the low energy moiré flat bands have a catalytic effect that at least partially triggers the QCP cascade boosting Tc at couplings also below the QCP.

Alongside Tc extracted from the linear gap equation, we also plot in Fig. 2a (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, Tc 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 Tc ≈ 4 K2,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. 2a. The superconducting amplitude predictably drops sharply at the two moiré band edges. In between, the amplitude is rather stable but has two local maxima 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. 1c.

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 \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x}\) and \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{y}\) of the linear gap equation, see Methods. This makes TBG a nematic superconductor54, 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.

Degenerate ground state

To further investigate the ground state symmetry at zero temperature, we introduce the parametrization \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}({{\Theta }},\varphi )=\parallel \hat{{{{{{{{\mathbf{\Delta }}}}}}}}}\parallel \left(\cos {{\Theta }}{\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x}+{e}^{i\varphi }\sin {{\Theta }}{\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{y}\right)\). Here, \(\parallel \hat{{{{{{{{\mathbf{\Delta }}}}}}}}}\parallel\) is the norm of the order parameter field, the nematic angle Θ controls the angle between the C2 nematic axis and the x-axis, and a relative complex phase, φ ≠ 0, breaks TRS. Here, the (Θ, φ)  [0, π/2] × [0, 2π] manifold is a (Bloch) sphere because the Θ = 0, π/2 lines collapse to points due to gauge invariance.

At Tc all the states \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}({{\Theta }},\varphi )\) 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. 2e the relative free energy of all possible \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}({{\Theta }},\varphi )\) 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. 2e) TRS breaking complex combinations \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}(\pi /4,\pm \pi /2)={\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x}\pm i{\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{y}\). The free energy maxima (i.e. smallest condensation energy) are instead attained around the rotationally symmetric (see Fig. 2e) TRS breaking complex combinations \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{\pm }=\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}(\pi /4,\pm \pi /2)={\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x}\pm i{\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{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. 2f. As a result, \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x}\) is one of the three gauge-inequivalent nematic ground states, all with the nematic C2-axis directed towards one of the next-nearest-neighbor AA regions. In contrast, \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{y}\) is one of the three least stable gauge-inequivalent nematic states that instead has the C2-axis directed towards one of the nearest-neighbor AA regions, as seen in Fig. 2b, c.

The amplitude of the energy splitting among the real-valued nematic states depends non-monotonically on the field norm (or equivalently on J or Tc), as seen in Fig. 2g, with notably large values, up to 10% of kBTc, around the experimentally observed Tc.

In contrast, the energy of the TRS breaking chiral state \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{\pm }\) relative to the nematic ground state is large around the experimentally observed critical temperatures and also increases monotonically with the field norm (or Tc 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 regime23,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, \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{\pm }\), just barely becomes the ground state (beyond the resolution used in Fig. 2g) 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 \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}({{\Theta }},\varphi )\) states, as revealed by the low-energy DOS in Fig. 2h–j at an experimentally relevant temperature. The ground state \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x}=\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}(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 \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{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 \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{\pm }\), 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 C2 nematic axis affecting both the condensation energy and gap structure. In particular, the real-valued 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 gap31,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 found56,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 C2 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 model23,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 nematicity35,36,37.

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 D6h 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 D6h irreducible representations: an extended s-wave and two d-waves, dxy 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: \({{{{{{{\mathbf{\Delta }}}}}}}}({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})=| {{{{{{{\mathbf{\Delta }}}}}}}}({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})| \big(\cos \tau ({{{{{{{{\boldsymbol{x}}}}}}}}}_{i}){{{{{{{{\boldsymbol{f}}}}}}}}}_{{d}_{{x}^{2}-{y}^{2}}}+\sin \tau ({{{{{{{{\boldsymbol{x}}}}}}}}}_{i}){{{{{{{{\boldsymbol{f}}}}}}}}}_{{d}_{xy}}\big)\), where the vector field \({{{{{{{\boldsymbol{\chi }}}}}}}}({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})=\cos \tau ({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})\hat{{{{{{{{\boldsymbol{x}}}}}}}}}+\sin \tau ({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})\hat{{{{{{{{\boldsymbol{y}}}}}}}}}\) captures the spatially varying d-wave orientation. As shown by the streamlines of the vector field χ(xi) in Fig. 2b, 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 C2 nematic axis orientation, while the two vortices stay close to the AA-region on opposite sides of the nematic axis, following its orientation. This atomic-scale 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 state23,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 normal-state Dirac mass term in a TBG continuum model61, which together might suggest an overall tendency for vortex formation within large moiré patterns.

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 C2 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 TBG40.

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–d, we plot the LDOS on sublattice A in the top layer along the corresponding four different line cuts shown in e–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–h form an progressively increasing angle against the nematic C2 axis (dashed line), and the corresponding LDOS in a–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–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 models23. 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.

Fig. 3: Nematic signatures in the local density of states (LDOS).
figure 3

ad Sublattice-resolved (sublattice A) top layer electronic LDOS nA of the nematic superconducting state contour plotted along four different real space cuts (black line shown in eh). The moiré cell center, x/LM = 0, is marked by dashed lines. eh Illustrations of the line cuts along which the LDOS is shown in ad, starting at the base and ending at the arrow tip. The cuts are shown on top of the superconducting order parameter field amplitude where the dashed line shows the nematic axis of the order. In the normal state, the cuts a and c are equivalent by symmetry, and so are cuts b and d. The normal state LDOS is therefore the same in each pair of cuts and also symmetric around the moiré cell center, but this symmetry is lifted by the nematic superconducting order in ad. il Contour plots of the sublattice LDOS polarization, i.e the LDOS difference between sublattice A and B, nA − nB, along the cuts of eh, where green shades indicate a larger LDOS on sublattice A, while red shades indicate a larger LDOS on sublattice B. The LDOS polarization is measured relative to the maximal LDOS, Max n, and each contour marks a 10% change relative to this maximal LDOS. In all figures, the coupling strength is J = 0.43t and the Fermi level is aligned with the upper van Hove singularity (conclusions are qualitatively unchanged for other J/t).

Interlayer π-locked Josephson coupling

The superconducting order parameters at both Tc and zero temperature always show a rigid interlayer π-shift in the superconducting phase and the atomic-scale nematic vectors in Fig. 2b 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 velocity21,27, and was also recently reported numerically using atomistic modeling24, albeit seemingly missing for a rescaled model23. This rigid π-shift suggests a strong and unconventional interlayer Josephson coupling that we here explore by manually tuning the interlayer phase of the self-consistent order parameter, as is standard in Josephson setups.

Adding a phase factor ΔeiϕΔ in the bottom layer, we plot in Fig. 4b 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.

Fig. 4: Interlayer Josephson coupling.
figure 4

a Color intensity plot of the density of states (DOS) as a function of interlayer phase difference ϕ, tuned by shifting the phase of the order parameter in the bottom layer starting from the the self-consistent solution found at ϕ = π. Away from ϕ = π energy gap and coherence peaks are monotonically reduced with the gap closing around ϕ = 0, 2π. b Superconducting condensation energy, i.e. energy of superconducting state relative to the normal state (N) as a function of interlayer phase ϕ. Pink region marks positive condensation energy, describing an unstable superconducting state. ce DOS for interlayer phases ϕ = 0, π/2, π, respectively, taken as vertical cuts in a. Here, the coupling constant is J = 0.43t and Fermi level at the upper van Hove singularity (conclusions are qualitatively unchanged for other J).

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. 4a, with specific spectra extracted in Fig. 4c–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 real-valued 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.

Conclusion

Using full-scale atomistic modeling of magic-angle TBG, we find a highly inhomogeneous and nematic d-wave 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 state24 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 evidence2,6,43 and theoretical work on TBG24, graphene31,44,55, and pπ-bonded organic molecules46. Similar interactions are also present in the cuprate superconductors45, with which TBG shares many phase diagram features12. Moreover, even if additional longer-range interactions were considered, results from the honeycomb lattice indicate that our results will remain qualitatively correct32,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 Tc:s already at very weak coupling J/t ~ 0.4, which clearly illustrates the strong tendency of the moiré flat bands to electronically order2,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 scenario6,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 weight39,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 models39. 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.

Finally, let us comment on the correlated insulating states existing at integer fillings1,2,3,4,5,6,7,8,9,10,11,12. In addition to the superconducting state, a complete phase diagram should be supplemented by these interspersing insulating states. 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 domes53. 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 TBG12, 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.

Methods

Moiré lattice structure

In the first layer of TBG, the two graphene unit vectors are \({{{{{{{{\boldsymbol{a}}}}}}}}}_{11}=a\hat{{{{{{{{\boldsymbol{x}}}}}}}}}\), \({{{{{{{{\boldsymbol{a}}}}}}}}}_{12}=(a/2)\hat{{{{{{{{\boldsymbol{x}}}}}}}}}+(3a/2)\hat{{{{{{{{\boldsymbol{y}}}}}}}}}\), while in the second layer, the lattice is rotated by the relative twist angle θ and spanned by the rotated unit vectors a2i = R(θ)a1i, 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 \(\{n{{{{{{{{\boldsymbol{a}}}}}}}}}_{l1}+m{{{{{{{{\boldsymbol{a}}}}}}}}}_{l2}+{\delta }_{SB}{{{{{{{{\boldsymbol{\eta }}}}}}}}}_{l}+({d}_{0}{\delta }_{l}/2)\hat{{{{{{{{\boldsymbol{z}}}}}}}}}\}\) for \(n,m\in {{{{{{{\mathcal{Z}}}}}}}}\), l = 1, 2, and SA, B, where \({{{{{{{{\boldsymbol{\eta }}}}}}}}}_{1}={a}_{{{{{{{{\rm{c}}}}}}}}}\hat{{{{{{{{\boldsymbol{y}}}}}}}}}\), η2 = R(θ)η1, \({a}_{{{{{{{{\rm{c}}}}}}}}}=a/\sqrt{3}\), \(\hat{{{{{{{{\boldsymbol{z}}}}}}}}}\) is the unit vector perpendicular to the layers, d0 = 3.35 Å, and δl = (−1)l. The two layers produces a periodic moiré interference pattern of period length \({L}_{{{{{{{{\rm{M}}}}}}}}}=a/(2\sin \theta /2)\).

We perform calculations using a large and periodically repeated moiré unit cell with the unrelaxed lattice positions65,66,67. This requires using commensurate twist angles, for which the lattices of the two layers periodically match up: n1a11 + n2a12 = m1a21 + m2a22, for integers mi, ni. Twist angles satisfying this condition are given by \(\cos \left(\theta \right)=(3{q}^{2}-{p}^{2})/(3{q}^{2}+{p}^{2})\) and parametrized by a relative prime integer pair (gcd(p, q) = 1, greatest common divisor), with \(p,q\in {\mathbb{N}}\), p ≥ q > 042. Specifically, the number of carbon atoms in the moiré unit cell is \({N}_{{{{{{{{\rm{c}}}}}}}}}(p,q)=4\gcd (3,p)\left({p}^{2}+3{q}^{2}\right)/\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 D3 point symmetry group. Due to the long wave length moiré pattern at small twist angles, this lattice also has an an approximate D6 symmetry that even becomes exact if the twist axis is instead taken through the center of a honeycomb lattice hexagon68.

Normal state Hamiltonian

We employ a standard tight-binding model including all carbon atoms: HN = Hintra + Hinter. The intralayer Hamiltonian is given by graphene:

$${H}_{{{{{{{{\rm{intra}}}}}}}}}=-\mu \mathop{\sum}\limits_{i,\sigma ,l}{c}_{i\sigma l}^{{{{\dagger}}} }{c}_{i\sigma l}-t\mathop{\sum}\limits_{{\langle i,j\rangle}\atop{\sigma ,l}}\left({c}_{i\sigma l}^{{{{\dagger}}} }{c}_{j\sigma l}+\,{{{{\mathrm{H}}}}}.{{{{\mathrm{c}}}}}\,.\right),$$
(1)

with in-plane nearest neighbor hopping t between the carbon pz electrons, created by the operator \({c}_{i,\sigma ,l}^{{{{\dagger}}} }\) 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 dependence22:

$${H}_{{{{{{{{\rm{inter}}}}}}}}}=-\mathop{\sum}\limits_{i,j}{t}_{ij}\left({c}_{i\sigma 1}^{{{{\dagger}}} }{c}_{j\sigma 2}+{{{{\mathrm{H}}}}}.{{{{\mathrm{c}}}}}\right)$$
(2)
$${t}_{ij}=-t{e}^{({a}_{{{{{{{{\rm{c}}}}}}}}}-{r}_{ij})/\lambda }{\sin }^{2}({\beta }_{ij}^{z})+{t}_{\sigma }{e}^{({d}_{0}-{r}_{ij})/\lambda }{\cos }^{2}{\beta }_{ij}^{z},$$
(3)

where rij is the displacement between the carbon sites i and j and \(\cos {\beta }_{ij}^{z}={\hat{{{{{{{{\boldsymbol{r}}}}}}}}}}_{ij}\cdot \hat{{{{{{{{\boldsymbol{z}}}}}}}}}\). The input parameters are fixed by matching to the electronic band structure single and AB-stacked bilayer graphene22: t = 2.7 eV, tσ = 0.48 eV, and λ = 0.184ac. To keep the Hamiltonian matrix sparse, we cutoff tij 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:

$${H}_{{{{{{{{\rm{SC}}}}}}}}}=\mathop{\sum}\limits_{\langle i,j\rangle }{{{\Delta }}}_{ij}\left({c}_{i\uparrow }^{{{{\dagger}}} }{c}_{j\downarrow }^{{{{\dagger}}} }-{c}_{i\downarrow }^{{{{\dagger}}} }{c}_{j\uparrow }^{{{{\dagger}}} }\right)+{{{{{{{\rm{H.c.}}}}}}}}+\frac{{\left|{{{\Delta }}}_{ij}\right|}^{2}}{J},$$
(4)

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:

$${{{\Delta }}}_{ij}=-J\langle {s}_{ij}\rangle$$
(5)
$$\delta {{{\Delta }}}_{ij}=\mathop{\sum}\limits_{\langle r,s\rangle }\left[-J\frac{\partial \langle {s}_{ij}\rangle }{\partial {{{\Delta }}}_{rs}}\right]\delta {{{\Delta }}}_{rs}=\mathop{\sum}\limits_{\langle r,s\rangle }{S}_{ij}^{rs}(T)\delta {{{\Delta }}}_{rs},$$
(6)

where \({s}_{ij}^{{{{\dagger}}} }={c}_{i\uparrow }^{{{{\dagger}}} }{c}_{j\downarrow }^{{{{\dagger}}} }-{c}_{i\downarrow }^{{{{\dagger}}} }{c}_{j\uparrow }^{{{{\dagger}}} }\). The non-linear self-consistency gap equation, Eq. (5), is equivalent to \(\partial F/\partial {\bar{\Delta }}_{ij}=0\), and thus finds the order parameters Δij which minimize the Free energy F. This equation is solved iteratively until convergence is reached. The linear gap equation, Eq. (6), only has solutions for the infinitesimal order field \(\delta \hat{{{{{{{{\mathbf{\Delta }}}}}}}}}\) (defined up to a complex constant) when the stability matrix \({S}_{ij}^{rs}(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 \({\partial }^{2}F/\partial {\bar{{{\Delta }}}}_{ij}\partial {{{\Delta }}}_{rs}=\partial \langle {s}_{ij}\rangle /\partial {{{\Delta }}}_{rs}+{\delta }_{ir}{\delta }_{js}/J\) changing signature with a zero eigenvalue in the normal state (\(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}=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}_{\beta }(H)={({e}^{\beta H}+1)}^{-1}\) using a rational pole expansion69,70. Specifically, we use the minimax rational approximation of J. E. Moussa that minimizes the uniform norm \(\epsilon =\mathop{\max }\nolimits_{z\in [-\beta {E}_{min},\infty ]}| {F}_{\beta }(z)-\mathop{\sum }\nolimits_{n = 1}^{{N}_{{{{{{{{\rm{p}}}}}}}}}}{R}_{n}/(\beta z-{P}_{n})|\) for Np poles at Pn with residues Rn, because of its rapid convergence at low temperatures71.

To illustrate the method, we first show that all single particle (anomalous) expectation values are elements of Fβ(HBdG(k)), within the Bogoliubov de-Gennes (BdG) formalism. In the 2Nc dimensional block Nambu-spinor basis \({X}_{{{{{{{{\boldsymbol{k}}}}}}}}}={(\{{c}_{{{{{{{{\boldsymbol{k}}}}}}}},\uparrow }\},\{{c}_{-{{{{{{{\boldsymbol{k}}}}}}}},\downarrow }^{{{{\dagger}}} }\})}^{{{{{{{{\rm{T}}}}}}}}}\), the complete Hamiltonian H = HN + HSC for the Nc carbon atoms and spins has the BdG bilinear form

$$H =\,\mathop{\sum} \limits_{{{\boldsymbol{k}}}} {X}_{{{\boldsymbol{k}}}}^{{\dagger}} \left( \begin{array}{cc} {H}_{{{\rm{N}}}} ({{{\boldsymbol{k}}}}) & {\Delta}({{{\boldsymbol{k}}}}) \\ {\Delta}^{{\dagger}}({{{\boldsymbol{k}}}}) & -{H}_{{{\rm{N}}}}^{{{\rm{T}}}}(-{{{\boldsymbol{k}}}}) \end{array} \right) {X}_{{{\boldsymbol{k}}}} + C\\ =\, {X}_{{{\boldsymbol{k}}}}^{{\dagger}} {H}_{{{\mathrm{BdG}}}}({{{\boldsymbol{k}}}}) {X}_{{{\boldsymbol{k}}}} + C$$
(7)

with the constant energy shift C = − μNc + ∑i, jΔij2/J. The BdG bilinear form is diagonalized by a unitary transformation with the block structure

$$U({{{{{{{\boldsymbol{k}}}}}}}})=\left(\begin{array}{cc}u({{{{{{{\boldsymbol{k}}}}}}}})&v({{{{{{{\boldsymbol{k}}}}}}}})\\ \bar{v}(-{{{{{{{\boldsymbol{k}}}}}}}})&\bar{u}(-{{{{{{{\boldsymbol{k}}}}}}}})\end{array}\right)$$
(8)
$$U{({{{{{{{\boldsymbol{k}}}}}}}})}^{{{{\dagger}}} }{H}_{{{{{\mathrm{BdG}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})U({{{{{{{\boldsymbol{k}}}}}}}})={{{{{{{\mathcal{E}}}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})=\,{{{{\mathrm{diag}}}}}\,(\{{E}_{{{{{{{{\boldsymbol{k}}}}}}}}}\},\{-{E}_{-{{{{{{{\boldsymbol{k}}}}}}}}}\}).$$
(9)

The accompanying canonical transformation defines the fermionic Bolgoliubov quasiparticles of opposite energy ± E,

$${Y}_{{{{{{{{\boldsymbol{k}}}}}}}}}={(\{{\gamma }_{{{{{{{{\boldsymbol{k}}}}}}}},\uparrow },\}\{{\gamma }_{-{{{{{{{\boldsymbol{k}}}}}}}},\downarrow }^{{{{\dagger}}} }\})}^{{{{{{{{\rm{T}}}}}}}}}\quad \quad {X}_{{{{{{{{\boldsymbol{k}}}}}}}}}=U({{{{{{{\boldsymbol{k}}}}}}}}){Y}_{{{{{{{{\boldsymbol{k}}}}}}}}}.$$
(10)

In this diagonal basis, the Fermi operator is

$${F}_{\beta }({H}_{{{{{\mathrm{BdG}}}}}}({{{{{{{\boldsymbol{k}}}}}}}}))=U({{{{{{{\boldsymbol{k}}}}}}}}){F}_{\beta }({{{{{{{\mathcal{E}}}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})){U}^{{{{\dagger}}} }({{{{{{{\boldsymbol{k}}}}}}}})=\left(\begin{array}{cc}{G}^{{ < }}({{{{{{{\boldsymbol{k}}}}}}}})&F({{{{{{{\boldsymbol{k}}}}}}}})\\ {F}^{{{{\dagger}}} }({{{{{{{\boldsymbol{k}}}}}}}})&{G}^{ { > }}({{{{{{{\boldsymbol{k}}}}}}}})\end{array}\right)$$
(11)

where the blocks are:

$${G}^{{ < }}({{{{{{{\boldsymbol{k}}}}}}}}) =\, u({{{{{{{\boldsymbol{k}}}}}}}}){F}_{\beta }({E}_{{{{{{{{\boldsymbol{k}}}}}}}}}){u}^{{{{\dagger}}} }({{{{{{{\boldsymbol{k}}}}}}}})+v({{{{{{{\boldsymbol{k}}}}}}}}){F}_{\beta }(-{E}_{-{{{{{{{\boldsymbol{k}}}}}}}}}){v}^{{{{\dagger}}} }({{{{{{{\boldsymbol{k}}}}}}}})\\ {G}^{{ > }}({{{{{{{\boldsymbol{k}}}}}}}}) =\, \bar{v}(-{{{{{{{\boldsymbol{k}}}}}}}}){F}_{\beta }({E}_{{{{{{{{\boldsymbol{k}}}}}}}}}){v}^{{{{{{{{\rm{T}}}}}}}}}(-{{{{{{{\boldsymbol{k}}}}}}}})+\bar{u}(-{{{{{{{\boldsymbol{k}}}}}}}}){F}_{\beta }(-{E}_{-{{{{{{{\boldsymbol{k}}}}}}}}}){u}^{{{{{{{{\rm{T}}}}}}}}}(-{{{{{{{\boldsymbol{k}}}}}}}})\\ F({{{{{{{\boldsymbol{k}}}}}}}}) =\, u({{{{{{{\boldsymbol{k}}}}}}}}){F}_{\beta }({E}_{{{{{{{{\boldsymbol{k}}}}}}}}}){v}^{{{{{{{{\rm{T}}}}}}}}}(-{{{{{{{\boldsymbol{k}}}}}}}})+v({{{{{{{\boldsymbol{k}}}}}}}}){F}_{\beta }(-{E}_{-{{{{{{{\boldsymbol{k}}}}}}}}}){u}^{{{{{{{{\rm{T}}}}}}}}}(-{{{{{{{\boldsymbol{k}}}}}}}})$$

From the transformations of Eq. (10), we find that all single-particle expectation values are indeed given by the elements of Fβ(HBdG(k)):

$$\langle {c}_{{{{{{{{\boldsymbol{k}}}}}}}}i}^{{{{\dagger}}} }{c}_{{{{{{{{\boldsymbol{k}}}}}}}}j}\rangle = \mathop{\sum}\limits_{\nu }\big[{u}_{j\nu }({{{{{{{\boldsymbol{k}}}}}}}}){F}_{\beta }({E}_{{{{{{{{\boldsymbol{k}}}}}}}}\nu }){u}_{\nu i}^{{{{\dagger}}} }({{{{{{{\boldsymbol{k}}}}}}}})\\ +{v}_{j\nu }({{{{{{{\boldsymbol{k}}}}}}}}){F}_{\beta }(-{E}_{-{{{{{{{\boldsymbol{k}}}}}}}}\nu }){v}_{\nu i}^{{{{\dagger}}} }({{{{{{{\boldsymbol{k}}}}}}}})\big]\\ = \,{{{{{{{{\boldsymbol{e}}}}}}}}}_{j}^{{{{\dagger}}} }{F}_{\beta }({H}_{{{{{\mathrm{BdG}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})){{{{{{{{\boldsymbol{e}}}}}}}}}_{i}={[{G}^{{ < }}({{{{{{{\boldsymbol{k}}}}}}}})]}_{ji}$$
(12)
$$\langle {c}_{-{{{{{{{\boldsymbol{k}}}}}}}}i}{c}_{{{{{{{{\boldsymbol{k}}}}}}}}j}\rangle = \, \mathop{\sum}\limits_{\nu }\big[{u}_{j\nu }({{{{{{{\boldsymbol{k}}}}}}}}){F}_{\beta }({E}_{{{{{{{{\boldsymbol{k}}}}}}}}\nu }){v}_{\nu i}^{{{{{{{{\rm{T}}}}}}}}}(-{{{{{{{\boldsymbol{k}}}}}}}})\\ +{v}_{j\nu }({{{{{{{\boldsymbol{k}}}}}}}}){F}_{\beta }(-{E}_{-{{{{{{{\boldsymbol{k}}}}}}}}\nu }){u}_{\nu i}^{{{{{{{{\rm{T}}}}}}}}}(-{{{{{{{\boldsymbol{k}}}}}}}})\big]\\ = \, {{{{{{{{\boldsymbol{e}}}}}}}}}_{j}^{{{{\dagger}}} }{F}_{\beta }({H}_{{{{{\mathrm{BdG}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})){{{{{{{{\boldsymbol{h}}}}}}}}}_{i}={[F({{{{{{{\boldsymbol{k}}}}}}}})]}_{ji}$$
(13)

where we introduce 2Nc dimensional basis vectors for the electron \({[{{{{{{{{\boldsymbol{e}}}}}}}}}_{i}]}_{j}={\delta }_{ij}\) and hole \({[{{{{{{{{\boldsymbol{h}}}}}}}}}_{i}]}_{j}={\delta }_{(i+{N}_{{{{{{{{\rm{c}}}}}}}}})j}\) blocks of the BdG Matrix, respectively.

Next, we apply the rational pole expansion. Specifically, the (anomalous) expectation values for qi = ei (qi = hi) are

$$\langle {c}_{-{{{{{{{\boldsymbol{k}}}}}}}}i}^{({{{\dagger}}} )}{c}_{{{{{{{{\boldsymbol{k}}}}}}}}j}\rangle =\, {{{{{{{{\boldsymbol{e}}}}}}}}}_{j}^{{{{\dagger}}} }{F}_{\beta }({H}_{{{{{\mathrm{BdG}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})){{{{{{{{\boldsymbol{q}}}}}}}}}_{i}\\ \approx\, \mathop{\sum }\limits_{n=1}^{{N}_{{{{{{{{\rm{p}}}}}}}}}}{R}_{n}{{{{{{{{\boldsymbol{e}}}}}}}}}_{j}^{{{{\dagger}}} }{\left[\beta {H}_{{{{{\mathrm{BdG}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})-{P}_{n}\right]}^{-1}{{{{{{{{\boldsymbol{q}}}}}}}}}_{i}\\ =\, \mathop{\sum }\limits_{n=1}^{{N}_{{{{{{{{\rm{p}}}}}}}}}}{R}_{n}{{{{{{{{\boldsymbol{e}}}}}}}}}_{j}^{{{{\dagger}}} }{{{\Pi }}}_{{{{{{{{\boldsymbol{k}}}}}}}}}^{n}{{{{{{{{\boldsymbol{q}}}}}}}}}_{i}=\mathop{\sum }\limits_{n=1}^{{N}_{{{{{{{{\rm{p}}}}}}}}}}{R}_{n}{{{{{{{{\boldsymbol{e}}}}}}}}}_{j}^{{{{\dagger}}} }{{{{{{{{\boldsymbol{p}}}}}}}}}_{i}^{n},$$
(14)

where we define the “propagated" vectors \({\left[\beta {H}_{{{{{\mathrm{BdG}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})-{P}_{n}\right]}^{-1}{{{{{{{{\boldsymbol{q}}}}}}}}}_{i}={{{\Pi }}}_{{{{{{{{\boldsymbol{k}}}}}}}}}^{n}{{{{{{{{\boldsymbol{q}}}}}}}}}_{i}={{{{{{{{\boldsymbol{p}}}}}}}}}_{i}^{n}\). Thus, by solving the set of linear equations for the propagated vectors

$$\left[\beta {H}_{{{{{\mathrm{BdG}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})-{P}_{n}\right]{{{{{{{{\boldsymbol{p}}}}}}}}}_{i}^{n}={{{{{{{{\boldsymbol{q}}}}}}}}}_{i},$$
(15)

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 HBdG, the static response of any expectation value, can also be computed using the rational pole expansion. Namely, using that \({\partial }_{\lambda }{A}_{\lambda }^{-1}=-{A}_{\lambda }^{-1}[{\partial }_{\lambda }{A}_{\lambda }]{A}_{\lambda }^{-1}\) for a matrix Aλ, the static response is,

$$\frac{\partial \langle {c}_{-{{{{{{{\boldsymbol{k}}}}}}}}i}^{({{{\dagger}}} )}{c}_{{{{{{{{\boldsymbol{k}}}}}}}}j}\rangle }{\partial \lambda }\approx -\mathop{\sum }\limits_{n=1}^{{N}_{{{{{{{{\rm{p}}}}}}}}}}\beta {R}_{n}{{{{{{{{\boldsymbol{e}}}}}}}}}_{j}^{{{{\dagger}}} }\left[{{{\Pi }}}_{{{{{{{{\boldsymbol{k}}}}}}}}}^{n}\frac{\partial {H}_{{{{{\mathrm{BdG}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})}{\partial \lambda }{{{\Pi }}}_{{{{{{{{\boldsymbol{k}}}}}}}}}^{n}\right]{{{{{{{{\boldsymbol{q}}}}}}}}}_{i}.$$
(16)

In particular, for derivatives with respect to the superconducting order parameter Δrs, only the off-diagonal blocks \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})\) are non-zero. In addition, HBdG(k) is block diagonal and contains only the normal state Hamiltonian when evaluated at \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}=0\), as is the case for the linear gap equation. As a consequence we find

$$\frac{\partial \langle {c}_{-{{{{{{{\boldsymbol{k}}}}}}}}i}{c}_{{{{{{{{\boldsymbol{k}}}}}}}}j}\rangle }{\partial {{{\Delta }}}_{rs}}\approx \mathop{\sum }\limits_{n=1}^{{N}_{{{{{{{{\rm{p}}}}}}}}}}\beta {R}_{n}{({{{{{{{{\boldsymbol{y}}}}}}}}}_{j}^{n})}^{{{{\dagger}}} }\left[\frac{\partial \hat{{{{{{{{\mathbf{\Delta }}}}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})}{\partial {{{\Delta }}}_{rs}}\right]{{{{{{{{\boldsymbol{x}}}}}}}}}_{i}^{n}$$
(17)

where \({{{{{{{{\boldsymbol{x}}}}}}}}}_{i}^{n}\) and \({{{{{{{{\boldsymbol{y}}}}}}}}}_{j}^{n}\) are Nc dimensional vectors satisfying the linear equations:

$$\left[\beta {H}_{{{{{{{{\rm{N}}}}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})-{P}_{n}\right]{{{{{{{{\boldsymbol{x}}}}}}}}}_{i}^{n}={{{{{{{\boldsymbol{i}}}}}}}}\quad \left[\beta {\bar{H}}_{{{{{{{{\rm{N}}}}}}}}}(-{{{{{{{\boldsymbol{k}}}}}}}})+{\bar{P}}_{n}\right]{{{{{{{{\boldsymbol{y}}}}}}}}}_{j}^{n}={{{{{{{\boldsymbol{j}}}}}}}}$$
(18)

where [i]s = δis and [j]s = δjs.

Together, Eqs. (14) and (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 method72, which, based on Lanczos iterations, takes advantages of the sparseness of HBdG, 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 orthogonality73. We have verified by comparing with direct diagonalization of HBdG 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 <10−8%. More specifically, to solve the linear gap equation we first compute the static response matrix \({{{\Gamma }}}_{ij}^{rs}(T)=-\partial \langle {s}_{ij}\rangle /\partial {{{\Delta }}}_{rs}\) at the fixed temperatures Tc = 6.25 × 2n10−5t ≈ 2n 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 PRIMME74. 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 nth (subordinate) order has the critical temperature of Tc. In Fig. 2a 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 operators \({\hat{P}}_{I}=({d}_{I}/| G| ){\sum }_{g\in G}{\chi }_{I}^{* }(g)\hat{{{\Gamma }}}(g)\) of the point group G = D3, where I index the Irreps of dimension dI with group characters χI and where \(\hat{{{\Gamma }}}(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 \(\hat{{{\Gamma }}}(g)\) of the D3 symmetry group and necessarily with the same symmetry as the whole eigenvector. In particular, this linear subspace is the direct sum of the D3 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 1. We then easily classify the moiré-scale symmetry of each eigenvector by projecting each subset of bonds on these Irrep basis vectors.

Table 1 Irreducible representation decomposition.

As previously noted, the moiré lattice model also has an approximate D6 symmetry, besides the exact D3 symmetry68. We have verified that also the leading superconducting order in Fig. 2a has an approximate D6 symmetry by calculating that it has a 99% weight in the E2 Irrep of D6. Here the weight \(\parallel {\hat{P}}_{I}\hat{{{\Delta }}}{\parallel }^{2}\!/\!\parallel \hat{{{\Delta }}}{\parallel }^{2}\) in the D6 Irrep I is calculated from the projection operators with the D6 geometric transformations \(\hat{{{\Gamma }}}(g)\) that interchanges all bonds of TBG with a zeroth order interpolation around the nearest hexagon center to the AA region.

We also extract, throughout the moiré unit cell, the local atomic-scale symmetry of the superconducting order. We do this by using the basis vectors of the Irreps of the graphene symmetry group D6h. Expressed on all nearest neighbor bonds these are \({{{{{{{{\boldsymbol{f}}}}}}}}}_{s}=(1,1,1)/\sqrt{3}\) for the identity representation A1g (s-wave symmetry), \({{{{{{{{\boldsymbol{f}}}}}}}}}_{{d}_{{x}^{2}-{y}^{2}}}=(2,-1,-1)/\sqrt{6}\) (\({d}_{{x}^{2}-{y}^{2}}\)-wave) and \({{{{{{{{\boldsymbol{f}}}}}}}}}_{{d}_{xy}}=(0,-1,1)/\sqrt{2}\) (dxy-wave) for the two-dimensional E2g representation31. The atomic-scale symmetry of the three-dimensional bond order parameter vector \({{{{{{{\mathbf{\Delta }}}}}}}}({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})=({{{\Delta }}}_{i{j}_{1}},{{{\Delta }}}_{i{j}_{2}},{{{\Delta }}}_{i{j}_{3}})\) at site xi is then extracted by projection on to these form factors.

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−7t) and with the same k-point sampling (ks = 4) as the linear equation Eq. (6). We find that the zero temperature self-consistent solutions \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{{{{{{{{\rm{SC}}}}}}}}}\) are always completely contained in the eigenspace Λ1 of the leading solutions of the linear equation. Formally, this is established by the projection \(\parallel {P}_{{{{\Lambda }}}_{1}}{\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{{{{{{{{\rm{SC}}}}}}}}}{\parallel }^{2}\!/\!\parallel {\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{{{{{{{{\rm{SC}}}}}}}}}{\parallel }^{2}\) being as high as 0.99993 for J = 0.43t and 0.99992 for J = 1.0t, where \(\parallel \hat{{{{{{{{\mathbf{\Delta }}}}}}}}}\parallel =\sqrt{{\sum }_{\langle i,j\rangle }| {{{\Delta }}}_{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. 1c, the 16 lowest eigenvalues of the normal state were computed evenly along the high symmetry cuts shown in Fig. 1a with PRIMME74.

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 2 for parameters of each figure. Furthermore, in Fig. 1d, 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. 1d 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. 1e. For Fig. 4a, 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.

Table 2 Attributes of the density of states (DOS) figures.

Local density of states

From the lowest eigenvalues together with eigenvectors computed with PRIMME, the normal state local density of states (LDOS) of Fig. 1g was computed using

$$n({{{{{{{{\boldsymbol{x}}}}}}}}}_{i},E)={N}_{{{{{{{{\rm{k}}}}}}}}}^{-1}\mathop{\sum}\limits_{\lambda ,{{{{{{{\boldsymbol{k}}}}}}}}}| {U}_{\lambda i}({{{{{{{\boldsymbol{k}}}}}}}}){| }^{2}\delta (E-{E}_{{{{{{{{\boldsymbol{k}}}}}}}}\lambda })$$
(19)

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:

$$n({{{{{{{{\boldsymbol{x}}}}}}}}}_{i},E)={N}_{{{{{{{{\rm{k}}}}}}}}}^{-1}\mathop{\sum}\limits_{\lambda ,{{{{{{{\boldsymbol{k}}}}}}}}}| {u}_{i\lambda }({{{{{{{\boldsymbol{k}}}}}}}}){| }^{2}\delta (E-{E}_{{{{{{{{\boldsymbol{k}}}}}}}}\lambda })+| {v}_{i\lambda }({{{{{{{\boldsymbol{k}}}}}}}}){| }^{2}\delta (E+{E}_{{{{{{{{\boldsymbol{k}}}}}}}}\lambda })$$
(20)

where uiλ(k) and viλ(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. 1f was obtained by first computing the particle occupation \({N}_{i}(\mu )={\sum }_{\sigma }\langle {c}_{i\sigma }^{{{{\dagger}}} }{c}_{i\sigma }\rangle\) 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(xi) = Ni(μ2) − Ni(μ1), is thus the charge density of the moiré bands and is also equivalent to the integrated LDOS: \({N}_{m}({{{{{{{\boldsymbol{x}}}}}}}})=\int\nolimits_{{\mu }_{1}}^{{\mu }_{2}}n({{{{{{{\boldsymbol{x}}}}}}}},E){{{{{{{\rm{d}}}}}}}}E\). 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. 1f.

Relative energy differences

The relative energy differences in Fig. 2e, f, were calculated by reinserting the linear combinations \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}({{\Theta }},\varphi )\) of the T = 1.25 × 10−4t eigenvectors \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x,y}\) in to the BdG Hamiltonian with same field amplitude \(\parallel \hat{{{{{{{{\mathbf{\Delta }}}}}}}}}\parallel\) as the zero temperature self-consistent solution of the corresponding coupling strength J(Tc) ≈ 0.43t. The zero temperature energy difference between the \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}({{\Theta }},\varphi )\) state and the ground state at \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}(0,0)\) were then computed as

$${{\Delta }}E({{\Theta }},\varphi )=\mathop{\sum}\limits_{n,{{{{{{{\boldsymbol{k}}}}}}}}}\left[\theta (-{E}_{{{{{{{{\boldsymbol{k}}}}}}}},{{\Theta }},\varphi }^{n}){E}_{{{{{{{{\boldsymbol{k}}}}}}}},{{\Theta }},\varphi }^{n}-\theta (-{E}_{{{{{{{{\boldsymbol{k}}}}}}}},0,0}^{n}){E}_{{{{{{{{\boldsymbol{k}}}}}}}},0,0}^{n}\right]$$
(21)

for the eigenenergy bands \({E}_{{{{{{{{\boldsymbol{k}}}}}}}},\theta ,\varphi }^{n}\) 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. 2g we calculated in the same way the relative energy differences for \(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}({{\Theta }},\varphi )\) but over a range of field norms \(\parallel \hat{{{{{{{{\mathbf{\Delta }}}}}}}}}\parallel\), i.e. not constrained to the self-consistent values at the appropriate coupling strength. Since the used basis vectors \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x,y}({T}_{{{{{{{{\rm{c}}}}}}}}})\) might change with Tc and thus field norm, this is technically an approximation. However, we find that basis vectors \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x,y}({T}_{{{{{{{{\rm{c}}}}}}}}})\) of the leading solution are remarkably stable, and that the resulting energy differences are insensitive to the \({\hat{{{{{{{{\mathbf{\Delta }}}}}}}}}}_{x,y}({T}_{{{{{{{{\rm{c}}}}}}}}})\) used for the variation, a result we checked by using both eigenvectors from T = 10−3t and T = 1.25 × 10−4t, as well as a higher (ks = 6) k-point sampling. To accurately capture the energy difference also in the very small \(\parallel \hat{{{{{{{{\mathbf{\Delta }}}}}}}}}\parallel\) regime, we computed the energy difference using a very high k-point sampling (ks = 20).

Condensation energies

The zero-temperature condensation energies in Fig. 2h–j, as well as the the maximum condensation energy of the self-consistent ϕ = π state in Fig. 4b, were computed using all the eigenvalues of the normal state and of the BdG Hamiltonian attained with Armadillo75,76. In terms of the eigenvalues En of Hamiltonian, the normal state energy is \({E}_{{{{{{{{\rm{N}}}}}}}}}={N}_{{{{{{{{\rm{k}}}}}}}}}^{-1}{\sum }_{n,{{{{{{{\boldsymbol{k}}}}}}}}}\theta (-{E}_{{{{{{{{\boldsymbol{k}}}}}}}}}^{n}){E}_{{{{{{{{\boldsymbol{k}}}}}}}}}^{n}\), while the energy in the superconducting state \(E(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}})={N}_{{{{{{{{\rm{k}}}}}}}}}^{-1}{\sum }_{n,{{{{{{{\boldsymbol{k}}}}}}}}}\theta (-{E}_{{{{{{{{\boldsymbol{k}}}}}}}}}^{n}(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}})){E}_{{{{{{{{\boldsymbol{k}}}}}}}}}^{n}(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}})-2\mu {N}_{{{{{{{{\rm{c}}}}}}}}}+{\sum }_{\langle i,j\rangle }\left|{{{\Delta }}}_{ij}{| }^{2}\right./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}_{{{{{{{{\rm{cond}}}}}}}}}={E}_{{{{{{{{\rm{N}}}}}}}}}-E(\hat{{{{{{{{\mathbf{\Delta }}}}}}}}})\).

Streamlines

The streamlines in Fig. 2b, 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 \({{{{{{{\boldsymbol{\chi }}}}}}}}({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})=\cos \tau ({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})\hat{{{{{{{{\boldsymbol{x}}}}}}}}}+\sin \tau ({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})\hat{{{{{{{{\boldsymbol{y}}}}}}}}}\) was the linearly interpolated from the angle \(\tau ({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})=\arctan ({{{{{{{\mathbf{\Delta }}}}}}}}({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})\cdot {{{{{{{{\boldsymbol{f}}}}}}}}}_{{d}_{{x}^{2}-{y}^{2}}}/{{{{{{{\mathbf{\Delta }}}}}}}}({{{{{{{{\boldsymbol{x}}}}}}}}}_{i})\cdot {{{{{{{{\boldsymbol{f}}}}}}}}}_{{d}_{xy}})\), with separate accounting for the appropriate quadrant for the angle.