Bypassing the lattice BCS–BEC crossover in strongly correlated superconductors through multiorbital physics

The collective and phase coherent condensation of electrons into bound Cooper pairs leads to the emergence of superconductivity. This macroscopic coherence enables dissipationless charge currents, perfect diamagnetism, fluxoid quantization, and technical applications^1,2 ranging from electromagnets in particle accelerators to quantum computing hardware. Often, superconducting (SC) functionality is controlled by the critical surface spanned by critical magnetic fields, currents, and temperatures which a SC condensate can tolerate. Fundamentally, these are determined by the characteristic length scales of a superconductor – the London penetration depth, λ_L, and the coherence length, ξ₀.

λ_L and ξ₀ quantify different aspects of the SC condensate: The penetration depth is the length associated with the mass term that the vector potential gains through the Anderson-Higgs mechanism³. In consequence, magnetic fields decay exponentially over a distance λ_L inside a superconductor. Through this, λ_L is connected to the energy cost of order parameter (OP) phase variations and hence the SC stiffness D_s. The coherence length, on the other hand, is the intrinsic length scale of OP amplitude variations and is associated with the amplitude Higgs mode. ξ₀ sets the scale below which amplitude and phase modes significantly couple such that spatial variations of the OP’s phase reduce its amplitude^3,4.

In addition to influencing the macroscopic properties of superconductors, λ_L and ξ₀ play an important role to understand strongly correlated superconductors, as is epitomized in the Uemura plot^5,6,7,8. For instance, the interplay of λ_L and ξ₀ impacts critical temperatures⁹, it is relevant for the pseudogap formation^10,11,12,13, it influences magneto-thermal transport properties like the Nernst effect^8,14,15, and it might underlie the light-enhancement of superconductivity^{15,16,17,18,19}. An important concept in this context is the BCS–BEC crossover phenomenology^{20,21,22,23,24}. It continuously connects the two limiting cases of weak-coupling Bardeen–Cooper–Schrieffer (BCS) superconductivity with weakly-bound and largely overlapping Cooper pairs to tightly-bound molecule-like pairs in the strong-coupling Bose–Einstein condensate (BEC) as the interaction strength or the density is varied (Fig. 1).

Evolution of the critical temperature T_c and Cooper pair size ξ_p in the BCS–BEC crossover (dark blue line) for (a) Fermi gases and (b) lattice systems. Both cases display a dome-shaped behavior of T_c in the intermediate crossover regime but behave qualitatively different in the strong coupling BEC phase: T_c remains finite in the Fermi gas (a)³¹ but approaches zero in the lattice case (b). During the crossover, ξ_p is reduced to the order of interparticle spacing r (a, b) and lattice constant a (b). For the lattice system, we contrast the evolution towards the multiorbital strong coupling phase (light blue line) discussed in this article. Here, the localization of pairs differently affects the bosonic hopping t_B, as drawn in the insets on the right. In the lattice BEC limit, bosonic hopping relies on a second-order process involving two fermionic hoppings, t_F, and a virtual intermediate state with broken pairs that is inhibited by the strong attraction g < 0. Consequently, this process and also ({T}_{{rm{c}}}propto {t}_{{rm{B}}}={t}_{{rm{F}}}^{2}/| g|) are quenched at large ∣g∣^21,24. In the multiorbital strong coupling case, the local coexistence of paired and unpaired electrons fluctuating between different orbitals enables bosonic hopping t_B as a first order process in t_F without any intermediate broken pair states. See Supplementary Note 7 for details on the involved states and energy diagram specific to the case of A₃C₆₀. A second temperature scale T* is drawn in both panels as its splitting from T_c (corresponding to the opening of a pseudogap) marks the beginning of the crossover regime²⁴.

Full size image

The BCS–BEC crossover has been studied in ultracold Fermi gases²³, low-density doped semiconductors²⁵, and is under debate for several unconventional superconductors^{6,8,24,26,27,28,29,30}. However, quasi-continuous systems, for instance Fermi gases, and strongly correlated superconducting solids show a crucially different behavior of how their SC properties, most importantly the critical temperature T_c, change towards the strong coupling BEC limit. While T_c converges to a constant temperature T_BEC for Fermi gases in a continuum³¹ (Fig. 1a), T_c can become arbitrarily small in strongly correlated lattice systems due to the quenching of kinetic energy (Fig. 1b). Since the movement of electron pairs, i.e., bosonic hopping, necessitates intermediate fermionic hopping, it becomes increasingly unfavorable for strong attractions. Thus, as Cooper pairs become localized on the scale of the lattice constant, the condensate’s stiffness and hence T_c are compromised^21,24. Figure 1 contrasts this generic BCS–BEC crossover picture for Fermi gases and correlated lattice systems in terms of the change of T_c and the pair size ξ_p as a function of pairing strength. Due to the decrease of T_c in the BCS and BEC limits, a prominent dome-shape of T_c can be expected in the crossover region for solid materials. Because of this, recent experimental efforts to increase T_c concentrate on stabilizing SC materials in this region^24,25,27,30.

