Intrinsic constraint on Tc for unconventional superconductivity

Despite the century-long pursuit of high-temperature superconductors, the possible existence of a theoretical upper limit to their transition temperature (T_c) under ambient pressure remains unsettled^1,2,3,4,5. Both mean-field and weak-coupling Eliashberg theories⁶ predict an artificial T_c that grows continuously with increasing pairing interaction, while experiments often find superconducting domes with maximum T_c near the phase boundaries of some long- or short-range orders associated with spin, charge, orbital, or structural degrees of freedom^{7,8,9,10,11,12,13,14,15,16,17,18}. The dome implies a dual role of the many-body interaction¹⁹, which may not only provide the pairing glue but also induce competing orders that constrain the maximum T_c. However, they are mostly external factors associated with instabilities of other channels. One may wonder if any intrinsic constraint on T_c may exist owing solely to the pairing instability.

Important lessons may be learned from cuprate high-temperature superconductors in the underdoped region, where strong pairing interactions relative to the renormalized effective quasiparticle hopping parameters favour short-range electron pairs²⁰ that in some literatures are thought to form already at high temperatures but only become superconducting when a (quasi-)long-range phase coherence is developed^21,22. This raises a few general questions: What is the true strong coupling limit of the pairing state? How is this strong coupling state related to the high-temperature superconductivity? Would it put any intrinsic constraint on the maximal value of T_c? Since the absolute magnitude of T_c is determined by certain basic energy scale, such as the pairing interaction J, the question of maximum T_c turns into the question of their maximum dimensionless ratio T_c/J. To address these important issues and gain insights into possible intrinsic constraints on T_c, we propose here to discard instabilities from all other channels such as magnetic or charge orders and focus only on the pairing instability, since all other instabilities are expected to compete with the superconductivity and further suppress T_c. Their contributions to the pairing can all be included phenomenologically in a pairing interaction term.

Theoretically, one may derive various ratios with respect to other measurable energy scales such as the Fermi energy and the superfluid density. However, it has been shown that these constraints may be violated in artificial models⁴, or even in real materials²³, thus preventing a useful bound for constraining T_c. To avoid such complication, instead of deriving a model-independent constraint, we first restrict ourselves to a minimal effective model that is most relevant in real correlated superconductors and includes only the quasiparticle hopping and nearest-neighbour spin-singlet pairing interaction. For the one-band model on a square lattice, we find a plaquette state in the strong-coupling limit that breaks both the translational and time-reversal symmetries and exhibits unusual spectral properties with a pseudogap or insulating-like normal state. The d-wave superconductivity emerges as the plaquettes melt and short-range electron pairs get mobilized to attain long-distance phase coherence at a reduced pairing interaction. A tentative phase diagram is then constructed where T_c/t reaches its maximum at the plaquette quantum critical point (QCP), resembling those often observed in experiments with other competing orders. This suggests some intrinsic constraints that prevent T_c from exhausting all kinetic or pairing energies in order to achieve a delicate balance between pairing and phase coherence. We then extend the calculations to more general models with either nearest-neighbour or onsite pairings, single layer or two-layer structures, intralayer or interlayer pairings, and obtain a maximum T_c/J ≈ 0.04−0.07. A close examination of existing experiments in known unconventional superconductors, including cuprate, iron-based, nickelate, and heavy fermion superconductors, seems to quite universally support the obtained ratio, indicating that these families, possibly except the infinite-layer nickelates, have almost reached their maximum T_c allowed by their respective spin exchange interactions. A room-temperature superconductor would then require a much larger pairing interaction beyond 400–700 meV within the current theoretical framework, which seems unrealistic from a single mechanism in correlated electron systems under ambient pressure. Our work therefore provides a useful criterion that may help to avoid futile efforts in exploring high-temperature superconductors along wrong directions. It also points out the necessity of new pairing mechanisms, possibly combining different pairing interactions, in order to achieve the room-temperature superconductivity.

Model and method

We start by first considering an effective one-band t–J type model on the square lattice, which will later be extended to more general cases (see Methods). The model Hamiltonian is written as:

where t_ij is the renormalized quasiparticle hopping parameter, μ is the chemical potential, and the pairing interaction is written in terms of the spin-singlet operator ({psi }_{ij}=frac{1}{sqrt{2}}({d}_{idownarrow }{d}_{juparrow }-{d}_{iuparrow }{d}_{jdownarrow })) between nearest-neighbour sites, where ({d}_{isigma }({d}_{isigma }^{dagger })) is the annihilation (creation) operator of the quasiparticles to be paired. For the exchange mechanism, the pairing interaction may be induced by the nearest-neighbour antiferromagnetic interaction and the attractive charge density interaction. As shown in Methods, J is then simply the nearest-neighbour exchange interaction, whose importance in cuprates has been justified in numerous experiments^{24,25,26,27,28}. To study the superconductivity, we introduce the complex auxiliary fields Δ_ij to decouple the pairing term²⁹:

To avoid the negative sign problem, we assume a static approximation, ({Delta }_{ij}(tau )to {Delta }_{ij}=| {Delta }_{ij}| {e}^{i{theta }_{ij}}), and employ the auxiliary field Monte Carlo approach^{30,31,32,33,34,35,36,37,38}. Following the standard procedure, we integrate out the fermionic degrees of freedom and simulate the final effective action only of the pairing fields by the Metropolis algorithm³⁸.

This method ignores dynamic fluctuations of the pairing fields but takes full consideration of their spatial and thermal fluctuations. We can investigate phase correlations of the pairing and determine T_c based on long-distance phase coherence rather than the BCS-type mean-field transition. The validity of our method in estimating the maximum T_c has been verified in the recently-discovered bilayer and trilayer nickelate superconductors^39,40 and by its consistency with the rigorous Quantum Monte Carlo simulations for the attractive Hubbard model². However, its applications in analyzing certain dynamical properties of the pairing fields are limited. To maximize T_c, we have also ignored all other instabilities outside the pairing channel. The effect of the Gutzwiller constraint is approximated by treating t_ij as free tuning parameters.

