The fate of $s_{\pm}$ superconductivity in the bilayer nickelate La$_3$Ni$_2$O$_7$

Introduction

The discovery of high-temperature superconductivity in the cuprates is a landmark achievement of physics [2, 3, 4] which has since spurred the community to consider other related oxides, such as the ruthenates [5, 6] and more recently, the nickelates. Anisimov and Rice Anisimov and Rice [7] originally reasoned that nickel, with one less electron than copper, might, in a strongly reducing environment, be pushed towards the $d^8/d^9$ configuration of copper yielding a new family of high temperature superconductors. The chemistry required to fulfill their proposal required three decades of further development in solid-state synthesis, culminating in the discovery of superconductivity in the infinite layer nickelates [8]. Further advances in synthesis and epitaxy over the past decade now leave us at a point where nickelate transition temperatures are comparable with the cuprates, and the discovery of ambient pressure materials opens up the potential of experiments at the same level of detail as in the cuprates [9]. This makes it an ideal time to revisit the theory.

A variety of theories have been put forth to explain how strong correlations between electrons drive pairing [3, 10]. A central question that each of these theories must answer is how electrons can accommodate the effects of strong Coulomb repulsion? In Bardeen Cooper Schrieffer (BCS) theory, the attractive interaction mediated by phonon exchange avoids the Coulomb interaction through the effect of retardation [11, 12]: in simple terms, the deformation of the lattice produced by an electron remains long after it has passed. However, in strongly interacting systems where the bare interaction between the electrons exceeds their kinetic energy, the pairing fermions avoid the Coulomb interaction by developing nodes in the pair wavefunction: nodes cost kinetic energy, but they pay off by keeping the electrons apart and minimizing the Coulomb cost.

In a weak-coupling treatment the development of nodes in the pair wavefunction can be seen from the zero-temperature gap equation,

where $\Lambda$ is a low energy cut-off. If we rewrite the gap equation in the form

