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

The study investigates whether and how nematic superconductivity emerges in magic-angle twisted bilayer graphene (TBG), a system with flat moiré bands and strong electronic correlations. Around the magic twist angle (~1.1°), TBG exhibits large density of states and correlated phases, including insulating states and superconductivity. While many theoretical works on reduced or continuum models predict a time-reversal symmetry breaking chiral d-wave state, recent proposals and experiments indicate possible nematic superconductivity. The authors aim to resolve the superconducting pairing symmetry and its spatial structure using full-scale atomistic modeling that can incorporate strong local interactions and the full moiré-scale inhomogeneity. The purpose is to clarify the ground state near experimentally relevant temperatures and fillings, its degeneracy and gap structure, and provide experimentally accessible signatures, thereby addressing discrepancies between previous theoretical predictions and experimental observations of nematicity.

Prior works using continuum and effective lattice models reproduced normal-state features but often favored chiral d + id superconductivity on hexagonal lattices due to symmetry-enforced degeneracy of d-wave components (e.g., in doped graphene). Several studies proposed chiral superconductivity in TBG, while others suggested nematic states stabilized via coupling to additional orders or nearly degenerate pairing channels. More recently, both rescaled and atomistic lattice models reported nematic phases in certain parameter regimes, albeit with differing properties. Experimentally, magnetotransport in magic-angle TBG has revealed intrinsic superconducting nematicity, and a pseudogap and strange metal behaviors have been observed above the superconducting domes. These findings motivate a careful, atomistically resolved analysis of superconductivity in TBG with local electronic interactions.

- System: Full-scale atomistic tight-binding model of commensurate-angle TBG including all carbon atoms in the moiré unit cell. Intralayer nearest-neighbor hopping t = 2.7 eV between p_z orbitals; interlayer hybridization modeled by exponentially decaying hopping with Koster–Slater angular dependence. Dense k-point sampling of the moiré Brillouin zone. - Superconducting pairing: Intralayer spin-singlet nearest-neighbor bond pairing with uniform coupling strength J, motivated by strong on-site repulsion and t–J physics relevant to graphene-like systems. Bond order parameters defined on every in-plane nearest-neighbor carbon bond forming an order parameter field Δ over the moiré cell. - Equations solved: Both (i) linearized gap equation near Tc to obtain critical temperatures and symmetry classification (irreps of D3: A1, A2, E), and (ii) full non-linear, zero-temperature self-consistent Bogoliubov–de Gennes (BdG) gap equation to determine the ground state and energy landscape beyond Tc. The order parameter is initialized generally (including complex configurations) and iterated to convergence. - Computational strategy: Use a rational pole (minimax) expansion of the Fermi operator to compute single-particle expectation values and static responses efficiently, avoiding full diagonalization for large systems. Solve resulting linear systems (Lanczos/MINRES) exploiting sparse BdG structure; parallelize over poles and k-points. Validate accuracy against direct diagonalization for larger angles (~≥1.5°). Kernel density estimation used for DOS/LDOS with specified Gaussian widths and k-grids. - Symmetry analysis: Global (moiré-scale) classification via D3 projection operators; local (atomic-scale) decomposition by projecting three-bond vectors onto D6h (or D3h) form factors: extended s-wave and two d-wave components (d_xy, d_x2−y2). Track local d-wave orientation vector r(x) to visualize nematic textures. - Parameter regimes: Magic-angle regime (~1.1–1.2°), fillings near the upper and lower van Hove singularities (|ν|≈1), tuning J/t (e.g., ~0.42–1.0) to match experimental Tc ~ 3–4 K, and scanning chemical potential across the moiré bands. Analyze interlayer phase locking by imposing a relative phase φ between layers and computing condensation energy vs φ.