For numerical calculations, we consider a 10 × 10 square lattice with periodic boundary conditions and include only the nearest-neighbour hopping t and the next-nearest-neighbour hopping ({t}^{{prime} }=-0.45t) as in cuprate high-temperature superconductors^41,42. The chemical potential is fixed to μ = −1.4t. The presented results have been examined and found qualitatively consistent for other values of the parameters or on a larger lattice. A twisted boundary condition is used for spectral calculations⁴³.

Theoretical phase diagram

Figure 1a shows a typical theoretical phase diagram for the one-band square lattice model, where we have intentionally plot T_c/t against t/J. A nonuniform plaquette state emerges at sufficiently strong paring interaction formed of 2 × 2 blocks induced by high-order pair hopping in the effective action of the pairing fields after integrating out the electron degrees of freedom. A typical pairing configuration of the plaquette state is given in the inset of Fig. 1. The paring amplitudes are relatively stronger on internal bonds of the 2 × 2 plaquettes and weaker on their links. As will be discussed later and in the Methods, the plaquette state simultaneously breaks the lattice translational symmetry and the time-reversal symmetry and coexists with a bond charge order, though the electron density remains uniform on all sites. Its transition temperature T_□ decreases with increasing t/J and diminishes at the QCP (t/J ≈ 0.27), where the plaquettes melt completely and uniform superconductivity emerges with a maximum T_c/t ≈ 0.08 for the chosen parameters (see Methods). Tuning the next-nearest-neighbour hopping and the chemical potential may slightly change the ratio and the location of the QCP, but does not alter the qualitative physics. Inside the plaquette state, the superconductivity is also spatially modulated and its T_c is greatly reduced as the pairing interaction further increases. The nonmonotonic evolution of T_c resembles typical phase diagrams observed in many unconventional superconductors with other competing orders such as long-range magnetism, charge density wave, or nematicity^{7,8,9,10,11,12,13,14,15,16}. However, the plaquette state reflects the internal instability in the pairing channel that constrains the magnitude of T_c for the d-wave superconductivity. Near the plaquette QCP, the superconductivity also breaks the time-reversal symmetry below T_BTRS. As t/J decreases, the T_BTRS transition line merges with the plaquette transition T_□, marking the simultaneous time-reversal symmetry breaking of the plaquette state. At high temperatures, the normal state exhibits pseudogap-like behaviour whose onset temperature T_p follows closely the variation of T_c or T_□^44,45 determined from the specific heat C_v or the temperature derivative of the quasiparticle density of states at the Fermi energy dN(0)/dT. As shown in Fig. 1b, we find peaks in the specific heat for all transitions at T_□, T_c, and T_BTRS, while in dN(0)/dT the feature at T_c is greatly suppressed for t/J < 0.27. Here and after, J is set as the energy unit if not explicitly noted.

a The T/t − t/J phase diagram, showing the onset temperature T_p of pseudogap-like behaviour determined from the suppression of the quasiparticle density of states at the Fermi energy N(0), the plaquette transition temperature T_□ from the pairing field amplitude distribution p(∣Δ∣) and the peaks in the specific heat C_v and the temperature derivative of the quasiparticle density of states dN(0)/dT, the superconducting transition temperature T_c from the long-distance phase coherence based on the phase mutual information and the properties of the Berezinskii–Kosterlitz–Thouless (BKT) transition for two-dimensional superconductivity, and the temperature T_BTRS for superconductivity with broken time-reversal symmetry from the deviation of the pairing field phases along x and y bonds attached to the same site. The inset shows a typical configuration of the pairing field inside the plaquette state at low temperatures, where the size of the symbols represents the amplitude ∣Δ_ij∣ and the colour denotes the sign of the phase θ_ij. Note that for small t/J, T_BTRS extends to the plaquette transition T_□, which breaks simultaneously the time-reversal symmetry and the translational symmetry and coexists with a bond charge order (see Methods). b Temperature evolution of C_v and dN(0)/dT for t/J = 0.15, 0.22, 0.30, showing peaks or shoulders at T_□, T_c, and T_BTRS.

Full size image

Plaquette state at strong coupling

The plaquette state and its phase transition may be seen in the joint distribution (p(| Delta {| }_{{boldsymbol{0}}}^{x},,| Delta {| }_{{boldsymbol{0}}}^{y})) of the paring amplitudes along the x and y directions attached to the same site 0 or the marginal distribution p(∣Δ∣) of the pairing field amplitudes on all bonds. As shown in Fig. 2a, (p(| Delta {| }_{{boldsymbol{0}}}^{x},,| Delta {| }_{{boldsymbol{0}}}^{y})) at low temperatures displays a four-point structure due to the nonuniform pairing configurations. As t/J increases, the four points gradually shrink into a single point, where the translational symmetry is recovered and the plaquette state melts into the uniform superconductivity. Correspondingly, the amplitude distribution p(∣Δ∣) also contains two peaks in the plaquette state. As shown in Fig. 2b for t/J = 0.15, these peaks get gradually broadened with increasing temperature and merge into a single peak above T_□.

a The joint distribution function (p(| Delta {| }_{{boldsymbol{0}}}^{x},,| Delta {| }_{{boldsymbol{0}}}^{y})) of the pairing field amplitudes (| {Delta }_{{boldsymbol{0}}}^{x}|) and (| {Delta }_{{boldsymbol{0}}}^{y}|) along x and y directions attached to the same site 0 for t/J = 0.15, 0.23, 0.30 at a very low temperature T/J = 0.0001. b Evolution of the marginal distribution p(∣Δ∣) of the pairing amplitude on all bonds with temperature for t/J = 0.15. c Comparison of the energy-momentum dependent spectral function and extracted dispersions (solid lines) at k_x/π = 0.44 at low and high temperatures for t/J = 0.15. The grey vertical lines mark the Fermi vector k_F that clearly differs from the wave vector k_G where the dispersions bend backwards.

Full size image

