Introduction

For over three decades, research on the cuprate superconductivity primarily focused on the underdoped and optimally doped region of the phase diagram. Here, it is now widely accepted that Tc is not set by the scale of Cooper pairing (as in BCS theory), but is instead largely determined by the onset of phase coherence (i.e., by the superfluid density)1,2. Phenomena such as the pseudogap, intertwined orders, and strange metal behavior remain the focus of considerable research today. In contrast, it is commonly believed that the physics of the overdoped cuprates is more conventional. For example, angle-resolved photoemission spectroscopy (ARPES) shows a large untruncated Fermi surface in the normal state with reasonably well-defined quasi-particle peaks, and a superconducting gap that decreases with increasing doping, more or less in tandem with Tc3. Moreover, in at least one material4, quantum oscillations, of the sort expected on the basis of band-theory, have been documented.

Thus, it was a surprise that recent penetration depth measurements5,6 on crystalline LSCO films suggest that the superconductivity in the overdoped cuprates is also limited by the onset of phase coherence. Consistent with this result, recent ARPES measurements7 of overdoped Bi2212 found spectroscopic evidence that Cooper pairs are already formed at temperatures about 30% higher than Tc. Adding to the puzzle, recent optical conductivity measurements8 showed that below Tc a large fraction of the normal-state Drude weight remains uncondensed. This is consistent with earlier specific heat measurements which show a T-linear term that persists to the lowest temperatures, T ≪ Tc, with a magnitude that is a substantial fraction of its normal-state value9,10. A possibly related observation11,12,13 from scanning tunneling microscopy (STM) is that, at least up to moderate levels of overdoping, a spectroscopic gap persists in isolated patches up to temperatures well above Tc, so that the normal-state electronic structure is suggestive of superconducting grains embedded in a normal metal matrix. Other than the STM results (for which the relevant data do not exist at very high overdoping), these phenomena become increasingly dramatic as the doped hole concentration, p, approaches the critical value, psmt, at which the superconductor-to-metal transition occurs at the overdoped end of the superconducting dome.

The primary goals of this paper are to present a simple theoretical model that captures what we believe to be the essence of the above phenomena, and to explain the cause of the superconductor–metal transition. We are aware that our model does not capture various quantitative aspects of the actual materials. In the next two paragraphs, we present the physical picture underlying this work.

A sketch of the Fermi surface of an overdoped cuprate is shown in Fig. 1; it is shown as being hole-like, although in some cuprates the Fermi surface passes through a Lifshitz transition at a doping concentration, pLif < psmt, in which case it would be electron-like14,15. The neighborhood of the van-Hove points—which we will refer to as the antinodal regions—is also the portion of the Fermi surface farthest from the gap nodes in the d-wave superconducting state and so is where the gap is largest. As we will discuss, the fact that the Fermi surface passes near the van-Hove point, meaning that the Fermi velocity is small in the antinodal regions, plays a significant role in the results we obtain; whether it is electron or hole-like is relatively less important.

Fig. 1: Model Fermi surface.
figure 1

Antinode to antinode scattering induced by the disorder potential is indicated by the yellow arrow.

Full size image

In considering the effects of disorder, the scattering between antinodal regions (indicated by the orange line in the figure) is particularly important, as it is pair-breaking. Such antinode to antinode scattering is apparent in STM quasi-particle interference measurements in both non-superconducting16 and superconducting17 overdoped Bi2201. In particular, in ref. 17, it is shown that at voltages corresponding to the antinodal gap energy, the Fourier transform of the local density of states exhibits a broad maximum at momentum q = (π, π). Moreover, from the coherence factor, it is inferred that this scattering occurs between momentum regions having opposite signs of the gap function17.

When the Cooper pair coherence length is comparable to the correlation length of the disorder potential, the prior mentioned pair-breaking causes the pair-field amplitude to be spatially heterogeneous18,19. The superfluid stiffness is large in regions with high pair-field amplitude, whereas the stiffness is low where the pair-field amplitude is small. This is reminiscent of a granular superconductor. The small stiffness in the inter-granular regions causes the averaged zero-temperature superfluid density to be low, hence superconducting phase fluctuations (both classical and quantum) are enhanced. In the metallic regions the Cooper pairing instability is inhibited since (a) the repulsive interactions between electrons force the average of the superconducting order parameter to be zero around the Fermi surface, and (b) a sign-changing order parameter is suppressed by disorder scattering. The unpaired electrons in the inter-granular regions give rise to a substantial uncondensed Drude component in the optical conductivity and to a residual T-linear term in the specific heat. A superconductor-to-metal transition occurs when the superconducting islands grow sufficiently sparse18.

