{'fraction_non_alphanumeric': 0.04211212886027848, 'fraction_numerical': 0.015426934747469019, 'mean_word_length': 4.185882352941176, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coherent Potential Approximation for `d - wave' Superconductivity in Disordered Systems

A Coherent Potential Approximation is developed for s-wave and d-wave superconductivity in disordered systems. We show that the CPA formalism reproduces the standard pair-breaking formula, the self-consistent Born Approximation and the self-consistent T-matrix approximation in the appropriate limits. We implement the theory and compute T_c for s-wave and d-wave pairing using an attractive nearest neighbor Hubbard model featuring both binary alloy disorder and a uniform distribution of scattering site potentials. We determine the density of states and examine its consequences for low temperature heat capacity. We find that our results are in qualitative agreement with measurements on Zn doped YBCO superconductors.

with measurements on Zn doped YBCO superconductors.

I. INTRODUCTION
A treatment of disorder is an essential part of the theory of superconductivity. After all, one must explain why impurity scattering does not cause resistance. Thus it is natural that as evidence for novel superconducting states multiplies the foundations of the subject, due mainly to Anderson 1 and Abrikosov and Gorkov 2,3 , are being re-examined. The experiments which stimulate most strongly the current revival of interest in the problem are those on the high temperature superconductors 4 , which are now universally regarded as 'd-wave superconductors' 5 , and those involving some of the heavy fermion systems which display signs of 'p-wave' pairing 6 . In what follows we wish to contribute to the theoretical discussion 7-14 of the issues raised by these very interesting developments.
The case of classic, 's-wave', superconductors is by now well understood. If the perturbation does not break time reversal symmetry and the coherence length is sufficiently long, so that the pairing potential ∆ does not fluctuate, the Anderson Theorem 1 guarantees that there is an absolute gap in the quasi-particle spectrum and the main effect of disorder is that the density of normal states in the gap equation is replaced by its average over configurations 15 . On the other hand if the perturbation breaks time reversal invariance, as is the case with paramagnetic impurities, the effect is more dramatic. For instance, the transition temperature T c is reduced from its clean limit value T c0 , according to the well known pair-breaking formula: where ψ(x) is the digamma function and ρ c = (2πτ T c ) −1 is a measure of the strength of the scattering and τ −1 is the scattering rate 2,7 .
By contrast in the case of superconductors whose Cooper pairs are of exotic 'p-wave' or 'd-wave' character even simple potential scattering, which does not break time reversal symmetry, causes pair breaking 7 . This fact was noted already in the early contributions to the field 17 , but has become a subject of intense scrutiny only recently [7][8][9][10][11][12][13][14] . Of particular interest are two dimensional models featuring 'd-wave' pairing as these may be relevant to experiments on high T c superconductors. Notably, for cuprates many experiments have explored the variation of T c , the density of states and other properties as a function of Ni and Zn substitutions on the copper sites 16,[18][19][20][21][22][23][24] or irradiation damage [25][26][27] . Although a wide variety of theoretical ideas [28][29][30][31][32] and methods [33][34][35][36][37][38] have been applied to interpret the experiments a comprehensive picture of the role of disorder is far from complete. On a more formal level, an intriguing problem arises from the observation of Gorkov and Kalugin 10 that the scattering in models where the order parameter has a line of zeros on the Fermi surface is highly singular and this may be a manifestation of interesting new physics. Indeed in 2-d, Nersesyan et al. 14 predict that the quasi particle density of states N(E) approaches zero, even in the disordered state, as power low, ∼ |E| α , with positive exponent α, instead of going to a finite value where γ measures strength of the interaction, as was found by Gorkov and Kalugin 10 . Another interesting and controversial issue is the relative importance of the self-consistent Born approximation (SCBA) and resonant scattering in the unitarity limit 11,12 . Our aim here is to explore the subject systematically on the bases of explicit calculations, albeit for a simple, extended Hubbard model with attractive interactions and site diagonal randomness only.
In short we will examine the problem of disordered unconventional superconductors making use of the coherent potential approximation (CPA). The CPA is the most reliable approximation developed for the theory of electronic structure of random metallic alloys in the normal state 39,40 . Notably it has been shown to be exact in both the weak and the strong scattering limits, and applicable to systems with low as well as high concentration of impurities. Significantly, the CPA reduces to the self-consistent Born approximation (SCBA) for weak scattering impurities, and agrees with the self-consistent T-matrix approximation (SCTA) results for strongly scattering impurities of low concentrations. Indeed it remains a good approximation in the unitarity limit of resonant scattering 11,12 . Finally, on account of the fact that it becomes exact as the number of nearest neighbors goes to infinity, the CPA is often referred to as a mean field theory of disorder 41 .
Given these desirable features it is clearly worthwhile to explore the consequences of the CPA for disordered superconductors. For the case of conventional s-wave pairing this has already been done, generating many useful results 42,43 . Apart from our earlier brief report 44 and the limited discussion in Ref. 45, the case of superconductors with Cooper pairs of d symmetry will be treated here, within CPA, for the first time.
We will describe fully our numerical method, demonstrate that in various limits our formalism reproduces many of the well known results for disordered superconductors, and examine in detail the phase diagram of the local and non-local attractive two-dimensional where c † iσ and c iσ are, respectively, the usual creation and annihilation operators for electrons on site i with spin σ, and the local charge operator isn i =n i↑ +n i↓ withn iσ = c † −iσc iσ . The chemical potential is µ, t ij are the hopping integrals (for i = j) and ǫ i is the local site energy.
The interaction term, U − ij, can be either be a local attractive interaction (U ii < 0) giving rise to s-wave pairing, or a non-local attractive interaction (U ij < 0 for i = j) giving rise to d-wave or extended s-wave pairing. Disorder is introduced into the problem by allowing the site energies, ǫ i , to vary randomly from site to site.
For computational convenience we shall take the hopping integral t ij to be non zero only when the sites i and j are nearest neighbors. The mean field pairing potential ∆ ij can either be local (i = j) or nearest neighbor non-local, see Fig. 1, and couples particle and hole amplitudes between sites i and j. Of course, the above equations are completed by the self-consistency condition that where η is a positive infinitesimal. To simplify matters we have assumed that the normal Hartree and exchange terms can be absorbed into the definitions of the chemical potential, µ, or the hopping integrals t ij . As usual, Eqs. (4) and (5) are to be solved subject to the requirement on the chemical potential that where n is the average number of electrons per unit cell. Clearly, the Greens function matrix G(i, j; ıω n ) determined by the above equations depends on the set of site energies {ε i }. Our task is to find the configurationally averaged Greens function matrix G(i, j; ıω n ) . Evidently, this is made much easier if we assume that the pairing potential does not fluctuate from configuration to configuration. As was argued by Gyorffy et al. 15 this is a good approximation when the T=0 coherence length ξ 0 is large. Thus our results will have to be treated with appropriate care when applied to superconductors with short coherence length such as superconducting cuprates.
Let us now proceed to deploy the CPA strategy for calculating the averaged Greens function matrix G(i, j; ıω n ) subject to the self consistency conditions: The first move in deriving the fundamental equations of the coherent potential approximation is to define a coherent medium Greens function matrix G c (i, j; ıω n ) by As will be clear latter G c (i, j; ıω n ) = G(i, j; ıω n ) and hence Σ 11 (ıω n ) and Σ 22 (ıω n ) are the diagonal components of the usual self-energy. Note that we did not introduce any off diagonal self-energies such as Σ 12 (ıω n ) and Σ 21 (ıω n ) because for the single site perturbations of our model they are zero. The next step is to consider the scattering of the quasi-particles, propagating according to G c (i, j; ıω n ) by the defects described by the potentials: where l labels one of the m different site energies we wish to consider.
In a straightforward application of the CPA principles, 40 Σ(ıω n ) and therefore G c (i, j; ıω n ) is determined by the condition that these defects do not scatter on the av- where and m l=1 c l = 1.
From Eqs. (11,12) it is now possible, in conjunction with equations (7-9), to calculate Σ(ıω n ) and G c (i, j; ıω n ). The numerical methodology for calculating G c (i, j; ıω n ) and Σ(ıω n ) closely follows that in Ref. 48 and is outlined in Appendix A.
A number of recent studies of superconductors with unconventional pairing suggest that the consequences of disorder depend sensitively on the models used to describe the randomness. 49,50 Thus we are going to investigate two different models. The first corresponds to binary alloy disorder, where m = 2. Namely we consider two types of sites with site energies ε 1 and ε 2 and concentrations of c and 1 − c respectively. The second model is described by a uniform distribution of site energies. Here we shall have in mind the limit . Consequently in Eq. (11) the sum l becomes the integral In the bimodal case, where m = 2, we can simplify Eq. (11) to find where |ε 1 − ε 2 | = δ.
In the second, uniform distribution, case we transform the sum in equation (11) into an integral so that After some straightforward algebra this leads to Thus our CPA calculations will consist of solving numerically either Eq. (14) for the bimodal distribution of the site energies, or Eq. (16) for the case of uniform distribution, to determine the self-energies Σ 11 (ıω n ) and Σ 22 (ıω n ).

III. PAIR BREAKING FORMULA IN CPA
We now relate the CPA formulae derived above to the usual results of disordered superconductors, corresponding to the well known pair breaking formula Eq. (1). As is well known 7 the pair breaking formula was first derived for magnetic impurities in s-wave superconductors 2 but it also applies in many other interesting circumstances such as our present concern, namely the case of non-magnetic impurities in d-wave superconductors 17 .
To derive it within the CPA let us start with the gap equation As a motivation for our argument we recall the method of Abrikosov and Gorkov 2 for solving the gap equation at T c for a clean superconductor. In that case, to find T c we linearize the analogue of Eq.(17) by approximating the off diagonal Greens function G c 12 as follows: where ξ q = ε q + µ, and for our tight binding model with a square lattice ε q = −2t(cos(q x ) + cos(q y )). Then, we note that the kernel of the linear integral equation for ∆( k) is a four term degenerate kernel: where η k = 2(cos(k x ) − cos(k y )) and γ k = 2(cos(k x ) + cos(k y )). Consequently, the general ∆( k) will be a linear superposition of η k , γ k , sin k x and sin k y . However, when the internal symmetry of the singlet Cooper pair is pure d-wave we may take ∆( k) to be of the form Then the condition for non-zero order parameter becomes Let us now define a d-wave weighted density of states and write the above condition, which determines the transition temperature T c0 , as where ω n = πT c0 (2n + 1). In the above equation the integral and the sum are divergent, so we need to introduce a cut-off, ω c n . In the usual way we assume the density of states N d (E) is slowly varying up to the cut off energy, so we will make the approximation N d (E) = N d (0). Then, considering that we can write and hence rewrite Eq. (23) as This is the BCS result for the superconducting transition temperature in the case of d-wave pairing 46 . It differs from the conventional result only in that the d-projected density of states N d (0) has replaced the usual full density of states N(0).
Let us now return to disordered superconductors and examine how the above well known argument is modified when the randomness is dealt with within the CPA. Using Eq. (9) it can be easily seen that instead of Eq. (18) we should use .
to linearize Eq. (17) at T c . Thus, noting that the condition which determines T c can be written as . (29) Now, at this point we need to know the form of Σ 11 (ıω n ) to progress any further. As a first approximation we assume that the most important component to the self-energy is the component at the Fermi energy E = E F = µ. Later on we will test the accuracy of this approximation by examining our numerical results for Σ 11 (ıω n ). For now, however, let us proceed by taking Evidently this leads to .
Again taking N d (E) outside of the integration as N d (0) and performing the integration over where again the sum is cut off, as in the clean limit, by ω c n . If we now add and subtract the terms corresponding to Σ 0 = 0 (the clean case) we find Clearly the term and the second sum is convergent. Hence the cutoff ω c n can be extended to infinity. As has been noted frequently before this infinite sum can be readily performed 2 and we find ln where Eqs. (34) and (35)  Evidently, the self energy at the Fermi energy, E − µ = 0, Σ 0 , gives a good description of the pair breaking parameter ρ c via Eq. (34). Also it is clear that, as expected, the Self-consistent Born Approximation Σ 0 ≡h τ = πδ 2 N(0) only works well in the weak scattering limit.

IV. ANALYTICAL FEATURES AND PREDICTIONS OF CPA EQUATIONS
In this section we examine various analytically accessible limits of the CPA formalism described above. Firstly, we demonstrate that Anderson's theorem is obeyed for s-wave superconductors and the CPA equations are consistent with the results of Abrikosov and Gorkov 2 . Secondly we show that for d-wave superconductors the quasi particle density of states at the Fermi energy, N(0), is non-zero in the presence of non magnetic disorder scattering and is consistent with the Gorkov-Kalugin formula, Eq. (2).

A. The Anderson's Theorem in the Coherent Potential Approximation
Formally, the CPA Eqs. (11) , ( 12) can be written in terms of renormalized Matsubara frequencies,ω n , pairing parameter,∆ k , and particle energiesξ k . This quantities are defined as followsω consequently and for the alloy type disorder with c = 0.5 and ε = ±δ/2, the self energy Σ 11 (ıω n ) which renormalizes ω n , ∆ k and ξ k , and is defined by Eq. (14) can be written as The alternative expression for Σ 11 (ıω n ) in the case of a uniform distribution of local potentials, −δ/2 < ε l < δ/2 is given in Eq. (16).
Note that in the case of a non isotropic, d-wave gap, ∆ k , is not renormalized if the disorder is diagonal both in site and Nambu space. This is different from the case of local s-wave pairing where, as can be readily shown, ∆ k is renormalized by the same factor asω n in Eq. (36). Thus for conventional superconductors, in contradiction to Eqs. (36) -(39), we find that the CPA yieldsω in agreement with Born approximation or Abrikosov-Gorkov theory 2 . As is widely is a case where the symmetry of ∆ k is of extended s-wave symmetry s * ∝ (cos(k x ) + cos(k y )) type. We shall encounter this interesting circumstance later on in this paper. As it turns out for a given gap parameter∆ k = ∆η k and in the limit of small disorder δ → 0 the CPA equations can be solved analytically. To affect the solution note that in Eq. (39) the major part of the summation is coming from the four singular points in the Brillouin zone where the denominator vanishes: Linearizing around these points and performing the summation over k analytically, we find that G c 11 (0) = where α = 2t∆π. Clearly in the limit |ImΣ 11 (0)| << 4∆ this leads to Moreover, when Σ 11 (ıω n ) is small compared to the band width we can rewrite Eq. (40), using the Self-Consistent Born Approximation as On substituting this result into Eq. (44) the later becomes an equation for G c 11 (0) which determines the density of quasi particle state N(0) via the formula N(0) = 1 π G c 11 (ıω n = 0).
Indeed we find and

V. LOCAL QUASI-PARTICLE DENSITY OF STATES CALCULATIONS
In this section we present results for the single site local quasi particle density of states calculations. As mentioned before we consider two types of disorder, (i) binary alloys where we have two on-site potentials randomly distributed throughout the lattice and (ii) a uniform distribution of random on-site potentials. For both types of disorder we have solved, numerically, Eq. (9) in conjunction with its appropriate self-consistency conditions, for the order parameter, ∆, the average number density, n, and the self-energy, Σ 11 (ıω n ).
The first situation we consider is a binary alloy of random on site energies ε 1 and ε 2 with equal concentration c = 1 − c = 0.5. The parameter we have chosen to use to describe the strength of the disorder is δ = ε 1 − ε 2 . Figure 5 Before turning to the problem of disordered d-wave superconductors for reference we have studied, briefly, the s-wave case. In short, we have introduced a site diagonal, local, attractive interaction with strength U, into the above model. As expected such interaction leads to conventional s-wave pairing. When implementing the CPA we have assumed that in this case the pairing order parameter ∆, which is now site diagonal does not vary from site to site. This is consistent with the Anderson theorem which was shown to be adequate for systems with a long coherence length. 15 Fig. 7(a) shows the s-wave quasi-particle density of states as calculated for various values of the binary alloy disorder strength δ. The results confirm that, for s-wave pairing, the gap is absolute and, whilst the edges may move, they remain well defined as required by Anderson's Theorem 1 . In Fig. 7(b) we have also plotted the self-consistent self-energy as a function of the quasi-particle energy, for the same disorder strengths used to obtain Fig 7(a). From this graph we can see that −ImΣ 11 (E) is zero inside the gap and hence the normal disorder does not act as a pair breaker, for s-wave superconductivity.
The quasi-particle density of states for d-wave superconductors is dramatically different from the above BCS spectrum in the s-wave case even in the clean limit. As it is well known it has the characteristic v-like dip shown in Fig. 8(a), for ε near the chemical potential (µ = 0). Upon introducing binary alloy disorder into the d-wave system one would naturally expect similarly different behavior from that described above and in Fig.   7. As is clear from the results reported in Fig. 8 this is indeed the case. Strikingly, as the analytic results of the previous section suggested, N(0) becomes non-zero for the slightest disorder. This implies gapless superconductivity in contrast to the gaped one in the s-wave case.
In Fig. 8(b) we have plotted the self-energies, corresponding to the quasi-particle density of states results presented in Fig. 8(a). This shows a complex evolution with disorder, for small δ the imaginary part of the self-energy, Σ(E), reflects the pure d-wave density of states (as expected in the SCBA limit) and hence Im[Σ(0)] ∼ 0. However increasing δ leads to a finite Σ(0), with a cusp like minimum in Im[Σ(E)] at E = 0. Increasing the disorder even further, to δ = 2.8t, the d-wave pairing is completely destroyed, and ImΣ(E) reverts to the normal system self-energy. In this case −ImΣ(E) is a maximum at E = 0, since the Fermi energy, E F = µ, was set exactly to 0.

VI. CRITICAL TEMPERATURE CALCULATIONS
To analyze the effect of the disorder upon T c we have solved the gap equation together with the CPA equations. Again for the sake of comparison we have developed the analogous theory for the conventional s-wave superconductors based on the site diagonal particleparticle interaction, of strength U 0 , mentioned briefly in a previous sections. In this case neglecting the spatial fluctuations of ∆ i , and linearizing the gap equation for the configurationally averaged single site pairing potential ∆ we find the condition for the temperature T c , below which ∆ = 0, to be and in the second extended s-wave case, T s c is given by where η k and γ k refer to s-wave and d-wave like 'harmonics' defined previously following Eq. 19.
In Fig. 10  To compare the effect the two types of disorder have upon T c we have, in Fig. 12, calculated T d c vs. δ/t at half filling for both binary alloy (solid line) and uniform (dashed line) disorder. From this graph we can see that the two types of disorder effect T c in almost exactly the same manner. If we rescale the graphs such that the standard deviation in the disorder strength is the same, for uniform disorder σ = δ 2 √ 3 and for binary disorder σ = δ 2 , then the two curves lie very close to one another. This prompts the conclusion, which may however be premature, 50 that T c is not very sensitive to the details of the fluctuation in the site energy.
Other interesting point to investigate concerns the relative robustness of T c against degradation by disorder in the cases of conventional s and d-wave pairing. 37 Fig. 13 shows the results for T c for inter-site d-wave and on site s-wave superconductors versus alloyed disorder strength δ, where n = 0.6. Here one can easily recognize a typical difference between these two superconducting states. Clearly in case of a d-wave superconducting state disorder acts as an effective pair breaker while for s-wave superconductors it decreases T c only slowly or not at all. An interesting physical consequence of this effect is that if both U ii and Uij (i = j) are attractive increasing disorder could lead to a transition from d-wave to s-wave pairing, as suggested by Abrikosov 37 .
Evidently, the binary alloy disorder has certain features not shared by the uniform distribution of site energies model. One of these is the possibility of band splitting. For a very large δ parameter the band is split into two sub-bands, as shown in Fig. 14. The on site s-wave superconducting state is created within one of sub-bands. For disorder parameters: c = 0.5 and ǫ 1 = −ǫ 2 = δ/2 we have not found any similar solution for d-wave, but for c=0.25 and n=0.5 a single sub-band d-wave solution is found. In Fig. 15 we show the very different behavior of T c for a small and a large non-local interaction parameter, U ij . For large enough U ij the d-wave T c does not go to zero when the band splits. The quasi-particle density of states corresponding to such a possibility is presented in Fig. 16 and, as in Fig. 14, we see the superconducting gap, which is now d-wave, has formed in one of the sub-bands.
To complete the discussion Figs. 17, 18 and 19 show the critical temperature plotted versus concentration c for various strength of binary alloy disorder δ with band fillings of n = 1 (Fig. 17), n = 0.7 (Fig. 18) and n = 0.1 (Fig. 19). and ImΣ 11 (ε) as disorder is added to the problem is also of general interest and, as it turns out, of lively controversy. [10][11][12][13][14]51,52 In this section we wish to examine those predictions of the CPA calculations which are relevant to these issues in particular. Using the methods outlined in the previous sections we calculate numerically the quasi-particle local density of states, N(ε), and the self-energy, ImΣ 11 (ε), and then investigate how these two quantities change with both temperature and disorder.
For simplicity we have studied the half filled band (µ = 0, n = 1). This is the band filling for which T c , in the d-wave case, is a maximum. We will examine the effects of alloyed disorder on the system and show that as we increase the amount of scattering the specific heat vs. temperature relation changes from T 2 , in the clean limit, to T in the limit of strong scattering. As is widely recognized 5 this behavior is the consequence of N(E) ∼ |E| changing to N(E) = constant and is consistent with experiments. 8 To get an impression of the form of N(E) and −ImΣ 11 (E) in the region of the chemical potential, µ, we fit them to the function a + b|E| c in the energy range −∆E < E < ∆E, where ∆E is small compared to the gap. Using the coefficients a, b and c to fit the curves N(E) and −ImΣ 11 (E) in the region of µ gives us a tool to analyze their functional form.
For example a tells us if the curves are finite at E = 0, b controls the gradient of the curve and c controls the curvature.
In Figs. 20(a-c) we have plotted these coefficients, for N(E), as a function of T , for different values of δ. We also included in Fig. 20(d)  Firstly we note that in Fig. 20(a) the parameter a(T ) tends to a finite limit as T → 0 and this limit increases more and more rapidly as the disorder described by δ, increases. Hence N(0) is finite in agreement with our earlier discussion in Sec. IV where we derived for low scattering and low temperatures the dependence of N(0) on δ (see Eq. (46)). It lends credit to the general consistency of our results that the a(T ) ∼ = N(0, T ) curves rise to their normal state values as T → T c . Interestingly, as can be seen in Fig. 20(b) the gradient of N(E), namely b, changes dramatically only near T c . In fact for the case where δ is small (squares) the gradient even changes sign. This means that we go from having a gap in the density of states, below T c to having a Van Hove singularity above T c . The curvature, represented by c and plotted in Fig. 20(c), shows that for low disorder (squares) as we increase the temperature we go from a curvature of c ≤ 1, ie. a cusp, to a curvature c > 1. Finally when the critical temperature is reached for each of the four different disorder types c goes to 2.
In the same manner, we now examine the corresponding coefficients, a ′ , b ′ and c ′ for −ImΣ 11 (E). In Figs. 21(a-c) we have plotted the calculated coefficients for −ImΣ 11 (E), at different values of δ as a function of the temperature. In Fig. 21(a) we can see that at low temperatures −ImΣ 11 (0) increases with disorder, again this agrees with our results in Sec.
IV where an analytical form for the dependence of −ImΣ 11 (0) upon δ was derived at T = 0, see Eq. (47). As we increase the temperature we can see the self-energy at the chemical potential also rises until the normal state value is reached at the critical temperature. In Fig. 21(b) we can see how the gradient changes from being large for large δ and small for small δ at low T , and at temperatures greater than T c to be large and negative for large δ and small and negative for small δ. Finally Fig. 21(c) shows how, c, the curvature goes from being almost linear, c = 1, at low temperatures, regardless of δ, to c = 2 above T c .
To illustrate the consequences of these results for the low temperature thermodynamic properties we have calculated the effect of disorder upon the specific heat. Using the above temperature dependent coefficients for the density of states near ε ∼ 0 we calculated the limiting behavior of the specific heat as T → 0. The results are shown in Fig. 22 for different strengths of disorder δ. In the case where δ is small (squares) we can see that the specific heat has a T 2 dependence (the plotted continuous solid line with no points) and for large δ the dependence upon temperature is linear, as expected. 6

VIII. CONCLUSIONS
We have compared and contrasted the effect of disorder on conventional s-and d-wave superconductors on the bases of an extended, negative U Hubbard model and a mean field, CPA, treatment of disorder. On the one hand we have derived many of the well known results, such as the pair breaking formula in Eq. (34) or that for the quasi-particle density of states N(0), Eq. (46). On the other hand we have solved the Gorkov-CPA equations Eqs. (14), (17) numerically and surveyed the salient features of their consequences by explicit calculations. The use of CPA in describing disordered d-wave superconductors is an advance in this very active field 4-14,49,50 because it allows us to avoid the usually delicate choice between methods, sets of diagrams, designed to deal with either weakly or strongly scattering perturbations. The CPA treats both kinds of problems equally accurately and it is known to provide a very credible interpolation between the two 40 . As an example where above feature of CPA may have a crucial role to play we recall the use of the resonant scatterer model in interpreting experimental data both on the cuprates 8,53 and some Heavy Fermion systems. 11,54 In short, note that in the impurity, c → 0, limit the self energy for the Greens function describing an electron moving on a lattice, is given by where the T -matrix is defined, in terms of the impurity 'potential' V imp , a real local shift in the site energy, as Evidently T (0, 0; ǫ) is complex number with an amplitude: and a phase, the phase shift φ(E), given by Now, observe that while for weak scatterers for a resonant scatterer in the unitary limit, defined by φ(E) = π 2 , where N(E) = − 2 π ImG(0, 0; E). It is the above dramatic difference in the dependence of ImΣ B and ImΣ R on the density of states N(E) that the cited authors rely on in interpreting the relevant experimental data. Evidently, since the individual scattering events described by the local T -matrices are always treated exactly in the CPA, the CPA describes weak scatterers and resonant scatterers equally well. Moreover, since it is a reliable approximation not only in the impurity limit, c ∼ 0, but also for arbitrary concentrations it deals with resonance scatterers even when Eq. (55) no longer holds. Thus CPA should be the preferred treatment for models with strong, even resonant scattering.
Although we have not specifically concentrated on this aspect of the method, the principle feature of Eqs. (54) and (55) , namely the dependence of Σ(E) on N(E), can be seen to be at work in our calculations of the previous section. To demonstrate this we have calculated Σ(0) and N(0) as functions of band filling, n, in the most interesting region n ∼ 1.0, and compared in Fig. 23 their relationships for weak and strong scatterers at c = 0.5. Clearly, for weak scatterers δ/t = 1, −ImΣ(0) is a more or less linear function of N(0) as in Eq.
(54), while for strong scatterers δ/t = 3, ImΣ(0) is inversely proportional to N(0) as in the resonant scatterers model described by Eq. (55). Thus we conclude that the CPA employed in the calculation gives a reliable account of disorder in both the weak and strong scattering regimes.
Having listed the above desirable properties of the CPA we hasten to emphasize that it is a 'mean-field' theory of disorder and, hence, does not describe such interesting phenomena as localization 53 even in the normal state. 54 Consequently, our result that N(0) is finite for the smallest amount of disorder can not be taken as evidence against the conclusion that N(E) ∼ |E| α of Nersesyan et al.. 14 As this originates from the divergence of the vertex corrections in perturbation series for Σ(E) we may conjecture that it has to do with where the indices i and j denote sites, while α and α ′ represent the particle or hole degree of freedom on each site. We denote particle degrees of freedom by α = 1 and hole degrees of freedom by α = 2. For example G c 1 2 (i, j; E) represents the Greens function between the particle degree of freedom on site i and the hole degree of freedom on site j at energy E.
To compute the Greens functions (A1) we can closely follow the method described by Martin and Annett. 48 Using this method we can transform the Hamiltonian to a block tridiagonal form where a n and b n are 2 × 2 matrices. Given this form for iα|E1 − H|jα ′ and expressing the Greens function as the Greens functions above can be evaluated as a matrix continued fraction so that Within Eqs. (A2) and (A4) we have a formally exact representation of the Greens functions. However in general both the tridiagonal representation of the Hamiltonian, and the matrix continued fraction (A2) will be infinite. In practice one can only calculate a finite number of terms in the continued fraction exactly. In the terminology of the recursion method it is necessary to terminate the continued fraction. 55,57,[59][60][61][62][63] If we were to calculate up to and including a n and b n and then simply set subsequent coefficients to zero then the Greens function would have 2n poles along the real axis. The density of states would then correspond to a set of 2n delta functions. In order to obtain accurate results it would be necessary to calculate a large number of exact levels, which would be expensive in terms of both computer time and memory.
As a more efficient alternative we choose to terminate the continued fraction using the extrapolation method, as used previously by Litak, Miller and Györffy. 57 We calculate the values for a n and b n exactly up to the first m coefficients using the recursion method. Then, noting the fact that the elements of the matrices a n and b n vary in a predictable manner, 57 we extrapolate the elements of the matrices for a further k iterations, where k is usually very much greater than m. This enables us to compute the various Greens functions accurately with relatively little computer time and memory.
We now have a method for calculating the Greens Function G c α α ′ (i, j; E). To calculate our system of equations self-consistently for a given site i, we need to know G c 1 2 (i, j 1 ; E), G c 1 2 (i, j 2 ; E), G c 1 2 (i, j 3 ; E) and G c 1 2 (i, j 4 ; E). If the interaction is non-local, else if the interaction is local then we only need to calculate G 1 2 (i, i; E), see Fig. 1. To calculate n we also need G c 1 1 (i, i; E).
In performing our self-consistency calculation we can make use of the fact that Σ(E) just acts as a renormalization to E →Ẽ (see Sec. IV). This means we do not need to recalculate a n and b n when we calculate the self-consistent self-energy. We first calculate the order parameter self-consistently for a given attractive interaction, then we calculate G c (i, i; E) and hence the self-energy self-consistently (remembering that in this self-consistent loop we do not need to recalculate a n and b n , since the self-energy does not affect them). Then we recalculate the superconducting order parameter self-consistently, incorporating the calculated self energy and repeat this procedure until both the order parameter and the self energy have reached stable self-consistent solutions.