Abstract
Twisted bilayer graphene (TBG) develops large moiré 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é length scales. The nematic state has a locally anisotropic real-valued d-wave pairing, with a nematic vector winding throughout the moiré pattern, and is three-fold degenerate. Although d-wave symmetric, the superconducting state has a full energy gap, which we tie to a π-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.
Similar content being viewed by others
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.
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.
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.
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.
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 S ∈ A, 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:
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:
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:
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 \({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
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β(HBdG(k)):
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
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
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,
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
where \({{{{{{{{\boldsymbol{x}}}}}}}}}_{i}^{n}\) and \({{{{{{{{\boldsymbol{y}}}}}}}}}_{j}^{n}\) are Nc dimensional vectors satisfying the linear equations:
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.
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.
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
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:
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
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.
Data availability
Data are available from the corresponding authors upon request.
Code availability
Computer program source codes used for all calculations and analysis are available for review from the corresponding authors upon request.
References
Cao, Y. et al. Correlated insulator behaviour at half-filling in magic-angle graphene superlattices. Nature 556, 80–84 (2018).
Cao, Y. et al. Unconventional superconductivity in magic-angle graphene superlattices. Nature 556, 43–50 (2018).
Li, G. et al. Observation of Van Hove singularities in twisted graphene layers. Nat. Phys. 6, 109–113 (2009).
Po, H. C., Zou, L., Vishwanath, A. & Senthil, T. Origin of mott insulating behavior and superconductivity in twisted bilayer graphene. Phys. Rev. X 8, 031089 (2018).
Lu, X. et al. Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene. Nature 574, 653–657 (2019).
Kerelsky, A. et al. Maximized electron interactions at the magic angle in twisted bilayer graphene. Nature 572, 95–100 (2019).
Jiang, Y. et al. Charge order and broken rotational symmetry in magic-angle twisted bilayer graphene. Nature 573, 91–95 (2019).
Xie, Y. et al. Spectroscopic signatures of many-body correlations in magic-angle twisted bilayer graphene. Nature 572, 101–105 (2019).
Choi, Y. et al. Electronic correlations in twisted bilayer graphene near the magic angle. Nat. Phys. 15, 1174–1180 (2019).
Sharpe, A. L. et al. Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene. Science 365, 605–608 (2019).
Balents, L., Dean, C. R., Efetov, D. K. & Young, A. F. Superconductivity and strong correlations in moiré flat bands. Nat. Phys. 16, 725–733 (2020).
Andrei, E. Y. & MacDonald, A. H. Graphene bilayers with a twist. Nat. Mater. 19, 1265–1275 (2020).
Yankowitz, M. et al. Tuning superconductivity in twisted bilayer graphene. Science 363, 1059–1064 (2019).
Polshyn, H. et al. Large linear-in-temperature resistivity in twisted bilayer graphene. Nat. Phys. 15, 1011–1016 (2019).
Cao, Y. et al. Strange metal in magic-angle graphene with near planckian dissipation. Phys. Rev. Lett. 124, 076801 (2020).
Keimer, B., Kivelson, S. A., Norman, M. R., Uchida, S. & Zaanen, J. From quantum matter to high-temperature superconductivity in copper oxides. Nature 518, 179–186 (2015).
Chubukov, A. Pairing mechanism in Fe-based superconductors. Annu. Rev. Condens. Matter Phys. 3, 57–92 (2012).
Peltonen, T. J., Ojajärvi, R. & Heikkilä, T. T. Mean-field theory for superconductivity in twisted bilayer graphene. Phys. Rev. B 98, 220504 (2018).
Wu, F., MacDonald, A. & Martin, I. Theory of phonon-mediated superconductivity in twisted bilayer graphene. Phys. Rev. Lett. 121, 257001 (2018).
Lian, B., Wang, Z. & Bernevig, B. A. Twisted bilayer graphene: a phonon-driven superconductor. Phys. Rev. Lett. 122, 257002 (2019).
Bistritzer, R. & MacDonald, A. H. Moiré bands in twisted double-layer graphene. Proc. Natl Acad. Sci. USA 108, 12233–12237 (2011).
Moon, P. & Koshino, M. Optical absorption in twisted bilayer graphene. Phys. Rev. B 87, 205404 (2013).
Su, Y. & Lin, S.-Z. Pairing symmetry and spontaneous vortex-antivortex lattice in superconducting twisted-bilayer graphene: Bogoliubov-de Gennes approach. Phys. Rev. B 98, 195101 (2018).
Fischer, A., Klebl, L., Honerkamp, C. & Kennes, D. M. Spin-fluctuation-induced pairing in twisted bilayer graphene. Phys. Rev. B 103, L041103 (2021).
Fidrysiak, M., Zegrodnik, M. & Spałek, J. Unconventional topological superconductivity and phase diagram for an effective two-orbital model as applied to twisted bilayer graphene. Phys. Rev. B 98, 085436 (2018).
Xu, C. & Balents, L. Topological superconductivity in twisted multilayer graphene. Phys. Rev. Lett. 121, 087001 (2018).
Wu, F. Topological chiral superconductivity with spontaneous vortices and supercurrent in twisted bilayer graphene. Phys. Rev. B 99, 195114 (2019).
Kennes, D. M., Lischner, J. & Karrasch, C. Strong correlations and d + id superconductivity in twisted bilayer graphene. Phys. Rev. B 98, 241407 (2018).
Liu, C.-C., Zhang, L.-D., Chen, W.-Q. & Yang, F. Chiral spin density wave and d + id superconductivity in the magic-angle-twisted bilayer graphene. Phys. Rev. Lett. 121, 217001 (2018).
Guo, H., Zhu, X., Feng, S. & Scalettar, R. T. Pairing symmetry of interacting fermions on a twisted bilayer graphene superlattice. Phys. Rev. B 97, 235453 (2018).
Black-Schaffer, A. M. & Doniach, S. Resonating valence bonds and mean-field d-wave superconductivity in graphite. Phys. Rev. B 75, 134512 (2007).
Nandkishore, R., Levitov, L. S. & Chubukov, A. V. Chiral superconductivity from repulsive interactions in doped graphene. Nat. Phys. 8, 158–163 (2012).
Black-Schaffer, A. M. & Honerkamp, C. Chiral d-wave superconductivity in doped graphene. J. Phys. Condens. Matter 26, 423201 (2014).
Venderbos, J. W. F. & Fernandes, R. M. Correlations and electronic order in a two-orbital honeycomb lattice model for twisted bilayer graphene. Phys. Rev. B 98, 245103 (2018).
Chichinadze, D. V., Classen, L. & Chubukov, A. V. Nematic superconductivity in twisted bilayer graphene. Phys. Rev. B 101, 224513 (2020).
Kozii, V., Isobe, H., Venderbos, J. W. F. & Fu, L. Nematic superconductivity stabilized by density wave fluctuations: Possible application to twisted bilayer graphene. Phys. Rev. B 99, 144507 (2019).
Dodaro, J. F., Kivelson, S. A., Schattner, Y., Sun, X. Q. & Wang, C. Phases of a phenomenological model of twisted bilayer graphene. Phys. Rev. B 98, 075154 (2018).
Wang, Y., Kang, J. & Fernandes, R. M. Topological and nematic superconductivity mediated by ferro-SU(4) fluctuations in twisted bilayer graphene. Phys. Rev. B 103, 024506 (2021).
Julku, A., Peltonen, T. J., Liang, L., Heikkilä, T. T. & Törmä, P. Superfluid weight and berezinskii-kosterlitz-thouless transition temperature of twisted bilayer graphene. Phys. Rev. B 101, 060505 (2020).
Cao, Y. et al. Nematicity and competing orders in superconducting magic-angle graphene. Science 372, 264–271 (2021).
Morell, E. S., Correa, J. D., Vargas, P., Pacheco, M. & Barticevic, Z. Flat bands in slightly twisted bilayer graphene: Tight-binding calculations. Phys. Rev. B 82, 121407 (2010).
Shallcross, S., Sharma, S., Kandelaki, E. & Pankratov, O. A. Electronic structure of turbostratic graphene. Phys. Rev. B 81, 165105 (2010).
Wu, F. & Sarma, S. D. Identification of superconducting pairing symmetry in twisted bilayer graphene using in-plane magnetic field and strain. Phys. Rev. B 99, 220507 (2019).
Wehling, T. O. et al. Strength of effective coulomb interactions in graphene and graphite. Phys. Rev. Lett. 106, 236805 (2011).
Lee, P. A., Nagaosa, N. & Wen, X.-G. Doping a Mott insulator: physics of high-temperature superconductivity. Rev. Mod. Phys. 78, 17–85 (2006).
Pauling, L. The Nature of the Chemical Bond and the Structure of Molecules and Crystals. (Cornell University Press, 1960).
Anderson, P. W. The resonating valence bond state in La2CuO4 and superconductivity. Science 235, 1196–1198 (1987).
Uchoa, B. & Neto, A. H. C. Superconducting states of pure and doped graphene. Phys. Rev. Lett. 98, 146801 (2007).
Kopnin, N. B. & Sonin, E. B. BCS superconductivity of dirac electrons in graphene layers. Phys. Rev. Lett. 100, 246808 (2008).
Heikkilä, T. T., Kopnin, N. B. & Volovik, G. E. Flat bands in topological media. JETP Lett. 94, 233–239 (2011).
Kopnin, N. B., Heikkilä, T. T. & Volovik, G. E. High-temperature surface superconductivity in topological flat-band systems. Phys. Rev. B 83, 220503 (2011).
Volovik, G. E. Flat band in topological matter. J. Supercond. Nov. Magn. 26, 2887–2890 (2013).
Löthman, T. & Black-Schaffer, A. M. Universal phase diagrams with superconducting domes for electronic flat bands. Phys. Rev. B 96, 064505 (2017).
Fu, L. Odd-parity topological superconductor with nematic order: application to cuxbi2se3. Phys. Rev. B 90, 100509 (2014).
Kiesel, M. L., Platt, C., Hanke, W., Abanin, D. A. & Thomale, R. Competing many-body instabilities and unconventional superconductivity in graphene. Phys. Rev. B 86, 020507 (2012).
Chubukov, A. V., Vafek, O. & Fernandes, R. M. Displacement and annihilation of Dirac gap nodes in d-wave iron-based superconductors. Phys. Rev. B 94, 174518 (2016).
Agterberg, D., Shishidou, T., O’Halloran, J., Brydon, P. & Weinert, M. Resilient nodeless d-wave superconductivity in monolayer FeSe. Phys. Rev. Lett. 119, 267001 (2017).
Pang, G. et al. Fully gapped d-wave superconductivity in CeCu2Si2. Proc. Natl Acad. Sci. USA 115, 5343–5347 (2018).
Nakayama, T., Shishidou, T. & Agterberg, D. F. Nodal topology in d-wave superconducting monolayer FeSe. Phys. Rev. B 98, 214503 (2018).
Nica, E. M. & Si, Q. Multiorbital singlet pairing and d + d superconductivity. npj Quant. Mater. 6, 3 (2021).
Zhang, L. Lowest-energy moiré band formed by dirac zero modes in twisted bilayer graphene. Sci. Bull. 64, 495–498 (2019).
Peotta, S. & Törmä, P. Superfluidity in topologically nontrivial flat bands. Nat. Commun. 6, 8944 (2015).
Hu, X., Hyart, T., Pikulin, D. I. & Rossi, E. Geometric and conventional contribution to the superfluid weight in twisted bilayer graphene. Phys. Rev. Lett. 123, 237002 (2019).
Xie, F., Song, Z., Lian, B. & Bernevig, B. A. Topology-bounded superfluid weight in twisted bilayer graphene. Phys. Rev. Lett. 124, 167002 (2020).
van Wijk, M. M., Schuring, A., Katsnelson, M. I. & Fasolino, A. Relaxation of moiré patterns for slightly misaligned identical lattices: graphene on graphite. 2D Materials 2, 034010 (2015).
Nam, N. N. T. & Koshino, M. Lattice relaxation and energy band modulation in twisted bilayer graphene. Phys. Rev. B 96, 075311 (2017).
Lucignano, P., Alfè, D., Cataudella, V., Ninno, D. & Cantele, G. Crucial role of atomic corrugation on the flat bands and energy gaps of twisted bilayer graphene at the magic angle θ ~ 1.08∘. Phys. Rev. B 99, 195419 (2019).
Zou, L., Po, H. C., Vishwanath, A. & Senthil, T. Band structure of twisted bilayer graphene: Emergent symmetries, commensurate approximants, and wannier obstructions. Phys. Rev. B 98, 085435 (2018).
Goedecker, S. Integral representation of the Fermi distribution and its applications in electronic-structure calculations. Phys. Rev. B 48, 17573–17575 (1993).
Goedecker, S. Linear scaling electronic structure methods. Rev. Mod. Phys. 71, 1085–1123 (1999).
Moussa, J. E. Minimax rational approximation of the Fermi-Dirac distribution. J. Chem. Phys. 145, 164108 (2016).
Paige, C. C. & Saunders, M. A. Solution of sparse indefinite systems of linear equations. SIAM J. Numer. Anal. 12, 617–629 (1975).
Cullum, J. K. & Willoughby, R. A. Lanczos Algorithms for Large Symmetric Eigenvalue Computations (Society for Industrial and Applied Mathematics, 2002).
Stathopoulos, A. & McCombs, J. R. PRIMME: PReconditioned Iterative MultiMethod Eigensolver: methods and software description. ACM Trans. Math. Softw. 37, 21:1–21:30 (2010).
Sanderson, C. & Curtin, R. Armadillo: a template-based c + + library for linear algebra. J. Open Source Softw. 1, 26 (2016).
Sanderson, C. & Curtin, R. A user-friendly hybrid sparse matrix class in c + + . In Mathematical Software – ICMS 2018, 422-430 (Springer International Publishing, 2018).
Acknowledgements
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.
Funding
Open access funding provided by Uppsala University.
Author information
Authors and Affiliations
Contributions
A.M.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.M.B.S. and T.L. wrote the manuscript with input from all authors.
Corresponding authors
Ethics declarations
Competing interests
The authors declare no competing interests.
Peer review
Peer review information
Communications Physics thanks the anonymous reviewers for their contribution to the peer review of this work.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Löthman, T., Schmidt, J., Parhizgar, F. et al. Nematic superconductivity in magic-angle twisted bilayer graphene from atomistic modeling. Commun Phys 5, 92 (2022). https://doi.org/10.1038/s42005-022-00860-z
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s42005-022-00860-z
This article is cited by
-
Unconventional superconductivity in magic-angle twisted trilayer graphene
npj Quantum Materials (2022)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.