then it is clear that if the interaction is predominantly repulsive ($V_{k,k'}>0$), then there is no solution for the gap function $\Delta_k$ that is positive everywhere on the Fermi surface. Moreover, to maximize the gap it is advantageous that $\operatorname{sgn}(V_{k,k'})\operatorname{sgn}(\Delta_{k'})<0$, i.e that the gap $\Delta_k$ must become negative in regions where $V_{k,k'}>0$ is repulsive, giving rise to an anisotropic gap with nodes. Famously, the gap $\Delta_k$ in superfluid Helium-3 has $p$-wave symmetry [13], while that of cuprate superconductors has $d$-wave symmetry (Figure 2). The presence of $d$-wave pairing in the cuprate superconductors has been confirmed by multiple experimental methods, such as corner-sharing SQUID experiments [14], flux quantization experiments [15], angle-resolved photoemission spectroscopy (ARPES) [16] and scanning tunnelling microscopy (STM) [17, 18, 19].

While weak-coupling theory provides a qualitative understanding of the development of gap anisotropy and the development of $d$-wave symmetry under the influence of antiferromagnetic interactions [20, 21, 22], in the cuprates the bare strength of interactions is comparable with their band-width, taking them out of the realm of weak-coupling. One approach is to try to renormalize the physics from weak to strong coupling [23]. The alternative attempts to incorporate strong coupling effects at the outset. A well-known strong-coupling model for cuprate superconductors is the $t$-$J$ model[24]

which describes electrons hopping with amplitude $-t$ between neighboring sites on a square lattice, with an additional antiferromagnetic coupling $J$ between two neighboring electrons. The antiferromagnetic interaction arises from integrating out the states of double occupancy from the Hubbard model, giving rise to a super-exchange coupling between sites of magnitude $J\sim t^2/U$, where $U$ is the onsite interaction energy of the Hubbard model. The $t$-$J$ model is defined on the subspace that excludes double occupancy on each site, spanned by the Hubbard operators $X_{\alpha\beta}=|\alpha\rangle\langle\beta|$, where $|\alpha\rangle,|\beta\rangle\in\{|0\rangle,|\uparrow\rangle,|\downarrow\rangle\}$. These non-canonical operators are neither boson nor fermion, but instead obey the Hubbard superalgebra $([X_{\alpha\beta},X_{\gamma\delta}]_{\pm}=\delta_{\beta\gamma}X_{\alpha\delta}\pm\delta_{\alpha\delta}X_{\gamma\beta})$[25].

Contrasting the nickelates and cuprates

In both the nickelates and the cuprates, the $d$-electrons are subject to a cubic crystal field from the neighboring oxygen atoms, leading to the $e_g$ shell $(3d_{x^2-y^2},3d_{z^2})$ being higher in energy. The physics of the both classes of material is dominated by the $e_g$ states, while the $t_{2g}$ orbitals are fully filled even at moderate hole doping and do not participate in the physics. Until only a few years ago, the superconducting $T_c$s in the nickelates were limited below 30 K. Amongst these early materials Nd$_{0.82}$Sr$_{0.18}$NiO$_2$ thin films [26] achieved the highest $T_c$s. In these materials, the Ni hosts an oxidation of +1, endowing it with an outer shell configuration identical to the cuprates ($3d^9$). In 2023, the Ruddlesden-Popper nickelate La$_3$Ni$_2$O$_7$ was synthesized with a $T_c$ of 80 K [1]. This system is markedly different from the cuprates - the Ni cation has an oxidation of +2.5, corresponding to an outer shell configuration of $3d^{7.5}$. Unlike the cuprates, which can be modeled with a single Fermi surface, La$_3$Ni$_2$O$_7$ hosts multiple Fermi surfaces, leading many groups to propose the possibility of $s_{\pm}$ pairing - an isotropic superconducting gap that switches sign on different Fermi pockets, as illustrated in Figure 3. Similar proposals have been explored for the iron pnictides.

Theories of pairing in La$_3$Ni$_2$O$_7$ range from weak-coupling treatments of spin- and charge-fluctuation-mediated pairing [27, 28] to strong-coupling, RVB-based descriptions [29, 30, 31, 32, 33, 34, 35]. Recent work has sharpened the competition between $s_{\pm}$ and $d$-wave pairing. A microscopic two-orbital $t$-$J$ derivation identifies exchange processes beyond the simplest single-orbital model [36], while cluster quantum Monte Carlo calculations find a leading interlayer $s_{\pm}$ instability in a two-orbital Hubbard–Hund model [37]. Constrained-path Monte Carlo instead identifies a crossover between $d$-wave and $s_{\pm}$ pairing controlled by Hund’s coupling and crystal-field splitting [38]. Other multiorbital calculations find a strong sensitivity to spin-fluctuation character and even possible odd-frequency pairing [39]. For ambient-pressure films, weak- and strong-coupling calculations find interlayer $s_{\pm}$ pairing even without a $\gamma$ Fermi-surface pocket [40]. Finally, functional-renormalization-group calculations show that nonlocal Coulomb interactions can either suppress or promote interlayer $s_{\pm}$ pairing, depending on their orbital structure [41].

These studies incorporate electronic correlations in substantially more detailed orbital models, but address a different question from the one posed here. We deliberately retain only a minimal bilayer $d_{x^2-y^2}$ sector, with the effects of the other orbitals encoded in effective couplings, in order to isolate the consequences of enforcing the local no-double-occupancy constraint on the pairing saddle point. In particular, a sign reversal between Fermi sheets does not by itself eliminate the on-site anomalous amplitude, whereas a symmetry-pure $d_{x^2-y^2}$ gap has a vanishing on-site projection in the tetragonal limit. Our $SU(2)$ construction enforces this constraint self-consistently, allowing us to quantify its effect on the competition between interlayer $s_{\pm}$ and intralayer $d$-wave pairing. Thus our result is not that multiorbital $s_{\pm}$ pairing is impossible, but that it must satisfy an additional local condition whose impact can be exposed most clearly in a minimal model.

The distinction between the two materials thus motivates the question: how can one compare interlayer $s_{\pm}$ and intralayer $d$-wave pairing on equal footing while excluding forbidden on-site pairs? We now introduce the strong-coupling gauge formulation that supplies the necessary constraints.

Strong Coupling Formulation

The comparison above exposes a limitation of conventional BCS mean-field theory: a sign-changing gap on different Fermi surfaces need not satisfy the local exclusion of double occupancy. In a projected $t$–$J$ description, this restriction must remain operative even after the electronic degrees of freedom are represented by spinons and bosons. Slave-boson fractionalization replaces the projection by local constraints, whose Lagrange multipliers appear as emergent gauge fields. The familiar $U(1)$ formulation enforces the local charge constraint, but a pairing mean field must also respect the particle–hole structure of the projected Hilbert space.

The $SU(2)$ formulation of Wen and Lee [42], developed further by Flint et al. [43, 44], makes this requirement explicit. In a gauge where the boson condensate points along the $b$ direction, $\langle b\rangle=\sqrt{x}$ and $\langle a\rangle=0$, the three local constraints become, at the uniform mean-field level,

Here $x$ is the hole density per site, $\ell$ labels the layer, and $f$ denotes a fermionic spinon. The first equation fixes the spinon density; the remaining two are conjugate conditions that eliminate an on-site spinon pair amplitude. Although they represent two equations, together they impose the real and imaginary parts of a single complex pairing constraint. Their associated Lagrange multiplier acts as a self-consistent Coulomb pseudopotential, penalizing mean-field pairing states that would otherwise violate the local projection.

This distinction is especially important for the bilayer nickelates. A $d_{x^2-y^2}$ gap has no on-site component by symmetry in the tetragonal limit, whereas a multiband $s_{\pm}$ state is not automatically protected by its sign change between Fermi sheets. The constraints can therefore alter the competition between intralayer and interlayer pairing. In Section II, we discuss the relevant solid-state chemistry of La$_3$Ni$_2$O$_7$ and motivate a minimal bilayer $t$–$J$ model. Section III derives the fractionalized operators, the local $SU(2)$ gauge structure, and the gauge-fixed constraints; Section IV shows how the corresponding Lagrange multipliers enter the bilayer saddle-point equations and pairing instability. This is followed by a discussion of our main results in § 1.4, followed by a discussion in § 1.5.

Solid State Chemistry of La$_3$Ni$_2$O$_7$ - Induced Superexchange and Minimal Model

In the cuprates, in-plane antiferromagnetic superexchange governs the motion of doped Zhang–Rice singlets [24]. In our model of bilayer La$_3$Ni$_2$O$_7$, the $3d_{x^2-y^2}$ orbitals support both in-plane exchange, generated by virtual hopping between neighboring Ni sites, and interlayer exchange, induced by Hund’s coupling to the $3d_{3z^2-r^2}$ orbitals (Figure 5). For the half-filled $3d_{x^2-y^2}$ orbitals, the latter is

where $H_{\mathrm{Hund}}=-J_H\mathbf{S}_{x^2-y^2}\cdot\mathbf{S}_{z^2}$ and $t_{\perp}$ denotes the interlayer $3d_{3z^2-r^2}$ hopping (see Appendix § 1.0.1 for a detailed derivation). Taking $J_H\simeq0.6$–$1$ eV and a $3d_{z^2}$ bonding–antibonding splitting $2t_{\perp}\simeq1.2$–$1.6$ eV gives $J_{\perp}^{\mathrm{phys}\ } \simeq0.05$–$0.2$ eV. For comparison, the in-plane exchange is $J_{\parallel}^{\mathrm{phys}}=4t_{\parallel}^2/U\simeq0.3$–$0.5$ eV. Thus $J_{\perp}^{\mathrm{phys}}$ is smaller but comparable to $J_{\parallel}^{\mathrm{phys}}$; both estimates depend on pressure and strain [45, 46].

With these considerations in mind, we adopt a projected bilayer $t$–$J$ model. Let $i$ label an in-plane unit cell, $\ell=1,2$ its two layers, and $X_{i\ell}^{ab}=|a\rangle\langle b|$ the Hubbard operators on the local three-state space $|0\rangle$, $|\uparrow\rangle$, $|\downarrow\rangle$. We use the Pauli-normalized local moment $\boldsymbol{\sigma}_{i\ell}=\sum_{ss'}X_{i\ell}^{ss'}\boldsymbol{\sigma}_{ss'}=2\mathbf{S}_{i\ell}$. Throughout the following symplectic decoupling, $J_{\parallel,\perp}$ denote the coefficients mul-

tiplying $\boldsymbol{\sigma}_i\cdot\boldsymbol{\sigma}_j$ (the factor conventionally written $1/N$ has been absorbed into $J$). Consequently $J_{\parallel,\perp}^{\mathrm{phys}}=4J_{\parallel,\perp}$ when $J^{\mathrm{phys}}$ multiplies $\mathbf{S}_i\cdot\mathbf{S}_j$. This choice gives the same coupling $J$ in the two symplectic bond channels. The real-space Hamiltonian is

The hopping matrices are Hermitian, $t_{ji}^{\parallel}=(t_{ij}^{\parallel})^*$; the first sum includes both hopping directions. The interlayer hopping need not be restricted to vertical rungs: its in-plane displacement dependence is required by the momentum-dependent bilayer splitting used below. Our convention is $t_{i,i\pm\hat{x}}^{\parallel}=t_{i,i\pm\hat{y}}^{\parallel}=t/2$, so that $\epsilon_{\parallel}(\mathbf{k})=-t(\cos k_x+\cos k_y)$. We choose $t_{ij}^{\perp}$ such that its Fourier transform obeys $-t_{\perp}(\mathbf{k})=t'(\cos k_x-\cos k_y)^2$. Explicitly, its nonzero Fourier coefficients are $t_0^{\perp}=-t'$, $t_{\pm2\hat{x}}^{\perp}=t_{\pm2\hat{y}}^{\perp}=-t'/4$, and $t_{\pm\hat{x}\pm\hat{y}}^{\perp}=t'/2$ (each of the four diagonal displacements). The superexchange estimates above use the conventional physical-spin coefficients $J_{\parallel,\perp}^{\mathrm{phys}}$; the $J$ entering the Pauli-normalized Hamiltonian and the Hubbard–Stratonovich action is one quarter of these values. Density-density pieces $-J_{ij}n_i n_j/4$, often included in the canonical $t$–$J$ Hamiltonian, are omitted here.

For comparison with the band structure and with the subsequent mean-field equations, Fourier transforming the projected hopping as a whole, rather than Fourier transforming its constituent slave particles, gives the single-particle Bloch Hamiltonian

where $\epsilon_{\parallel}=-t(\cos k_x+\cos k_y)$ and $\epsilon_{\perp}=t'(\cos k_x-\cos k_y)^2$. The real-space form, Eq. (6), is the starting point for the fractionalization below.

This replicates the qualitative features of the $\alpha$- and $\beta$-pockets on the Fermi surface as seen in the ARPES data [47, 48]. In our perspective of the strong coupling model, since only the $3d_{x^2-y^2}$ orbitals participate in the formation of the RVB state, the system is at a filling of $1/4$ (one electron for every two orbitals at each site), and is consequently at a hole doping of $0.5$ — zero doping corresponds to two electrons per site (one electron per orbital) that are localized, and then half a hole is added to each nickel ion. This effective model therefore treats La$_3$Ni$_2$O$_7$ as a heavily hole-doped cuprate system.

Wen–Lee path integral formulation of the model

We combine the in-plane position and layer label into a single site index $j\equiv(\mathbf{R}_j,\ell)$, with $\ell=1,2$. A bond $\langle ij\rangle$ may therefore be either an in-plane nearest-neighbor bond or an interlayer rung; its exchange coupling is $J_{ij}=J_{\parallel}$ or $J_{\perp}$, respectively. Following Wen and Lee and the symplectic formulation of Flint et al. [42, 44], we represent the projected electron by a fermionic spinon and a two-component charge boson. Let us define

The physical projected electron is the gauge-singlet composite

Note that the two bosons do not describe independent holes, they are required to encode the full SU(2) gauge freedom. The Pauli-normalized spin density used in Eq. (6) becomes, in the projected spinon representation,

The three generators of the local $SU(2)$ redundancy may be expressed as

The longitudinal constraint fixes the charge, while the transverse constraints remove unphysical on-site spinon pairing. All three are needed to recover the physical projected Hilbert space.

In imaginary time the constraints can be incorporated into the covariant derivatives of the doublets. With $\boldsymbol{\tau}$ the Pauli matrices acting on the gauge (Nambu) indices, define $w_j=\mathbf{w}_j\cdot\boldsymbol{\tau}$. One can then, in terms of these fractionalized fields, express the path integral as

Here $H_{tJ}$ denotes Eq. (6) expressed through the fractionalized fields. The factor $1/2$ avoids double counting the two Nambu-related spinors; the precise normalization and constant shift of the field are fixed by the constraints in Eqs. (11). In particular, $w_j=\vec{w}_j\cdot\vec{\tau}$ is a traceless gauge connection that couples to the isospin of the $f$-electrons, imposing the constraints $\psi^\dagger\vec{\tau}\psi+B^\dagger\vec{\tau}B=0$. Under a site- and time-dependent $g_j(\tau)\in SU(2)$,

the composite $\psi^\dagger_{j\sigma}B_j$ and the full action are invariant. In a convention where the Lagrange multipliers are real, the Euclidean temporal connection is understood with the corresponding factor of $i$; we absorb this into $w_j$ here.

To decouple the spins in the exchange term, we use the symplectic completeness relationship of $Sp(2N)$ [44]. For even spin-component number $N$, define $\epsilon_{\sigma\sigma'}=\tilde{\sigma}\delta_{\sigma,-\sigma'}$ and generators $S^a$ normalized by $\operatorname{Tr}(S^aS^b)=2\delta^{ab}$. The symplectic projector $P_{S}A=(A-\epsilon A^T\epsilon^T)/2$ gives

At the physical two-component point the generators are Pauli matrices. We do not take a one-component spin representation: instead, we absorb the conventional prefactor $1/N$ into the exchange constant $J_{ij}$, so no explicit $N$ appears in the physical model. Inserting Eq. (16) into $J_{ij}\boldsymbol{\sigma}_i\cdot\boldsymbol{\sigma}_j$ and normal ordering operators on different sites gives the symplectic exchange form

If we introduce oriented hopping and pairing operators

so that Eq. (17) becomes

where the displayed operator order follows the normal-ordering convention in Eq. (17). The two channels have the same coefficient $J_{ij}$, with $J_{ij}=J_{\parallel}$ on in-plane bonds and $J_{\perp}$ on the bonds between the bilayer sites. For each bond we introduce complex Hubbard–Stratonovich fields $V_{ij}$ and $\Delta_{ij}$ conjugate to $D_{ji}$ and $P_{ji}$. Following standard procedure, completing the square in the hopping channel gives

and the pairing channel analogously gives $\lvert\Delta_{ij}\rvert^2/J_{ij}-\Delta_{ij}^{*}P_{ji}-\Delta_{ij}P_{ji}^{\dagger}$. Varying the auxiliary fields near the saddle point yields $V_{ij}=J_{ij}\langle D_{ji}\rangle$ and $\Delta_{ij}=J_{ij}\langle P_{ji}\rangle$. In Nambu form these channels assemble into the link matrix

The link transformation ensures that the quadratic spinon term is gauge covariant.

For a doped homogeneous saddle point, the boson is assumed to condense, $\langle B_j\rangle\ne0$. A local $SU(2)$ transformation rotates this doublet into the $b$ direction; we choose its remaining phase to be real:

With $j=(\mathbf{R}_j,\ell)$, the projected hopping becomes

The remaining mean-field constraints are, with the combined site label restored,

For the real, time-reversal-symmetric ansatz, the transverse components of the temporal field reduce to an onsite pair multiplier $w$, while its longitudinal component gives $\lambda$. The gauge-fixed constraint Hamiltonian is

Combining the hopping, exchange links and constraints, the full real-space saddle-point Lagrangian can be summarized as

The sums over $\langle ij\rangle$ include the in-plane bonds and interlayer rungs, with their respective couplings. We use equation (28) as the real-space starting point for Section IV, where we transform the link fields to momentum space and derive the bilayer Bogoliubov–de Gennes spectrum and saddle-point equations. While the longitudinal Lagrange multiplier $\lambda$ constrains the spinon filling to the doping $x$, The on-site field $w$ enforces the vanishing of the net anomalous spinon amplitude and thus acts as a strong-coupling Coulomb pair pseudopotential.

Methods: bilayer saddle point and pairing instability

We now specialize the real-space saddle-point Lagrangian, Eq. (28), to homogeneous intralayer and interlayer link fields. The projected hopping and the Hubbard–Stratonovich decoupling have already been derived in Section III; here we fix the momentum-space conventions and obtain the bilayer quasiparticle spectrum.

With $j=(\mathbf{R}_j,\ell)$, the projected hopping becomes

The uniform constraints are those in Eq. (27).

Let $V_x=V_y=V_{\parallel}/2$ in each layer, $V_z=V_{\perp}$ on each rung, and real in-plane pair fields $\Delta_x,\Delta_y$. Their Fourier transform is

where $\Delta_s=\Delta_x+\Delta_y$, $\Delta_d=\Delta_x-\Delta_y$, $\gamma_s=\cos k_x+\cos k_y$ and $\gamma_d=\cos k_x-\cos k_y$. If $N_s$ denotes the number of bilayer unit cells, the two layers together give the bond-field cost $2N_s(\Delta_x^2+\Delta_y^2)/J_{\parallel}=N_s(\Delta_s^2+\Delta_d^2)/J_{\parallel}$ in this normalization. The exchange constants are Pauli-normalized, as in Section III; bond counting fixes their relation to the momentum-space amplitudes. We compare separate symmetry sectors, $a=s$ or $a=d$, rather than assuming a state that spontaneously mixes them. The $s$ sector allows both $\Delta_s$ and a rung singlet $\Delta_z$, producing bonding/antibonding gaps $\Delta_s\gamma_s\pm\Delta_z$; the $d$ sector has $\Delta_d\gamma_d$, which is orthogonal to the on-site $s$-wave constraint by lattice symmetry.

To take care of the Nambu and interlayer degrees of freedom and to allow ourselves to write the final mean-ends_mid=1 field Bogoliubov-de Gennes (BdG) Hamiltonian with cleaner notation, we introduce the four-component spinor

spanned by the Pauli operators $\boldsymbol{\alpha}$ and $\boldsymbol{\tau}$. $\boldsymbol{\alpha}$ acts on the layer degree of freedom on the four-component space,

and

Similarly, the Nambu isospin degrees of freedom are acted on by $\boldsymbol{\tau}$,

With these operators defined, we can rewrite the mean-field Hamiltonian as shown below:

Here $N_s$ is the number of bilayer unit cells, so that there are $2N_s$ physical nickel sites. The constraint term $2N_s\lambda x$ is required by the chosen Nambu convention.

The free energy for the above Hamiltonian is calculated by finding the eigenvalues of the BdG Hamiltonian, and can be written as

where the two Bogoliubov quasiparticle bands (bonding and antibonding) have the energies

The energies $E_{\parallel}$, $E_{\perp}$ are a combination of the true dispersion on the bands combined with the dispersion generated from decoupling the superexchange into the particle-hole channel.

The solution to this mean-field free energy lies at the saddle point. By extremizing the free energy with respect to $\lambda$, $V_{\parallel}$, $V_{\perp}$, $w$, $\Delta_a$ and $\Delta_z$ respectively, we obtain the following self-consistency solutions.

where the integral $\int_{\mathbf{k}}=(2\pi)^{-2}\int_{\mathrm{BZ}}(\mathrm{d}^2 k)=\frac{1}{N_s}\sum_{\mathbf{k}}$ is the average over the Brillouin zone (discrete or continuous). At the critical temperature, the pairing amplitudes and the pair repulsion constraint vanish, so $E_{\pm}(\mathbf{k})=\lvert\alpha_{\pm}(\mathbf{k})\rvert$ (not $\alpha_{\pm}$ without an absolute value). In expressions containing $\alpha_{\pm}/E_{\pm}$, the combination is $\operatorname{sgn}(\alpha_{\pm})$, so that $\tanh(E_{\pm}/2T)\alpha_{\pm}/E_{\pm}=\tanh(\alpha_{\pm}/2T)$. The first three equations are then simplified to

where $\beta_c$ is the critical inverse temperature, and the second three equations are replaced with the condition $\det\chi=0$, where $\chi$ is the susceptibility matrix.

The elements of the susceptibility matrix are:

These self-consistency conditions are solved simultaneously to find the critical temperature $T_c$. At the superconducting transition, the determinant of the susceptibility $\det\{\chi\}$ changes sign, corresponding to a sign change in one of the eigenvalues of $\chi$. By finding the eigenvector corresponding to this unstable eigenvalue, we can obtain the relative proportions of extended-s ($\Delta_s(c_x+c_y)$) to interlayer ($\Delta_z$) pairing, allowing us to plot the gap at the spinon Fermi surface.

Results

By using the Wen-Lee fractionalization scheme and numerically minimizing the free energy by solving the coupled self-consistency equations, we have found the critical temperature $T_c$ for different values of doping. We present a full sweep for the hole doped regime for completeness; however, as stated before, the bilayer system La$_3$Ni$_2$O$_7$ corresponds to a doping of $x=0.5$. Our results are summarized in Figure 4. In both cases, we have chosen values

of $J_{\perp}$ and $J_{\parallel}$ much smaller than the bandwidth $W$ (in our case, $J_{\parallel}/W\sim1/8$), but the qualitative takeaways from the results are by and large insensitive to this fact. In the weak $J_{\perp}$ regime, we chose the value of $J_{\perp}=\frac{1}{2}J_{\parallel}$. In this regime, we see that for small doping, the $T_c$ for $d$-wave order far exceeds $s_{\pm}$ order. In the low doping regime, there is no effect of the Coulomb pseudopotential. This is to be expected, since at half-filling the Coulomb constraint is satisfied automatically. Upon hole doping, the $T_c$ for both $d$-wave and $s_{\pm}$ fall, but at a doping of 0.25 a sharp upturn is noticed in the $s_{\pm}$ $T_c$. This sharp upturn coincides with the pairing gap acquiring the same sign on both Fermi surfaces (Figure 6a), leading to a mean-field solution with no nodes, a higher free energy and therefore a higher $T_c$. At even larger doping ($x\approx0.75$), $s_{\pm}$ manages to surpass $d$-wave in $T_c$. However, once the Coulomb pseudopotential constraint is added to the solution, the $s_{\pm}$ $T_c$ is suppressed significantly, never crossing the $d$-wave $T_c$. By imposing the Coulomb constraint, the pair function is forced to change signs between Fermi surfaces (Figure 6b).

We had argued previously in Sec. II that the induced interlayer superexchange between the $3d_{x^2-y^2}$ orbitals is expected to be smaller than the in-plane superexchange. However, we consider for the sake of completeness the case where $J_{\perp}$ is of comparable magnitude to $J_{\parallel}$ — namely, $J_{\perp}=J_{\parallel}$ (Figure 4, right). In this case, the competition between the two orders is stronger. For low doping, the $T_c$ for $d$-wave is nearly identical to $s_{\pm}$. This is not that surprising since the interlayer and intralayer couplings are the same strength. This coincidence continues down to moderate doping ($x\sim0.3$), where $d$-wave wins out over $s_{\pm}$. Note that in both this and the previous plot, the $T_c$ for $d$-wave follows identical trajectories. This is because band renormalization effects on the $d$-wave pairing state are minimal for weaker interlayer coupling. At half-doping, same as before, the $s_{\pm}$ $T_c$ starts to rise, emerging as the dominant state for $x\gtrsim0.7$. Upon the addition of Coulomb repulsion, the $T_c$ for $s_{\pm}$ is suppressed below $d$-wave for all doping beyond $x\sim0.7$.

The addition of the Coulomb pseudopotential at the mean-field level, by suppressing onsite pairing, always reduces the $T_c$ of the $s_{\pm}$ pairing order. One way to think about this that is particularly elucidating is to look at the susceptibility matrix $\chi$ defined in (48). We can rewrite this as

where $\chi_0$ is the bare susceptibility without the Coulomb psuedopotential constraint. At $T = T_c$, $\det\chi$ vanishes, and using block matrix properties, one can rewrite this vanishing condition as $\chi_{ww} \det(\widetilde{\chi}_0) = 0$, where

We observe that, according to (55), the addition of the Coulomb psuedopotential creates an effective susceptibility $\widetilde{\chi}_0$ that is the original susceptibility with a positive semidefinite correction. This correction reduces the eigenvalues of the susceptibility matrix, thereby reducing $T_c$ for $s_{\pm}$ pairing. For $d$-wave pairing, the correction to the bare susceptibility is 0, since the $d$-wave symmetry is completely orthogonal to the Coulomb pseudopotential as well as the dimer pair order $\Delta_z$.

Discussion and Outlook

In this work, we have investigated the influence of strong local Coulomb repulsion on the competition between $d$-wave and $s_{\pm}$ superconductivity in bilayer nickelates. Using the SU(2) fractionalization of the bilayer $t$-$J$ model, we have shown that the suppression of on-site pairing introduces a constraint on the superconducting order parameter that is absent in conventional mean-field treatments. This constraint suppresses the $s_{\pm}$ pairing susceptibility through a positive-semidefinite correction, while leaving the $d$-wave channel unaffected by symmetry.

The Coulomb constraint does not, however, prohibit $s_{\pm}$ superconductivity. Strong interlayer exchange $J_{\perp}$ enhances interlayer pairing and can partially overcome the Coulomb suppression. The resulting competition is controlled by the balance between intralayer and interlayer superexchange, subject to the requirement that the condensate avoid local double occupancy. For the parameters studied here, the constraint favors $d$-wave superconductivity over much of the phase diagram, including the nominal doping $x = 0.5$ of La$_3$Ni$_2$O$_7$.

An important motivation for studying the nickelates is the possibility of identifying the essential ingredients of high-temperature superconductivity in both the cuprates and nickelates. Both families share closely related square-planar transition-metal–oxygen structures, strongly correlated $d_{x^2-y^2}$ orbitals, and charge fluctuations involving $d^8$ and $d^9$ configurations. Both also exhibit strange-metal behavior. Given the rarity of high-temperature superconductivity, these similarities provide a compelling motivation to seek a common microscopic description rather than separate, material-specific explanations.

Several directions for future work emerge from this perspective.

• Effective low-energy models. In the cuprates, Zhang–Rice singlet formation provides a route to an effective single-band $t$-$J$ model. In the bilayer nickelates, strong interlayer hybridization of the $d_{z^2}$ orbitals may play an analogous role by splitting these states away from the Fermi energy. Establishing whether this leaves a predominantly $d_{x^2-y^2}$ low-energy electronic structure is an important experimental and theoretical challenge. Comparison with infinite-layer nickelates, where the orbital hierarchy arises without strong bilayer hybridization, may help determine whether these materials admit a common low-energy description.

• Experimental tests of pairing symmetry. Recent ARPES measurements on ambient-pressure (La,Pr,Sm)La$_3$Ni$_2$O$_7$ films [49] reveal nodeless gaps on the $\beta$ and $\gamma$ Fermi-surface pockets, challenging the $d$-wave state favored by our calculations. The additional $\gamma'$ pocket, neglected in our minimal model, may modify the competition between pairing channels by allowing cancellations of the local anomalous amplitude between different bands. Phase-sensitive measurements, together with comparisons between ambient-pressure films and pressurized bulk materials, could help establish whether the observed superconductivity is consistent with an $s_{\pm}$ state satisfying the Coulomb constraint.

• Beyond mean-field theory. Multiorbital calculations incorporating the $d_{z^2}$ states and Hund’s coupling are needed to establish the robustness of our results. Cluster DMFT and related approaches could also examine how the Coulomb constraint is realized when the anomalous self-energy acquires momentum, orbital, and frequency dependence. The central question is whether these additional degrees of freedom qualitatively alter the competition between $d$-wave and $s_{\pm}$ pairing, or can be integrated out to recover a simpler effective theory.

We conclude by noting that the discovery of ambient-pressure nickelate superconductivity opens an exceptional opportunity to compare the cuprates and nickelates through common experimental and theoretical probes. It transforms the cuprate problem from the study of a unique family of materials into a comparative search for universal principles. The challenge is not to construct separate explanations for each superconductor, but to identify the minimal physics underlying strange-metal behavior and high-temperature superconductivity in a unified family of materials—the cupra-nickelates.

Acknowledgements

The authors would like to thank Liam L. H. Lau, Indra Gankhuyag, Po-Jui Chen and Andreas Gleis for helpful discussions. This work was supported by the Office of Basic Energy Sciences, Material Sciences and Engineering Division, U.S. Department of Energy (DOE) under Contract DE-FG02-99ER45790.

Calculation for the Hund's induced superexchange in the $3d_{x^2-y^2}$ orbitals of La$_3$Ni$_2$O$_7$

Consider a two-site model (site 1 and 2) with two orbitals on each site labelled by $a$ and $b$ (in our case, the $3d_{z^2}$ and $3d_{x^2-y^2}$ orbitals in Ni). Assume that the $a$ orbitals have a strong overlap, so the dominant energy scale in the problem is the hopping amplitude $t$.

Now assume there is a weaker Hund’s coupling between the $a$- and $b$-orbitals on each site.

We want to find the induced superexchange between the $b$-orbitals from virtual spin-flip fluctuations from the Hund’s exchange. We first recast $H_{aa}$ in terms of bonding and antibonding orbitals.

which gives us

The ground state of the system is the bonding-orbital singlet $\lvert S_B\rangle\equiv B^\dagger_\uparrow B^\dagger_\downarrow\lvert0\rangle$. An excitation to the antibonding orbital costs energy $2t$. Now we rewrite the onsite spin operators for the $a$- orbital in terms of the bonding and antibonding orbitals.

Similarly, the other spin operator is defined as

Here, $\mathbf{T}$ is the spin-flip operator. The action of the Hund’s coupling on $\lvert S_B\rangle$ is clear - both $\mathbf{S}_B$ and $\mathbf{S}_A$ evaluate to zero, so the only remaining contribution is:

The action of the spin-flip operator on the bonding orbital singlet can be calculated component by component. We can ignore the first terms in $\mathbf{T}$ since they include a B creation operator, which acts trivially on the ground state singlet.

z-component:

x-component:

y-component:

Each of these three states are orthogonal triplets with the same energy above the ground state ($2t$). The effective coupling between the $b$-orbitals is

The cross terms vanish ($\langle T_xT_y\rangle=\langle T_zT_y\rangle=\langle T_xT_z\rangle=0$), and the only contributions left are

This leaves us with

We have dropped irrelevant constant terms in the last line. This gives us an effective effective AFM superexchange between the $b$-orbitals of $J_b=\frac{J_H^2}{8t}$.