Twisted bilayer graphene (TBG), formed by two graphene layers rotated at a small twist angle (approximately 1.1°), exhibits a rich phase diagram featuring correlated insulating states and superconductivity. The emergence of superconductivity within superconducting domes flanking integer flat band fillings (ν = ±2) and across the moiré bands presents a significant puzzle. While the critical temperature (Tc ≈ 3 K) is relatively low, the ratio Tc/TF ≈ 0.1 surpasses the weak coupling regime, suggesting strong local electronic interactions may play a role, alongside electron-phonon coupling. The large moiré pattern (∼104 carbon atoms) and complex band topology have greatly hindered experimental and theoretical investigations into the nature of superconductivity in TBG. Continuum models successfully replicate the normal-state band structure but struggle to capture the strong local interactions likely involved in the superconducting mechanism. Previous studies, many relying on effective lattice models, have suggested a time-reversal symmetry (TRS)-breaking chiral d-wave state, mirroring observations in heavily doped graphene. However, recent experimental and theoretical work has hinted at the presence of nematic superconductivity in TBG, raising questions about its prevalence and observable characteristics. This study employs full-scale atomistic modeling to accurately address the superconducting behavior in TBG, considering all carbon atoms, comprehensive k-point sampling, and realistic electron interactions consistent with experimental findings in TBG and graphene-like systems, as well as known interactions in cuprates. By solving the mean-field self-consistent and linear gap equations at experimentally relevant temperatures, the study aims to comprehensively characterize the superconducting state.

The study of superconductivity in TBG has been actively pursued, with many studies focusing on the interplay of correlated insulating states and superconducting domes. The large scale of the moiré pattern and complex band structure posed challenges for theoretical modeling. Continuum models successfully capture the normal state band structure, but struggle to capture the strong correlations. Effective lattice models, while computationally tractable, must balance accuracy with orbital proliferation. Many proposed a time-reversal symmetry (TRS) breaking chiral d-wave state. This state is considered a natural outcome given the transformation of the cuprate d-wave state into dxy and dx²-y² wave states with a symmetry-enforced degeneracy at Γ on lattices with three- or six-fold symmetry, such as graphene and TBG. In the honeycomb lattice the chiral d + idxy wave becomes fully gapped, which energetically favors it over any real-valued nematic d-wave state which always has a nodal spectrum. However, very recent research has proposed alternative nematic superconducting states for TBG. Some of these proposals rely on higher-order coupling to additional coexisting orders or pairing channels to energetically favor the nematic state over the chiral state. Others, utilizing rescaled or atomistic lattice models, directly observe a nematic phase, but with seemingly varying properties. The recent experimental observation of intrinsic nematicity in TBG's superconducting state through magnetotransport measurements strengthens the focus on the role of nematic superconductivity.

This work uses full-scale atomistic tight-binding modeling to accurately capture superconductivity in TBG. The model includes all carbon atoms in a large periodic moiré unit cell, employs a dense k-point sampling of the moiré Brillouin zone (MBZ), and incorporates electron interactions consistent with experimental findings in TBG and graphene-like systems. The out-of-plane pz carbon orbitals hybridize in-plane through nearest-neighbor hopping (t = 2.7 eV), while interlayer hybridization is modeled using exponentially decaying hopping with Koster-Slater angular dependence. Superconducting pairing is modeled using spin-singlet order parameters on in-plane nearest-neighbor carbon bonds with uniform coupling strength J. The superconducting order parameter field Δ is determined by solving both the full non-linear, self-consistent, and linearized gap equations. The linear gap equation, valid near Tc, helps determine the critical temperature and symmetry of the superconducting state. The non-linear gap equation, solved iteratively at zero temperature, accounts for non-linear contributions and allows for the exploration of symmetry breaking. The rational pole expansion method is used to efficiently solve these gap equations due to the large system size of TBG. This method calculates the Fermi operator using minimax rational approximation, enabling fast computation of expectation values and static responses necessary for solving the gap equations. To analyze the superconducting state, the study parametrizes the order parameter field as a function of nematic angle (θ) and complex phase (φ), enabling exploration of different superconducting states (including chiral and real-valued nematic states). The relative free energy of these states is calculated to determine the ground state at zero temperature. The atomic-scale nematicity is characterized by introducing a three-dimensional local order vector and projecting it onto the irreducible representations of the D6h point group of the graphene honeycomb lattice. Finally, the study examines the impact of interlayer Josephson coupling on the superconducting gap structure by tuning the interlayer phase.