Note that the normal-state transport is dominated by the nodal quasiparticles, which are known to be less affected by impurity scattering20,21. In the rest of the paper, we present results corroborating the physical picture presented above.

Results

The model

The model we use to describe the superconducting state contains hopping terms, interaction terms, and disorder potential terms. The Hamiltonian is

$$H=-\sum _{i,j,\sigma }{t}_{ij}({c}_{i\sigma }^{\dagger }{c}_{j,\sigma }+h.c.)+\sum _{i,\sigma }({w}_{i}-\mu ){c}_{i\sigma }^{\dagger }{c}_{i,\sigma }+{H}_{{\rm{int}}}$$
(1)

where tij is the hopping integral between sites i and j on a square lattice which we take (to produce a cuprate-like Fermi surface) to be tij = 1 between nearest-neighbor sites, tij = −0.35 between second-neighbor sites, and tij = 0 for all further neighbors. To compare different pairing symmetries, we consider two different forms of Hint: (1) As a model of a d-wave superconductor (relevant to the cuprates) we adopt a model with a nearest-neighbor antiferromagnetic Heisenberg exchange interaction, Hint = J∑〈ij〉Si ⋅ Sj. (2) As a model of an s-wave superconductor, we consider an attractive Hubbard interaction, \({H}_{{\rm{int}}}=-U{\sum }_{i}{c}_{i\uparrow }^{\dagger }{c}_{i\uparrow }{c}_{i\downarrow }^{\dagger }{c}_{i\downarrow }\). We fix the strength of the pairing interactions to J = 0.8 and U = 1.35, respectively, so that the two superconducting gap scales in the absence of disorder are approximately the same.

In treating the problem with disorder, we consider a finite system of size 40 × 40 and, unless otherwise indicated, assume periodic boundary conditions. The random potentials, wj, represent the effects of disorder: on a randomly chosen fraction nimp of sites we set wj = w > 0, with wj = 0 on all other sites. Modeling disorder as point-like impurities is a simplification. In reality, the potential produced by the dopants can extend over multiple unit cells. The point-like impurity is considered because it can cause the large momentum transfer of antinode to antinode scattering seen in STM16,17, which is a key ingredient in our theoretical framework. In the main text, we report results for disorder strength, w = 1, but in Supplementary Note 2 we include results for other values—the main qualitative results do not depend sensitively on the value of w. We repeat this with multiple different impurity configurations in order to compute the configuration averages of physical observables; typically, we average over 64 distinct impurity configurations but we average over 128 configurations when the impurity concentration is large and the superconducting pairing is highly inhomogeneous.

We solve the model by self-consistent BCS mean-field approximation. However, the disorder scattering is treated exactly. Since H lacks translation invariance, the self-consistency equations need to be solved numerically. Since our calculation does not capture the thermal and quantum fluctuations, we focus primarily on zero temperature and on doping sufficiently away from the quantum critical doping of the superconductor–metal transition.

The carrier concentration is controlled by the chemical potential μ and the impurity concentration nimp. Sometimes we fix μ while changing nimp to reach the desired carrier (hole) concentration. In the case where we want to study the effect of disorder at a fixed carrier density, we tune μ while varying nimp to achieve the desired carrier density. To avoid the complex issues of the pseudogap, intertwined orders, and strange metals, we take our lowest carrier concentration to be slightly larger than optimal doping. Thus except in section “Role of the antinodal dispersion”, the impurity concentration is measured relative to that present at optimal doping. For example, in Fig. 2 we show the hole concentration p as a function of nimp at fixed μ with w = 1.

Fig. 2: Plot of the density of doped holes as a function of impurity concentration.
figure 2

The disorder strength is fixed with w = 1.

Full size image

The mean-field solution

The local value of the gap parameter, Δij, that enters the mean-field equations, which we will refer to as the pair field, is given by the product of the pairing interaction times the expectation value of the pair annihilation operator. In the s-wave case, Δ is site diagonal, Δj ≡ Δjj = U〈cj↑cj↓〉 while for the d-wave case, Δij = J〈ci↑cj↓ + cj↑ci↓〉, where i, j are any pair of nearest-neighbor sites. The self-consistently computed values of the pair field for two impurity configurations with different doped hole concentrations are shown in Fig. 3.

