text:
# Pairing fluctuation theory of high 𝑇_𝑐 superconductivity in the presence of nonmagnetic impurities

## I Introduction

The pseudogap phenomena in high $`T_c`$ superconductors have been a great challenge to condensed matter physicists since over a decade ago. These phenomena manifestly contradict BCS theory by, e.g, presenting a pseudo excitation gap in single particle excitation spectrum. Yet the origin of the pseudogap and, in general, the mechanism of the superconductivity are still not clear. Many theories have been proposed, which fall into two classes, based on whether the pseudogap has a superconducting origin. Some authors propose that the pseudogap may not be related to the superconductivity; instead, it is associated with another ordered state, such as the antiferromagnetism related resonating valence bond (RVB) state,Anderson (1987) d-density waveChakravarty et al. (2001) and spin density wave order.Demler et al. (2001) On the other hand, many others believe that the pseudogap has the same origin as the superconductivity, such as the phase fluctuation scenario of Emery and KivelsonEmery and Kivelson (1995) and the various precursor superconductivity scenarios.Randeria (1995); Maly et al. (1999a, b); Janko et al. (1997); Chen et al. (1998); Kosztin et al. (1998) Previously, Chen and coworkers have worked out, within the precursor conductivity school, a pairing fluctuation theoryChen et al. (1998, 1999, 2000a) which enables one to calculate quantitatively physical quantities such as the phase diagram, the superfluid density, etc. for a clean system. In this theory, two-particle and one-particle Green’s functions are treated on an equal footing, and equations are solved self-consistently. Finite center-of-mass momentum pair excitations become important as the pairing interaction becomes strong, and lead to a pseudogap in the excitation spectrum. In this context, these authors have been able to obtain a phase diagram and calculate the superfluid density, in (semi)quantitative agreement with experiment.

However, to fully apply this theory to the cuprates, we need to extend it to impurity cases, since impurities are present even in the cleanest samples of the high $`T_c`$ materials, such as the optimally doped YBa2Cu3O7-δ (YBCO) single crystals. In addition, this is necessary in order to understand the finite frequency conductivity issue. Furthermore, study of how various physical quantities respond to impurity scattering may help to reveal the underlying mechanism of the superconductivity. For example, it can be used to determine whether the pseudogap has a superconducting origin.Kruis et al. (2001) Particularly, it is important to address how $`T_c`$ and the pseudogap itself vary with impurity scattering, especially in the underdoped regime. To this end, one needs to go beyond BCS theory and include the pseudogap as an intrinsic part of the theory. Due to the complexity and technical difficulties of this problem, there has been virtually no work in the field on this important problem.

Among all physical quantities, the density of states (DOS) $`N(\omega )`$ close to the Fermi level ($`\omega =0`$) is probably most sensitive to the impurities. Yet different authors have yielded contradictory results in this regard. BCS-based impurity $`T`$ matrix calculations predict a finite DOS at $`\omega =0`$,Hirschfeld and Vollhardt (1986); Hirschfeld et al. (1988) which has been used to explain the crossover from $`T`$ to $`T^2`$ power law for the low temperature superfluid density.Hirschfeld and Goldenfeld (1993) Nonperturbative approaches have also been studied and have yielded different results. Senthil and FisherSenthil and Fisher (1999) find that DOS vanishes according to universal power laws, Pépin and LeePépin and Lee (2001) predict that $`N(\omega )`$ diverges as $`\omega 0`$, assuming a strict particle-hole symmetry, and ZieglerZiegler et al. (1998) and coworkers’ calculation shows a rigorous lower bound on $`N(\omega )`$. Recently, Atkinson et al.Atkinson et al. (2000a) try to resolve these contradictions by fine-tuning the details of the disorders. Nevertheless, all these calculations are based on BCS theory and cannot include the pseudogap in a self-consistent fashion, and thus can only be applied to the low $`T`$ limit in the underdoped cuprates. Therefore, it is necessary to extend the BCS-based calculations on impurity issues to include the pseudogap self-consistently.

In this paper, we extend the pairing fluctuation theory from clean to impurity cases. Both the impurity scattering and particle-particle scattering $`T`$ matrices are incorporated and treated self-consistently. This goes far beyond the usual self-consistent (impurity) $`T`$-matrix calculations at the BCS level by, e.g., Hirschfeld and others.Hirschfeld and Vollhardt (1986); Hirschfeld et al. (1988) In this context, we study the evolution of $`T_c`$ and various gap parameters as a function of the coupling strength, the impurity concentration, the hole doping concentration and the impurity scattering strength. In addition, we study not only the Born and the unitary limits, but also at intermediate scattering strength. We find that the real part of the frequency renormalization can never be set to zero, the chemical potential adjusts itself with the impurity level. As a consequence, the positive and negative strong scattering limits do not meet. The residue density of states at the Fermi level is generally finite at finite impurity concentrations, in agreement with what has been observed experimentally. Both $`T_c`$ and the total excitation gap decrease with increasing impurity level, as one may naively expect. Born scatterers are about twice as effective as unitary scatters in suppressing $`T_c`$ and the gap. In the unitary limit, the zero temperature superfluid density decreases faster with $`n_i`$ when $`n_i`$ is still small, whereas in the Born limit, it is the opposite. At given $`T<T_c`$, both the order parameter and the total gap are suppressed, but the pseudogap is not. Finally, incoherent pair excitations contribute an additional $`T^{3/2}`$ term to the low $`T`$ temperature dependence of the superfluid density, robust against impurity scattering.

In the next section, we first review the theory in a clean system, and then present a theory at the Abrikosov-Gor’kov level. Finally, we generalize it to include the full impurity $`T`$-matrix, in addition to the particle-particle scattering $`T`$-matrix, in the treatment, and obtain a set of three equations to solve for $`\mu `$, $`T_c`$ and various gaps. In Sec. III, we present numerical solutions to these equations. We first study the effects of impurity scattering on the density of states, then study the effects on $`T_c`$ and the pseudogap at $`T_c`$, followed by calculations of the effects on the gaps and the superfluid density below $`T_c`$. Finally, we discuss some related issues, and conclude our paper.

## II Theoretical Formalism

The excitation gap forms as a consequence of Cooper pairing in BCS theory, while the superconductivity requires the formation of the zero-momentum Cooper pair condensate. As these two occurs at the same temperature in BCS theory, one natural way to extend BCS theory is to allow pair formation at a higher temperature ($`T^{}`$) and the Bose condensation of the pairs at a lower temperature ($`T_c`$). Therefore, these pairs are phase incoherent at $`T>T_c`$, leading to a pseudogap without superconductivity. This can nicely explain the existence of the pseudogap in the cuprate superconductors. Precursor superconductivity scenarios, e.g., the present theory, provides a natural extension of this kind. At weak coupling, the contribution of incoherent pairs is negligible and one thus recovers BCS theory, with $`T^{}=T_c`$. As the coupling strength increases, incoherent pair excitations become progressively more important, and $`T^{}`$ can be much higher than $`T_c`$, as found in the underdoped cuprates. In general, both fermionic Bogoliubov quasiparticles and bosonic pair excitations coexist at finite $`T<T_c`$.

### II.1 Review of the theory in a clean system

The cuprates can be modeled as a system of fermions which have an anisotropic lattice dispersion $`ϵ_𝐤=2t_{}(2\mathrm{cos}k_x\mathrm{cos}k_y)+2t_{}(1\mathrm{cos}k_z)\mu `$, with an effective, short range pairing interaction $`V_{𝐤,𝐤^{}}=g\phi _𝐤\phi _𝐤^{}`$, where $`g<0`$. Here $`t_{}`$ and $`t_{}`$ are the in-plane and out-of-plane hopping integrals, respectively, and $`\mu `$ is the fermionic chemical potential. For the cuprates, $`t_{}t_{}`$. The Hamiltonian is given by

$`^0`$ $`=`$ $`{\displaystyle \underset{𝐤\sigma }{}}ϵ_𝐤c_{𝐤\sigma }^{}c_{𝐤\sigma }^{}`$ (1)
$`+{\displaystyle \underset{\mathrm{𝐤𝐤}^{}𝐪}{}}V_{𝐤,𝐤^{}}c_{𝐤+𝐪/2}^{}c_{𝐤+𝐪/2}^{}c_{𝐤^{}+𝐪/2}^{}c_{𝐤^{}+𝐪/2}^{},`$
The pairing symmetry is given by $`\phi _𝐤=1`$ and $`(\mathrm{cos}k_x\mathrm{cos}k_y)`$ for $`s`$\- and $`d`$-wave, respectively. Here and in what follows, we use the superscript “0” for quantities in the clean system, to be consistent with the notations for the impurity dressed counterpart below. For brevity, we use a four-momentum notation: $`K=(𝐤,i\omega ),_K=T_{𝐤,\omega }`$, etc.

To focus on the superconductivity, we consider only the pairing channel, following early work by Kadanoff and Martin;Kadanoff and Martin (1961) the self-energy is given by multiple particle-particle scattering. The infinite series of the equations of motion are truncated at the three-particle level $`G_3`$, and $`G_3`$ is then factorized into single- ($`G`$) and two-particle ($`G_2`$) Green’s functions. The final result is given by the Dyson’s equations for the single particle propagator \[Refer to Refs. (Chen et al., 1998, 2000a) for details\]

$`\mathrm{\Sigma }^0(K)`$ $`=`$ $`G_{0}^{0}{}_{}{}^{1}(K)G_{}^{0}{}_{}{}^{1}(K)`$ (2a)
$`=`$ $`{\displaystyle \underset{Q}{}}t^0(Q)G_0^0(QK)\phi _{𝐤𝐪/2}^2,`$
and the $`T`$ matrix (or pair propagator)
$$t^0(Q)=t_{sc}^0(Q)+t_{pg}^0(Q),$$
(2b)
with
$$t_{sc}^0(Q)=\frac{\mathrm{\Delta }_{sc}^2}{T}\delta (Q),$$
(2c)
where $`\mathrm{\Delta }_{sc}0`$ at $`TT_c`$, and
$$t_{pg}^0(Q)=\frac{g}{1+g\chi ^0(Q)},$$
(2d)
where $`\mathrm{\Delta }_{sc}`$ is the superconducting order parameter, $`G_0^0(K)=1/(i\omega ϵ_𝐤)`$ is the bare propagator. and
$$\chi ^0(Q)=\underset{K}{}G^0(K)G_0^0(QK)\phi _{𝐤𝐪/2}^2$$
(2e)
is the pair susceptibility. This result can be represented diagrammatically by Fig. 1. The single and double lines denote the bare and full Green’s functions, respectively, and the wiggly double lines denote the pair propagator.

The superconducting instability is given by the Thouless criterion

$$1+g\chi ^0(0)=0,(TT_c),$$
(3)
which leads to the approximation