In this work, we demonstrate how multiorbital effects^32,33 can enhance superconductivity beyond the expectations of the lattice BCS–BEC crossover phenomenology (as contrasted in Fig. 1b) with a model inspired by alkali-doped fullerides (A₃C₆₀ with A = K, Rb, Cs). The material family of A₃C₆₀ hosts exotic s-wave superconductivity of critical temperatures up to T_c = 38 K, being the highest temperatures among molecular superconductors^34,35,36,37, and they possibly reach photo-induced SC at even higher temperatures^15,17,18,38. In order to theoretically characterize the SC state, the knowledge of the intrinsic SC length scales is essential. While BCS theory and Eliashberg theory provide a microscopic description of λ_L and ξ₀ for weakly correlated materials^39,40,41, their validity is unclear for superconductors with strong electron correlations. To the best of our knowledge, ξ₀ is generally not known from theory in strongly correlated materials. Only approaches to determine λ_L exist where an approximate, microscopic assessment of the SC stiffness from locally exact theories has been established, albeit neglecting vertex corrections^{42,43,44,45,46,47}.

For this reason, we introduce a novel theoretical framework to microscopically access λ_L and ξ₀ from the tolerance of SC pairing to spatial OP variations. Central to our approach are calculations in the superconducting state under a constraint of finite-momentum pairing (FMP). Via Nambu–Gor’kov Green functions we get direct access to the superconducting OP and the depairing current j_dp which in turn yield ξ₀ and λ_L. The FMP constraint is the SC analog to planar spin spirals applied to magnetic systems^48,49,50. As in magnetism, a generalized Bloch theorem holds that allows us to consider FMP without supercells; see Supplementary Note 2 for a proof. As a result, our approach can be easily embedded in microscopic theories and ab initio approaches to tackle material-realistic calculations.

In this work, we implement FMP in Dynamical Mean-Field theory (DMFT)⁵¹ to treat strongly correlated superconductivity. In DMFT, the interacting many-body problem is solved by self-consistently mapping the lattice model onto a local impurity problem. By this, local correlations are treated exactly. We apply the FMP-constrained DMFT to A₃C₆₀ where we find ξ₀ and λ_L in line with experiment for model parameter ranges derived from ab initio estimates, validating our approach. For enhanced pairing interaction, we then reveal a multiorbital strong coupling SC state with minimal ξ₀ on the order of only 2 – 3 lattice constants, but with robust stiffness D_s and high T_c which increases with the pair interaction strength without a dome shape. This strong coupling SC state is distinct to the lattice BCS–BEC crossover phenomenology showing promising routes to optimize superconductors with higher critical temperatures. We discuss this possibility in-depth after the presentation of our results.

The remaining paper is organized as follows: First, we motivate how the FMP constraint is linked to λ_L and ξ₀ from phenomenological Ginzburg–Landau theory after which we summarize the microscopic approach from Green function methods. Technical details are available in the Supplementary Information as cross-referenced at relevant points. Subsequently, we discuss our results for the multiorbital model of A₃C₆₀. Readers primarily interested in the analysis of these results might skip the first part on the physical motivation for obtaining superconducting length scales from the FMP constraint.

Extracting superconducting length scales from the constraint of finite-momentum pairing

In most SC materials, Cooper pairs do not carry a finite center-of-mass momentum q = 0. Yet, in presence of external fields, coexisting magnetism, or even spontaneously SC states with FMP, i.e., q ≠ 0, might arise^{52,53,54,55,56} as originally conjectured in Fulde–Ferrel–Larkin–Ovchinnikov (FFLO) theory^57,58,59. Here, we introduce a method to access the characteristic SC length scales ξ₀ and λ_L of strongly correlated materials through the calculation of FMP states with a manually constrained pair center-of-mass momentum q^{4,40,60,61,62}. For this, we enforce the OP to be of the FMP form Ψ_q(r) = ∣Ψ_q∣e^iqr corresponding to FF-type pairing⁵⁷, which is to be differentiated from pair density waves with amplitude modulations^58,59,63. We contrast the OP for zero and finite momentum in real and momentum space in the top panel of Fig. 2. For FMP, the OP’s phase is a helix winding around the direction of q, while the OP for zero-momentum pairing is simply a constant.