At sufficiently low temperatures, the plaquette state may develop long-distance phase coherence to form the superconductivity but exhibits unusual spectral features due to the nonuniform spatial distribution of the pairing amplitudes. As shown in Fig. 2c for t/J = 0.15, its momentum-energy dependent spectral function at negative energies splits into two sets of dispersions. One dispersion resembles that of uniform superconductivity, but its back-bending vector k_G differs consistently from the Fermi vector k_F, which has also been observed experimentally for possible pair density wave (PDW) state^46,47,48. At high temperatures, the two dispersions recombine into a single curve pointing upwards even in the normal state. The gap indicates a pseudogap or insulating-like phase due to the large nearest-neighbour pairing interaction. This suggests that the normal state may also undergo a metal-insulator transition as t/J decrease, a phenomenon observed in cuprate superconductors under high pressure but unexplained⁴⁹. At intermediate temperature T_c < T < T_□, the superconducting phase coherence is lost and the plaquette state with broken time-reversal symmetry is in a sense similar to the fermionic quadrupling phase proposed earlier in experiment^50,51.

Time-reversal symmetry breaking

The time-reversal symmetry breaking may be seen from the probabilistic distribution p(δθ_xy) of the phase difference (delta {theta }_{xy}={theta }_{{boldsymbol{0}}}^{x}-{theta }_{{boldsymbol{0}}}^{y}) of the pairing fields along the x and y directions. The results are shown in Fig. 3a for three different values of t/J at very low temperature. For small t/J = 0.15 in the plaquette state, the existence of multiple peaks mark the phase difference on different bonds. These peaks develop into a two-peak structure at higher temperatures and then merge into a single peak at δθ_xy = π above T_□ (see Methods). For large t/J = 1.0, there exists a single maximum around δθ_xy/π = 1, which signals the uniform d-wave superconductivity with opposite sign of the pairing field along the x and y directions. Quite unexpectedly, for t/J = 0.3, we still have a single peak but its position deviates from δθ_xy/π = 1. To see such a variation more clearly, Fig. 3b plots the average deviation (∣〈δθ_xy〉∣) and the smallest deviation ((| delta {theta }_{xy}^{min }|)) of the peak positions. While ∣〈δθ_xy〉∣ evolves nonmonotonically and reaches a minimum at the plaquette QCP, (| delta {theta }_{xy}^{min }|) keeps increasing with t/J. Interestingly, the two quantities become equal beyond the plaquette QCP but only approach π at a much larger t/J ≈ 0.5.

a The probabilistic distribution of the phase difference along x and y directions (p(delta {theta }_{xy})=p({theta }_{{boldsymbol{0}}}^{x}-{theta }_{{boldsymbol{0}}}^{y})) for different values of t/J. b Evolution of the average phase difference ∣〈δθ_xy〉∣ and the minimum phase difference (| delta {theta }_{xy}^{min }|) determined by the peak positions as functions of t/J. The vertical line marks the plaquette QCP at t/J = 0.27. c Comparison of the gap function Δ(ϕ) near nodal and antinodal directions as functions of t/J, determined by the position of the positive-energy peak in the spectral function A(ϕ, ω). The temperature is T/J = 0.0001. d Temperature evolution of p(δθ_xy) at t/J = 0.3 in the intermediate phase. e The angle-dependent gap function Δ(ϕ) for different temperatures at t/J = 0.3, showing the evolution from a full gap at low temperatures to a partial gap at high temperatures.

Full size image

Under time-reversal operation, the phase of the pairing field changes sign so that δθ_xy → − δθ_xy (mod 2π). Thus, the deviation of the peak position from π around t/J = 0.15 indicates that the time-reversal symmetry is broken inside the whole plaquette state. For t/J = 0.3, it marks an intermediate region of uniform superconductivity that breaks the time-reversal symmetry, with the gap function ({Delta }_{{boldsymbol{k}}}propto cos ({{boldsymbol{k}}}_{x})+{e}^{-idelta {theta }_{xy}}cos ({{boldsymbol{k}}}_{y})propto {Delta }_{{boldsymbol{k}}}^{d}-icot frac{delta {theta }_{xy}}{2}{Delta }_{{boldsymbol{k}}}^{s}), representing d + is pairing with a nodeless gap. Here ({Delta }_{{boldsymbol{k}}}^{d}=cos ({{boldsymbol{k}}}_{x})-cos ({{boldsymbol{k}}}_{y})) is the d-wave component and ({Delta }_{{boldsymbol{k}}}^{s}=cos ({{boldsymbol{k}}}_{x})+cos ({{boldsymbol{k}}}_{y})) denotes an extended s-wave component from the nearest-neighbour pairing interaction. The onsite pairing is not included due to the strong Coulomb repulsion. We have therefore a two-stage transition from the plaquette to the uniform d-wave superconductivity, with an intermediate region that recovers the translational symmetry but still breaks the time-reversal symmetry. Similar d + is pairing may have been found under certain conditions in twisted double-layer cuprates⁵² and infinite-layer nickelates^53,54. In the latter case, it arises from the interplay of Kondo and superexchange interactions⁵⁵. Here it is associated with the quasiparticle hopping, (ito i+hat{x}to i+hat{x}+hat{y}to i+hat{y}to i). Integrating out the electron degrees of freedom leads to a term like ({rm{Re}}({Delta }_{i,i+hat{x}}{Delta }_{i+hat{y},i+hat{x}+hat{y}}{Delta }_{i,i+hat{y}}^{* }{Delta }_{i+hat{x},i+hat{x}+hat{y}}^{* })to {rm{Re}}({Delta }_{x}^{2}{Delta }_{y}^{* 2})propto cos (2delta {theta }_{xy})), while the second order hopping process such as (i+hat{x}to ito i+hat{y}) contributes a term ({rm{Re}}({Delta }_{x}{Delta }_{y}^{* })propto cos delta {theta }_{xy}). Their combined free energy may be minimized at δθ_xy away from 0 and π⁵⁶. Thus, time-reversal symmetry breaking represents an intrinsic tendency of the superconductivity with nearest-neighbour pairing at strong coupling, where the normal state is no longer a Fermi liquid.