Fig. 3: The real-space distribution of the pair field.
figure 3

The upper two panels are for d-wave pairing (where the pair fields lie on nearest-neighbor bonds); the lower two panels are for s-wave pairing (where the pair fields are on-site). The size of the symbols, namely the thickness of the bonds in (a) and (b), and the size of the dots in (c) and (d), represents the magnitude of the pair field whereas the color (red positive, blue negative) the sign. The left and right columns correspond to two different impurity concentrations. The magnitude of pair field in panel a ranges from 0.0003 to 0.103 while that in (b) ranges from 0.000005 to 0.1008. In c and d, the magnitude of the pair field on each site ranges from 0.076 to 0.18 and 0.073 to 0.20, respectively.

Full size image

In the s-wave case, we find that Δj has a uniform sign and a magnitude that is weakly dependent on position. Moreover, it does not depend on the doped hole concentration strongly. The red symbols in Fig. 4a show the configuration averaged value of ∣Δi∣, as a function of p.

Fig. 4: The doping dependence of the spatial averaged zero-temperature magnitude of the pair field and the mean-field value of the superfluid density.
figure 4

The band structures used in these plots are the same, the only difference is the pairing interaction. The red symbols represent s-wave pairing while the black symbols represent the d-wave pairing. In the band structures used in constructing these plots, the Fermi surface crosses van-Hove point, namely Lifshitz transition occurs, at the doping level pLif ≈ 0.33.

Full size image

In the d-wave case, Δij has a magnitude that varies significantly as a function of position and which is strongly doping dependent. It also reflects the d-wave symmetry of the uniform state from which it descends in that, with minor exceptions, Δij is positive on bonds oriented in the \(\hat{x}\) direction, and negative on bonds in the \(\hat{y}\) direction. The black symbols in Fig. 4b show the configuration average of ∣Δij∣ as a function of p. Notice that it drops dramatically with increasing p, but then has a long tail with small magnitude that extends to high values of p.

Note that the band structures used in the two cases are the same; only the pairing interaction is different. The dramatic contrast between the two cases manifests Anderson’s theorem for the s-wave, and pair-breaking by the scalar disorder for the d-wave.

Specifically, in the d-wave case, pair-breaking induced by the antinode to antinode scattering tends to strongly suppress superconducting pairing. A consequence of this is that when disorder is strong, the pair-field amplitude becomes granular (heterogeneous) with significant pairing occurring only in rare regions where disorder is weak19. This can be seen clearly in Fig. 3b. The superconducting order parameter on different grains are connected by effective SNS (superconductor–metal–superconductor) Josephson junctions which, as suggested in ref. 22, can vary randomly in sign due to the sign-changing superconducting order parameter. This can cause frustration in the superconducting phase coherence, and ultimately, as we will see shortly, to spontaneous time-reversal-symmetry breaking and the existence of local super-current loops.

The net superfluid density is computed from the standard Kubo formula, Supplementary Eqs. (2) and (3). In Fig. 4b we compare the p dependence of the T = 0 mean-field superfluid density for the s-wave and d-wave pairing cases. Notice that for the d-wave case, the superfluid density drops considerably more rapidly than does the pair-field amplitude. (This can be seen more quantitatively in Supplementary Fig. 1, where the d-wave case is shown on a log-linear scale.) This implies that phase fluctuation effects, beyond the mean-field treatment, must inevitably become large in this range of doping.

We have also computed the T dependence of the mean-field superfluid density. For the d-wave case, the results are shown in Supplemental Fig. 4. At low T (where we see a T-linear decrease), the results may be physically meaningful, but at higher temperatures, thermal-phase fluctuations, which are ignored in our mean-field treatment, must certainly play a role in the vanishing of the superfluid density as T → Tc.

Equilibrium current loops

A subtle but remarkable feature of the mean-field solution in the highly overdoped regime is shown in Fig. 5. Here, the blue-colored bonds represent ∣Δij∣ while the red arrows represent equilibrium supercurrents. Here, the current operator on bond \(\left\langle ij\right\rangle\) is given by

$${J}_{ij}=i{t}_{ij}\sum _{\sigma }\langle {c}_{i\sigma }^{\dagger }{c}_{j\sigma }-h.c.\rangle$$
(2)