The top panel insets sketch the position and momentum space representation of the order parameter (OP) Ψ_q(r) = ∣Ψ_q∣e^iqr in the zero-momentum (left, q = 0) and finite-momentum pairing states (right, q > 0). The main panels show (a) the momentum dependence of the OP modulus and (b) the supercurrent density j_q = ∣j_q∣ as function of Cooper pair momentum q = ∣q∣ in reciprocal lattice units (r.l.u.). Gray lines indicate the points of extracting ξ and j_dp (see text). The data shown are results for the A₃C₆₀ model (c.f. Eq. (7)) with lattice constant a and interaction parameters U/W = 1.4, J/W = − 0.04 evaluated at different temperatures T (color coded; see color bar with white triangular markers).

Full size image

Before turning to the implementation in microscopic approaches, we motivate how the FMP constraint relates to λ_L and ξ₀. The Ginzburg–Landau (GL) framework provides an intuitive picture to this connection^4,40,61 which we summarize here and discuss in detail in Supplementary Note 1, including a derivation of the following equations. The GL low-order expansion of the free energy density f_GL in terms of the FMP-constrained OP reads

with q = ∣q∣. α, b, and m* are the material and temperature dependent GL parameters. Here, the temperature-dependent (GL) correlation length ξ(T) appears as the natural length scale of the amplitude mode (∝∣α∣) and the kinetic energy term

with its zero-temperature value ξ₀ being the coherence length⁴. The stationary point of Eq. (1) shows that the q-dependent OP amplitude (leftvert {Psi }_{{boldsymbol{q}}}rightvert =leftvert {Psi }_{0}rightvert sqrt{1-{xi }^{2}{q}^{2}}) (T-dependence suppressed) decreases with increasing momentum q. For large enough q, SC order breaks down (i.e., Ψ_q → 0) as the kinetic energy from phase modulations becomes comparable to the gain in energy from pairing. The length scale associated with this breakdown is ξ(T) and can, therefore, be inferred from the q-dependent OP suppression. We employ, here, the criterion (xi (T)=1/(sqrt{2}| {boldsymbol{Q}}|)) with Q such that (| {Psi }_{{boldsymbol{Q}}}(T)/{Psi }_{0}(T)| =1/sqrt{2}) for fixed T (see Supplementary Note 5-A for more information).

The finite center-of-mass momentum of the Cooper pairs entails a charge supercurrent j_q ∝ ∣Ψ_q∣²q, c.f. Eq. (8) in Supplementary Note 1. This current density is a non-monotonous function of q with a maximum called depairing current density j_dp. It provides a theoretical upper bound to critical current densities, j_c, measured in experiment. We note that careful design of SC samples is necessary for j_c reaching j_dp as its value crucially depends on sample geometry and defect densities^40,61,64. As the current is related to the London penetration depth λ_L through the second London equation ({mu }_{0}{boldsymbol{j}}=-{lambda }_{{rm{L}}}^{-2}{boldsymbol{A}}), we can determine λ_L from j_dp in GL theory via

with the magnetic flux quantum Φ₀ = h/2e. The temperature dependence with the quartic power stated here is empirical^39,40 and we find that it describes our calculations better than the linearized GL expectation as discussed in Supplementary Note 5-B. Note that taking the q → 0 limit in our approach is related to linear-response-based methods to calculate λ_L or, equivalently, the stiffness D_s^47,62.

The GL analysis shows that the OP suppression and supercurrent induced by the FMP constraint connect to ξ₀ and λ_L. In a microscopic description, we acquire the OP and supercurrent density from the Nambu–Gor’kov Green function

where ({psi }_{{boldsymbol{k}},{boldsymbol{q}}}^{dagger }=({c}_{{boldsymbol{k}}+frac{{boldsymbol{q}}}{2}uparrow }^{dagger },,{c}_{-{boldsymbol{k}}+frac{{boldsymbol{q}}}{2}downarrow })) (orbital indices suppressed) are Nambu spinors that carry an additional dependence on q due to the FMP constraint. G (F) denotes the normal (anomalous) Green function component for electrons (G) and holes ((bar{G})) in imaginary time τ. For s-wave superconductivity as in A₃C₆₀^65,66,67,68, we use the orbital-diagonal, local anomalous Green function as the OP

which is the same for all orbitals α. This allows us to work with a single-component OP. The current density can be calculated via (c.f. Eq. (26) in the Supplementary Information)