To further confirm the two-stage transition, Fig. 3c plots the gap function Δ(ϕ) with t/J near the nodal and antinodal directions in the momentum space deduced from the spectral function. The gap near the antinode is always finite, but varies nonmonotonically with a maximum at the plaquette QCP t/J = 0.27, in good correspondence with the maximum T_c. By contrast, the gap near the nodal direction decreases continuously and only diminishes at t/J ≈ 0.5, confirming a full gap for 0.27 ≤ t/J ≤ 0.5 consistent with the above phase analysis. The transition temperature T_BTRS of the d + is phase may also be extracted from the temperature evolution of p(δθ_xy). As shown in Fig. 3d for t/J = 0.3, the peak in p(δθ_xy) gets broadened and moves gradually to δθ_xy = π as the temperature increases across T_BTRS. The angle-dependent gap functions are given in Fig. 3e, showing a fully gapped d + is pairing state and a nodal d-wave pairing state below and above T_BTRS, respectively. Note that the higher-temperature d-wave gap contains a finite gapless region on the Fermi surface, which has also been observed previously in some experiments⁵⁷.

Superconducting phase coherence

The superconducting transition is determined from the phase mutual information ({I}_{{boldsymbol{R}}}^{x/y}) of the pairing fields as well as the vortex number n_v (see Methods)^38,39. Figure 4a shows the semilog plot of the phase mutual information between two bonds of the largest distance R = (5, 5) for t/J = 0.15, 0.22, 0.30 on the 10 × 10 lattice. We find a slope change at low temperature, marking the establishment of long-distance phase coherence of the pairing fields. The slope change at higher temperature is associated with the onset of the spatial phase correlation, which has a temperature scale in rough agreement with T_p for t/J > 0.27 in Fig. 1 and is therefore responsible for the pseudogap above the superconducting T_c.

a The mutual information between two pairing field phases ({theta }_{(0,0)}^{x/y}) and ({theta }_{(5,5)}^{x/y}) of the distance (5, 5) and the normalized numerical derivatives of the vortex number dn_v/dT with temperature for t/J = 0.15, 0.22, 0.30, respectively. The vertical lines show the extracted T_c. b The distribution p(δθ₁) for different hopping at T/J = 0.0001, where δθ₁ is the phase difference between nearest-neighbour pairing fields (delta {theta }_{1}={theta }_{{boldsymbol{0}}}^{x/y}-{theta }_{left.(1,0)/(0,1)right)}^{x/y}). c The peak position in p(δθ₁) and the inverse of the fluctuation ({rm{std}}(delta {theta }_{i})=sqrt{langle {(delta {theta }_{i})}^{2}rangle }), where δθ_i is the phase difference between two nearest-neighbour (i = 1) or next-nearest-neighbour (i = 2) bonds along x or y directions. The two behave similarly for uniform superconductivity but differ in the plaquette state.

Full size image

The low-temperature transition coincides with the peak position of dn_v/dT also plotted in Fig. 4a. The maximum of dn_v/dT implies a rapid development of the vortex number n_v with increasing temperature, which is a characteristic feature of the Berezinskii–Kosterlitz–Thouless (BKT) transition for two-dimensional superconductivity^58,59,60. We thus identify this transition as the superconducting transition. The value of T_c is examined for other lattice size and found to vary only slightly, confirming the robustness of our qualitative conclusions.

The final phase diagram is already discussed in Fig. 1a, showing nonmonotonic variation of T_c/t with t/J and a maximum at the plaquette QCP. The overall evolution of T_c may be understood from the phase difference of the pairing fields on neighbouring bonds. Figure 4b plots the probabilistic distribution p(δθ₁) of (delta {theta }_{1}={theta }_{{boldsymbol{0}}}^{x/y}-{theta }_{(1,0)/(0,1)}^{x/y}). We find two symmetric peaks around zero in the plaquette state and a single peak in the uniform superconducting state. Interestingly, as shown in Fig. 4c, while the peak position (| delta {theta }_{1}^{max }|) decreases gradually and diminishes above t/J = 0.27, the inverse of its fluctuation, as well as that between next-nearest-neighbour bonds, also varies nonmontonically with t/J and exhibits a maximum near the plaquette QCP, in good correspondence with the evolution of T_c/t. This coincidence is unexpected at first glance but easy to understand, since a smaller fluctuation of δθ₁ around zero indicates a larger phase stiffness of the pairing fields on neighbouring bonds, thus favoring larger superfluid density and T_c. Theoretically, this is usually described by the free energy²⁹, (F=frac{{rho }_{{rm{s}}}}{2}{int}_{x}{(delta theta )}^{2}), such that the phase fluctuation 〈(δθ)²〉 is inversely related to the superfluid density ρ_s. This explains our observed correlation between the fluctuation of the phase difference and the magnitude of T_c in Fig. 4c.

Discussion on the plaquette state

The plaquette state may also have other exotic properties detectable in experiments. For example, pairing field modulation may affect local spin susceptibility^61,62 and cause some spin resonance mode⁶³. In fact, the plaquette state shares many similarities with the supersolid phase realized in dipolar cold atoms^{64,65,66,67,68}. Both break translational symmetry and U(1) phase symmetry at zero temperature. Similar to the plaquette state, the microscopic configurations of supersolid consist of weakly connected droplets. Both occupy an intermediate region of their respective phase diagram: the plaquette state occurs between the uniform superconductivity and a disordered phase of coexisting plaquettes and dimers for extremely large pairing interaction, while the supersolid exists between the superfluid phase and an incoherent droplet solid. Given these similarities, one may anticipate that vortices may exist in the supersolid phase, while two modes with different dispersions for some dynamic structure factor observed in supersolid⁶⁵ may also emerge in the plaquette state.

Though the Bose–Einstein condensation (BEC)⁶⁹ has traditionally been argued to be the strong coupling limit of the superconductivity, our results question its naive extension to unconventional superconductor that cannot be described by the local s-wave pairing with onsite attractive interaction. In particular, for unconventional superconductors with strong onsite Coulomb repulsion such as the cuprates, onsite s-wave pairing is generally unfavored and the pairs tend to occupy different sites. As a result, short-range pairing emerges for nearest-neighbour spin exchange interaction and, at strong coupling, may cause plaquette states that break the lattice translational symmetry. This differs from the local two-particle bound state typical of the BEC. On the other hand, our derived plaquette state does share some similarities with the BEC, such as the U-shaped density of states near the Fermi energy, the flat dispersion around k_x = 0, and the pseudogap in the normal state at high temperatures.