$`\mathrm{\Sigma }_{pg}^0(K)`$ $`=`$ $`{\displaystyle \underset{Q}{}}t_{pg}^0(Q)G_0^0(QK)\phi _{𝐤𝐪/2}^2`$
$``$ $`\mathrm{\Delta }_{pg}^2G_0^0(K)\phi _𝐤^2,`$
where the pseudogap is defined by

$$\mathrm{\Delta }_{pg}^2=\underset{Q0}{}t_{pg}^0(Q).$$
(5)
As a consequence, the self-energy takes the standard BCS form

$$\mathrm{\Sigma }(K)=\mathrm{\Delta }^2G_0^0(K)\phi _𝐤^2=\mathrm{\Delta }_𝐤^2G_0^0(K);$$
(6)
where $`\mathrm{\Delta }^2=\mathrm{\Delta }_{sc}^2+\mathrm{\Delta }_{pg}^2`$, and $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }\phi _𝐤`$. In this way, the full Green’s function $`G^0(K)`$ also takes the standard BCS form, with the quasiparticle dispersion given by $`E_𝐤=\sqrt{ϵ_k^2+\mathrm{\Delta }^2\phi _𝐤^2}`$. So does the excitation gap equation

$$1+g\underset{𝐤}{}\frac{12f(E_𝐤)}{2E_𝐤}\phi _𝐤^2=0.$$
(7)
We emphasize that although this equation is formally identical to its BCS counterpart, the $`\mathrm{\Delta }`$ here can no longer be interpretted as the order parameter as $`\mathrm{\Delta }_{pg}0`$ in general. For self-consistency, we have the fermion number constraint

$$n=2\underset{K}{}G^0(K)=2\underset{𝐤}{}\left[v_𝐤^2+\frac{ϵ_𝐤}{E_𝐤}f(E_𝐤)\right],$$
(8)
The gap equation Eq. (7), the number equation Eq. (8), and the pseudogap parametrization Eq. (5) form a complete set, and can be used to solve self-consistently for $`T_c`$, $`\mu (T_c)`$, and $`\mathrm{\Delta }(T_c)`$ by setting $`\mathrm{\Delta }_{sc}=0`$, or $`\mu (T)`$, $`\mathrm{\Delta }(T)`$, and $`\mathrm{\Delta }_{pg}(T)`$ at given $`T<T_c`$. Here $`v_𝐤^2=\frac{1}{2}(1ϵ_𝐤/E_𝐤)`$, as in BCS.

### II.2 Impurity scattering at the Abrikosov-Gor’kov level

For simplicity, we restrict ourselves to nonmagnetic, elastic, isotropic $`s`$-wave scattering. At the same time, we will keep the derivation as general as possible. In the presence of impurities of concentration $`n_i`$, the Hamiltonian is given by

$$=^0+_I,$$
(9)
where in the real space the impurity term is given by

$$_I=\underset{i}{}𝑑\text{x}u(\text{x}\text{x}_i)\psi ^{}(\text{x})\psi (\text{x}),$$
(10)
with $`u(\text{x})=u\delta (\text{x})`$ for isotropic $`s`$-wave scattering.

To address the impurity scattering, we begin at the Abrikosov-Gor’kov (AG) level,Abrikosov and Gor’kov (1959, 1960); AGD  which is a good approximation in the Born limit. Following AG, we include all possible configurations of impurity dressing, but excluding bridging diagrams like Fig. 2(a), crossing diagrams like Fig. 2(b), and higher order terms like Fig. 2(c). The dashed lines denote impurity scattering, and the crosses denote the impurity vertices. We clarity, in most diagrams, we do not draw the fermion propagation arrows. It is understood, however, that they change direction at and only at a pairing vertex. As in subsec. II.1, we use plain double lines to denote the Green’s function ($`G^0`$) fully dressed by the pairing interaction but without impurity scattering, i.e.,

$$G^0(K)=\frac{1}{i\omega ϵ_𝐤\mathrm{\Sigma }^0(K)},$$
(11)
where

$$\mathrm{\Sigma }^0(K)=\underset{Q}{}t(Q)G_0^0(QK)\phi _{𝐤𝐪/2}^2.$$
(12)
However, since we will address the impurity dressing of the pairing vertex or, equivalently, the pair susceptibility $`\chi (Q)`$, we assume the pair propagator $`t(Q)`$ in the above equation is already dressed with impurity scattering, with

$$t(Q)=t_{sc}(Q)+t_{pg}(Q),$$
(13)
$$t_{sc}(Q)=\frac{\mathrm{\Delta }_{sc}^2}{T}\delta (Q),$$
(14)
and

$$\mathrm{\Delta }_{pg}^2=\underset{Q0}{}t_{pg}(Q)=\underset{Q0}{}\frac{g}{1+g\chi (Q)},$$
(15)
as in the clean case. The shaded double lines denote impurity-dressed full Green’s function $`G`$, and “shaded” single lines denote impurity-dressed bare Green’s function (which we call $`\widehat{G}_0`$), i.e,

$$\widehat{G}_0(K)=\frac{1}{i\omega ϵ_𝐤\overline{\widehat{G}}_{0,\omega }},$$
(16)
where the bar denotes impurity average $`\overline{\widehat{G}}_{0,\omega }=n_i_𝐤^{}|u(\text{k}\text{k}^{})|^2\widehat{G}_0(K^{})`$. We use open circles to denote bare pairing vertex $`\gamma (K|Q)=\phi _{𝐤+𝐪/2}`$, and shaded circles full pairing vertex $`\mathrm{\Gamma }(K|Q)`$, where $`Q`$ is the pair four-momentum. To obtain the Feynman diagrams for the impurity dressed full $`G`$, we first expand the pairing self-energy diagram as an infinite series which contains only bare single particle Green’s function and pair propagators, and then insert all possible impurity scattering on the single particle propagators at the AG level. We assume that the pair propagators are always self-consistently dressed by the impurity scattering. After regrouping all non-impurity dressed lines on the left, the final result for the diagrams is shown in Fig. 3(a). Here following AG, the subdiagrams inside the two impurity legs are assumed to be self-consistently dressed by impurity scattering. To make direct comparison with the BCS case easier, we present the corresponding diagrams for the BCS case in Appendix A. The first term on the right hand side (RHS) of Fig. 3(a) contains all diagrams without impurity dressing (except via the pair propagators). The second term corresponds to the third term of the first equation in Fig. 18. The third term corresponds to the last term, the last term to the second. The fourth and the fifth together correspond to the fourth term \[see Fig. 4(a)\]. The fifth term in Fig. 3(a) arises since the two impurity legs can cross two separate pairing self-energy dressing parts; it can be eliminated using the equality shown in Fig. 3(b). Here the shaded elliptical region denotes self-consistent impurity dressing of the double pairing vertex structure inside the two impurity scattering legs, as shown in Fig. 4(a). It is worth pointing out that these diagrams reduce to their BCS counterpart if one removes the pairing propagators. The Dyson’s equation $`\widehat{G}_0^1(K)=G_{0}^{0}{}_{}{}^{1}(K)\overline{\widehat{G}}_{0,\omega }`$ can then be used to eliminate the fourth term in Fig. 3(a). We now obtain the greatly simplified diagrams for $`G`$ as shown in Fig. 3(c), which can be further reduced into Fig. 3(e), upon defining a reduced pairing vertex $`\mathrm{\Gamma }_0(K|Q)`$ (shaded triangles) as shown in Fig. 3(d),

$`\mathrm{\Gamma }_0(K|Q)`$ $`=`$ $`\phi _{𝐤𝐪/2}`$ (17)
$`+n_i{\displaystyle \underset{𝐤^{}}{}}|u(\text{k}\text{k}^{})|^2\widehat{G}_0(QK^{})\mathrm{\Gamma }(K^{}|Q)G(K^{}).`$
One can then read off the impurity-induced quasiparticle self-energy $`\mathrm{\Sigma }^{}(K)=G_{}^{0}{}_{}{}^{1}(K)G^1(K)`$ immediately, as shown in Fig. 3(e),

$$\mathrm{\Sigma }^{}(K)=\overline{G}_\omega +\mathrm{\Sigma }(K)\mathrm{\Sigma }_0(K),$$
(18)
where the impurity average

$$\overline{G}_\omega =n_i\underset{𝐤^{}}{}|u(\text{k}\text{k}^{})|^2G(K^{}),$$
(19)
and the “full” self-energy

$$\mathrm{\Sigma }(K)=\underset{Q}{}\mathrm{\Gamma }_0(K|Q)t(Q)\mathrm{\Gamma }(K|Q)\widehat{G}_0(QK).$$
(20)
Therefore, we have finally

$`G^1(K)`$ $`=`$ $`G_{}^{0}{}_{}{}^{1}(K)\mathrm{\Sigma }^{}(K)=i\omega ϵ_𝐤\overline{G}_\omega \mathrm{\Sigma }(K)`$ (21)
$`=`$ $`i\stackrel{~}{\omega }ϵ_𝐤\mathrm{\Sigma }(K)=G_0^1(K)\mathrm{\Sigma }(K),`$
where we have defined the renormalized frequency $`i\stackrel{~}{\omega }=i\omega \overline{G}_\omega `$ and the “bare” Green’s function

$$G_0(K)=\frac{1}{i\stackrel{~}{\omega }ϵ_𝐤}.$$
(22)
Note here $`G_0(K)\widehat{G}_0(K)`$.

Now we deal further with the pairing vertex $`\mathrm{\Gamma }(K|Q)`$. First, we notice that the impurity-dressed double-vertex structure in Fig. 3(a) can be simplified as shown in Fig. 4(a) using the approximation

$$\underset{Q}{}t(Q)f(K,Q)\left[\underset{Q}{}t(Q)\right]f(K,0),$$
(23)
where $`f(K,Q)`$ is an arbitrary slow-varying function of $`Q`$. This is due to the fact that $`t(Q)`$ diverges as $`Q0`$ at $`TT_c`$. The Dyson’s equation for the Green’s function $`G^0(K)`$ in a clean system is also used in getting the second line of Fig. 4(a). Therefore, we have approximately the impurity-dressed pairing vertex as shown in Fig. 4(b), which implies the equality shown in Fig. 4(c). Fig. 4(c) can be written as

$`\widehat{G}_0(K)\mathrm{\Gamma }(K|Q)G(QK)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }_0(K|Q)G(QK)}{\widehat{G}_0^1(K)\overline{G}_\omega +\overline{\widehat{G}}_{0,\omega }}}`$ (24)
$`=G_0(K)\mathrm{\Gamma }_0(K|Q)G(QK).`$
This result demonstrates the following important relationship:

$$\widehat{G}_0(K)\mathrm{\Gamma }(K|Q)=G_0(K)\mathrm{\Gamma }_0(K|Q),(\text{small }Q).$$
(25)
Using this relationship, now the self-energy can be simplified as follows:

$`\mathrm{\Sigma }(K)`$ $`=`$ $`{\displaystyle \underset{Q}{}}\mathrm{\Gamma }_0^2(K|Q)t(Q)G_0(QK)`$ (26)
$``$ $`\mathrm{\Delta }^2\mathrm{\Gamma }_0^2(K)G_0(K)=\stackrel{~}{\mathrm{\Delta }}_𝐤^2G_0(K),`$
where $`\mathrm{\Gamma }_0(K)\mathrm{\Gamma }(K|Q=0)`$ and $`\stackrel{~}{\mathrm{\Delta }}_𝐤=\mathrm{\Delta }\mathrm{\Gamma }_0(K)`$.

Finally, the impurity dressing of each rung \[i.e., $`\chi (Q)`$\] of the particle-particle scattering ladder diagrams is topologically identical to the impurity dressing of the pairing vertex and the two associated single particle lines. And summing up all the ladders gives the pairing $`T`$ matrix. Therefore, the pair susceptibility becomes

$`\chi (Q)`$ $`=`$ $`{\displaystyle \underset{K}{}}\mathrm{\Gamma }(K|Q)\widehat{G}_0(QK)G(K)\phi _{𝐤𝐪/2}`$ (27)
$`{\displaystyle \underset{K}{}}\mathrm{\Gamma }_0(K|Q)G_0(QK)G(K)\phi _{𝐤𝐪/2}.`$
And the gap equation is given by

$$1+g\chi (0)=0=1+g\underset{K}{}\mathrm{\Gamma }_0(K)G(K)G_0(K)\phi _𝐤.$$
(28)
This result can be easily verified to be consistent with the self-consistency condition. Define formally the generalized Gor’kov $`F`$ function:

$$F^{}(K)\mathrm{\Delta }\widehat{G}_0(K)\mathrm{\Gamma }(K)G(K).$$
(29)
Using Eq. (25), we have

$$F^{}(K)\mathrm{\Delta }\mathrm{\Gamma }_0(K)G_0(K)G(K)=\stackrel{~}{\mathrm{\Delta }}_𝐤G_0(K)G(K).$$
(30)
The formal difference between this $`F^{}`$ and that in BCS is that $`\phi _𝐤`$ is now replaced by the renormalized vertex $`\mathrm{\Gamma }_0(K)`$. One immediately sees that the condition

$$\mathrm{\Delta }_𝐤g\underset{K^{}}{}\phi _𝐤\phi _𝐤^{}F^{}(K^{})$$
(31)
is consistent with the gap equation Eq. (28). However, it should be emphasized that the $`F`$ function so defined does not vanish above $`T_c`$ in the pseudogap regime, different from the BCS case.

### II.3 Impurity scattering beyond the AG level

In this subsection, we include both the impurity scattering $`T`$-matrix with the particle-particle scattering $`T`$-matrix, and, thus, go beyond the AG level. We notice that if one replaces the second-order impurity scattering subdiagrams at the AG level with the corresponding impurity $`T`$-matrices, as shown in Fig. 5, the derivation for $`G(K)`$ goes through formally without modification. Now we only need to determine the impurity $`T`$-matrices $`T_\omega `$ and $`T_\mathrm{\Delta }^{}`$ (as well as their complex conjugate) in terms of their AG-level counterpart, $`\overline{G}_\omega `$ and $`\overline{F}_\omega ^{}=n_i_𝐤^{}|u(𝐤𝐤^{})|^2F^{}(K^{})`$, respectively. In other words, except that $`i\stackrel{~}{\omega }`$ and $`\stackrel{~}{\mathrm{\Delta }}_𝐤`$ now have different expressions, everything else remains the same in terms of $`i\stackrel{~}{\omega }`$ and $`\stackrel{~}{\mathrm{\Delta }}_𝐤`$ (as well as their complex conjugate), just as in the BCS case (see Appendix A).

The Feynman diagrams for $`T_\omega `$ and $`T_\mathrm{\Delta }^{}`$ are shown in Fig. 6. To obtain the second line for $`T_\mathrm{\Delta }^{}`$, we make use of the approximation in Fig. 4(a) to convert the left part of the second and the third term on the first line to the full $`G`$. This result is direct analogy with its BCS counterpart as shown in Fig. 20. One can now write down the equations for $`T_\omega `$ and $`T_\mathrm{\Delta }^{}`$ without difficulty.

$`T_\omega (𝐤,𝐤^{})`$ $`=`$ $`u(𝐤,𝐤^{})+{\displaystyle \underset{𝐤^{\prime \prime }}{}}u(𝐤,𝐤^{\prime \prime })G(K^{\prime \prime })T_\omega (𝐤^{\prime \prime },𝐤^{})`$ (32)
$`+{\displaystyle \underset{𝐤^{\prime \prime }}{}}{\displaystyle \underset{Q}{}}u(𝐤,𝐤^{\prime \prime })G(K^{\prime \prime })\widehat{G}_0(QK^{\prime \prime })`$
$`\times \mathrm{\Gamma }(K^{\prime \prime }|Q)t(Q)T_\mathrm{\Delta }^{}(K^{\prime \prime }Q,K^{}|Q),`$
and

$`T_\mathrm{\Delta }^{}(KQ,K^{}|Q)={\displaystyle \underset{𝐤^{\prime \prime }}{}}u(𝐤^{\prime \prime }𝐪,𝐤𝐪)`$
$`\times \widehat{G}_0(QK^{\prime \prime })G(K^{\prime \prime })\mathrm{\Gamma }(K^{\prime \prime }|Q)T_\omega (𝐤^{\prime \prime },𝐤^{})`$
$`+{\displaystyle \underset{𝐤^{\prime \prime }}{}}u(𝐤^{\prime \prime }𝐪,𝐤𝐪)G(QK^{\prime \prime })T_\mathrm{\Delta }^{}(K^{\prime \prime }Q,K^{}|Q).`$
Note here $`T_\mathrm{\Delta }^{}`$ does not contain the factor $`\mathrm{\Delta }`$, unlike its BCS counterpart. It has the same dimension as $`\mathrm{\Gamma }`$. Both $`\mathrm{\Gamma }`$ and $`\mathrm{\Gamma }_0`$ now contain the full impurity $`T`$-matrix beyond the AG level, and the vertex relation Eq. (25) remains valid. The new expression for $`\mathrm{\Gamma }_0`$ is given by

$$\mathrm{\Gamma }_0(K|Q)=\phi _{𝐤𝐪/2}+n_i\underset{𝐤^{}}{}T_\mathrm{\Delta }^{}(KQ,K^{}|Q).$$
(34)
So far, we have kept the derivation for a generic elastic scattering $`u(𝐤,𝐤^{})`$. For isotropic $`s`$-wave scattering, $`u(𝐤,𝐤^{})=u`$. In this case, $`T_\omega `$ and $`T_\mathrm{\Delta }^{}`$ are independent of $`𝐤`$ and $`𝐤^{}`$. Neglecting the momentum dependence, we obtain

$$T_\mathrm{\Delta }^{}(\omega \mathrm{\Omega },\omega )=\frac{u_𝐤G_0(QK)\mathrm{\Gamma }_0(K|Q)G(K)}{1u_𝐤G(QK)}T_\omega ,$$
(35a)
and
$`T_\omega `$ $`=`$ $`u+u{\displaystyle \underset{𝐤}{}}G(K)T_\omega `$ (35b)
$`+u{\displaystyle \underset{Q}{}}t(Q)\left[{\displaystyle \underset{𝐤}{}}G(K)\mathrm{\Gamma }_0(K|Q)G_0(QK)\right]`$
$`\times T_\mathrm{\Delta }^{}(\omega \mathrm{\Omega },\omega ).`$
Use has been made of the vertex relation Eq. (25). Here again, we need to make use of the approximation Eq. (23). Defining $`\overline{\overline{G}}_\omega =_𝐤G(K)`$ and

$$\overline{\overline{F}}_\omega ^{}=\underset{𝐤}{}F^{}(K)=\underset{𝐤}{}\mathrm{\Delta }\mathrm{\Gamma }_0(K)G_0(K)G(K),$$
(36)
we obtain

$$T_\omega =\frac{u(1u\overline{\overline{G}}_\omega )}{\left(1u\overline{\overline{G}}_\omega \right)\left(1u\overline{\overline{G}}_\omega \right)+u^2\overline{\overline{F}}_\omega \overline{\overline{F}}_\omega ^{}},$$
(37a)
and
$`T_\mathrm{\Delta }^{}(\omega \mathrm{\Omega },\omega )`$ $`=`$ $`{\displaystyle \frac{u^2_𝐤G_0(QK)\mathrm{\Gamma }_0(K|Q)G(K)}{\left(1u\overline{\overline{G}}_\omega \right)\left(1u\overline{\overline{G}}_\omega \right)+u^2\overline{\overline{F}}_\omega \overline{\overline{F}}_\omega ^{}}}`$ (37b)
$`\times {\displaystyle \frac{1u\overline{\overline{G}}_\omega }{1u\overline{\overline{G}}_{\mathrm{\Omega }\omega }}}.`$
Letting $`\mathrm{\Omega }0`$, the last equation becomes
$$\mathrm{\Delta }T_\mathrm{\Delta }^{}(\omega )=\frac{u^2\overline{\overline{F}}_\omega ^{}}{\left(1u\overline{\overline{G}}_\omega \right)\left(1u\overline{\overline{G}}_\omega \right)+u^2\overline{\overline{F}}_\omega \overline{\overline{F}}_\omega ^{}}.$$
(37c)
The frequency and gap renormalizations are given by

$$i\stackrel{~}{\omega }=i\omega \mathrm{\Sigma }_\omega ,i\underset{¯}{\overset{~}{\omega }}=i\omega \mathrm{\Sigma }_\omega ,$$
(38a)
$$\stackrel{~}{\mathrm{\Delta }}_𝐤=\mathrm{\Delta }_𝐤+\mathrm{\Sigma }_\mathrm{\Delta },\stackrel{~}{\mathrm{\Delta }}_𝐤^{}=\mathrm{\Delta }_𝐤^{}+\mathrm{\Sigma }_\mathrm{\Delta }^{}$$
(38b)
where $`\mathrm{\Sigma }_\omega =n_iT_\omega `$ and $`\mathrm{\Sigma }_\mathrm{\Delta }=n_i\mathrm{\Delta }T_\mathrm{\Delta }`$. Here $`\underset{¯}{\overset{~}{\omega }}=\stackrel{~}{(\omega )}`$. The expression for $`\chi (Q)`$ remains the same as in previous subsection.

For $`d`$-wave, $`T_\mathrm{\Delta }^{}=T_\mathrm{\Delta }=0`$, and $`\stackrel{~}{\mathrm{\Delta }}_𝐤=\mathrm{\Delta }_𝐤`$. Then Eq. (37a) is greatly simplified,

$$T_\omega =\frac{u}{1u\overline{\overline{G}}_\omega }.$$
(39)
The full Green’s function is given by

$$G(K)=\frac{i\underset{¯}{\overset{~}{\omega }}ϵ_𝐤}{(i\stackrel{~}{\omega }ϵ_𝐤)(i\underset{¯}{\overset{~}{\omega }}ϵ_𝐤)+\mathrm{\Delta }_𝐤^{}\mathrm{\Delta }_𝐤}.$$
(40)
Due to the approximation Eq. (23), we are able to bring the final result Eqs. (39) and (40) into the BCS form. It is easy to show that they are equivalent to the more familiar form in Nambu formalism, as used in Ref. Hirschfeld et al., 1988. Define

$$T_{A\omega }=\frac{1}{2}(T_\omega T_\omega ),T_{S\omega }=\frac{1}{2}(T_\omega +T_\omega ),$$
(41)
and similarly for $`\mathrm{\Sigma }_\omega `$ and $`\overline{\overline{G}}_\omega `$. Here the subscript “A” and “S” denote antisymmetric and symmetric part, respectively. Further define

$$i\stackrel{~}{\omega }_A=i\omega \mathrm{\Sigma }_{A\omega },\stackrel{~}{ϵ}_K=ϵ_𝐤+\mathrm{\Sigma }_{S\omega },$$
(42)
then we obtain (with $`\mathrm{\Delta }^{}=\mathrm{\Delta }`$)

$$G(K)=\frac{i\stackrel{~}{\omega }_A+\stackrel{~}{ϵ}_K}{(i\stackrel{~}{\omega }_A)^2\stackrel{~}{ϵ}_K^2\mathrm{\Delta }_𝐤^2},$$
(43)
and

$`T_{A\omega }`$ $`=`$ $`{\displaystyle \frac{u^2\overline{\overline{G}}_{A\omega }}{(1u\overline{\overline{G}}_{S\omega })^2u^2\overline{\overline{G}}_{A\omega }^2}},`$ (44a)
$`T_{S\omega }`$ $`=`$ $`{\displaystyle \frac{u(1\overline{\overline{G}}_{S\omega })}{(1u\overline{\overline{G}}_{S\omega })^2u^2\overline{\overline{G}}_{A\omega }^2}},`$ (44b)
where

$`\overline{\overline{G}}_{A\omega }`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{i\stackrel{~}{\omega }_A}{(i\stackrel{~}{\omega }_A)^2\stackrel{~}{ϵ}_K^2\mathrm{\Delta }_𝐤^2}},`$ (45a)
$`\overline{\overline{G}}_{S\omega }`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\stackrel{~}{ϵ}_K}{(i\stackrel{~}{\omega }_A)^2\stackrel{~}{ϵ}_K^2\mathrm{\Delta }_𝐤^2}},`$ (45b)
are the antisymmetric and symmetric parts of $`\overline{\overline{G}}_\omega `$, respectively. It is evident that $`T_A`$ and $`T_S`$ correspond to $`T_0`$ and $`T_3`$, respectively, in the Nambu formalism in Ref. Hirschfeld et al., 1988, (and similarly for $`G_A`$ and $`G_S`$).

It should be emphasized, however, that unless $`u=0`$ or $`u=\pm \mathrm{}`$, the symmetric part of the impurity $`T`$-matrix, $`T_{S\omega }`$, can never be set to zero, even if one could in principle have $`\overline{\overline{G}}_{S\omega }=0`$. This means that $`\stackrel{~}{ϵ}_K`$ will always acquire a non-trivial, frequency dependent renormalization, $`\mathrm{\Sigma }_{S\omega }`$. While this renormalization is small for weak coupling BCS superconductors, it is expected to be significant for the cuprate superconductors.

## III Numerical solutions for $`d`$-wave superconductors

### III.1 Analytical continuation and equations to solve

Since there is no explicit pairing vertex renormalization for $`d`$-wave superconductors, i.e., $`\mathrm{\Gamma }_0(K)=\phi _𝐤`$ or $`\stackrel{~}{\mathrm{\Delta }}_𝐤=\mathrm{\Delta }_𝐤`$, a major part of the numerics is to calculate the frequency renormalization. Everything else will follow straightforwardly.

Numerical calculations can be done in the real frequencies, after proper analytical continuation. Since $`T_{S\omega }0`$, $`\mathrm{\Sigma }_\omega `$ and $`\mathrm{\Sigma }_\omega `$ are independent of each other. To obtain the frequency renormalization $`\mathrm{\Sigma }_\omega `$, one has to solve a set of four equations for $`\overline{\overline{G}}_\omega `$, $`\overline{\overline{G}}_\omega `$, $`\mathrm{\Sigma }_\omega `$, and $`\mathrm{\Sigma }_\omega `$ self-consistently for given $`\omega `$. Because $`i\underset{¯}{\overset{~}{\omega }}i\stackrel{~}{\omega }`$, and one needs to analytically continue both simultaneously, the analytical continuation must be done carefully. For $`n>0`$, $`i\stackrel{~}{\omega }_n\omega _+^R=\omega _++i\mathrm{\Sigma }_+^{\prime \prime }`$, and $`i\underset{¯}{\overset{~}{\omega }}_n\omega _{}^A=\omega _{}i\mathrm{\Sigma }_{}^{\prime \prime }`$. For $`n^{}=n<0`$, $`i\stackrel{~}{\omega }_n^{}\omega _{}^R=\omega _{}+i\mathrm{\Sigma }_{}^{\prime \prime }`$ and $`i\underset{¯}{\overset{~}{\omega }}_n^{}\omega _+^A=\omega _+i\mathrm{\Sigma }_+^{\prime \prime }`$. Here $`\omega _\pm =\pm \omega \mathrm{\Sigma }_\pm ^{}`$, and we choose $`\omega >0`$ and $`\mathrm{\Sigma }_\pm ^{\prime \prime }>0`$. Then we obtain four equations as follows:

$`\overline{\overline{G}}_{\omega >0}^R`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\omega _{}i\mathrm{\Sigma }_{}^{\prime \prime }ϵ_𝐤}{(\omega _++i\mathrm{\Sigma }_+^{\prime \prime }ϵ_𝐤)(\omega _{}i\mathrm{\Sigma }_{}^{\prime \prime }ϵ_𝐤)+\mathrm{\Delta }_𝐤^2}}`$
$`\overline{\overline{G}}_{\omega <0}^R`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\omega _+i\mathrm{\Sigma }_+^{\prime \prime }ϵ_𝐤}{(\omega _{}+i\mathrm{\Sigma }_{}^{\prime \prime }ϵ_𝐤)(\omega _+i\mathrm{\Sigma }_+^{\prime \prime }ϵ_𝐤)+\mathrm{\Delta }_𝐤^2}}`$
$`\mathrm{\Sigma }_{\omega >0}^R`$ $`=`$ $`{\displaystyle \frac{n_iu}{1u\overline{\overline{G}}_\omega ^R}}=\mathrm{\Sigma }_+^{}i\mathrm{\Sigma }_+^{\prime \prime }`$
$`\mathrm{\Sigma }_{\omega >0}^R`$ $`=`$ $`{\displaystyle \frac{n_iu}{1u\overline{\overline{G}}_\omega ^R}}=\mathrm{\Sigma }_{}^{}i\mathrm{\Sigma }_{}^{\prime \prime }`$ (46)
These equations are solved self-consistently for $`(\mathrm{\Sigma }_\pm ^{},\mathrm{\Sigma }_\pm ^{\prime \prime })`$, as well as the real and imaginary parts of $`\overline{\overline{G}}_{\pm \omega }`$, as a function of $`\omega `$. No Kramers-Kronig relations are invoked in these numerical calculations. Note in real numerics, we subtract $`\mathrm{\Sigma }_{\omega =0}^{}`$ from $`\mathrm{\Sigma }_\omega ^R`$ so that $`\mathrm{\Sigma }_{\omega =0}^{}=0`$. This subtraction is compensated by a constant shift in the chemical potential $`\mu `$.

Having solved the frequency renormalization $`\mathrm{\Sigma }_\omega `$, one can evaluate the pair susceptibility in the gap equation Eq. (28),

$`\chi (0)={\displaystyle \underset{K}{}}{\displaystyle \frac{\phi _𝐤^2}{(i\stackrel{~}{\omega }ϵ_𝐤)(i\underset{¯}{\overset{~}{\omega }}ϵ_𝐤)+\mathrm{\Delta }_𝐤^2}}`$
$`=\text{Im}{\displaystyle \underset{𝐤}{}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{d\omega }{\pi }}{\displaystyle \frac{[12f(\omega )]\phi _𝐤^2}{(\omega _++i\mathrm{\Sigma }_+^{\prime \prime }ϵ_𝐤)(\omega _{}i\mathrm{\Sigma }_{}^{\prime \prime }ϵ_𝐤)+\mathrm{\Delta }_𝐤^2}},`$
where $`f(x)`$ is the Fermi distribution function. It is easy to check that the gap equation Eq. (28) does reduce to its clean counterpart Eq. (7) as $`n_i0`$.

The particle number equation becomes

$$n=2\underset{𝐤}{}_{\mathrm{}}^{\mathrm{}}\frac{d\omega }{2\pi }A(𝐤,\omega )f(\omega )=2\mathrm{Im}_{\mathrm{}}^{\mathrm{}}\frac{d\omega }{\pi }\overline{\overline{G}}_\omega ^Rf(\omega ).$$
(48)
The real and imaginary parts of $`\chi (Q)`$ are given respectively by

$`\chi ^{}(\mathrm{\Omega }+i0^+,𝐪)`$ (49a)
$`=\text{Im}{\displaystyle \underset{𝐤}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}\{G^R(\omega ,𝐤)G_0^R(\mathrm{\Omega }\omega ,𝐪𝐤)`$
$`\times \left[f(\omega \mathrm{\Omega })f(\omega )\right]`$
$`+G^R(\omega ,𝐤)G_0^A(\mathrm{\Omega }\omega ,𝐪𝐤)`$
$`\times [1f(\omega )f(\omega \mathrm{\Omega })]\}\phi _{𝐤𝐪/2}^2,`$
and
$`\chi ^{\prime \prime }(\mathrm{\Omega }+i0^+,𝐪)`$ (49b)
$`={\displaystyle \underset{𝐤}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}\mathrm{Im}G^R(\omega ,𝐤)A_0(\mathrm{\Omega }\omega ,𝐪𝐤)`$
$`\times \left[f(\mathrm{\Omega }\omega )f(\omega )\right]\phi _{𝐤𝐪/2}^2,`$
where $`A_0(\omega ,𝐤)=2\text{Im}G_0^R(\omega ,𝐤)`$ is the “bare” spectral function.

The pseudogap is evaluated via Eq. (15). To this end, we expand the inverse $`T`$-matrix to the order of $`\mathrm{\Omega }`$ and $`q^2`$ via a (lengthy but straightforward) Taylor expansion

$`t_{pg}^1(\mathrm{\Omega }+i0^+,𝐪)`$ $`=`$ $`\chi (\mathrm{\Omega }+i0^+,𝐪)\chi (0,\mathrm{𝟎})`$ (50)
$`=`$ $`(a_0^{}+ia_0^{\prime \prime })\mathrm{\Omega }+b^{}q^2+ia_1^{\prime \prime }\mathrm{\Omega }^2.`$
Here the imaginary part $`b^{\prime \prime }`$ vanishes. The term $`b^{}q^2`$ should be understood as $`b_{}^{}q_{}^2+b_{}^{}q_{}^2`$ for a quasi-two dimensional square lattice. We keep the imaginary part up to the $`\mathrm{\Omega }^2`$ order. Substituting Eq. (50) into Eq. (15), we obtain

$$\mathrm{\Delta }_{pg}^2=\underset{𝐪}{}_{\mathrm{}}^{\mathrm{}}\frac{d\mathrm{\Omega }}{\pi }\frac{(a_0^{\prime \prime }+a_1^{\prime \prime }\mathrm{\Omega })\mathrm{\Omega }}{(a_0^{}\mathrm{\Omega }+b^{}q^2)^2+(a_0^{\prime \prime }+a_1^{\prime \prime }\mathrm{\Omega })^2\mathrm{\Omega }^2}b(\mathrm{\Omega }),$$
(51)
where $`b(x)`$ is the Bose distribution function.

Numerical solutions confirm that the coefficients $`a_0^{}`$ and $`b^{}`$ have very weak $`T`$ dependence at low $`T`$. Therefore, in a three-dimensional (3D) system, we have $`\mathrm{\Delta }_{pg}^2T^{3/2}`$ at low $`T`$. As the system dimensionality approaches two, the exponent decreases from $`3/2`$ to $`1`$. However, in most physical systems, e.g., the cuprates, this exponent is close to $`3/2`$. The product $`n_pa_0^{}\mathrm{\Delta }_{pg}^2`$ roughly measures the density of incoherent pairs.

The gap equation (28) \[together with Eq. (LABEL:eq:chi-int)\], the fermion number equation (48), and the pseudogap equation (51) form a closed set, which will be solved self-consistently for $`T_c`$, $`\mu `$, and gaps at and below $`T_c`$. For given parameters $`T_c`$, $`\mu `$, and gaps, we can calculate the frequency renormalization, $`\mathrm{\Sigma }_\omega `$, and then solve the three equations. A equation solver is then used to search for the solution for these parameters. The momentum sum is carried out using integrals. Very densely populated data points of $`\omega `$ are used automatically where $`\mathrm{\Sigma }_\omega `$ and/or $`\overline{\overline{G}}_\omega `$ change sharply. A smooth, parabolic interpolation scheme is used in the integration with respect to $`\omega `$. The relative error of the solutions is less than $`1.0\times 10^5`$, and the equations are satisfied with a relative error on both sides less than $`1.0\times 10^7`$. In this way, our numerical results are much more precise than those calculated on a finite size lattice.

The solutions of these equations can be used to calculate the superfluid density $`n_s/m`$. Without giving much details, we state here that the impurity dressing of the current vertex for the (short-coherence-length) cuprate superconductors does not lead to considerable contributions for $`s`$-wave isotropic scattering and the long wavelength $`q0`$ limit. The expression for the in-plane $`n_s/m`$ is given formally by the formula for a clean system (before the Matsubara summation is carried out), as in Ref. Kosztin et al., 2000,

$$\frac{n_s}{m}=\frac{n}{m}+P(0),$$
(52)
where the in-plane current-current correlation function $`P_{ij}(Q)=P(Q)\delta _{ij}`$ can be simply derived from Eqs. (31) and (32) of Ref. Kosztin et al., 2000. The result is given by

$`P(Q)`$ $`=`$ $`{\displaystyle \underset{K}{}}G(K)G(KQ)\{[1+(\mathrm{\Delta }_{sc}^2\mathrm{\Delta }_{pg}^2)\phi _𝐤\phi _{𝐤𝐪}`$ (53)
$`\times G_0(K)G_0(QK)\left]\right({\displaystyle \frac{ϵ_{𝐤𝐪/2}}{𝐤_{}}})^2`$
$`\mathrm{\Delta }_{pg}^2G_0(K){\displaystyle \frac{ϵ_{𝐤𝐪/2}}{𝐤_{}}}{\displaystyle \frac{\phi _𝐤^2}{𝐤_{}}}\}.`$
Using spectral representation and after lengthy but straightforward derivation, we obtain

$`{\displaystyle \frac{n_s}{m}}`$ $`=`$ $`4\mathrm{\Delta }_{sc}^2{\displaystyle \underset{𝐤}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}\mathrm{Im}[\left(\stackrel{~}{F}^A(\omega ,𝐤)\right)^2(\stackrel{}{}ϵ_𝐤)^2\phi _𝐤^2`$ (54)
$`+{\displaystyle \frac{1}{2}}G^A(\omega ,𝐤)\stackrel{~}{F}^A(\omega ,𝐤)\stackrel{}{}ϵ_𝐤\stackrel{}{}\phi _𝐤^2]f(\omega ),`$
where

$$\stackrel{~}{F}(K)=G(K)G_0(K)=\frac{1}{(i\stackrel{~}{\omega }ϵ_𝐤)(i\underset{¯}{\overset{~}{\omega }}ϵ_𝐤)+\mathrm{\Delta }_𝐤^2},$$
(55)
which differs from $`F(K)`$ by a factor $`\mathrm{\Delta }_𝐤`$.

As in the clean system, Eq. (54) differs from its BCS counterpart $`(n_s/m)_{BCS}`$ only by the overall prefactor, $`\mathrm{\Delta }_{sc}^2`$,

$$\frac{n_s}{m}=\frac{\mathrm{\Delta }_{sc}^2}{\mathrm{\Delta }^2}\left(\frac{n_s}{m}\right)_{BCS}.$$
(56)
For $`d`$-wave superconductors, $`(n_s/m)_{BCS}ABT`$ or $`ABT^2`$ at very low $`T`$, depending on whether the system is clean or dirty. Bearing in mind that $`\mathrm{\Delta }_{sc}^2/\mathrm{\Delta }^2=1\mathrm{\Delta }_{pg}^2/\mathrm{\Delta }^21C^{}T^{3/2}`$, we have

$`{\displaystyle \frac{n_s}{m}}`$ $``$ $`ABTCT^{3/2},\text{(clean)},`$ (57)
$``$ $`ABT^2CT^{3/2},\text{(dirty)},`$
at very low $`T`$. In other words, pair excitations lead to a new $`T^{3/2}`$ term in the low $`T`$ superfluid density.

### III.2 Renormalization of the frequency by impurity scattering and the density of states

Except in the Born limit, for $`u4t_{}`$, impurity scattering usually introduces a sharp resonance close to the Fermi level in the frequency renormalization for a $`d`$-wave superconductor. In addition, it induces an impurity band outside the main particle band. Both the low energy resonance and the high energy impurity band arise from the vanishing of the real part of the denominator of the impurity $`T`$-matrix, Eq. (39), while the imaginary part is small. In Fig. 7 we show an example of (a) the retarded, impurity induced renormalization of the frequency, $`\mathrm{\Sigma }_\omega ^R=\stackrel{~}{\omega }^R\omega `$ and (b) the corresponding impurity average of the single-particle Green’s function, $`\overline{\overline{G}}_\omega ^R`$, which is related to the density of states by $`N(\omega )=2\text{Im}\overline{\overline{G}}_\omega ^R`$. The curves for $`\mathrm{\Sigma }_\omega ^R`$ corresponding to the impurity band are replotted in the inset, to show the strong renormalization of the frequency inside the impurity band. (The curves for $`\mathrm{\Sigma }_\omega ^{}`$ has been offset so that $`\mathrm{\Sigma }_{\omega =0}^{}=0`$. This offset is compensated by a shift in the chemical potential $`\mu `$.) The van Hove singularity and its mirror image via particle-hole mixing are clearly seen in the density of states, and are also reflected in $`\mathrm{\Sigma }_\omega ^R`$ as the small kinks in Fig. 7(a). The real part $`\mathrm{\Sigma }_\omega ^{}`$ is usually neglected in most non-self-consistent calculations.Hirschfeld et al. (1988) It is evident, however, that it has a very rich structure, and, in general, cannot be set to zero in any self-consistent calculations. This conclusion holds even in the presence of exact particle-hole symmetry, as can be easily told from Eq. (39).

For $`u<0`$, the low energy resonance in $`\text{Im}\mathrm{\Sigma }_\omega ^R`$ will appear on the positive energy side in Fig. 7(a). Regardless of the sign of $`u`$, the resonance peak will become sharper as $`n_i`$ decreases and as $`|u|`$ increases. For larger $`|u|`$, the resonant frequency will be closer to $`\omega =0`$, where $`\text{Im}\overline{\overline{G}}_\omega ^R`$ is small because of the $`d`$-wave symmetry; A smaller $`n_i`$ further reduces $`\text{Im}\overline{\overline{G}}_\omega ^R`$. Both factors help minimize the imaginary part of the denominator of Eq. (39) and, thus, lead to a stronger resonance. It should be emphasized that a resonance peak does not show up in $`\text{Im}\overline{\overline{G}}_\omega ^R`$ since the resonance in $`\text{Im}\mathrm{\Sigma }_\omega ^R`$ requires that $`\text{Im}\overline{\overline{G}}_\omega ^R`$ be small at the resonance point.

The location of the impurity band is sensitive to the sign and strength of impurity scattering. For $`u<0`$, the impurity band on the negative energy (left) side of the plot. As $`|u|`$ gets smaller, the impurity band merges with the main band; as $`|u|`$ gets larger, it moves farther away, with a much stronger renormalization of $`\omega `$. For large $`|u|`$, the spectral weight under the impurity band in Fig. 7(b) is given by $`2n_i`$, and the weight in the main band is reduced to $`2(1n_i)`$. This leads to a dramatic chemical potential shift as a function of the impurity concentration $`n_i`$ (as well as $`u`$). For $`u<0`$, the impurity band will always be filled, so that increasing $`n_i`$ pushes the system farther away from the particle-hole symmetry. For $`u>0`$, on the contrary, the impurity band is empty, and the system becomes more particle-hole symmetrical as $`n_i`$ increases from 0, and reaches the particle-hole symmetry (in the main band) at $`n_i=1n`$. This fact implies that Pépin and Lee’s assumption of an exact particle-hole symmetry is not justified so that their prediction of a diverging $`N(\omega )`$ as $`\omega 0`$ is unlikely to be observed experimentally. It should be mentioned that the appearance of the impurity band has not been shown in the literature, largely because most authors concentrate on the low energy part of the spectra only, and do not solve for the full spectrum of the renormalization of $`\omega `$ self-consistently.

In the Born limit, only the product $`\gamma =n_iu^2`$ is a meaningful parameter, not $`n_i`$ or $`u`$ separately. In Fig. 8, we plot (a) the frequency renormalization $`\text{Im}\mathrm{\Sigma }_\omega ^R`$ and (b) the impurity average of Green’s function $`\text{Im}\overline{\overline{G}}_\omega ^R`$ as a function of $`\omega `$ for various values of $`u`$ but with a fixed $`\gamma =0.02`$. (Note: when $`u`$ is small, this requires an unphysically large $`n_i`$.) In the Born limit, these two quantities are identical up to a constant coefficient. As $`u`$ increases, a resonance develops in $`\text{Im}\mathrm{\Sigma }_\omega ^R`$ at small $`\omega `$. $`\text{Im}\mathrm{\Sigma }_\omega ^R`$ and $`\text{Im}\overline{\overline{G}}_\omega ^R`$ become very different. And a impurity band develops gradually ($`u=0.2`$ and $`u=0.5`$), until it splits from the main band ($`u=1`$). At fixed $`\gamma `$, the Born limit is more effective in filling in the DOS in the mid-range of $`\omega `$ within the gap and smearing out the coherence quasiparticle peaks, whereas the large $`u`$ limit is more effective in filling in the DOS in the vicinity of $`\omega =0`$ but keeping the quasiparticle peaks largely unchanged. In addition, the main band becomes narrower at large $`u`$ than that in the clean system or the Born limit, so that part of the spectral weight has now been transferred to the impurity band. Also note that the DOS at $`\omega =0`$ is essentially zero in Fig. 8(b) because $`n_i`$ is very small when $`u`$ becomes large for the current choice of $`\gamma =0.02`$.

The effects of the scattering strength $`u`$ and the impurity density $`n_i`$ on the DOS are shown in Fig. 9(a) and (b), respectively. For the effect of $`u`$ in Fig. 9(a), we choose an intermediate $`n_i=0.02`$. And for the effect of $`n_i`$, we focus on the unitary limit, and choose $`u=1000`$. There is a dip in the DOS at $`\omega =0`$ in both the small $`u`$ and small $`n_i`$ cases, mimicking a fractional power law dependence on $`\omega `$. At higher $`n_i`$ and higher $`u`$, the DOS is filled in mainly at small $`\omega `$.

Shown in Fig. 10 are the residue DOS at the Fermi level, $`N(0)`$, as a function of (a) the impurity concentration $`n_i`$ for different values of $`u`$ from the unitary limit $`u=1000`$ through $`u=10`$, 5, 3, 2, down to 1, and (b) as a function of the scattering strength $`u`$ for $`n_i=0.005`$, 0.01, 0.02, 0.03, and 0.05. Figure 10(a) indicates that below certain “critical” value of $`n_i`$, $`N(0)`$ remains essentially zero. This behavior is also implied by the presence of the dip at small $`n_i`$ in Fig. 9(b). The “critical” value for $`n_i`$ in Fig. 10(a) is clearly scattering strength dependent. The smaller $`u`$, the larger this value. A replot (not shown) of these curves in terms of $`\mathrm{log}_{10}N(0)`$ as a function of $`1/n_i`$ reveals that for small $`n_i`$, $`N(0)`$ vanishes exponentially as $`e^{A/n_i}`$, where $`A`$ is a constant. For comparison, we also show in Fig. 10(a) the Born limit as a function of $`\gamma (=n_i)`$. As one may expect, the Born limit is rather different from the rest, since it is equivalent to a very small $`u<<1`$ and unphysically large $`n_i`$. A similar activation behavior of $`N(0)`$ as a function of $`u`$ is seen in Fig. 10(b), where the “critical” value for $`u`$ is strongly $`n_i`$ dependent. The asymmetry between positive and negative $`u`$ reflects the particle-hole asymmetry at $`\mu =0.9`$. It should be noted that it is not realistic to vary $`u`$ continuously in experiment.

An earlier experiment by Ishida et al.Ishida et al. (1993) suggests that $`N(0)`$ varies as $`n_i^{1/2}`$. In our calculations, however, $`N(0)`$ does not follow a simple power law as a function of $`n_i`$. The curve for $`u=1000`$ in Fig. 10(a) fits perfectly with $`an_i^\alpha b`$, with $`\alpha 0.175`$, for $`0.002<n_i<0.05`$. The damping of the zero frequency (not shown), $`\text{Im}\mathrm{\Sigma }_{\omega =0}^R`$, also fits this functional form very well, with $`\alpha 0.61`$. The exponents are different for different values of $`u`$. Our calculation for the $`n_i`$ dependence of $`N(0)`$ is consistent with the result of FehrenbacherFehrenbacher (1996) in that $`N(0,n_i)`$ it is strongly $`u`$ dependent. However, it does not seem likely that the simple power law $`N(0)n_i^{1/2}`$ may be obtained in an accurate experimental measurement. Further experiments are needed to double check this relationship between $`N(0)`$ and $`n_i`$.

From Figs. 7-10, we conclude that for very small $`n_i`$ and $`u`$, the zero frequency DOS $`N(0)`$ is exponentially small. At high $`n_i`$, $`N(0)`$ is finite in both the Born and the unitary limit. However, for certain intermediate values of $`n_i`$ and $`u`$ \[e.g., $`u1`$ in Fig. 9(a), and $`n_i0.002`$ in Fig. 9(b)\], $`N(\omega )`$ vanishes with (very small) $`\omega `$ according to some fractional power $`\omega ^\alpha `$ where $`\alpha <1`$. We see neither the universal power laws for $`N(\omega )`$ predicted by Senthil and Fisher,Senthil and Fisher (1999) nor the divergent DOS predicted by Pépin and LeePépin and Lee (2001) and others.Bocquet et al. (2000)

When the values of $`u`$ and $`n_i`$ are such that $`N(\omega )\omega ^\alpha `$ with $`\alpha <1`$, one may expect to see a fractional low $`T`$ power law in the superfluid density. However, such a power law is not robust as it is sensitive to the impurity density for a given type of impurity. The situation with a negative $`u`$ is similar to Fig. 9.

For given chemical potential $`\mu `$ and the total excitation gap $`\mathrm{\Delta }`$, the calculation of the frequency renormalization with impurities does not necessarily involve the concept of the pseudogap. It is essentially the “self-consistent” impurity $`T`$-matrix treatment by Hirschfeld et al.Hirschfeld et al. (1988) except that we now have to solve for both the real and imaginary parts of $`\mathrm{\Sigma }_\omega ^R`$ simultaneously in a self-consistent fashion. Our numerical results agree with existing calculations in the literature.

Finally, we emphasize the difference between the self-consistent impurity $`T`$-matrix treatment of the one-impurity problem Kruis et al. (2001); Atkinson et al. (2000b) and the current many-impurity averaging. For the former case, the impurity $`T`$-matrix will be given by Eq. (37a) but with $`\overline{\overline{G}}_\omega `$ replaced by $`\overline{\overline{G}}_\omega ^0`$, i.e., the impurity average of the clean $`G^0(K)`$. As a consequence, the position of the poles of $`T_\omega `$ is independent of the renormalized DOS, and, therefore, a resonance peak may exist at low $`\omega `$ in the DOS,Kruis et al. (2001) whereas it cannot in the current many-impurity case.

### III.3 Effects of the impurity scattering on $`T_c`$ and the pseudogap

In this section, we study the influence of impurity scattering on the behavior of $`T_c`$ and the pseudogap $`\mathrm{\Delta }_{pg}`$ as a function of the coupling strength as well as the hole doping concentration.

First, we study the effect of the scattering strength $`u`$ and whether it is repulsive ($`u>0`$) or attractive ($`u<0`$). In Fig. 11, we plot $`T_c`$ as a function of $`u`$, for a pseudogapped $`d`$-wave superconductor with $`n_i=0.05`$. The corresponding curves for $`\mu `$ and $`\mathrm{\Delta }=\mathrm{\Delta }_{pg}(T_c)`$ are shown in the upper left and upper right insets, respectively. All three quantities, $`T_c`$, $`\mu `$, and $`\mathrm{\Delta }`$, vary with $`u`$. For either sign of $`u`$, both $`T_c`$ and $`\mathrm{\Delta }`$ are suppressed by increasing $`|u|`$. It should be emphasized that the chemical potential $`\mu `$ in the two (large $`\pm u`$) unitary limits does not meet, nor does $`T_c`$ or $`\mathrm{\Delta }`$. This is because a (large) negative $`u`$ creates a filled impurity band below the main band, and is effective in bringing down the chemical potential, whereas a positive $`u`$ creates an empty impurity band above the main band, and tends to raise $`\mu `$ toward the particle-hole symmetrical point, $`\mu =1`$. This result cannot and has not been observed in previous, non-self-consistent calculations where the real part of frequency renormalization is set to zero.

In Fig. 12, we compare the effect of the impurity concentration for different scattering strengths: the Born limit, both unitary limits ($`u=\pm 100`$), and intermediate $`u=1`$. Both (b) $`T_c`$ and (a) $`\mathrm{\Delta }_{pg}(T_c)`$ are suppressed by increasing impurity density. This is natural in a model where the pseudogap originates from incoherent pair excitations. As will be seen below, $`\mathrm{\Delta }_{pg}(T_c)`$ is suppressed mainly because $`T_c`$ is lowered. Except in the Born limit, the chemical potential $`\mu `$ is fairly sensitive to $`n_i`$, as shown in the inset. It is clear that the scattering in the Born limit is the most effective in suppressing $`T_c`$. In comparison with experiment,Franz et al. (1997) calculations at the AG level (i.e., the Born limit) tend to overestimate the $`T_c`$ suppression by as much as a factor of 2. This is in good agreement with the current result in the unitary limit. At large $`n_i`$ for large positive $`u`$, the system is driven to the particle-hole symmetrical point, where the effective pair mass changes sign. It is usually hard to suppress $`T_c`$ by pairing at the particle-hole symmetrical point, as indicated by the solid curve in the lower panel. In fact, exactly at this point, the linear $`\mathrm{\Omega }`$ term $`a_0^{}`$ in the inverse $`T`$ matrix expansion vanishes, so that one needs to go beyond the current approximation and expand up to the $`\mathrm{\Omega }^2`$ term. Large negative $`u=100`$ is more effective in suppressing $`T_c`$ and $`\mathrm{\Delta }`$ than intermediate negative $`u=1`$, in agreement with Fig. 11 and the DOS shown in Fig. 9(a) and Fig. 10.

There is enough evidence that zinc impurities are attractive scatterer for electrons in the cuprates.Fehrenbacher (1996) Therefore, we concentrate ourselves on negative $`u`$ scattering in the unitary limit. Plotted in Fig. 13 are $`T_c`$ (main figure), $`\mu `$ (lower inset), and $`\mathrm{\Delta }`$ as a function of $`g`$ for increasing $`n_i=0,0.02,0.05,0.1`$ with $`u=100`$. Also plotted for comparison are the results assuming the Born limit with $`\gamma =0.05`$. Clearly, the Born limit is more effective in suppressing $`T_c`$ at relatively weak coupling, $`g/4t_{}0.75`$, consistent with Fig. 12. Both $`T_c`$ and $`\mathrm{\Delta }`$ are suppressed continuously with $`n_i`$. However, it should be noted that a larger $`n_i`$ helps $`T_c`$ to survive a larger $`g/4t_{}`$. This is mainly because the filled impurity band at large $`n_i`$ pushes the system far away from particle-hole symmetry (see $`\mu `$ in the lower inset), and reduces the effective fermion density, so that the pair mobility is enhanced and the pair mass does not diverge until a larger $`g/4t_{}`$ is reached.

To make contact with the cuprates, we use the non-double occupancy condition associated with the Mott insulator physics, as in Ref. Chen et al., 1998, so that the effective hopping integral is reduced to $`t_{}(x)t_0x`$, where $`x=1n`$ is the hole concentration, and $`t_0=0.6`$ eV is the hopping integral in the absence of the on-site Coulomb repulsion. We assume $`g/4t_0=0.047`$, which is $`x`$ independent. Then we can compute $`T_c`$, $`\mu `$, and $`\mathrm{\Delta }=\mathrm{\Delta }_{pg}(T_c)`$ as a function of $`x`$. The result for $`T_c`$ (main figure) and $`\mathrm{\Delta }`$ (inset) is shown in Fig. 14 for the clean system and $`n_i=0.02`$ and $`0.05`$. In the overdoped regime, $`T_c`$, as well as the small $`\mathrm{\Delta }`$, are strongly suppressed by impurities. This provides a natural explanation for the experimental observation that $`T_c`$ vanishes abruptly at large $`x`$; it is well-known that high crystallinity, clean samples are not available in the extreme overdoped regime. On the other hand, neither $`T_c`$ nor $`\mathrm{\Delta }`$ is strongly suppressed in the highly underdoped regime, where the gap is too large. At this point, experimental data in this extreme underdoped regime are still not available. Our result about the suppression of $`T_c`$ and $`\mathrm{\Delta }`$ in the less strongly underdoped regime ($`x>0.1`$) are in agreement with experimental observationsWilliams et al. (1995, 1998) and other calculations.Franz et al. (1997)

It should be noted, however, that in our simple model, we do not consider the fact that disorder or impurities may reduce the dimensionality of the electron motion and thus suppress $`T_c`$. Furthermore, since induced local spin and Kondo effects have been observed near zinc sites in both zinc-doped YBCOMacFarlane et al. (2000); Bobroff et al. (1999) and zinc-doped Bi2Sr2-xLaxCuO6+δ,Hanaki et al. (2001) this raises an important question whether zinc can be treated as a nonmagnetic impurity.

### III.4 Gaps and superfluid density below $`T_c`$ in the presence of nonmagnetic impurities

In this subsection, we study the effect of nonmagnetic impurities on the behavior of the excitation gap $`\mathrm{\Delta }`$, the order parameter $`\mathrm{\Delta }_{sc}`$, and the pseudogap $`\mathrm{\Delta }_{pg}`$ as well as the chemical potential $`\mu `$ as a function of temperature in the superconducting phase. The numerical solutions for these quantities are then used to study the temperature dependence of the superfluid density $`n_s/m`$ at various impurity levels. We concentrate on the unitary limit, which is regarded as relevant to the cuprates. To be specific, we use $`n=0.85`$, $`t_{}/t_{}=0.01`$, and $`u=100`$ in the calculations presented below.

We first study the impurity effect in the BCS case, without the complication of the pseudogap. Plotted in Fig. 15 are the superfluid density $`n_s/m`$ (main figure) and the corresponding gap $`\mathrm{\Delta }`$ in a $`d`$-wave BCS superconductor as a function of the reduced temperature $`T/T_c^0`$ for the clean system (solid curves), impurity density $`n_i=0.02`$ (dotted), and $`n_i=0.05`$ (dashed) at $`g/4t_{}=0.3`$. Here $`T_c^0=0.0416`$ is the $`T_c`$ in the clean case. As expected, both $`T_c`$ and $`\mathrm{\Delta }(T)`$, as well as $`n_s/m`$, are suppressed by impurity scattering. In agreement with experiment, the low $`T`$ normal fluid density is linear in $`T`$ in the clean case, and becomes quadratic in the two dirty cases.

Now we add pseudogap for the underdoped cuprates. We show in Fig. 16 the temperature dependence of $`n_s/m`$ (main figure) and various gaps (inset) in a $`d`$-wave pseudogapped superconductor for impurity concentration $`n_i=0`$ (clean, solid curve), 0.02 (dotted), and 0.05 (dashed) in the unitary limit at $`g/4t_{}=0.5`$. As the order parameter develops below $`T_c`$, the pseudogap decreases with decreasing $`T`$. This reflects the fact that the pseudogap in the present model is a measure of the density of thermally excited pair excitations. The total gap $`\mathrm{\Delta }`$, the order parameter $`\mathrm{\Delta }_{sc}`$, and the superfluid density $`n_s/m`$ are suppressed by increasing $`n_i`$, similar to the BCS case above. However, at given $`T<T_c`$, the pseudogap $`\mathrm{\Delta }_{pg}`$ remains roughly unchanged. Furthermore, the low $`T`$ power law for the superfluid density is different from the BCS case, as predicted in Eq. (57). It is now given by $`T+T^{3/2}`$ and $`T^2+T^{3/2}`$ for the clean and dirty cases, respectively. Due to the presence of the $`T^{3/2}`$ term, the low $`T`$ portion of the curves for $`n_i=0.02`$ and 0.05 are clearly not as flat as in Fig. 15. Nevertheless, it may be difficult to distinguish experimentally $`T^2+T^{3/2}`$ from a pure $`T^2`$ power law. This $`T^{3/2}`$ contribution of the pair excitations has been used to explain successfullyChen et al. (1998) the quasi-universal behavior of the normalized superfluid density $`n_s(T)/n_s(0)`$ as a function of $`T/T_c`$. We have also foundChen et al. (2000b); Carrington et al. (1999) preliminary experimental support for this $`T^{3/2}`$ term in low $`T`$ penetration depth measurement in the cuprates as well as organic superconductors. Systematic experiments are needed to further verify this power law prediction.

A careful look at the values of the zero temperature superfluid density $`(n_s/m)_0`$ for different values of $`n_i`$ in both Fig. 15 and Fig. 16 suggests that in the unitary limit, $`(n_s/m)_0`$ drops faster with $`n_i`$ when $`n_i`$ is still small. This is manifested in a systematic study of $`(n_s/m)_0`$ as a function of $`n_i`$, as shown in Fig. 17, with $`u=100`$ (solid curve). This behavior has been observed experimentally.Bernhard et al. (1996) Also plotted in the inset is the corresponding zero temperature gap $`\mathrm{\Delta }_0`$ versus $`n_i`$. Clearly, the slope $`d\frac{n_s}{m}/dn_i`$ is much steeper as $`n_i`$ approaches zero, very different from the behavior of $`\mathrm{\Delta }_0`$. This demonstrates that $`n_s/m`$ is influenced more by the DOS than by the gap size. A very small amount of impurities may strongly suppress $`(n_s/m)_0`$. This conclusion is significant in data analysis of the penetration depth measurement, especially when the $`T=0`$ value of penetration depth is not measured directly.Har  For comparison, we also plot the corresponding curves in the Born limit. While the gap is suppressed faster, in contrast to unitary case, the slope $`d\frac{n_s}{m}/dn_i`$ is smaller for smaller $`n_i`$.

## IV Discussion

In Sec. II, we have used the approximation Eq. (23) to bring the single-particle self-energy and thus the gap equation into a BCS-like form. This approximation derives from the divergence of the $`T`$ matrix $`t_{pg}(Q)`$ as $`Q0`$, which is the pairing instability condition. The spirit of this approximation is to “put” the incoherent, excited pairs into the condensate, by setting $`Q=0`$. The contribution of these pseudo Cooper pairs to the single-particle excitation gap is calculated via Eqs. (15) and (51), weighted by the Bose function. Therefore, the incoherent pairs and the zero momentum condensate are not distinguished from each other in terms of the single particle self-energy, as they add up to a total excitation gap. However, they are distinct when phase sensitive quantities are involved, e.g., in the calculation of $`T_c`$ and of the superfluid density.

With this approximation, there is a close analogy between the Feynman diagrams in the current pairing fluctuation theory in Sec. II and its BCS counterpart in Appendix A. When the finite momentum pair propagators are removed (or “pushed into the condensate”) from Fig. 3 through Fig. 6, these diagrams will become their BCS counterpart in Fig. 18 through Fig. 20. (The diagram for BCS pairing vertex is not shown in Appendix A).

This approximation is in general good when the gap is large in comparison with $`T`$, and when the impurity concentration is low. When the gap is small, the contribution of the incoherent pair excitations is usually small, and does not have a strong effect on $`T_c`$, so that $`T_c`$ is roughly determined by its BCS mean-field solution. When $`n_i`$ is large, the fermionic frequency renormalization is strong, and the pair dispersion also becomes highly damped. In this case, approximation Eq. (23) may not be quantitatively accurate.

Even without the complication of impurities, incoherent pairs are not expected to deplete completely the spectral weight within the two quasiparticle peaks of the spectral function.Chen et al. (2001); Maly et al. (1999a) This, however, cannot be captured by the approximation Eq. (23). Unfortunately, we still do not know yet how to solve the Dyson’s equations without this approximation due to technical difficulties.

Another simplification comes from the $`d`$-wave symmetry of the cuprate superconductors under study. Although we have kept the theoretical formalism general for both $`s`$\- and $`d`$-wave in Sec. II, the pairing vertex renormalization drops out when we finally carry out numerical calculations for $`d`$-wave superconductors. For $`s`$-wave superconductors, one would have to include self-consistently one more complex equation for the renormalization of $`\mathrm{\Delta }_𝐤`$, when solving for the renormalization of $`\omega `$. And the equations (37a) and (37c) also look much more complicated than Eq. (39). Nevertheless, since there is no node in the excitation gap for $`s`$-wave, the numerics is expected to run faster.

It is well-known that for $`d`$-wave superconductors, the Anderson’s theoremAnderson (1959) breaks down.Gor’kov  For Anderson’s theorem to hold, it requires that the gap and the frequency are renormalized in exactly the same fashion. This condition can be satisfied (approximately) only in weak coupling, isotropic BCS $`s`$-wave superconductors, for which the real part of the frequency renormalization is negligible. Since the frequency $`\omega `$ is a scalar, this condition is violated when the gap $`\mathrm{\Delta }_𝐤`$ have any anisotropic dependence on k. Furthermore, when the gap is considerably large in comparison with the band width so that the upper limit of the energy integral cannot be extended to infinity, this condition will not be satisfied, either. In both cases, $`T_c`$ will be suppressed.

## V Conclusions

In this paper, we extend the pairing fluction theory to superconductors in the presence of non-magnetic impurities. Both the pairing and impurity $`T`$-matrices are included and treated self-consistently. We obtain a set of three equations for ($`T_c`$, $`\mu `$, $`\mathrm{\Delta }(T_c)`$) or ($`\mu `$, $`\mathrm{\Delta }_{sc}`$, $`\mathrm{\Delta }_{pg}`$) at $`T<T_c`$, with the complex equations for the frequency renormalization. In consequence, we are able to study the impurity effects on $`T_c`$, the order parameter, and the pseudogap. In particular, we carry out calculations for $`d`$-wave superconductors and apply to the cuprate superconductors. Instead of studying the physical quantities with all possible combinations of the parameters $`n_i`$, $`u`$, $`g`$, and $`n`$, we mainly concentrate on the negative $`u`$ unitary limit, which is regarded as relevant to the zinc impurities in the cuprates.Fehrenbacher (1996)

Calculations show that in addition to the low energy resonance in the imaginary part of the renormalized frequency, a considerably large $`|u|`$ leads to a separate impurity band, with a spectral weight $`2n_i`$. The real part of the frequency renormalization, in general, cannot be set to zero in a self-consistent calculation. The chemical potential varies with the impurity concentration, so that the assumption of exact particle-hole symmetry is not justified when one studies the impurity effects. One consequence of this chemical potential shift is that the repulsive and attractive unitary scattering limits do not meet as has been widely assumed in the non-self-consistent treatment in the literature. Unitary scatterers fill in the DOS mostly in the small $`\omega `$ region, whereas Born scatterers do in essentially the whole range within the gap. At small $`n_i`$ and/or small $`u`$, there is a dip at $`\omega =0`$ in the DOS, so that $`N(\omega )`$ vanishes as a fractional power of $`\omega `$, which may in turn contribute a fractional power law for the low $`T`$ temperature dependence of the penetration depth.

Both $`T_c`$ and the pseudogap $`\mathrm{\Delta }_{pg}(T_c)`$ are suppressed by impurities. In this respect, Born scatterers are about twice as effective as unitary scatterer. Treating zinc impurities as unitary scatterers explains why the actual $`T_c`$ suppression is only half that predicted by calculations at the AG level (i.e., in the Born limit). In the overdoped regime, the gap is small, and therefore the superconductivity can be easily destroyed by a small amount of impurities. In contrast, it takes a larger amount of impurities to destroy the large excitation gap in the underdoped regime.

The reason $`\mathrm{\Delta }_{pg}(T_c)`$ is suppressed is mainly because $`T_c`$ is suppressed. In fact, for a given $`T<T_c`$, the pseudogap remains roughly unchanged (actually it increases slightly). The suppression of the total excitation gap arises from the suppression of the order parameter. The density of incoherent pairs, as measured by $`n_p`$, slightly increases for not-so-large $`n_i`$. This supports the notion that nonmagnetic impurities do not mainly break incoherent pairs. Instead, they scatter the Cooper pairs out of the condensate.Not

Our self-consistent calculations show that in the unitary limit, the low $`T`$ superfluid density is quadratic in $`T`$ in a BCS $`d`$-wave superconductor, in agreement with existing calculations and experiment. Strong pair excitations add an additional $`T^{3/2}`$ term, with preliminary experimental support. As a function of increasing $`n_i`$, the zero $`T`$ superfluid density decreases faster at first for unitary scatterers, whereas the opposite holds for scattering in the Born limit. The former behavior is in agreement with experiment.

###### Acknowledgements.

We would like to thank A. V. Balatsky, P. J. Hirschfeld, and K. Levin for useful discussions. The numerics was in part carried out on the computing facilities of the Department of Engineering, the Florida State University. This work is supported by the State of Florida.

## Appendix A Impurity dressing for BCS theory at the Abrikosov-Gor’kov level

In this appendix, we present the impurity dressing for a BCS superconductor, following Abrikosov-Gor’kovAbrikosov and Gor’kov (1959); AGD , but in a more general form, namely, we do not assume $`\overline{G}_\omega =\overline{G}_\omega `$. This will make it easier to understand the current theory in the presence of strong pairing correlations, as there is a strong similarity between the impurity dressing diagrams for both BCS theory and the pairing fluctuation theory.

For a pure BCS superconductor, we have the Gor’kov equations,

$$G_{0}^{0}{}_{}{}^{1}(K)G^0(K)=1\mathrm{\Delta }_𝐤F_{}^{0}{}_{}{}^{}(K),$$
(58a)
$$G_{0}^{0}{}_{}{}^{1}(K)F_{}^{0}{}_{}{}^{}(K)=\mathrm{\Delta }_𝐤^{}G^0(K).$$
(58b)
At the AG level, the relationship between the impurity dressed Green’s functions $`G`$ and $`F`$ is represented by the Feynman diagrams shown in Fig. 18. Define the impurity average $`\overline{G}_\omega `$ as in Eq. (19), and

$$\overline{F}_\omega ^{}=n_i\underset{𝐤^{}}{}|u(\text{k}\text{k}^{})|^2F^{}(K^{}),$$
(59)
as well as their complex conjugate. Note Fig. 18 is actually Fig. 105 in Ref. AGD, . Without giving details, we give the result following AG:

$$(i\omega ϵ_𝐤\overline{G}_\omega )G(K)+(\mathrm{\Delta }_𝐤+\overline{F}_\omega )F^{}(K)=1,$$
(60a)
$$(i\omega +ϵ_𝐤+\overline{G}_\omega )F^{}(K)+(\mathrm{\Delta }_𝐤^{}+\overline{F}_\omega ^{})G(K)=0.$$
(60b)
Define $`i\stackrel{~}{\omega }=i\omega \overline{G}_\omega `$, $`i\underset{¯}{\overset{~}{\omega }}=i\omega \overline{G}_\omega `$, $`\stackrel{~}{\mathrm{\Delta }}_𝐤=\mathrm{\Delta }_𝐤+\overline{F}_\omega `$, and $`\stackrel{~}{\mathrm{\Delta }}_𝐤^{}=\mathrm{\Delta }_𝐤^{}+\overline{F}_\omega ^{}`$. Then we obtain

$$G(K)=\frac{i\underset{¯}{\overset{~}{\omega }}ϵ_𝐤}{(i\stackrel{~}{\omega }ϵ_𝐤)(i\underset{¯}{\overset{~}{\omega }}ϵ_𝐤)+\stackrel{~}{\mathrm{\Delta }}_𝐤^{}\stackrel{~}{\mathrm{\Delta }}_𝐤},$$
(61a)
$$F^{}(K)=\frac{\stackrel{~}{\mathrm{\Delta }}_𝐤^{}}{(i\stackrel{~}{\omega }ϵ_𝐤)(i\underset{¯}{\overset{~}{\omega }}ϵ_𝐤)+\stackrel{~}{\mathrm{\Delta }}_𝐤^{}\stackrel{~}{\mathrm{\Delta }}_𝐤}.$$
(61b)
For $`d`$-wave, the first equation becomes Eq. (40). Note $`G(K)`$ is no longer symmetrical in $`\omega `$ in general as a consequence of impurity scattering, but $`F(K)`$ still is, since $`F(K)`$ involves $`\pm \omega `$ pairs.

The above result can be easily extended to self-consistent impurity $`T`$-matrix calculations, by replacing the AG-level impurity scattering with the self-consistent impurity $`T`$-matrices, as shown in Fig. 19. The relationship between the regular and anomalous impurity $`T`$ matrices $`T_\omega `$ and $`T_\mathrm{\Delta }^{}`$ are shown in Fig. 20. One can easily write down the corresponding equations, as follows.

$$T_\omega =u+u\overline{\overline{G}}_\omega T_\omega u\overline{\overline{F}}_\omega T_\mathrm{\Delta }^{},$$
(62a)
$$T_\mathrm{\Delta }^{}=u\overline{\overline{F}}_\omega ^{}T_\omega +u\overline{\overline{G}}_\omega T_\mathrm{\Delta }^{}.$$
(62b)
Finally, one has

$$T_\omega =\frac{u(1u\overline{\overline{G}}_\omega )}{\left(1u\overline{\overline{G}}_\omega \right)\left(1u\overline{\overline{G}}_\omega \right)+u^2\overline{\overline{F}}_\omega \overline{\overline{F}}_\omega ^{}},$$
(63a)
$$T_\mathrm{\Delta }^{}(\omega )=\frac{u^2\overline{\overline{F}}_\omega ^{}}{\left(1u\overline{\overline{G}}_\omega \right)\left(1u\overline{\overline{G}}_\omega \right)+u^2\overline{\overline{F}}_\omega \overline{\overline{F}}_\omega ^{}},$$
(63b)
where $`\overline{\overline{G}}=_𝐤G(K)`$, $`\overline{\overline{F}}^{}=_𝐤F^{}(K)`$, and similarly for their complex conjugate. Note these two equations are formally identical to Eqs. (37a) and (37c), except that the current $`T_\mathrm{\Delta }^{}`$ contains the factor $`\mathrm{\Delta }`$ already.

Now with the new definition $`i\stackrel{~}{\omega }=i\omega \mathrm{\Sigma }_\omega `$, $`i\underset{¯}{\overset{~}{\omega }}=i\omega \mathrm{\Sigma }_\omega `$, $`\stackrel{~}{\mathrm{\Delta }}_𝐤=\mathrm{\Delta }_𝐤+\mathrm{\Sigma }_\mathrm{\Delta }`$, and $`\stackrel{~}{\mathrm{\Delta }}_𝐤^{}=\mathrm{\Delta }_𝐤^{}+\mathrm{\Sigma }_\mathrm{\Delta }^{}`$, as well as $`\mathrm{\Sigma }_\omega =n_iT_\omega `$ and $`\mathrm{\Sigma }_\mathrm{\Delta }^{}=n_iT_\mathrm{\Delta }^{}`$, Eqs. (61b) for $`G`$ and $`F`$ remain valid.