where (hslash underline{{boldsymbol{v}}}={nabla }_{{boldsymbol{k}}}underline{h}({boldsymbol{k}})) is the group velocity obtained from the one-electron Hamiltonian (underline{h}({boldsymbol{k}})) and the trace runs over the orbital indices of (underline{{boldsymbol{v}}}) and ({underline{G}}_{{boldsymbol{q}}}). Underlined quantities indicate matrices in orbital space, N_k is the number of momentum points, and e is the elementary charge. See Methods and Supplementary Note 4 for details on the DMFT-based implementation and Supplementary Notes 3 for a derivation and discussion of Eq. (6).

The bottom panels of Fig. 2 show an example of our DMFT calculations which illustrates the q-dependence of the OP amplitude and current density for different temperatures T. Throughout the paper, we choose q parallel to a reciprocal lattice vector q = qb₁. We find a monotonous suppression of the OP with increasing q. The supercurrent initially grows linearly with q, reaches its maximum j_dp and then collapses upon further increase of q. Thus, both ∣Ψ_q∣ and j_q behave qualitatively as expected from the GL description. For decreasing temperature, the point where the OP gets significantly suppressed moves towards larger momenta q (smaller length scales ξ(T)), while j_dp(T) increases. We indicate the points where we extract ξ(T) and j_dp(T) with gray circles connected by dashed lines.

Superconducting coherence in alkali-doped fullerides

We apply FMP superconductivity to study a degenerate three-orbital model H = H_kin + H_int, where H_kin is the kinetic energy and the electron-electron interaction is described by a local Kanamori–Hubbard interaction⁶⁹

with total number (hat{N}), spin (hat{{boldsymbol{S}}}), and angular momentum operator (hat{{boldsymbol{L}}}). The independent interactions are the intraorbital Hubbard term U and Hund exchange J, as we use the SU(2) × SO(3) symmetric parametrization. This model is often discussed in the context of Hund’s metal physics relevant to, e.g., transition-metal oxides like ruthenates with partially filled t_2g shells^69,70. In the special case of negative exchange energy J < 0, it has been introduced to explain superconductivity in A₃C₆₀ materials^65,66,67,68. In fullerides, exotic s-wave superconductivity exists in proximity to a Mott-insulating (MI)^35,36,71 and a Jahn–Teller metallic phase^36,72,73,74. The influence of strong correlation effects and inverted Hund’s coupling were shown to be essential for the SC pairing^{65,66,67,68,75} utilizing orbital fluctuations⁷⁶ in a Suhl–Kondo mechanism⁷². In the weakly correlated regime of U → 0 and small J, the system follows BCS phenomenology⁶⁶.

The inversion of J seems unusual from the standpoint of atomic physics, where it dictates the filling of atomic shells via Hund’s rules. In A₃C₆₀, a negative J is induced by the electronic system coupling to intramolecular Jahn–Teller phonon modes^68,75,77,78. As a result, Hund’s rules are inverted such that states which minimize first spin S and then angular momentum L are energetically most favorable, see the second term of Eq. (7) for J < 0.

We connect the model to A₃C₆₀ by using an ab initio derived model for the kinetic energy⁷⁹

where t_αγ(R_ij) is the hopping amplitude between half-filled t_1u orbitals (α, γ = 1, 2, 3) of C₆₀ molecules on sites connected by lattice vector R_ij. We take the bandwidth W as the unit of energy (W ≈ 0.3 – 0.5 eV for Cs to K based A₃C₆₀)^68,79, see Supplemental Note 4-D for further details. For the interaction, we take the first principles’ estimates of fixed J/W = − 0.04 and varying U/W^{65,66,67,68,73} to emulate unit cell volumes as resulting from the size of different alkali dopants. We solve this Hamiltonian using DMFT, which explicitly takes into account superconducting order. Through the momentum dependence of ∣Ψ_q(T)∣ and j_q(T), ξ(T) and λ_L(T) can be extracted as discussed in the previous subsection. We show the derived temperature dependence of ∣Ψ_q=0(T)∣, ξ(T), and λ_L(T) for different U/W in Fig. 3.

The temperature dependence of (a) the zero-momentum order parameter ∣Ψ₀∣, (b) the correlation length ξ, and (c) the London penetration depth λ_L are shown for different ratios U/W. The data were obtained for fixed J/W = − 0.04 as estimated from ab initio data. Length scales are given in units of the lattice constant a. Fits to extract critical temperatures T_c and zero-temperature values ξ₀ and λ_L,0 are shown by dashed lines. The region of uncertainty to fit T_c is indicated by shaded regions.

Full size image