Our proposed plaquette state is different from the widely-studied PDW state^{46,47,70,71,72}, even though both exhibit real-space modulation of the pairing fields. While the PDW may generally lead to a charge density wave, the plaquette state breaks simultaneously the time-reversal symmetry and the translational symmetry and coexists with a bond charge order¹⁹. The PDW is by far only found experimentally in superconducting region^73,74,75 and might arise theoretically from the interplay of magnetism and superconductivity⁷⁶, while the plaquette state proposed here represents an intrinsic pairing instability at strong coupling.

Relevance to the cuprates

We note again that t is the effective hopping parameter of the quasiparticles, which should already take account of the Gutzwiller constraint. It is a small number proportional to the hole doping in underdoped cuprates, and reaches about 100−200 meV in overdoped cuprates as indicated by ARPES measurements⁷⁷. By contrast, the exchange interaction J is almost doping-independent and roughly 100−200 meV as revealed by RIXS experiments^24,25,26. Hence, our choice of the t/J range is reasonable according to these experiments. Numerically, our results cover the strong to weak coupling regions of the superconducting pairing, and provide useful information on the maximum T_c and the plaquette instability. Many of our findings are in good correspondence with experimental observations^{48,49,78,79,80}. In particular, the plaquette structure might be closely related to the 4 × 4 structure observed in recent STM measurements on underdoped cuprates^78,79, suggesting that Cooper pairs may first form within these local structures and only achieve long-range phase coherence as the doping reaches a certain threshold^49,78. For overdoped cuprates, a pseudogap feature, possibly driven by phase fluctuations, has also been observed within a narrow region above T_c⁷⁷. Such overall correspondence supports potential relevance of our work to the cuprate physics and highlights the importance of a strong-coupling real-space perspective in exploring high-temperature superconductivity.

Constraint on maximum T
_c/J

Another important observation of our calculations is that the superconductivity may intrinsically be suppressed for sufficiently strong pairing interaction even without considering competing orders from other channels. Thus, T_c is constrained from both sides of strong and weak pairing interactions. It is then sensible to study the ratio T_c/J to have a feeling about the maximum T_c allowed by the pairing interaction J^27,28,81. For the one-band square lattice model discussed so far, we find a maximum ratio T_c/J ≈ 0.04. We have also tested other parameters and find that tuning the next-nearest-neighbour hopping ({t}^{{prime} }) and the chemical potential μ can only slight improve this ratio. Specifically, at half-filling with ({t}^{{prime} }=0) and μ = 0 near the van Hove singularity, the maximum T_c/J may be enhanced to 0.045. Motivated by the possible importance of apical oxygen on T_c⁸², we have also studied a model with an extra conduction layer, and find the maximum T_c/J may be at most enhanced to about 0.06 for certain special (nearest-neighbour) interlayer hopping. On the other hand, local interlayer hopping is found to suppress this maximum ratio. Taking t ≈ 100−200 meV from the angle-resolved photoemission spectroscopy (ARPES) and the specific heat analysis^77,83 and J ≈ 100−190 meV from the resonant inelastic X-ray scattering measurement (RIXS)^81,84, these ratios yield the highest T_c to be 100–130 K, consistent with the reported ({T}_{c}^{max }=)97 K for single-layer and 135 K for multi-layer cuprate superconductors under ambient pressure^82,85.

To further explore the above idea, we extend our calculations to other variations of the minimum effective model, covering nearest-neighbour or onsite pairings, single or multi-layer structures, and intralayer or interlayer pairings (see Methods). It is important to note that our models do not depend on fine details of the microscopic pairing mechanism, as long as the effective pairing interaction and the low-energy Hamiltonians remain the same. Figure 5a shows the variations of T_c/J versus t/J in these models, where J is the local attractive Hubbard interaction for onsite pairing, and the interlayer superexchange interaction for interlayer pairing as discussed previously for La₃Ni₂O₇ under high pressure^39,86. We see all curves behave nonmonotonically with the pairing interaction, although they may have different strong-coupling limit (e.g. BEC for onsite pairing and preformed local interlayer pairing for bilayer nickelates), with the maximum T_c/J lying within the interval from 0.04 to 0.07. Notably, for the attractive Hubbard model, our simulations yield consistent results compared with previous quantum Monte Carlo simulations (open down-pointing triangles)², which reinforces the reliability of our approach, and introducing an additional conduction layer gives a similar maximum ratio^87,88. Note that we have ignored long distance pairing since it is typically weaker than onsite or nearest-neighbour ones for reaching the maximum T_c. All other instabilities are also ignored to maximize the pairing instability. Hence our phase diagram is not the full phase diagram with all possible ground states of a physical model, but a phase diagram that intentionally exaggerates the superconductivity and other possible instabilities in the pairing channel, so that the derived T_c/J could be a better estimate of its potential upper limit.

a Typical results of T_c/J as functions of t/J for several effective models with interlayer, intralayer onsite, or intralayer nn (nearest-neighbour) pairing interactions. For simplicity, only the nearest-neighbour hopping t is considered and the chemical potential is set to μ = 0. For nearest-neighbour pairing in a one-layer model (away from the van Hove singularity), introducing an extra layer of conduction electrons with nearest-neighbour interlayer hopping (t_p = 0.7t) is found to enhance the maximum T_c/J. Also compared is a typical result of the attractive Hubbard model from previous quantum Monte Carlo simulations (open down-pointing triangles) away from the half-filling². b Collection of experimental T_c/J ratios for a number of cuprate, nickelate, iron-based, and heavy fermion superconductors, where J are estimated from their respective spin interactions. The shaded area marks the region T_c/J ≈ 0.04−0.07. The large error bar exceeding this region comes from bulk FeSe as discussed in the main text. All error bars come from the experimental uncertainty of J as given by the original literatures listed in Table 1. The points circled by a square box mark those uncertain data from CePd₂Si₂, SmO_1−xF_xFeAs, LaO_1−xF_xFeAs, and La₃Ni₂O₇ discussed in the main text.