We stress that we have not added any time-reversal symmetry-breaking perturbation. The currents form loops, thus satisfying the continuity equations, and are manifestations of the spontaneous breaking of time-reversal symmetry expected from the random-in-sign Josephson couplings that emerge when superconducting islands are small and sparse. The patterns of currents are analogous to those in an XY spin-glass, with near-degeneracies associated with reversing the local currents around localized loops in different regions of the system. This near degeneracy is reminiscent to the existence of two-level centers in a glass. They can give rise to orbital paramagnetism. However, since all these effects occur where the superfluid density is small, there may be important qualitative changes in this behavior when the effect of thermal and quantum phase fluctuations are included. For instance, the near degeneracy of the two-level centers can be lifted by tunneling processes in which the local currents reverse direction. To the extent that they survive fluctuational effects, these spontaneous current loops are a qualitatively significant result of the combination of d-wave pairing and scalar disorder.

Fig. 5: Equilibrium current loops.
figure 5

The blue-colored bonds represent the absolute value of the d-wave pair field. The red arrows represent the spontaneously generated super-current. The thickness of the arrow denotes the magnitude of the current. The impurity concentration is 0.19 in (a) and 0.31 in (b). No detectable current exists (smaller than 10−12) in (a) while they are quite apparent for (b). c Zoom-in view of the lower-left corner indicated by the black dashed line in panel (b). For clarity, only when the magnitude of current is >0.001 (which is ~2/10 of the maximum current value) do we plot a dark red arrow. For smaller current values we use pink arrows to represent them.

Full size image

Role of the antinodal dispersion

Because the d-wave gap is maximal in the antinodal region of the Brillouin zone (BZ), many features of the mean-field solution depend sensitively on the band structure in this region. Here we show the effects of the flat antinodal dispersion in overdoped cuprates. Such effects have been emphasized in the recent photoemission work7. Due to the enhanced density of states, we find that the existence of the flat dispersion amplifies the disorder-induced antinodal scattering, and hence enhances pairing heterogeneity. Moreover, it leads to more rapid suppression of the zero-temperature superfluid density with increasing p.

In Fig. 6, we compare results for two band structures: one with a flat dispersion near the antinodes (as is generic in the cuprates), and the other in which the antinodal dispersion is relatively steep. The band-structure parameters are chosen such that the doping level is fixed at 23% and the Fermi energy (3.875t) is the same for the two different band structures (see Supplementary Note 3). The most significant difference between the flat and steep bands is the existence/absence of a flat dispersion along the BZ boundary as shown in the relevant part of the BZ in Fig. 6a. The spatial averaged zero-temperature pair-field amplitudes and superfluid densities as a function of impurity concentration are shown in panels b and c. Note that in these figures, when the impurity concentration is varied, we tune the chemical potential so that the hole concentration is unchanged. It is apparent that the suppression of the pair-field amplitude and of the superfluid density are considerably more rapid in the case of a flat band. This result is consistent with recent ARPES results7 and the interpretation therein. We attribute the more rapid suppression of the pair-field amplitude and the superfluid density to the larger density of states associated with the flat band, and the concomitant enhancement of the pair-breaking antinode to antinode scattering.

Fig. 6: The comparison between flat and steep bands.
figure 6

The spatially averaged zero-temperature d-wave pair-field amplitude (b) and superfluid phase stiffness (c) computed with the two different band structures shown in (a). Results for the partially flat band are shown as the black squares and for the steep band by the red.

Full size image

The specific heat and optical conductivity

A feature of the inhomogeneous state is that there remains a large density of gapless quasi-particle states arising from the approximately normal metallic regions between the superconducting grains. This is reflected in a residual T-linear contribution to the specific heat and to the ω → 0 optical conductivity that survives even as T → 0. In Fig. 7, we plot the ratio between the low-temperature specific heat coefficient γ ≡ c/T and the corresponding value in the normal state as a function of doping concentration. Notably, when the doping concentration is high, the ratio approaches one. In Supplementary Fig. 3, we plot the real part of the optical conductivity as a function of frequency at T = 0, which shows that a large portion of normal-state Drude weight is uncondensed. In Supplementary Note 4, we discuss some discrepancies in the frequency dependence of our result when compared with the experimental data of ref. 8.

Fig. 7: Doping dependence of specific heat coefficient γ1 normalized by the corresponding normal-state value γ1N.
figure 7