Close to the transition point, the OP vanishes and the critical temperature T_c can be extracted from ∣Ψ₀(T)∣² ∝ T − T_c. We find that T_c increases with U, contrary to the expectation that a repulsive interaction should be detrimental to electron pairing. This behavior is well understood in the picture of strongly correlated superconductivity: As the correlations quench the mobility of carriers, the effective pairing interaction ~ J/ZW increases due to a reduction of the quasiparticle weight Z^65,66,67. The trend of increasing T_c is broken by a first-order SC to MI phase transition for critical U ~ 2 W which is indicated by a dotted line. Upon approaching the MI phase, the magnitude of ∣Ψ₀∣ behaves in a dome-like shape. The T_c values of 0.8 – 1.4 × 10⁻² W that we obtain from DMFT correspond to 49 – 85 K which is on the order of but quantitatively higher than the experimentally observed values. The reason for this is that we approximate the interaction to be instantaneous as well as that we neglect disorder effects and non-local correlations which reduce T_c^67,68. The same effects could mitigate the dominance of the MI phase which prevents us from observing the experimental T_c -dome^36,68.

Turning to the correlation length, we observe that, away from T_c, ξ(T) is strongly reduced to only a few lattice constants (a ~ 14.2 – 14.5 Å) by increasing U, i.e., pairing becomes very localized. At the same time, λ_L(T) is enlarged. Hence, the condensate becomes much softer as there is a reduction of the SC stiffness ({D}_{{rm{s}}}propto {lambda }_{{rm{L}}}^{-2}) upon increasing U. Fitting Eqs. (2) and (3) to our data, we obtain zero-temperature values of ξ₀ = 3 – 10 nm and λ_L,0 = 80 – 120 nm. Comparing our results with experimental values of ξ₀ ~ 2 – 4.5 nm and λ_L ~ 200 – 800 nm^6,7,37,80,81, we see an almost quantitative match for ξ₀ and a qualitative match for λ_L. Both experiment and theory consistently classify A₃C₆₀ as type II superconductors (λ_L ≫ ξ₀)^37,81.

We speculate that disorder and spontaneous orbital-symmetry breaking⁷⁴ in the vicinity of the Mott state could lead to a further reduction of ξ₀ as well as an increase of λ_L beyond what is found here for the pure system. This could bring our calculations with minimal ξ₀ =3 nm closer to the experimental minimal coherence length of 2 nm revealed by measurements of large upper critical fields reaching up to a maximal H_c2 = 90 T in Ref. ⁸⁰ using H_c2 = Φ₀/(2πξ₀) with the flux quant Φ₀.

Circumvention of the lattice BCS–BEC crossover upon boosting inverted Hund’s coupling

The inverted Hund’s coupling is crucial for superconductivity in A₃C₆₀. This premise motivates us to explore in Fig. 4 the nature of the SC state in the interaction (U, J) phase space for J < 0 beyond ab initio estimates.

a Critical temperature T_c, zero temperature penetration depth λ_L,0, inverse square root of quasiparticle weight Z, coherence length ξ₀, and the statistical weight ({rho }_{leftvert {phi }_{0}rightrangle }^{N}) of the local lowest energy states (leftvert {phi }_{0}rightrangle) of the N = 2, 3, 4 particle sectors obeying inverted Hund’s rules as a function of U and J. Gray dots show original data points used for interpolation and the dashed line indicates a region where the proximity to the Mott state leads to a suppression of the superconducting state. There is no data point at the charge degeneracy line U = 2∣J∣ in the lower left corner as marked by black ‘canceling lines. b Distribution of statistical weights ({rho }_{leftvert {phi }_{n}rightrangle }^{N}) at four different interaction values U and J. (See Supplementary Note 7 for a listing of the eigenstates (leftvert {phi }_{n}rightrangle) and their respective local eigenenergies ({E}_{n}^{N})). Red (blue) bars denote the density matrix weight of ground states in the N = 3 (N = 2, 4) particle sector, the sum of which is plotted in the last panel of (a). c Exemplary depiction of representative lowest inverted Hund’s rule eigenstates. A delocalized doublon (electron pair) fluctuates between different orbitals due to correlated pair hopping J.

Full size image

As long as ∣J∣ < U/2, we find that strengthening the negative Hund’s coupling enhances the SC critical temperature with an increase up to T_c ≈ 5 × 10⁻² W, i.e., by a factor of seven compared to the ab initio motivated case of J/W = − 0.04. There is, however, a change in the role that U plays in the formation of superconductivity. While U was supportive for small magnitudes ∣J∣ ≲ 0.05 W, it increasingly becomes unfavorable for ∣J∣ > 0.05 W where T_c is reduced with increasing U. The effect is largest close to the MI phase where superconductivity is strongly suppressed. We indicate this proximity region by a dashed line (c.f. Supplementary Note 5-C).