Full size image

To see if the above constraint may indeed apply in real materials, Fig. 5b and Table 1 collect the data for a number of well-known unconventional superconductors^{24,81,84,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110}. The spin energy scale in cuprate, iron-based, and nickelate superconductors have been determined mainly by the spin wave fitting in inelastic neutron scattering (INS) or RIXS experiments^{24,25,26,106,111}, where J has been found to vary only slightly with doping. It differs from the renormalized one due to the feedback effect observed in low-energy measurements by INS¹¹² and two-magnon extraction in Raman spectra¹¹³. In iron-based superconductors such as CaFe₂As₂, SrFe₂As₂, BaFe₂As₂, and NaFeAs, the ratios T_c/J are less than 0.063^{98,99,100,101,102,103,104,105,106,107}, where the value of J is extracted from the reported SJ by taking the effective spin size S = 0.69 for SrFe₂As₂ and S = 1/2 for all others except for FeSe following the literatures^89,90,91. The large error bar exceeding the shaded area in Fig. 5b comes from the bulk FeSe (T_c = 8 K), for which neutron scattering measurements reported the ratio T/J = 0.86 ± 0.35 at T = 110 K¹⁰⁷. Unfortunately, we do not find the data for FeSe films, whose high T_c might involve contributions from the interface. To the best of our knowledge, there is also no exact estimate of J for the 1111 systems. It has been reported that SmOFeAs adopts an intermediate spin dispersion between those of NaFeAs and BaFe₂As₂¹¹⁴. Assuming that the spin interaction is not sensitive to the doping, as observed in BaFe_2−xNi_xAs₂ and NaFe_1−xCo_xAs^106,111, we might roughly estimate J ~ 80−118.4 meV for SmO_1−xF_xFeAs and thus obtain a maximum ratio T_c/J ≈ 0.040−0.059 given its maximum T_c= 55 K¹¹⁵. While for LaO_1−xF_xFeAs, experiments only indicate an overall magnitude of SJ ~ 40 meV along different directions¹¹⁶, which yields T_c/J ~ 0.046 with its maximum T_c = 43 K using S = 1/2¹¹⁷. Both fall within our proposed range.

The infinite-layer nickelate superconductors have a small maximum ratio of about 0.026, which indicates the potential to reach a higher T_c^{92,93,94,95,96,97,118}. RIXS measurements¹¹⁹ on the high-pressure high-temperature bilayer nickelate superconductor La₃Ni₂O₇ reported an interlayer spin interaction strength (J) of about 140 meV assuming its spin size S = 1/2, which also seems to be confirmed by inelastic neutron measurements¹²⁰. Although these measurements were performed under ambient pressure, it gives a rough estimate of the magnitude of J. If we naively apply this value to the high pressure region where the superconductivity was reported with ({T}_{c}^{max }approx 80) K, we find ({T}_{c}^{max }/Japprox 0.05) for the bilayer nickelate superconductors, which agrees well with our previous Monte Carlo simulations³⁹. Recently, superconductivity has been reported also in the trilayer nickelate superconductor La₄Ni₃O₁₀ under high pressure, albeit with a much smaller ({T}_{c}^{max }approx 30) K¹²¹. It has been proposed theoretically that competition and frustration of interlayer pairing between the inner layer and two outer layers may lead to strong superconducting fluctuations and thus reduce the maximum ratio of T_c/J to 0.02−0.03⁴⁰. This, together with layer imbalance and the possibly smaller interlayer J, may explain the much reduced ({T}_{c}^{max }) in the trilayer nickelate compared to those in the bilayer ones.

By contrast, the cuprate high-temperature superconductors have the highest ({T}_{c}^{max }) in the trilayer structure, and their overall ({T}_{c}^{max }/J) ratios can reach up to 0.067, as observed in HgBa₂CaCu₂O_6+δ, YBa₂Cu₄O₈, YBa₂Cu₃O_6+δ, NdBa₂Cu₃O_6+δ, Tl₂Ba₂CuO_6+δ, and Bi_2+xSr_2−xCa₂Cu₃O_10+δ⁸¹. This opposite tendency reflects an intrinsic distinction in the pairing mechanisms between multilayer nickelate and cuprate superconductors^39,40,86,122. In heavy-fermion superconductors such as CeCoIn₅ or PuCoGa₅, systematic measurements of J are lacking. We therefore estimate the spin interaction energy from the coherence temperature scale, namely the Ruderman–Kittle–Kasuya–Yosida (RKKY) scale, and find the highest T_c/J to be about 0.046^108,109,110. To the best of our knowledge, a spin wave fitting has only been applied to CePd₂Si₂ and yields J = 0.61 meV under ambient pressure¹²³. Combining naively this value with its T_c = 0.43 K at 3 GPa gives the ratio T_c/J ≈ 0.061, in good alignment with our suggested constraint.

Despite vast complexities across all these different families of unconventional superconductors that are far beyond our simplified models, none of their maximum T_c/J ratios exceeds the proposed range 0.04−0.07, suggesting that our calculations indeed capture some essence of the fundamental physics of unconventional superconductivity. While this range of the maximum T_c/J ratio is also supported by data collection in previous studies^24,124,125, one should keep in mind that our calculations only cover a small part of the parameter and model space, and it is not completely clear how other complex factors, such as the dimensionality and multi-orbitals, might affect the ratio. The maximum ratio proposed here seems to more represent a practical upper limit for some quite generic situations in real superconductors. This consistency urges for a more rigorous theoretical understanding.

