Interlayer pairing in bilayer nickelates

npj | quantum materials Article Published in partnership with Nanjing University https://doi.org/10.1038/s41535-026-00849-9 Interlayer pairing in bilayer nickelates Check for updates 1234567890():,; 1234567890():,; Thomas A. Maier1 , Peter Doak1, Ling-Fang Lin2, Yang Zhang2,3, Adriana Moreo2,3 & Elbio Dagotto2,3 The discovery of Tc ~ 80 K superconductivity in pressurized La3Ni2O7 has launched a new platform to study high-temperature superconductivity. Using non-perturbative dynamic cluster approximation quantum Monte Carlo calculations, we characterize the magnetic and superconducting pairing behavior of a realistic bilayer two-orbital Hubbard-Hund model of this system that describes the relevant Ni eg states with physically relevant interaction strengths. We find a leading s± superconducting instability in this model at a temperature T ~ 100 K close to the experimentally observed Tc. Analyzing the orbital and spatial structure of the effective pairing interaction giving rise to this state reveals that the interaction predominantly acts between local interlayer pairs of the d 3z2
r 2 orbital. By correlating the strength of the interaction with that of the magnetic spin fluctuations we show that it is driven by strong interlayer spinfluctuations arising from the d 3z2
r 2 orbital. These results provide first-time non-perturbative evidence supporting the picture that a simple single-orbital bilayer Hubbard model for the Ni d 3z2
r 2 orbital provides an excellent low-energy effective description of the superconducting behavior of La3Ni2O7. Superconductivity in pressurized La3Ni2O71,2 has been widely addressed in bilayer two-orbital Hubbard models that account for the eg manifold (d x2
y2 and d 3z2
r2 orbitals) of the Ni-d states near the Fermi level in these systems3–10. Due to the complexity of the electronic structure, most of these studies have used perturbative, either weak-coupling4–9,11–14 or strongcoupling10,15–17 approaches. Depending on details in the model parameters12 and the type of the approximation, these studies have found s±, d x2
y2 -, and dxy-wave superconducting states. Optical studies, however, show evidence that La3Ni2O7 is characterized by moderately strong electronic correlations18,19 in the intermediate coupling regime. Consistent with this, recent density functional theory and constrained random-phase approximation (RPA) calculations indeed find that the Hubbard U interaction on the Ni eg orbitals is approximately of the same size as the electronic bandwidth20. This raises the question of whether perturbative weak- or strong-coupling approaches can accurately characterize the nature of pairing in these systems, and underscores the need for non-perturbative methods to address this question. While numerical methods such as density matrix renormalization group and tensor networks21–25, auxiliary-field Monte Carlo26, and cluster dynamical meanfield theory (CDMFT)27,28 have been used to provide a non-pertubative picture, they are typically based on reduced effective t-J models, simplified interactions, or quasi-one dimensional lattice geometries. Here, by using state-of-the-art dynamical cluster approximation (DCA) quantum Monte Carlo29 calculations for a realistic bilayer two-orbital model on a two-dimensional lattice6 with physically relevant interaction parameters, we examine what this model tells us about the pairing mechanism in the bilayer La3Ni2O7 compound. We find a leading s± superconducting instability in this model at a temperature T ~ 100 K close to the experimentally observed Tc ~ 80 K. By analyzing the spatial and orbital structure of the effective pairing interaction, we show that this instability arises from interlayer electron pairing on neighboring sites in the top and bottom layers, primarily in the d 3z2
r2 orbital. We then analyze the spin correlations in the model and correlate their strength with that of the s± pairing interaction to demonstrate that the pairing is primarily driven by interlayer spin fluctuations arising from the d 3z2
r2 orbital. Our results provide first-time nonperturbative evidence supporting the picture7 that a single-orbital bilayer Hubbard model30 for the d 3z2
r2 orbital provides an appropriate low-energy effective description of the pairing behavior of La3Ni2O7. Results The model we consider was introduced in ref. 6 and is illustrated in Fig. 1a. Its Hamiltonian is given in the “Methods” section. It includes the two Ni-3d orbitals, d x2
y2 and d 3z2
r2 , that have been found to account for the lowenergy electronic structure of the bilayer nickelates. These two orbitals are located on a two-dimensional (2D) square bilayer lattice, and the model includes nearest neighbor hopping, both intra- and inter-layer. It also includes intra-orbital Coulomb U, inter-orbital Coulomb U 0 , and Hund’s rule coupling J interactions, which we set to U = 3 eV, U 0 ¼ 2 eV, and J = 0.5 eV (see also “Methods” section). Pair-field susceptibility We start by discussing results for the pair-field susceptibility defined as Z Pα ðTÞ ¼ β 0 dτhT τ Δα ðτÞΔyα ð0Þi ð1Þ 1 Computational Sciences and Engineering Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA. 2Department of Physics and Astronomy, University of e-mail: Tennessee, Knoxville, TN, USA. 3Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA. npj Quantum Materials | (2026)11:19 1 Article https://doi.org/10.1038/s41535-026-00849-9 Fig. 1 | Illustration of the bilayer two-orbital model and its leading pairing states. The two orbitals, dx2
y2 (blue) and d 3z2
r2 (red), with nearest neighbor hopping parameters txx = − 0.515, tzz = − 0.110, txz = 0.243, and t ? zz ¼ 0:666 taken from ref. 6 for 25 GPa pressure are shown in panel (a). The noninteracting Fermi surface consistst of α, β, and γ sheets as illustrated in (b). The superconducting order parameter form factors listed in Table 1 for s± (b), d x2
y2 (c), and dxy (d) states are illustrated by red (positive) and blue (negative) colors. Table 1 | Singlet order parameter form factors used for the calculation of the pair-field susceptibility in Eq. (1) with the pair operator 1 X ‘‘0 Δyα ¼ pffiffiffiffi g α ðkÞcyk‘" cy
k‘0 # : N k;‘‘0 ð2Þ Here we have used the momentum space Fourier representation with wavevector k ¼ ðkx ; ky ; kz Þ with kz = 0 or π representing bonding and antibonding pffiffiffiffi combinations, respectively, of the two layers, P 0 cyk‘σ ¼ 1= N i cyi‘σ eikri , and g ‘‘ α ðkÞ is a symmetry form-factor. The singlet form factors we will use are listed in Table 1 and illustrated in Fig. 1b, c, and d. They correspond to local inter-layer s± pairs and in-plane d x2
y2 and dxy pairs. We only include intra-orbital components in these form factors, which we expect to dominate, and confirm in the next section through an unbiased analysis that the leading states indeed have substantial intra-orbital pairing character and therefore are captured by these susceptibilit (...truncated)