The impact of U on the SC state can be understood from the U-dependence of the London penetration depth: λ_L grows monotonously with U and reaches its maximum close to the MI phase. Hence, the condensate is softest in the region where Mott physics is important and it becomes stiffer at smaller U. We find that this fits to the behavior of the effective bandwidth W_eff = ZW ∝ D_s where the quasiparticle weight Z is suppressed upon approaching the Mott phase. The behavior of Z shown in Fig. 4a confirms the qualitative connection ({lambda }_{{rm{L}}}propto {D}_{{rm{s}}}^{-1/2}propto 1/sqrt{Z}) for ∣J∣ > 0.05 W. The J-dependence of λ_L is much weaker than the U-dependence, as can be seen in Fig. 4a and the corresponding line-cuts in Fig. 5.

Critical temperature T_c, coherence length ξ₀, and stiffness ({D}_{{rm{s}}}propto {lambda }_{{rm{L}},0}^{-2}) when approaching the multiorbital strong coupling state as function of ∣J∣/W at fixed U/W = 0.8 (a) and when approaching the Mott insulating state as function of U/W at fixed J/W = − 0.04 (b). To indicate a general trend, each quantity is normalized to its maximal value within the line cut. The phenomenology of the SC regimes resembles the BCS to lattice BEC crossover when approaching the Mott insulator (b), which is distinct from the multiorbital strong coupling case (a) (c.f. summary listed in Table 1). Note that the T_c at U/W = 2 in (b) corresponds to the transition from a Mott insulating state instead of a metallic phase, as visually differentiated by a dashed connecting line.

Full size image

ξ₀, in contrast, depends strongly on J. By just slightly increasing ∣J∣ above the ab initio estimate of ∣J∣, the SC state becomes strongly localized with a short coherence length on the order of 2 – 3 a. Remarkably, the small value of ξ₀ is independent of U and thus the proximity to the MI phase. The localization of the condensate with ξ₀ on the order of the lattice constant is reminiscent of a crossover to the BEC-type SC state. However, the dome-shaped behavior, characteristic of the lattice BCS–BEC crossover^27,29, with decreasing T_c in the strong coupling limit is notably absent here. Instead, T_c still grows inside the plateau of minimal ξ₀ when increasing the effective pairing strength proportional to ∣J∣ for fixed U/W (c.f. Fig. 5a). Only by diagonally traversing the (U, J) phase space, it is possible to suppress T_c inside the short ξ₀ plateau with a dome structure as shown, e.g., in Ref. ⁷³.

The reason for this circumvention of the lattice BCS–BEC phenomenology can be understood from an analysis of the local density matrix weights ({rho }_{leftvert {phi }_{n}rightrangle }), where (leftvert {phi }_{n}rightrangle) refers to the eigenstates of the local Hamiltonian of our DMFT auxiliary impurity problem. We show ({rho }_{leftvert {phi }_{n}rightrangle }) of four different points in the interaction phase space in Fig. 4b. In the region of short ξ₀, the local density matrix is dominated by only eight states (red and blue bars) given by the “inverted Hund’s rule” ground states (leftvert {phi }_{0}rightrangle) of the charge sectors with N = 2, 3, 4 particles that are sketched in Fig. 4c. This can be seen in the last panel of Fig. 4a as the total weight of these eight states approaches one when entering the plateau of short ξ₀.

By increasing ∣J∣, the system is driven into a strong coupling phase where local singlets are formed as Cooper pair precursors⁷⁷ while electronic hopping is not inhibited. On the contrary, hopping between the N = 2, 3, 4 ground states is even facilitated via large negative J. It does so by affecting two different energy scales (c.f. Supplementary Note 7): Enhancing J < 0 reduces the atomic gap ({Delta }_{{rm{at}}}={E}_{0}^{N = 4}+{E}_{0}^{N = 2}-2{E}_{0}^{N = 3}=U-2| J|)⁶⁹ relevant to charge excitations and thereby supports hopping. A higher negative Hund’s exchange simultaneously increases the energy ΔE = 2∣J∣ necessary to break up the orbital singlets within a fixed charge sector. As a result, unpaired electrons in the N = 3 state become more itinerant while the local Cooper pair binding strength increases. Since the hopping is not reduced, the SC stiffness is not compromised by larger ∣J∣. This two-faced role or Janus effect of negative Hund’s exchange, that localizes Cooper pairs but delocalizes electrons, can be understood as a competition of a Mott and a charge disproportionated insulator giving way for a mixed-valence metallic state in between⁸².