The atomistic model reveals a highly inhomogeneous, real-valued nematic d-wave superconducting state that dominates the phase diagram at experimentally observed temperatures. The nematicity is present on both atomic and moiré length scales. The order parameter field exhibits a significant amplitude variation throughout the moiré unit cell, with enhanced amplitudes in the AA regions. The nematic order breaks the C3 rotational symmetry of the normal state. The critical temperature (Tc) exhibits a non-linear dependence on the coupling strength J/t, indicating contributions from dispersive bands and a catalytic effect from the low-energy moiré flat bands. At zero temperature, the self-consistent solutions are always real-valued linear combinations of the two leading degenerate solutions (Δx and Δy) of the linear gap equation, confirming the nematic nature of the superconductivity. The ground state is threefold degenerate, with the nematic C2 axis oriented toward one of the next-nearest-neighbor AA regions. The energy splitting among the real-valued nematic states is sixfold symmetric as a function of the nematic angle, with large values (up to 10% of kBTc) around the experimentally observed Tc. The TRS-breaking chiral state is found to be energetically unfavorable at experimentally relevant temperatures. The real-valued nematic ground state possesses a full energy gap, a feature tied to strong π-phase Josephson locking between graphene layers. This full gap is unexpected for a d-wave state and is in sharp contrast to the nodal spectra traditionally associated with such states. The atomic-scale nematicity shows a strong spatial variation, exhibiting two vortex-antivortex pairs outside the central AA region. Nematicity is directly detectable in the sublattice-resolved LDOS, showing significant asymmetries around the moiré cell center in the superconducting state but not in the normal state. This provides a measurable signal for nematic superconductivity in TBG. The π-phase interlayer Josephson coupling is crucial for the existence of the full superconducting gap. The strong interlayer coupling and the resulting π-shift prevent any nodal behavior in the superconducting spectrum, stabilizing the nematic state.

The findings of this study address the open question of the nature of superconductivity in magic-angle TBG. The observation of a nematic d-wave state with a full energy gap, driven by a π-phase interlayer Josephson coupling, significantly advances our understanding of this system. The importance of atomistic modeling is clearly demonstrated, as simpler models fail to capture the inhomogeneous nature of the superconducting state and its intricate features. The threefold degenerate ground state and its dependence on the C2 nematic axis orientation highlight the unique role of the moiré pattern in shaping the superconducting properties. The full energy gap, despite the d-wave symmetry, is an unexpected result and points to the necessity of considering multiorbital or multiband effects. The predicted sublattice-resolved LDOS anisotropy provides a direct experimental signature of the nematic superconductivity, offering a path for verifying the theoretical predictions. The comparison with high-temperature cuprate superconductors suggests a potential universality in the role of local electronic interactions, despite the very different resulting superconducting states.

This research reveals a highly inhomogeneous nematic d-wave superconducting state in magic-angle TBG at experimentally relevant temperatures. The work highlights the critical role of atomistic modeling in capturing the complex superconducting state and provides strong evidence for the existence of a full-gap nematic superconductor in TBG. The predicted sublattice-resolved LDOS anisotropy offers a crucial experimental test for future investigations. Future research should explore the precise nature of the electronic interactions in greater detail and investigate the effects of disorder and other external perturbations on the nematic superconducting state. The findings could stimulate research into similar phenomena in other moiré systems.

The model relies on the assumption of intralayer spin-singlet bond interactions. While this choice is supported by experimental evidence and theoretical work, variations in the specific interactions could potentially influence the detailed features of the superconducting state. The computational cost of full-scale atomistic modeling limits the size of the simulated systems and the range of explored parameters. Furthermore, the impact of structural imperfections and disorder, which are always present in real TBG samples, is not explicitly included in this model. The absence of detailed lattice relaxation effects could affect the quantitative aspects of the results, although it is likely that the qualitative conclusions remain unchanged.