- Leading instability and Tc: Across a wide range of J/t, the highest Tc solution belongs to the two-dimensional E irrep, spanned by two real order parameter fields Δx and Δy. Superconductivity is significantly enhanced in TBG versus pristine graphene (graphene QCP Jc/t ≈ 1.76), with finite Tc even at weak coupling; for weak J, Tc scales approximately linearly with J/t due to flat bands. Tuning J ≈ 0.42–0.43 t yields Tc ≈ 4 K, consistent with experiment. - Nematic ground state: Zero-temperature self-consistent solutions are purely real and lie entirely in the E-manifold as linear combinations of Δx and Δy, indicating spontaneous rotational symmetry breaking (nematicity) on the moiré scale. The order amplitude is enhanced in AA regions along the nematic C2 axis. The superconducting amplitude vs filling shows maxima near the upper and lower van Hove singularities and a minimum near half-filling, reflecting the normal-state DOS. - Degeneracy and energy landscape: Parameterizing Δ(θ, φ) = ||Δ||(cosθ Δx + e^{iφ} sinθ Δy), the free-energy landscape at T=0 exhibits minima for real combinations (φ=0), forming three gauge-inequivalent nematic ground states (three-fold degeneracy) with distinct C2-axis orientations. TRS-breaking chiral combinations Δx ± iΔy are maxima (least favorable). The energy splitting among nematic orientations has a six-fold structure in θ, with magnitude up to ~10% of kB Tc near experimental Tc; the chiral state is substantially higher in energy in this regime and only approaches competitiveness at vanishingly small order amplitudes. - Gap structure: The most favorable nematic state is fully gapped, producing large condensation energy with sharp coherence peaks; less favorable nematic and chiral states show reduced or partially filled gaps with closer coherence peaks and reduced condensation energies. This nodeless d-wave behavior contrasts with doped graphene, where nematic d-wave is nodal and chiral d + id is fully gapped. - Atomic-scale nematic texture: The local d-wave composition is almost purely d-wave (<0.025% s-wave). The d-wave orientation field r(x) aligns with the moiré-scale nematic axis in the AA region and forms two robust vortex–antivortex pairs outside AA; antivortices are pinned to AB and BA centers, while vortices lie near AA on opposite sides of the nematic axis. This pattern is time-reversal invariant (no supercurrents) and robust across fillings and J/t. - LDOS signatures: Sublattice-resolved LDOS along cuts at increasing angles to the nematic axis shows pronounced asymmetry about the moiré center, mainly via shifts of coherence peaks. A–B sublattices exhibit opposite shifts (LDOS polarization), causing anisotropy to cancel in combined LDOS. The effect can reach nearly half of the maximal moiré LDOS, providing a clear prediction for sublattice-resolved STS/STM. - Interlayer π-phase Josephson locking: The superconducting phases in the two layers are rigidly π-shifted in both linear and self-consistent solutions, reversing the local nematic vectors between layers. Condensation energy vs interlayer phase φ shows that superconductivity is stable only in a narrow phase window centered at φ=π; for much of φ-space the superconducting state is unstable (higher energy than normal). The π-shift generates the full superconducting gap; as φ→0 the gap closes. This π-locking is an intrinsic consequence of the moiré band structure. - Comparison and requirements: Accurate description requires atomistic spatial resolution and dense k-point sampling over the MBZ; sampling only the MBZ center yields an artificial chiral ground state at experimental temperatures. The results indicate that similar local interactions can yield markedly different superconducting states in TBG versus cuprates.

The results establish that magic-angle TBG hosts a nematic superconducting state at and above experimentally relevant temperatures and fillings. The spontaneous breaking of C3 rotational symmetry arises intrinsically from the superconducting E-channel without invoking additional orders, resolving tensions between chiral predictions from simplified models and recent experimental reports of nematicity. The unexpected fully gapped spectrum of a real d-wave state is attributed to a strong, moiré-induced interlayer π-phase Josephson coupling, highlighting the role of multiband/moiré physics absent in monolayer graphene and simple rescaled models. The three-fold degenerate nematic orientations produce distinct condensation energies and gap structures, and generate robust atomic-scale vortex–antivortex textures in the pairing field. Clear, testable signatures are predicted in sublattice-resolved LDOS anisotropy, implying that advanced STS/STM with sublattice resolution can directly detect superconducting nematicity. The strong dependence on k-point sampling and the necessity of atomistic modeling underscore that coarse-grained or MBZ-center-only approaches can qualitatively misidentify the ground state (e.g., erroneously favoring chiral order). Overall, the findings demonstrate that local, spin-singlet bond interactions are sufficient to produce nematic superconductivity in TBG, distinct from cuprates where chiral d + id is typically favored on the honeycomb lattice.

Full-scale atomistic modeling of magic-angle TBG reveals a highly inhomogeneous, real-valued d-wave nematic superconducting state with a full energy gap near the experimentally observed Tc. The leading instability is the twofold E irrep with three-fold degenerate nematic ground states. A strong interlayer π-phase locking, intrinsic to the moiré band structure, is responsible for the full gap and stabilizes nematic order over the chiral state at realistic temperatures and couplings. Nematic superconductivity produces pronounced, experimentally accessible sublattice-resolved LDOS anisotropies and robust atomic-scale vortex–antivortex textures in the pairing field. Accurate resolution of the superconducting state requires atomistic spatial detail and dense MBZ sampling; oversimplified sampling can yield incorrect chiral ground states. Future research directions include: refining the microscopic interaction model beyond uniform intralayer bond coupling (e.g., incorporating longer-range, interlayer, or phonon-mediated interactions); exploring strain/disorder effects and their control of nematic orientation; studying thermal/quantum fluctuations between the three nematic domains; and extending to other twisted multilayer systems or varying twist angles to map robustness of π-locking and gap nodelessness.

The main uncertainty is the precise form of electronic interactions in TBG; the model assumes uniform intralayer spin-singlet nearest-neighbor bond pairing with a single coupling J and neglects explicit electron-phonon and interlayer pairing channels. Results depend on detailed atomistic modeling and dense k-point sampling; restricted sampling (e.g., MBZ center only) can qualitatively change the inferred ground state. While the study examines relevant twist angles and fillings near VHS, material-specific variations (relaxation, strain, disorder, substrate effects) and exact parameter values for interlayer hopping may affect quantitative predictions. Finite-temperature fluctuation effects beyond mean-field and domain dynamics among the three nematic orientations are not included.