We finally comment on the T_c/t ratio widely used in previous literatures. Unlike T_c/J, we find the maximum T_c/t depends more sensitively on models and may reach 0.29, 0.15, 0.105 upon tuning the hopping parameters or the chemical potential for interlayer, intralayer onsite, intralayer nearest-neighbour (nn) pairings, respectively. For the attractive Hubbard model, our derived maximum T_c/t ≈ 0.15 is close to the quantum Monte Carlo result of 0.17, which also confirms the validity of our estimate². In general, our calculations across these different models yield the maximum T_c/t ratio within the range of 0.1−0.3 at an optimal t/J value of around 0.2, somewhat different from those for the maximum T_c/J ratio. The relatively lager variation of the T_c/t ratio may be ascribed to the fact that the long-range phase coherence determining T_c relies heavily on the cooperative hopping of paired electrons and may hence differ greatly for different pairing configurations and lattice geometries beyond the simple hopping parameters. The quasiparticle hopping is also strongly renormalized by correlation effects, which makes it difficult to estimate in practice. It is for these reasons that we have chosen to treat t as a tuning parameter and focus on the T_c/J ratio that might be better compared with experiment.

Route to room temperature superconductivity?

It is important to emphasize again that the above agreement by no means implies that all these superconductors, including hole-doped cuprates, are fully described by the specified pairing mechanisms in our simplified models. There is also no rigorous theoretical proof for a maximum T_c in unconventional superconductors⁴. Nevertheless, if we take the above constraint seriously, achieving room temperature superconductivity seems unlikely under ambient pressure within the current theoretical framework. For T_c to reach 300 K, we need a pairing interaction of the order 400−700 meV, which is twice higher than the spin exchange interaction in cuprates and seems unrealistic based on our survey of existing correlated materials. Moreover, the maximum T_c/J is only realized at an optimal ratio of t/J, thus also requiring a larger quasiparticle hopping t, a situation that seems to only occur under pressure. This is contrary to the weak-coupling BCS theory which predicts a higher T_c for a larger density of states (smaller t), while in the strong coupling limit, the pairing strength is sufficiently large and one requires sufficient kinetic energy, hence a larger t, to achieve the phase coherence. For instance, in the attractive Hubbard model with an infinite U, T_c is determined by a small fraction of t. Consequently, for room-temperature superconductivity, t must reach the order of hundred meV, which is not favoured in flat-band systems¹²⁶.

It is therefore imperative to explore alternative avenues to enhance the ratio under ambient pressure. It has been noticed that three-layer cuprate superconductors have the highest T_c. One may therefore speculate that multi-layer may promote T_c. Indeed, the maximum T_c increases from 97 K in the single-layer HgBa₂CuO_4+δ to 127 K in the two-layer HgBa₂CaCu₂O_6+δ and 135 K in the three-layer HgBa₂Ca₂Cu₃O_9+δ^82,85. However, the ratio T_c/J ≈ 0.062 seems to remain unchanged and the increase of T_c seems to come purely from the increase of J⁸¹. On the other hand, the maximum T_c/J does increase from 0.021 in the single-layer Bi₂Sr_2−xLa_xCuO_6+δ to 0.058 in the three-layer Bi_2+xSr_2−xCa₂Cu₃O_10+δ in Bi-systems⁸¹, but the latter still lies within our proposed range, implying that increasing the number of layers from Bi2201 to Bi2223 only helps to tune the optimal conditions for maximizing T_c/J, while the constraint itself is not touched. We have also examined the effect of additional local interlayer hopping and find that it actually reduces the maximum T_c/J. Additionally, one may follow the studies of FeSe films^127,128 and consider to improve T_c by introducing phonons, but this seems empirically at most to provide an increase of around 40 K, given the limited characteristic phonon frequencies under ambient pressure^129,130. A larger spin interaction occurs for the Hund’s rule coupling inside an atom. However, it is not clear if intra-atomic inter-orbital pairing may support a high T_c due to their very different orbital characters of the paired electrons.

Taken together, the known unconventional superconductor families seem to have almost exhausted their potentials in reaching the highest T_c allowed by their spin exchange interactions. As a result, room-temperature superconductivity at ambient pressure seems unlikely to arise from a single pairing mechanism within the current theoretical framework, unless one could find a way to substantially enhance the exchange interaction. Our results may not only help rule out some evidently wrong directions^131,132, but also point out the necessity of exploring alternative approaches to achieve room-temperature superconductors at ambient pressure^{1,86,88,110,133,134,135,136,137,138,139,140,141}. It encourages the possibility of incorporating different pairing mechanisms^{8,142,143,144,145,146}, including but not limited to magnetic, charge, orbital fluctuations, or excitons, bipolarons, etc, to improve the overall effective pairing interaction, for which FeSe films may be a good example^147,148,149. Our derived ratio provides a tentative guide for future material exploration of novel high-temperature superconductors. Theoretically, by utilising J from newly developed methods¹⁵⁰ and effective hopping t from strongly correlated calculations¹⁵¹, an approximate estimate of the upper limit of T_c may be predicted for the selection of promising candidates. Experimentally, estimating J from RIXS, INS, or other state-of-the-art techniques in newly discovered materials may also help identify their potential in reaching the desired T_c. Last but not least, understanding unconventional superconductivity from a real-space, strong-coupling perspective may already provide an operational and more practical avenue for material design compared to the momentum-space, weak-coupling approach.

Pairing interaction

For onsite pairing in the attractive Hubbard model, we have immediately

where ψ_i = d_i↓d_i↑ and J = U.

For the Hubbard model with a large Coulomb repulsion U, the pairing interaction J is given by the exchange interaction between nearest-neighbour sites. To see this, we may follow the standard derivation to first project out the double occupancy and map the Hubbard model to an effective low-energy model with the following interaction term:

where ({{boldsymbol{s}}}_{i}={sum }_{alpha ,beta }{d}_{ialpha }^{dagger }frac{{{boldsymbol{sigma }}}_{alpha beta }}{2}{d}_{ibeta }) is the spin density and ({n}_{i}={sum }_{sigma }{d}_{isigma }^{dagger }{d}_{isigma }) is the charge density.

To study the superconductivity, we introduce the spin-singlet (({psi }_{ij}^{{rm{S}}})) and spin-triplet (({{boldsymbol{psi }}}_{ij}^{{rm{T}}})) operators:

where ({psi }_{ij}^{{rm{S}}}) and ({{boldsymbol{psi }}}_{ij}^{{rm{T}}}) satisfy ({psi }_{ji}^{{rm{S}}}={psi }_{ij}^{{rm{S}}}) and ({{boldsymbol{psi }}}_{ji}^{{rm{T}}}=-{{boldsymbol{psi }}}_{ij}^{{rm{T}}}).