The zero-temperature specific heat coefficient is extracted by fitting C(T)/T using C(T)/T = γ1 + γ2T and extrapolating to T = 0. The fitting is illustrated in the inset, where the left panel is for the d-wave SC state while the right corresponds to the normal state.

Full size image

Discussion

There are several aspects of the superconductor-to-metal transition in overdoped cuprates that have been the subject of various recent theoretical studies.

One issue concerns the microscopic mechanism for the superconductor-to-insulator transition. While on the underdoped side of the superconducting dome, the gap scale appears to remain large even as Tc → 0, on the overdoped side the gap (measured in various ways) decreases significantly as Tc → 0. Since it is generally believed that spin-fluctuations are a dominant contributor to the d-wave pairing, and since signatures of incipient antiferromagnetic order become increasingly weak with increasing p, it is reasonable to associate the drop in Δ with a weakening of the pairing interaction. That such a trend occurs in a Hubbard model with a band structure suitable for the cuprates has recently been shown in ref. 23. On the other hand, whether the short-range antiferromagnetic correlation drops sufficiently strongly to account for the demise of the superconducting phase is under debate24. In the present study, we showed (at mean-field level) that even holding the strength of the pairing interaction fixed, a strong drop in the pairing scale can be accounted for simply as a consequence of an increased density of random scattering centers. We consider it likely that both effects play a role in the overdoped cuprates.

Another issue concerns how disorder is treated. In recent theoretical studies of overdoped cuprates23,25,26, the effects of disorder are treated in an effective medium approximation, in which macroscopic fluctuations in the local impurity configurations are averaged out, and the superconducting state is homogeneous. These studies ignore any self-organized granularity. Nonetheless, with suitable choices of parameters, they have been shown to produce phenomenologically reasonable results. However, we feel that the experimentally observed5 quantitative similarity between Tc and Tθ (i.e., the T = 0 superfluid density expressed2 in temperature units), implicates the reduced superfluid density as the cause of the superconductor-to-metal transition. The emergent granularity that we have found provides a theoretically sound origin for such a reduced superfluid density.

The current work, emphasizing the interplay of d-wave pairing and disorder as the cause of a superconductor-to-metal transition, is qualitatively different from older studies of the superconductor-to-insulator transition in the case of an s-wave SC27,28, despite the appearance of self-organized granularity in both cases. In particular, the residual resistivity at the end of SC dome is roughly 10 μΩ-cm, i.e., two orders of magnitude smaller than the quantum of resistivity5 per Cu–O plane. Although our model and calculation are similar to some previous works of mean-field calculation on dirty d-wave superconductors29,30,31,32,33,34, the issues we address, namely the disorder-driven superconductor-to-metal transition, have not been sharply articulated previously.

Thus, we propose the key to understanding the essence of overdoped cuprates is the combined effect of disorder and d-wave pairing. These features, especially when combined with relatively flat bands in the antinodal regions of the BZ, lead to a self-organized granular superconducting state. As an additional consequence of the d-wave pairing, the Josephson coupling between the superconducting islands is generically frustrated in the strong disorder limit. As a consequence, there are spontaneous current loops and associated local breaking of time-reversal symmetry. There are several possible observable signatures of this22; for instance, we find that under conditions in which the mean-field zero-temperature superfluid density is significantly suppressed by disorder, the superfluid density can be an increasing function of an applied field for small fields. In addition, due to the formation of local current loops, the system can exhibit a paramagnetic Meissner effect35, namely, a paramagnetic response upon applying a small magnetic field. Such effects are beyond the reach of theories which treat the disorder in effective medium approximation.

Finally, we reiterate that the model, and the mean-field treatment of it we have discussed, are simplified compared to any actual cuprate. Indeed, it follows from the present results that thermal and quantum fluctuations, not included explicitly in our treatment so far, are inevitably important in the vicinity of superconductor-to-metal transition. Thus, while our work is a step toward understanding the basic physics of superconductor-to-metal transition in overdoped cuprates, many important issues are still unsettled.

Methods

Self-consistency mean-field calculation

The results are obtained by numerically self-consistency mean-field calculation. The system size in our calculation is 40 × 40. We average over distinct disorder configurations to compute the physical observables. Typically, the number of disorder configurations is 64, but we average over 128 configurations when the impurity concentration is large and the superconducting pairing is highly inhomogeneous.