Correspondingly, as the superconducting state at −J/W > 0.05 relies on direct transitions between the local inverted Hund’s ground states from filling N = 3 to N = 2 and N = 4, the local Hubbard repulsion U has to fulfill two requirements for superconductivity with appreciable critical temperatures. First, significant occupation of the N = 2 and 4 states at half-filling requires that U is not too large. Otherwise, the system turns Mott insulating (around U/W ≳ 1.6 for enhanced ∣J∣) and the SC phase stiffness is reduced upon increasing U towards the Mott limit. At the same time, a significant amount of statistical weight of the N = 3 states demands that U must not be too small either. We find that Δ_at > 0 and thus U > 2∣J∣ is necessary for robust SC pairing. At Δ_at < 0 there is a predominance of N = 2 and 4 states, which couple kinetically only in second-order processes and which are susceptible to charge disproportionation⁸². Our analysis shows that a sweet spot for a robust SC state with high phase stiffness and large T_c exists for ∣J∣ approaching U/2.

We summarize the overall change from a weak-coupling BCS state to the multiorbital strong coupling SC state characterized by T_c, ξ₀, and ({D}_{{rm{s}}}propto {lambda }_{{rm{L}}}^{-2}) upon enhancing ∣J∣ in Fig. 5a and contrast it with the lattice BCS–BEC phenomenology in Table 1. Behavior as in the lattice BCS–BEC crossover can be found upon approaching the MI state (Fig. 5b), where strong local repulsion U decreases ξ₀ and superconductivity compromises stiffness. As discussed above, we cannot observe a proper T_c dome that is seen in the experimental phase diagram due to the MI state dominating over the SC state in our calculations for fixed J/W = − 0.04^36,67,68. However, an additional analysis of the SC gap Δ and coupling strength in Supplementary Note 6 indicates that our calculations are indeed in the vicinity of the crossover region. Diagonally traversing the (U, J) interaction plane by increasing U for J = − γU (γ ≤ 0.5) leads to a similar conclusion, although a T_c dome can emerge more easily⁷³ as charge excitations controlled by Δ_at are suppressed at a slower pace.

Thus, two distinct localized SC states exist in the multiorbital model A₃C₆₀ – one facilitated by local Hubbard repulsion U and the other by enhanced inverted Hund’s coupling ∣J∣. One might speculate under which conditions THz driving or more generally photoexcitation^15,17,18,38 could enhance ∣J∣ and steer A₃C₆₀ into this high-T_c and short-ξ₀ strong coupling region, e.g., via quasiparticle trapping or displacesive meta-stability⁸³. Further experimental characterization via observables susceptible to changes in ξ₀ and λ_L, like critical fields, currents (see Supplementary Note 1-B), or thermoelectric and thermal transport coefficients^8,14, can lash down the possibility of this scenario. Twisted bilayer graphene might be another platform to host the mechanism proposed in this work as an inverted Hund’s pairing similar to that in A₃C₆₀ is currently under discussion^84,85,86.

On general grounds, the bypassing of the usual lattice BCS–BEC scenario via multiorbital physics is promising for optimization of superconducting materials to achieve higher critical currents or temperatures. Generally, limits of accessible T_c are unknown with so far only a rigorous boundary existing for two-dimensional (2D) systems^87,88. An empirical upper bound to T_c emerges from the Uemura classification^5,6,7,8 which compares T_c to the Fermi temperature T_F = E_F/k_B ∝ D_s: the temperature ({T}_{{rm{BEC}}}^{3{rm{D}}}={T}_{{rm{F}}}/4.6) of a three-dimensional (3D) non-interacting BEC. Most superconducting materials, including cuprate materials with the highest known T_c at ambient conditions, however, are not even close to this boundary; typically, T_c/T_BEC = 0.1 – 0.2 for unconventional superconductors. Notable exceptions are monolayer FeSe⁸⁹ (T_c/T_BEC = 0.43), twisted graphitic systems^26,84 (T_c/T_BEC ~ 0.37 – 0.57), and 2D semiconductors at low carrier densities²⁵ (T_c/T_BEC ~ 0.36 – 0.56) which all reach close to the 2D boundary ({T}_{{rm{c,lim}}}^{2{rm{D}}}/{T}_{{rm{BEC}}}^{3{rm{D}}}approx 0.575)⁸⁷. Only ultracold Fermi gases can be tuned very close to the optimal T_BEC^23,24.

