Generalization of Anderson's Theorem for Disordered Superconductors
John F. Dodaro Steven A. Kivelson




#1⟨#1⟩
⟨#|1⟨#1|
|#⟩1|#1⟩
⟨#|1⟩#2⟨#1|.#2⟩
#1#2^#1_#2
#1|#1|
#1(#1)
#1[#1]
#1{#1}
Re
Im
sgn
 mod 
#1Tr[#1]
#1Det(#1)

    [ ∂>̄⇒→↑↓Phys. Rev. BPhys. Rev. Lett.2FathbbI kΤFathcalTFathcalEFathcalKFathbbZ#1c#1#1{cc #1. Cliff#1 Pf#1mattergaugeΔ_K⃗Δ_-K⃗Δ_K⃗^'Δ_-K⃗^'k̨̨⃗Q⃗q⃗K⃗P⃗r̊⃗̊0̣Δ_0#1Eq. (<ref>)#1Fig. <ref>Department of Physics, Stanford University, Stanford, CA 94305-4060, USA
        We show that at the level of BCS mean-field theory, the superconducting T_c is always increased in the presence of disorder, regardless of order parameter symmetry, disorder strength, and spatial dimension.  This result reflects the physics of rare events – formally analogous to the problem of Lifshitz tails in disordered semiconductors – and arises from considerations of spatially inhomogeneous solutions of the gap equation.  So long as the clean-limit superconducting coherence length, ξ_0, is large compared to disorder correlation length, a, when fluctuations about mean-field theory are considered, the effects of such rare events are small (typically exponentially in [ξ_0/a]^d);  however, when this ratio is ∼ 1, these considerations are important.  The linearized gap equation is solved numerically for various disorder ensembles to illustrate this general principle.
         Generalization of Anderson's Theorem for Disordered Superconductors
            John F. Dodaro and Steven A. Kivelson
            
        ====================================================================
        
        
        
        
        Most theoretical analyses of the effect of disorder on the superconducting transition rest on two separate assumptions:  the pairing instability is treated in the context of a generalized BCS mean-field theory, and the disorder is treated either perturbatively or in an effective medium approximation.  The primary focus of this paper is the demonstration that there is additional universal structure to the mean-field solution of this problem when the effect of disorder is treated exactly.  In particular, the presence of rare regions – which are neglected in effective medium approximations – leads to the unintuitive conclusion that disorder always increases the mean-field T_c, independent of the nature of the disorder and whether we are considering a conventional (s-wave) or unconventional (p-wave or d-wave) superconductor.  We demonstrate this explicitly with numerical solutions of the linearized gap equation for various disorder ensembles.
        
        There is a separate issue of whether this improved analysis of the BCS equations is physically significant.  When the regions that support a local superconducting order parameter at mean-field level are sufficiently rare, the inclusion of fluctuation effects beyond mean-field theory reveals the inferred high transition temperatures to be artifacts of mean-field theory;   in this case, an incorrect solution of the mean-field equations can give a more physically sensible solution.  However, when the T=0 superconducting correlation length, ξ_0, is comparable to the disorder correlation length, a, locally superconducting solutions of the mean-field equations can lead to significant consequences – such systems exhibit a broad fluctuational regime in which a form of self-organized granularity arises.  We stress this is unavoidable in superconductors with correlation lengths comparable to the lattice constant, even if the disorder is entirely homogeneous as in a “perfect” substitutional alloy.
        
        
        
        
        § THE MODEL
        
        We consider the tight-binding model
        
        H=-∑_a,b t_ab c_a^† c_b + ∑_ab b^' a^'V_ab b^' a^' c_a^† c_b^† c_ b^' c_a^'
        where a=(r⃗,σ) is an index specifying the lattice position, r⃗, the spin polarization, σ, and possibly – if relevant – an orbital index, and c^†_a creates an electron in single-particle state a.   We assume that both t and V are short-range in space, but they need not be translationally invariant.  Thus, to fully specify the problem we need to specify the ensemble which determines the probability of different “configurations,” i.e. realizations of H.  We will imagine that the ensemble has statistical translational symmetry, so that the configuration-averaged version
        H = - ∑_r⃗,r⃗^',σt(r⃗-r⃗^ ') c_r⃗,σ^† c_r⃗,σ 
        + ∑_r⃗,r⃗^'V(r⃗ -r⃗^ ')n̂(r⃗) n̂(r⃗^ ')
        
        is translationally invariant.  (For simplicity, we have ignored spin-orbit coupling and assumed a density-density interaction, n̂(r⃗) = ∑_σ c^†_r⃗σc_r⃗σ, although the same considerations apply to more general models.)  Then when we ask about the effect of disorder on T_c, we are comparing the case in which we replace H by H (translation symmetry restored) with the full problem using the configuration-dependent H.  We will always focus the analysis on systems in the thermodynamic limit, and will assume that the disorder ensemble is sufficiently well-behaved that the system properties are self-averaging in this limit.  Since in general V includes effective interactions obtained by integrating out high energy degrees of freedom, we consider disorder in V as well as single-particle disorder t_aa.  (See, for example, Ref. nunner_2005, andersen_2006)
        
        
        
        
          §.§.§ BCS mean-field theory
         The interacting problem can be treated at mean-field level by introducing the trial Hamiltonian
        
         H_tr=-∑_a,bτ_ab c_a^† c_b + ∑_ab[Δ_ab c_a^† c_b^† + H.C.]
        
        where the parameters are chosen so as to minimize the variational free energy
        
        F_var[τ,Δ] = F_tr + ⟨ H-H_tr⟩_tr
        where the expectation values are taken with respect to the Gibbs ensemble of H_tr.  From the condition that F_var is stationary with respect to variations of the parameters entering H_tr, we obtain the self-consistency conditions 
        Δ_ab = ∑_ a^' b^'V_ab b^' a^'⟨ c_ b^' c_ a^'⟩
        and a corresponding self-consistency conditions for τ_ab.  
          
        
        Still at mean-field level, the superconducting state is characterized by a broken symmetry, and thus occurs whenever the lowest free energy solution of these equations has a non-vanishing value of Δ_ab for any ab;  the mean-field T_c is  the upper bound on temperatures at which such a solution exists.   If the mean-field transition is continuous, T_c can be identified by studying solutions of the linearized gap equation. 
        Specifically, the pair field expectation value to linear order is given by
        
        
        
        
        
        
        ⟨ c_ac_b⟩ =∂ F_tr/∂Δ^⋆_ab= - χ_ab, a^' b^' (T) Δ_a^' b^' +  O(|Δ|^2Δ^⋆)
        
        where the notation is such that repeated indices are summed over, and χ is the pair-field susceptibility computed in the normal state (i.e. in the thermal ensemble of H_tr with Δ_ab=0). 
        It is easy to see that χ is a Hermitian, positive semidefinite matrix 
         so the matrix square-root  χ^1/2
        
        
        
        exists and is also Hermitian.
        
        With this in mind, we can expand the self-consistency equation (<ref>) to linear order.  Specifically,    consider the solutions of the eigenvalue equation
        λ_n(T)  F^(n)_ab = 
        M_ab, b^' a^'(T)   F^(n)_b^' a^'
        where M is the real Hermitian matrix 
        
        M_ab,cd≡ -
        χ^1/2_ab,a^' b^'(T) V_a^' b^',c^' d^' χ^1/2_ c^'  d^', cd(T)
        
        
        
        
        
        
        
        
        Ordering the eigenstates λ_n≥λ_n+1, the mean-field T_c is then identified as the temperature at which λ_max(T_c) ≡λ_0(T_c)=1.  If the mean-field transition is first order, the transition temperature must always be greater than that inferred in this way; thus, λ_max(T)=1 is a lower-bound estimate of the mean-field T_c.
        
        
        
        
        
        
        § QUALITATIVE ASPECTS OF THE SOLUTION
        
        
        The linearized gap equation can be solved numerically on moderately large systems for any given configuration, as discussed below.  Analytic solutions are difficult, but a few general observations allow for qualitative statements under broad circumstances.  
        
        We consider a short-ranged interaction V_a_1a_2,a_3a_4 which vanishes sufficiently rapidly so as to be negligible for |r⃗_i-r⃗_j| ≫ 1 for any i≠ j (with lattice constant = 1).  However, the physics of metals is reflected in the behavior of χ, which becomes increasingly long range as T→ 0.  Specifically, at T=0, in either a clean metal or a diffusive metal, χ_aa,a^' a^'∼ |R⃗|^-d for |R⃗| ≫ 1; 
          the fact that ∑_R⃗χ_(0⃗↑),(0⃗↓),(R⃗↓),(R⃗↑) is logarithmically divergent is  a direct manifestation of the Cooper instability of a Fermi liquid – and the fact (Anderson's theorem) that it persists even in the presence of  disorder.<cit.>
        
        Were we studying quantum phase transitions at T=0, the long-range character of χ would make a statistical analysis of M particularly subtle  (see Ref. spivak_2008).  However, as we are studying finite temperature transitions, there is always a finite length scale, ℓ(T), beyond which χ falls exponentially.  Specifically, in a metal, ℓ(T) is determined by a thermal coherence length, L_T, that diverges with a power of T as T→ 0:  as L_T = ħ v_F/T in a clean metal and as L_T = √(ħ D/T) in a diffusive metal where v_F and D are, respectively, the Fermi velocity and the electron diffusion constant.  If the disorder is sufficiently strong to produce localization (which is to say any disorder in 2d), the same considerations apply so long as L_T<ξ_loc where ξ_loc is the localization length, but at low enough temperatures, L_T saturates to ξ_loc.  In other words, ℓ(T) ≈min[L_T,ξ_loc].  Thus, M_a_1a_2,a_3a_4 is short-ranged in the sense that it vanishes exponentially when |r⃗_i-r⃗_j| ≫ L_T for any i≠ j.[Even at finite T,  the fact that ℓ(T) ≫ 1 is an essential feature of the theory of superconductivity.  In all cases of interest, the interaction strength in the appropriate pairing channel is always small compared to E_F.  Thus, even the largest single matrix elements in M are small in magnitude – it is only the fact the range of M diverges in a non-integrable fashion that gives rise to the Cooper instability in this limit.]
        
        
          §.§.§ Lifshitz tails and a  “theorem”
        
        The problem of determining the spectrum of M is thus structurally similar to the problem of a non-interacting quantum particle moving in a random potential.  We can view M as a temperature-dependent tight-binding model with moderate but finite range random hopping terms extracted from a complicated disorder ensemble; the higher the temperature, the more local is M.  The generic structure of the solutions is familiar, and therefore without further analysis we can infer various properties of the present problem on the basis of well-established results from the random potential problem.  We will adopt the terminology that anything that can be established on the basis of this precise analogy is “proven,” although (in common with “Anderson's theorem”) no proof in the mathematician's sense of propositions and lemmas will be attempted.  
        
        On this basis, we conclude that there is a T-dependent “density of states,”ρ(λ;T), for the eigenvalues of M, and that this distribution is self averaging, ρ(λ) = ρ(λ).  Where ρ is large, we expect solutions to be delocalized or at most weakly localized on exponentially long length scales.  But likewise, there are universal “Lifshitz tails” to the distribution <cit.> that extend to values of  λ well in excess of the mean.  In the tails, the eigenstates are strongly localized and are associated with unusual rare regions in which the configuration is exceptionally conducive to superconductivity.  The structure of this tail in the density of states depends on details of the disorder.  In many circumstances its asymptotic form, in the regime in which ρ(λ;T) is approaching zero, can be estimated by solving an optimization problem <cit.>, which in turn can be related to a replica symmetry breaking instanton solution of an appropriately replicated effective field theory. <cit.>  
        
        In the present problem, as in the problem originally analyzed by Lifshitz <cit.>, we expect that the distribution is bounded, i.e. that ρ(λ;T) ≠ 0 only within a finite range λ_min(T)< λ < λ_max(T).  Consider all possible ways the system could be organized within the allowed ensemble, including all possible concentrations and local (possibly highly ordered) arrangements of impurity atoms; λ_max(T) is then associated with the optimal configuration.  
        
        Of course, the concentration of regions of the sample that happen to well approximate this optimal configuration is extremely small.  The more precisely we require the optimal configuration, and the larger the region in question, the smaller the concentration of such regions; however, in the thermodynamic limit, for any given size and required precision, the probability of finding such a region is non-zero, and hence ρ(λ) >0  so long as λ < λ_max. 
        
          Thus, the mean-field critical temperature is bounded below by T_c,max, which is defined implicitly as the solution of λ_max(T_c,max) = 1. In other words, the mean-field transition of a disordered system is determined by  the highest possible T_c of any system that is allowed by the physics!  This is generically larger than T_c of the average Hamiltonian, H.  This is the  new “theorem.”
        
        
          §.§.§ Significance of states in the tails
        
        The eigenstates with λ close to λ_max are localized in rare regions with a local configuration peculiarly conducive to superconductivity.  For temperatures slightly smaller than T_c,max, these regions can be treated as isolated superconducting “puddles.”  The magnitude of the order parameter Δ^(n) on the puddle n– which depends on the character of the non-linear terms in Eq. <ref> and so cannot be directly computed from the solutions of the linearized equation we have studied here – can nonetheless be estimated to depend on λ_n as 
         Δ^(n)(T) ∼Δ_0 √(λ_n(T) -1). 
         
        So long as the puddles are both rare and uncoupled in this range of temperatures,  their effect on the global properties of the system will generally be relatively weak.  However, the presence of an increasing concentration n_sc of locally superconducting regions with decreasing temperature could be detected in various ways:  local spectroscopy will see the opening of what should look like a superconducting gap in any region in which λ_n(T) >1.  
        
        It is only when n_sc(T) ℓ^d(T) ceases to be small that the Josephson coupling between such superconducting puddles will start to become significant, so that the more familiar signatures of superconductivity, such as a significant drop in the resistivity and substantial diamagnetism, will appear.  Global phase coherence –i.e. the true T_c– is a still lower temperature where the puddles effectively begin to overlap.  Even though the electronic structure is that of a statistically homogeneous metal for T>T_c,max and a correspondingly homogeneous superconductor for T_c ≪ T, it is effectively  an inhomogeneous mixture of superconducting puddles embedded in a normal metal matrix for T_c ≲ T < T_c,max.  In weakly coupled superconductors, in which ξ_0 is large compared to the correlation length of the disorder, this regime is parametrically small (by a factor more or less equivalent to the usual Ginzburg parameter).  However, for superconductors with a short correlation length, this regime is not small, in which case significant deviations from mean-field behavior are to be expected, and T_c is determined more by phase ordering than by (local) gap formation. 
        
        
        
        
          §.§.§ Approximate analytic considerations
        
        
        
        
        
        
        
        
        To get a feeling for the nature of the states in the Lifshitz tail, we make a variational ansatz for a localized eigenstate,
        
        F_a_1a_2 =  f_σ_1,σ_2(r⃗_1-r⃗_2) F(r⃗_1+r⃗_2/2)
        
        with the normalization condition ∑_r⃗, σ, σ^'|f_σ,σ^'(r⃗)|^2 = 1.
        
        
        
         
        This ansatz assumes that there is a preferred form of the pair-wave-function, with the only variational freedom associated with the local magnitude, F(R⃗).  The simplified eigenvalue equation is then
        ∑_r⃗^' M(r⃗,r⃗^ ') F(r⃗^ ') = λ F(r⃗)
        
        
        
        
        
        
        Moreover, if we further assume that F is a slowly varying function over the range of M (i.e.ℓ) then we can treat r⃗ as a continuous variable and expand the eigenvalue equation in a gradient expansion as
        [ γ(r⃗;T) + (1/2)∇_μ W_μ,ν(r⃗;T)∇_ν + ⋯] F(r⃗) = λ F(r⃗).
        
        Finally, again invoking the central limit theorem, we can approximate γ(r⃗;T) ≈γ(T) + δγ(r⃗) where δγ(r⃗) is a local gaussian random variable with variance δγ(r⃗)δγ(r⃗^ ')=δ(r⃗-r⃗^ ')σ^2, and W_μ,ν(r⃗;T) ≈W_μ,ν(r⃗;T)  = δ_μ,ν(1/m).  We moreover expect m>0 reflecting the fact that the superconducting susceptibility generally favors spatially uniform pair-wave-functions over oscillatory ones, and on dimensional grounds we expect m ∼ℓ(T)^-2.  
        
        This leaves us precisely with the problem of Lifshitz tails for a particle moving in a random potential, and  asymptotic forms of ρ(λ) have been derived in several classic references. <cit.>  For d<4, the essential features of the results can be derived rather simply as follows:  assuming we are interested in localized solutions with λ > γ, we look to find a region of size R in which the average value of δγ is sufficiently positive to produce the desired eigenvalue.  This requires that the constraint δγ =(λ - γ) -β m^-1 R^-2 be satisfied where the precise value of β depends on the shape of the assumed puddle – it has its smallest value if the puddle is spherical.  The probability of finding such a puddle ∼exp[-S] where S=β^' R^d (δγ)^2/[2σ^2] and β^' R^d is the volume of the puddle.  Minimizing S with respect to δγ subject to the above constraint, yields
        
        S = ασ^-2ℓ^d  [λ - γ(T)]^(4-d)/2
        and R = α^' ℓ [λ -γ(T)]^-1/2 where α and α^' are numbers of order 1 (that can be expressed in terms of β, β^', m_0, and d.).   Finally, if we assume that γ(T) can be approximated as a linear function of T as γ(T) ≈ 1 - γ_0(T-T_0)/T_0, then this yields a concentration of superconducting regions
        
        n_sc(T)  ∼      ρ(λ=1,T) 
              ∼exp{- A (ξ_0/a)^d (T-T_0/T_0)^(4-d)/2}
        where T_0 is the average T_c, ξ_0 is the corresponding (T=0) superconducting coherence length, a is the correlation length of the disorder, and A∼γ_0/σ^2 is a constant that depends on the strength of the disorder and other  features we have swept under the rug.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        § SOME NUMERICALLY SOLVED EXAMPLES
        
        
        To see how the general considerations outlined above play out explicitly, we have computed ρ(λ) by numerical solution of Eq. <ref> on a square lattice with periodic boundary conditions and  sizes L× L with L between 24 and 50.  Typically we average our results over 600 distinct disorder configurations.  For the “uniform” problem we take H to be a nearest-neighbor tight-binding model on a square lattice with hopping matrix t=1.  The chemical potential is chosen such that n = 0.8 electrons per site.  To obtain smooth curves for ρ(λ), a gaussian broadening of the levels has been introduced with δλ = 0.01.
        
        We define different disorder ensembles as follows:  1)  To represent potential disorder, we add a random single-site energy, t_(r⃗,σ),(r⃗^',σ^') = ϵ(r⃗) δ_σ,σ^'δ_r⃗,r⃗^';  when we take ϵ from a square distribution of width W we will refer to it as “Born scattering”, while the binary distribution ϵ = 100 t̅ or ϵ=0 with probabilities n_imp and 1-n_imp, respectively, will be referred to as “unitary scatterers”.  2) To represent the pairing interaction, we assume an on-site V_0(r⃗) and nearest-neighbor V_1(r⃗) density-density interaction.  The case in which we take V_0 < 0 and V_1 = 0 we will refer to as “s-wave pairing,” while the case with  V_0 > 0 and V_1 < 0 will be called “d-wave.” We consider both uniform (V_0 and V_1 independent of r⃗) and inhomogeneous pairing interactions, where for instance to represent “negative U centers” we take V_0= -  U with probability x and V_0=0 with probability 1-x.  We work in units of hopping t = 1 and take values of U = 3, 4.  Since λ scales with U, we scale the x-axis by |U| in the figures below such that the distribution is independent of our chosen value.
        
        
        
          §.§.§ Born scattering and s-wave pairing
        
        The blue-dotted curve in Fig. <ref> shows ρ(λ) for the case of uniform s-wave pairing (V_0=-U) and Born scattering with W=3 at temperature T = 0.05.  As we plot λ in units of U, the results apply for arbitrary values of U.  For comparison, the dotted vertical lines represent the maximum value of λ for the disorder free model (W=0):  λ_max(W=0,T=0.05) / |U| = 0.74.  Since the results are represented at fixed T, the meaning of the dotted curve is that for the disorder free system described by H, the critical strength of U such that T_c > T for any U > U_c(T) is determined by U_c(T) = t λ_max(W=0,T).  Manifestly, at T=0.05, the tail of the distribution in the presence of disorder is easily seen to extend well above the λ_max, meaning that disorder enhances T_c as expected.  
        
          §.§.§ Born scattering plus negative U centers
        
        The red-solid and yellow-dashed curves in Fig. <ref> show ρ(λ) for a concentration x=1/2 of negative U centers in the absence and presence of Born scattering (W=0 and W=3), respectively.  Since the average strength of the attractive interaction is V_0 = - U/2, we have also, for comparison, shown as the vertical dashed line the maximum value of λ one would obtain for a disorder free system with this average value of the effective interaction.  The integral of the distribution is normalized to unity, though half of the eigenvalues are exactly zero (not shown) when x = 1/2.  The tails of the distribution extend above the dashed line showing that disorder enhances T_c.
        
        There is additional structure to the distribution with no Born scattering.  For example, we have checked that the separated peak at λ / |U| ∼ 0.5 appears for system sizes between 24 ≤ L ≤ 50 with over 1000 disorder configurations suggesting the peak is not due to insufficient statistics.  Such peaks in the distribution are associated with localized solutions associated with certain distinct sorts of local structures.  Similar peaks appear for x=1/10 where they can be associated with clusters of 6–8 negative U centers.  
         
        
        
          §.§.§ d-wave pairing with Unitary and Born disorder
        
        Disorder-averaged ρ(λ) at a temperature T = 0.05, electron density n=0.8, and uniform interactions of a sort that favor d-wave pairing, V_0 = +U > 0 and V_1 = -1.5 V_0, are shown in Fig. <ref>.  Here the blue-solid line is for Born scattering with W = 3, while the red-dashed line is for a concentration n_imp = 10% of unitary scatters.  For comparison, the dashed and dotted vertical lines correspond, respectively, to the maximum and minimum eigenvalues of the disorder free reference system (W=0 and n_imp=0), λ_max(W=0,n_imp=0,T=0.05) / |U|= 1.83 and λ_min(W=0,n_imp=0,T=0.05) / |U|= -0.77.  
        
        Notice that although ρ(λ) is very small by the time λ is comparable to λ_max, in both the Born and the unitary case the distribution extends beyond this point.  Thus, even in the case of “d-wave” pairing, the mean-field T_c is enhanced by disorder, as promised.  As in the previous example, there is additional structure to the distribution superimposed on a generally rapid decay at large λ.  For instance, we have checked that the multiple peaks seen in the Born case at the largest values of λ are not much affected by system size nor with the change in the number of disorder configurations over which the average is taken.  Here, the corresponding eigenstates are moderately well localized, and the peaks represent structures associated with particular classes of local disorder configurations. 
        
        If, as assumed, V_0>0, the solutions with negative λ are of no physical significance.  However, the same calculation applies to the case of V_0<0, where the model corresponds to an attractive on-site interaction with a weaker repulsive nearest-neighbor repulsion.  In this case, the role of positive and negative λ are interchanged.  Of course, in this case, the pairing is s-wave, and not surprisingly the results are quite reminiscent of these for the negative U problem.  Here, we can see an even more vivid example of local structures favoring highly localized large λ solutions;  the eigenstates corresponding to the peak in the distribution for the unitary scatterers at λ / |U| ∼ -2.5 are typically localized in a radius of only a few lattice sites.
         
        
        
        
        
        
        § RELATION TO OTHER WORK
        
        Needless to say, the effect of disorder on the superconducting transition temperature,  the character of superconducting fluctuations, and the structure of the superconducting state have been the subject of extensive theoretical investigation for many decades.  Since for the most part, these studies had in mind large coherence-length superconductors, the effects of inhomogeneous pairing were (rightly) largely neglected in these studies.  However, more recently, especially since the discovery of cuprate high temperature superconductivity, a number of studies have highlighted various circumstances in which inhomogeneous pairing correlations play a significant role in superconductors with short correlation lengths  – see, for example, Refs. podolsky_2005, yukalov_1995, yukalov_2004, garcia_2014, scalettar_2006, trivedi_1998, trivedi_2001, nunner_2005, andersen_2006.   These effects are particularly amplified near a T=0 (quantum) superconductor to metal<cit.> or superconductor to insulator<cit.> transition.  Circumstances in which disorder enhances T_c have likewise been discussed previously.<cit.>  For instance, in a semimetal with weak attractive interactions, T_c is zero in the absence of disorder, but non-zero in the presence of disorder.<cit.>  However, as far as we know, that universal Lifshitz tails give rise to a broad distribution of local mean-field T_c's, and the implied importance of thermal (classical) phase fluctuations for the thermal transition in short-coherence length superconductors, has not been sharply articulated previously.
        
        
        
        
        
        
        
        
        
        
        
        
        
        § POSSIBLE RELEVANCE TO THE CUPRATES
        
        
        Of course, the best known short-coherence length superconductors are the cuprates in which, depending somewhat on the range of doping and the method used to make the estimates, the superconducting coherence length is thought to vary from a few lattice constants up to perhaps a dozen lattice constants.  Moreover, the fact that these materials are highly quasi-2D, means that the large factor of (ξ_0/a)^d in the exponent of Eq. <ref> has d=2 rather than d=3 which enhances the range of temperatures over which the present considerations are relevant.  With that in mind, we note that there are several features of the body of experimental observations that indeed may reflect a degree of self-organized granularity in the superconducting phenomena.
        
        Inhomogeneity in the gap has been observed in scanning tunneling microscopy measurements of peaks in the local density of states in BSCCO.<cit.>
         Significantly, it has been observed<cit.> that regions that locally have large values of the gap at low T have a local gap onset temperature that is larger than average, roughly in proportion to the local value of the low T gap.  This is suggestive of the sort of local variations in the pair-wave-function discussed above.  Moreover, this correlates with the observation that T_c is typically determined by the superfluid stiffness, rather than the pairing scale.<cit.>  On the other hand, the above mentioned results primarily concern the underdoped cuprates, where even the “normal” state above T_c deviates dramatically from the usual Fermi liquid metal phase on which BCS theory is based.  Below the critical temperature, however, mean-field d-wave pairing inhomogeneity can still reproduce a relatively homogeneous LDOS at sufficiently low energy<cit.> such that the proposed granularity above the critical temperature does not directly conflict with local spectroscopy in BSCCO observing a homogeneous superconducting state.<cit.> 
         
        Recently, interesting experiments have begun to systematically probe the overdoped regime, where the normal state is at least more nearly Fermi-liquid like.  Somewhat unexpectedly, it is found that the same relation between T_c and the superfluid density persists, even as the quantum critical doping at which T_c vanishes is approached.<cit.>  Moreover, there is evidence both from optics<cit.> and specific heat<cit.> that even deep in the superconducting state, a large density of apparently normal metallic quasiparticles persists, with a density that approaches that of the normal state as T_c → 0.   We therefore consider it likely that these experiments reflect the sort of self-organized granularity that we have described here, though clearly local probe measurements are necessary to confirm this.  In making this suggestion, we wish to stress that what we are talking about is an intrinsic feature of short-coherence-length superconductors in the presence of statistically homogeneous disorder;  it has nothing to do with any large scale structural or chemical inhomogeneities that are typically meant when one discusses “sample inhomogeneities.”[A  superficially very different interpretation of these same experiments has been presented in Refs. broun_2017 and hirschfeld_2018. Here, neglecting the role of inhomogeneous pairing correlations, the effective medium theory of disordered d-wave BCS superconductors has been deployed and shown to rather naturally account for many of the salient features  of the experiments below T_c  using plausible model parameters.  It is unclear to what extent the two approaches are truly at odds.  Both approaches start from the perspective that BCS mean-field theory is a plausible first approximation, and that the effects of disorder play an essential role even in the cleanest samples.  It is likely that   the induced self-organized inhomogeneity on which we have focused  is most significant above and in the vicinity of T ≈ T_c;  below T_c, the familiar proximity effect tends to produce a more uniform electronic structure even in the presence of inhomogeneous pairing interactions.  However, close enough to the quantum critical point at which T_c→ 0, a variety of “phase-sensitive” effects could give rise to unambiguous macroscopic signatures of self-organized granularity, even in the superconducting state.<cit.>]  We acknowledge useful discussions with D. Broun, J. C. Seamus Davis, P. Hirschfeld, and N. Trivedi.  This work was supported in part by the Department of Energy, Office of Basic Energy Sciences, under contract no. DE-AC02-76SF00515 at Stanford.
        
        99nunner_2005
        T. S. Nunner et al.95, 177003 (2005).
        
        andersen_2006
        B. M. Andersen et al.74, 060501 (2006).
        
        anderson_1959
        P. W. Anderson, J. Phys. Chem. Solids 11, 26 (1959).
        
        AG_1959
        A. A. Abrikosov and L. P. Gor'kov, JETP 9, 220 (1959).
        
        spivak_2008
        B. Spivak, P. Oreto, and S. A. Kivelson, 77, 214523 (2008).
        
        
        lifshitz_1965
        I. M. Lifshitz, Sov. Phys. Usp. 7, 549 (1965).
        
        halperin_1966
        B. I. Halperin and M. Lax, Phys. Rev. 148, 722 (1966).
        
        zittartz_1966
        J. Zittartz and J. S. Langer, Phys. Rev. 148, 742 (1966).
        
        cardy
        J. L. Cardy, J. Phys. C 11, L321 (1978).
        
        yaida
        S. Yaida, 93, 075120 (2016).
        
        
        podolsky_2005
        I. Martin, D. Podolsky, and S. A. Kivelson, 72, 060502 (2005).
        
        yukalov_1995
        A. J. Coleman, E. P. Yukalova, V. I. Yukalov, Physica C 243, 76 (1995).
        
        yukalov_2004
        V. I. Yukalov and E. P. Yukalova, 70, 224516 (2004).
        
        garcia_2014
        J. Mayoh and A. M. Garcia-Garcia, 90, 134513 (2014).
        
        scalettar_2006
        K. Aryanpour et al.73, 104518 (2006).
        
        trivedi_1998
        A. Ghosal, M. Randeria, and N. Trivedi, 81, 3940 (1998).
        
        trivedi_2001
        A. Ghosal, M. Randeria, and N. Trivedi, 65, 014501 (2001).
        
        larkin_1998
        M. V. Feigelman and A. I. Larkin, Chem. Phys. 235, 107 (1998).
        
        trivedi_2011
        K. Bouadim et al., Nature Physics 7, 884 (2011).
        
        garcia_2015
        J. Mayoh and A. M. Garcia-Garcia, 92, 174526 (2015).
        
        feigelman_2007
        M. V. Feigel'man et al.98, 027001 (2007).
        
        mirlin_2012 
        I. S. Burmistrov, I. V. Gornyi, and A. D. Mirlin, 108, 017002 (2012)
        
        mirlin_2015
        I. S. Burmistrov, I. V. Gornyi, and A. D. Mirlin, 92, 014506 (2015).
        
        franz_1997
        M. Franz et al.56, 7882 (1997).
        
        walker_1998
        M. E. Zhitomirsky and M. B. Walker, 80, 5413 (1998).
        
        gastiasoro_2017
        M. N. Gastiasoro and B. M. Andersen, “Enhancing Superconductivity by Disorder,” arXiv:1712.02656 (2017).
        
        romer_2017
        A. T. Romer, P. J. Hirschfeld, and B. M. Andersen, “Boosting T_c with Disorder in Spin-Fluctuation Mediated Unconventional Superconductors,” arXiv:1712.07914 (2017).
        
        
        
        sondhi_2013
        R. Nandkishore et al.87, 174511 (2013).
        
        sondhi_2014
        I. D. Potirniche et al.90, 094516 (2014).
        
        
        howald_2001
        C. Howald, P. Fournier, and A. Kapitulnik, 64, 100504 (2001).
        
        mcelroy_2005  
        K. McElroy et al. Science 309, 1048 (2005).
        
        yazdani_2007
        K. K. Gomes et al. Nature 447, 569 (2007).
        
        emerykivelson
        V. J. Emery and S. A. Kivelson, Nature 374, 434 (1995).
        
        vishik_2012
        I. M. Vishik et al., Proc. Natl. Acad. Sci. 109, 18332 (2012).
        
        zhong_2018
        Y. G. Zhong et al., “Continuous doping of a cuprate surface: new insights from in-situ ARPES,” arXiv:1805.06450 (2018).
        
        fang_2006
        A. C. Fang et al.96, 017007 (2006).
        
        seamus_2012
        K. Fujita et al. J. Phys. Soc. Jap. 81, 011005 (2012).
        
        seamus_2015
        K. Fujita et al. (2015) “Spectroscopic Imaging STM: Atomic-Scale Visualization of Electronic Structure and Symmetry in Underdoped Cuprates,” In: Avella A., Mancini F. (eds) Strongly Correlated Systems. Springer Series in Solid-State Sciences, 180. Springer, Berlin, Heidelberg
        
        
        
        
        bozovic_2016
        I. Bozovic et al. Nature 536, 309 (2016).
        
        armitage_2018
        F. Mahmood, X. He, I. Bozovic, and N. P. Armitage, “Locating the missing superconducting electrons in overdoped cuprates,” arXiv:1802.02101 (2018).
        
        hhwen_2004
        H. H. Wen et al.70, 214505 (2004).
        
        broun_2017
        N. R. Lee-Hone, J. S. Dodge, and D. M. Broun, 96, 024501 (2017);  Erratum-ibid. 97, 219903 (2018).
        
        hirschfeld_2018
        N. R. Lee-Hone et al.“Optical conductivity of overdoped cuprate superconductors: application to LSCO,” arXiv:1802.10198 (2018).
        
        
        
        
        
        
        
        
        spivakandme  S. A. Kivelson and B. Spivak,  92, 184502 (2015).
        
         
         Supplementary Material
        
        
        
        § LINEARIZED GAP EQUATION
        
        
        The spectrum of M is computed numerically for different disorder realizations to extract the effect of rare events on mean-field T_c.  The disordered problem is diagonalized H = ∑_n E_n γ^†_n γ_n using c_r σ = ∑_n U_r σ ,  n ^* γ_n and Green's function
        
            ⟨ c_a σ(τ) c^†_b λ⟩ = ∑_n  U^*_a σ , n  U_b λ , n  e^ - τ E_n  ( 1 - f_n )
        
        where the columns of U_r σ , n are eigenvectors and f_n = (1+e^β E_n)^-1.  For singlet operator defined as ψ_a b = (c_a ↓ c_b ↑ + c_b ↓ c_a ↑ ) / 2, the pair susceptibility is
        
            χ_a b ; c d = ∫_0^β d τ ⟨ψ_a b(τ) ψ^†_c d⟩  
            = 1/2∫_0^β d τ ⟨  c_a ↓(τ) c^†_c ↓⟩⟨ c_b ↑(τ) c^†_d ↑⟩  + ⟨  c_a ↑(τ) c^†_d ↑⟩⟨  c_b ↓(τ) c^†_c ↓⟩
        
            = 1/2∑_ n , m = 1 ^ 2 N_x N_y ( 1 - f_n - f_m / E_n + E_m ) [ U^*_a ↓ , n  U_ c ↓ , n  U^*_b ↑ , m  U_ d ↑ , m  + U^*_a ↑ , n  U_d ↑ , n  U^*_b ↓ , m  U_ c ↓ , m ]
        
        for sites a,b,c,d = 1, ⋯ , N_x N_y.  The interaction term has the form
        
        H_int = ∑_r V_0(r) c^†_r ↑ c^†_r ↓ c_r ↓ c_r ↑ + ∑_⟨ r , r' ⟩ V_1(r) n_r n_r'  
        where n_r = ∑_σ c^†_r σ c_r σ, and the inhomogeneous on-site interaction V_0(r) and nearest neighbor interaction V_1(r) are chosen to favor either s-wave (V_0 < 0 and V_1 = 0) or d-wave (V_0 > 0 and V_1 < 0) solutions.  The simplified eigenvalue equation (<ref>) is solved for (5 N_x N_y) × (5 N_x N_y) matrix M_r , r + δ ; r' , r' + δ'  =  V_r,r+δχ_r , r+δ ; r' , r' + δ' for δ, δ' = { (0,0) , (± 1, ± 1) } and interactions V_r,r = V_0(r) and V_r,r+δ = V_1(r).  
        
        
        
        
        
        
        § S-WAVE PAIRING AT HIGHER TEMPERATURE
        
        
        The disorder-averaged distribution ρ(λ) at a higher temperature T = t is shown in Fig. <ref> for the same disorder ensembles shown in Fig. <ref>.  The maximum value of λ for the disorder free model (W = 0) is λ_max(W=0,T=1)/|U| = 0.2.  For the high temperature case, the range of M has dropped to being on the order of one lattice site, the fall-off of distribution is sufficiently sharp that the probability has dropped below the limits accessible to our numerical results before the clean-limit value is reached.  The increasingly fast drop of the tails of the distribution with increasing temperature is to be expected from the above considerations;  our failure to see the tail of the distribution at high temperatures extending above λ_max is, we believe, an artifact but one that could only be corrected by keeping a much larger number of disorder configurations.  
        For the case with negative U centers (x=1/2), the entire observed distribution lies above the dashed line. 
        
        
        
        
        § INVERSE PARTICIPATION RATIO
        
        
        The degree to which the gap solution is localized can be quantified using the inverse participation ratio I_λ for a given eigenvalue λ defined as
        
        I_λ = ⟨ |F_λ(r⃗)|^4 ⟩/⟨ |F_λ(r⃗)|^2 ⟩^2
        where ⟨⋯⟩ denotes a spatial average and F_λ(r⃗) is the eigenvector of the linearized BCS equation.  The inverse participation ratio has the property that I_λ∼ 1 / L^d → 0 as linear system size L →∞ for extended states, and I_λ∼ 1/ξ_loc^d for states localized within a characteristic length ξ_loc as L →∞.  The inverse participation ratio therefore provides a length scale 1/√(I_λ) which gives an estimate n_sc / I_λ of the degree to which superconducting solutions are isolated from each other in two dimensions.
        
         Negative U centers: The inverse participation ratio I_λ decreases with increasing λ meaning the solutions with the largest eigenvalues are most delocalized.  The characteristic length scale 1/√(I_λ) increases smoothly to roughly L/2 at a value λ∼ 1.8 at which point there is a jump in both λ and 1/√(I_λ) (visible as a peak at λ / V_0 ≈ 0.5).  
        
         Unitary scattering with on-site interaction:  Different types of disorder are more conducive to local enhancements of the density of states and therefore larger eigenvalues of M.  The results for V_0 = -4 and V_1 = 0 in the presence of pairing inhomogeneity x = 1/2 and a concentration n_imp = 10% of randomly located unitary scatterers of strength V_imp = 100 t are shown in Fig. <ref> at a temperature T = 0.05 t.  The inverse participation ratio for each of the eigenvectors is plotted (red circles) on top of the tail end of the ρ(λ) distribution (blue).  The vertical dashed line (black) corresponds to the maximum eigenvalue λ_max,0 for the clean system at T = 0.05 t.  Above λ≈ 1 in the exponential tail of the distribution, the inverse participation ratio remains small corresponding to relatively extended states; however, above clean system maximum λ_max,0≈ 3, the solutions with the largest λ show an increase in IPR.  The maximum eigenvalue solutions with λ / V_0 ∼ 2.3 have I_λ = 1 meaning the gap eigenvector is localized on a single site.  From visual inspection of the solution and local density of states (see Supplementary Material, "LDOS Enhancement"), it can be seen that these solutions are associated with clusters of multiple unitary scatterers in close proximity.  There appear to be two peaks near λ≈ 5 and λ≈ 9 that can correspond to impurity bands with IPR = 1/2 (states localized on two sites) and IPR = 1 (states localized on single site) respectively, which would hybridize at sufficiently large concentration.  For λ < 1 (not shown), the IPR spans a range from zero to one; however, since the density of states ρ(λ) is large, the IPR loses meaning due to mixing in the self-consistent mean-field solution.  The distribution at higher temperature T = 1 (not shown) similarly exhibits two peaks above the clean λ_max,0.      
        
        
        
        
        
        § LDOS ENHANCEMENT
        
        
        From visual inspection of the gap solution and local density of states, it can be seen that there is an enhancement due to rare microscopic impurity configurations.  The chemical potential, local density of states, and on-site pairing solutions F_λ(r⃗) with n_imp = 10% unitary scatterers for largest (λ_1 = 9.7) and second largest (λ_2 = 3.4) eigenvalues for a particular disorder realization at temperature T = 0.05 t are shown in Fig. <ref>.  The on-site interaction is taken as V_0 = -4 and nearest neighbor interaction V_1 = 0.  The peak in the local density of states is associated with a rare impurity cluster, and the largest eigenvalue solution F_λ_1(r) has I_λ_1 = 1 corresponding to non-zero value on a single site.  The second largest eigenvalue solution F_λ_2(r⃗) for this disorder realization has I_λ_2 = 0.48 and can be seen to have finite value over a localized region.
        
        
        
        
        § POTENTIAL DISORDER WITH 
        
        NEAREST NEIGHBOR INTERACTION
        
        
        The results for random potential and interaction disorder with nearest-neighbor interaction show qualitatively similar behavior to the case V_1 = 0:  the distributions decay exponentially with the eigenvalues in the tail exceeding the uniform clean case, and increasing temperature narrows the distribution.  In contrast to the previous example, the negative U centers with V_1 ≠ 0 at high temperature T = 1 show two separated peaks:  the peak with larger eigenvalues corresponds to solutions localized in clusters where the interaction is non-zero such that the effective interaction of the puddle is U.  The peak with smaller eigenvalues (absent for V_1 = 0) corresponds to solutions localized on “checkerboard” patterns where the interaction averages to U/2 - the eigenvalues associated with this second peak are roughly half that of the larger-eigenvalue peak (though the lower end of this second peak has a range of effective interactions less than U).  This characterization is confirmed quantitatively by defining an effective gap V_1 = ∑_r V_1(r) |Δ_x(r)|^2 / ∑_r' |Δ_x(r')|^2 of the puddle which shows V_1≈ V_1 / 2 for states in the smaller-eigenvalue peak.   
        
        
        
        
        
        
        
        
        
        § LINEARIZED BCS EQ.
        
        
        The spectrum of M is computed numerically for different disorder realizations extract the effect of rare events on mean-field T_c.  The disordered problem is diagonalized as H = ∑_n, σ E_n γ^†_n σγ_n σ using c_r σ = ∑_n U_r n^* γ_n σ where the columns of U_r n are eigenvectors.  For pair operator defined as ψ_a σ_1 , b σ_2 = c_a σ_1 c_b σ_2 and Green's function
        
            G_a σ_1 , b σ_2(τ) = ⟨ c_a σ_1(τ) c^†_b σ_2⟩
            
            = δ_σ_1 σ_2∑_n=1^N_x N_y U^*_a n U_n b (1-f_n) e^- τ E_n
        
        the pair susceptibility χ including spin indices is
        
            χ_b' σ_2' , a' σ_1' , c σ_3 , d σ_4 = ∫_0^β d τ ⟨ψ_b' σ_2' , a σ_1'(τ) ψ^†_c σ_3 , d σ_4⟩  
            =  ∫_0^β d τ ⟨  c_b' σ_2'(τ) c^†_c σ_3⟩⟨ c_a' σ_1'(τ) c^†_d σ_4⟩ - ⟨  c_a' σ_1'(τ) c^†_c σ_3⟩⟨ c_b' σ_2'(τ) c^†_d σ_4⟩
        
            =  ∑_n,m=1^N_x N_y( 1-f_n-f_m/E_n+E_m) [ δ_σ_2' σ_3δ_σ_1' σ_4 U^*_b' n U_n c U^*_a' m U_m d  -   δ_σ_1' σ_3δ_σ_2' σ_4 U^*_a' n U_n c U^*_b' m U_m d]
        
        on sites a',b',c,d = 1, ⋯ , N_x N_y.  The interaction term chosen to favor s/d-wave solutions has the form
        
        H_int = ∑_r V_0(r) c^†_r ↑ c^†_r ↓ c_r ↓ c_r ↑ + ∑_⟨ r , r' ⟩ V_1(r) n_r n_r'  
        where n_r = ∑_σ c^†_r σ c_r σ is the site occupancy, V_0(a) is the inhomogeneous on-site interaction, and V_1(r) is the inhomogeneous nearest-neighbor interaction.  Denoting σ = - σ, the corresponding interaction matrix is
        
            V_a σ_1 ,  b σ_2 , b' σ_2' , a' σ_1 ' = ( δ_σ_1 ' σ_1δ_a' aδ_σ_2' σ_2δ_b' b ) ×
            [ V_0(a)/2δ_abδ_σ_1 σ_2 + V_1(a) (δ_σ_1 σ_2 + δ_σ_1 σ_2 ) ∑_⟨Δ⟩δ_b,a+Δ]
        
        where ⟨Δ⟩ are nearest neighbors.  The matrix M is then
        
            M_a σ_1 , b σ_2 , c σ_3 , d σ_4 = V_0(a)/2δ_a bδ_σ_1 σ_2 ( δ_σ_1 σ_3δ_σ_1 σ_4 - δ_σ_1 σ_3δ_σ_1 σ_4 ) ∑_n,m( 1-f_n-f_m/E_n+E_m) U^*_a n U_n c U^*_a m U_m d
        
            + V_1(a) ( ∑_⟨Δ⟩δ_b,a+Δ) ∑_n,m( 1-f_n-f_m/E_n+E_m) [ δ_σ_1 σ_2δ_σ_1 σ_3δ_σ_1 σ_4( U^*_b n U_n c U^*_a m U_m d - U^*_a n U_n c U^*_b m U_m d)
        
            +  δ_σ_1 σ_2(  δ_σ_1 σ_3δ_σ_1 σ_4  U^*_b n U_n c U^*_a m U_m d  - δ_σ_1 σ_3δ_σ_1 σ_4  U^*_a n U_n c U^*_b m U_m d) ] ]