The above interaction can then be rewritten as:

where the spin-triplet term cancels. This is the Hamiltonian (with J = J_ex) used in our simulations for the d-wave superconductivity. For cuprates, J is given by the superexchange mechanism. Note that to maximize T_c, we have included both the nearest-neighbour antiferromagnetic spin-density interaction, Js_i ⋅ s_j, and the nearest-neighbour attractive charge-density interaction, −V^cn_in_j. While only the former is included in many works, both seem to be supported by recent experiments^152,153. For spin-fluctuation interaction^8,81, V(q)s_q ⋅ s_−q, the dominant contribution has a similar form in real space. These justify our choice of the phenomenological pairing term.

Similarly, we have for interlayer pairing:

where ({{boldsymbol{s}}}_{ai}={sum }_{alpha ,beta }{d}_{aialpha }^{dagger }frac{{{boldsymbol{sigma }}}_{alpha ,beta }}{2}{d}_{ajbeta }), ({n}_{i}={sum }_{sigma }{d}_{aisigma }^{dagger }{d}_{aisigma }), and ({psi }_{12i}=frac{1}{sqrt{2}}({d}_{1idownarrow }{d}_{2iuparrow }-{d}_{1iuparrow }{d}_{2idownarrow })). Again, the pairing interaction is J = J_ex.

Mutual information and vortex number

The superconducting transition temperature is always determined by the superconducting phase coherence from the phase mutual information defined as^38,39:

where (p({theta }_{{boldsymbol{0}}}^{x/y}),,p({theta }_{{boldsymbol{R}}}^{x/y})) is the marginal distribution of the pairing field phase on two bonds with a distance R, and p(θ₀, θ_R) is their joint probabilistic distribution. For onsite or interlayer pairing, ({theta }_{{boldsymbol{R}}}^{x/y}) simplifies to θ_R.

The vortex number is calculated using

where w_i is the winding number for ({theta }_{i}to {theta }_{i+hat{x}}to {theta }_{i+hat{x}+hat{y}}to {theta }_{i+hat{y}}to {theta }_{i}) with the phase θ_i of ({Delta }_{i}=({Delta }_{i,i+hat{x}}+{Delta }_{i,i-hat{x}}-{Delta }_{i,i+hat{y}}-{Delta }_{i,i-hat{y}})/4) for nearest-neighbour pairing and 〈〉 denotes the statistic average over all pairing configurations. For onsite or interlayer pairing, θ_i is the phase of the pairing field at site i.

We find that T_c determined from the phase mutual information agrees well with that estimated from the numerical derivatives of the vortex number (frac{d{n}_{{rm{v}}}}{dT}) as well as the BKT transition temperature obtained from the superfluid stiffness¹⁵⁴, which confirms the validity of our method in calculating T_c for these models. However, the superfluid stiffness is computationally more expensive.

Time-reversal symmetry breaking and coexisting bond charge order of the plaquette state

Figure 6a, b plots the probabilistic distributions of the phase difference δθ_xy along x and y directions for t/J = 0.20 and 0.15, where the plaquette state develops at low temperatures. We see two-peak structures below roughly T = 0.06 for t/J = 0.20 and T = 0.09 for t/J = 0.15, both of which coincide with T_□ determined in Fig. 1a from the distribution of the pairing amplitude and other quantities. As discussed in the main text, the deviation from δθ_xy = π marks the time-reversal symmetry breaking. Figure 6c, d further compare the distributions of the site charge density ({n}_{i}={sum }_{sigma }langle {d}_{isigma }^{dagger }{d}_{isigma }rangle) and the bond charge density ({n}_{ij}={sum }_{sigma }langle {d}_{isigma }^{dagger }{d}_{jsigma }+{d}_{jsigma }^{dagger }{d}_{isigma }rangle) between two nearest-neighbour sites. While the former always exhibits a single peak, the latter also develops a two-peak structure below T_□. These indicate a uniform charge distribution on all sites but spatial modulation of the bond charge density. We thus conclude that the plaquette state below T_□ breaks simultaneously the time-reversal symmetry and the translational symmetry and coexists with a bond charge order.

a, b Probabilistic distributions of the phase difference δθ_xy along x and y direction for t/J = 0.20 and 0.15, respectively. c, d Comparison of the distribution of the site electron density n_i and the bond charge density n_ij for t/J = 0.15.

Full size image

Other models

To derive the T_c/J constraint, we extend the simplest one-band model to the following variations:

(1) A two-layer model with intralayer nearest-neighbour pairing and interlayer hopping:

where the subscript a = 1, 2 represents the layer index, ({psi }_{aij}=frac{1}{sqrt{2}}({d}_{aidownarrow }{d}_{ajuparrow }-{d}_{aiuparrow }{d}_{ajdownarrow })), and t_p denotes the local interlayer hopping.

(2) A two-layer model with an extra conduction layer motivated by the possible importance of apex oxygens in cuprates:

where t_p denotes local (i = j) or nearest-neighbour ((j=ipm hat{x}) or (ipm hat{y})) interlayer hopping.

(3) A single-layer model with onsite pairing interaction as in the attractive Hubbard model:

where ψ_i = d_i↓d_i↑ and J is given by the local attractive Hubbard interaction.

(4) A two-layer model with onsite pairing in one layer and an extra conduction layer:

where ψ_i = d_i↓d_i↑ and t_p denotes local (i = j) or nearest-neighbour ((j=ipm hat{x}) or (ipm hat{y})) interlayer hopping.

(5) A two-layer model with interlayer pairing:

where ({psi }_{12i}=frac{1}{sqrt{2}}({d}_{1idownarrow }{d}_{2iuparrow }-{d}_{1iuparrow }{d}_{2idownarrow })).

Models (1) and (2) are constructed to reflect the effects of interlayer hopping and apical oxygen in cuprate superconductors. Model (2) may also be applied to the infinite-layer nickelate superconductors. Models (3) and (4) deal with onsite pairing with local attractive interaction. And model (5) is motivated by the bilayer nickelate superconductor.