The empirical limitations of T_c in most unconventional superconductors with rather high densities make sense from the standpoint and constraints of the lattice BCS–BEC crossover. In this picture, high densities seem favorable for reaching high T_c as electrons can reach high intrinsic energy scales. Yet, lattice effects and their negative impact on the kinetic energy of superconducting carriers become also more pronounced, preventing T_c to reach T_BEC. Possible routes to evade these constraints include quantum geometric and hybridization related band structure effects^{90,91,92,93,94,95,96,97,98} as well as the multiorbital Hund’s interaction effects triggering substantial localized stiffness uncovered here.

We emphasize that the framework to calculate the coherence length ξ₀ and the London penetration depth λ_L introduced in this work can be implemented in any Green function or density functional based approach to superconductivity without significant increase of the numerical complexity. Thus, our work opens the gate for “in silico” superconducting materials’ optimization targeting not only T_c but also ξ₀ and λ_L. On a more fundamental level, availability of ξ₀ and λ_L rather than T_c alone can provide more constraints on possible pairing mechanisms through more rigorous theory-experiment comparisons, particularly in the domain of superconductors with strong electronic correlations.

Dynamical mean-field theory under the constraint of finite-momentum pairing

We study the multiorbital interacting model H_kin + H_int of A₃C₆₀ using Dynamical Mean-Field theory (DMFT) in Nambu space under the constraint of finite-momentum pairing (FMP). In DMFT, the lattice model is mapped onto a single Anderson impurity problem with a self-consistent electronic bath that can be solved numerically exactly. The self-energy becomes purely local Σ_ij(iω_n) = δ_ijΣ(iω_n) capturing all local correlation effects. To this end, we have to solve the following set of self-consistent equations

in DMFT⁵¹. Here, the local Green function G_loc is obtained from the lattice Green function G_k (first line) in order to construct the impurity problem by calculating the Weiss field G_W (second line). Solving the impurity problem yields the impurity Green function G_imp from which the self-energy Σ is derived (third line). Underlined quantities denote matrices with respect to orbital indices (α), ({mathbb{1}}) is the unit matrix in orbital space, ω_n = (2n + 1)πT denote Matsubara frequencies, ({h}_{alpha gamma }({boldsymbol{k}})={sum }_{j}{t}_{alpha gamma }({{boldsymbol{R}}}_{j}){{rm{e}}}^{i{boldsymbol{k}}{{boldsymbol{R}}}_{j}}) is the Fourier transform of the hopping matrix in Eq. (8), μ indicates the chemical potential, and N_k is the number of k-points in the momentum mesh. Convergence of the self-consistency problem is reached when the equality ({underline{G}}_{{rm{loc}}}(i{omega }_{n})={underline{G}}_{{rm{imp}}}(i{omega }_{n})) holds.

We can study superconductivity directly in the symmetry-broken phase by extending the formalism to Nambu–Gor’kov space. The Nambu–Gor’kov function under the constraint of FMP (c.f. Eq. (4)) on Matsubara frequencies

then takes the place of the lattice Green function ({underline{G}}_{{boldsymbol{k}}}(i{omega }_{n})mapsto {underline{{mathcal{G}}}}_{{boldsymbol{q}}}(i{omega }_{n},{boldsymbol{k}})) in the self-consistency cycle in Eq. (9). In addition to the normal component Σ^N ≡ Σ, the self-energy gains an anomalous matrix element Σ^AN for which the gauge is chosen such that it is real-valued.

We use a 35 × 35 × 35 k-mesh and 43200 Matsubara frequencies to set up the lattice Green function in the DMFT loop. In order to solve the local impurity problem, we use a continuous-time quantum Monte Carlo (CT-QMC) solver⁹⁹ based on the strong coupling expansion in the hybridization function (CT-HYB)¹⁰⁰. Details on the implementation can be found in Refs. ^67,68. Depending on calculation parameters (T, U, J) and proximity to the superconducting transition, we perform between 2.4 × 10⁶ up to 19.2 × 10⁶ Monte Carlo sweeps and use a Legendre expansion with 50 up to 80 basis functions. Some calculations very close to the onset of the superconducting phase transition (depending on T or q) needed >200 DMFT iterations until convergence. We average 10 or more converged DMFT iterations and calculate the mean value and standard deviation in order to estimate the uncertainty of the order parameter originating from the statistical noise of the QMC simulation.

Further details on code implementation are given in Supplementary Note 4. It entails a comparison of finite-momentum pairing to the zero-momentum case⁵¹ in the Nambu–Gor’kov formalism, an explanation on the readjustment of the chemical potential μ to fix the electron filling, and a simplification of the Matsubara sum for calculating j_q in Eq. (6).