{'fraction_non_alphanumeric': 0.03373654792692083, 'fraction_numerical': 0.010678234754317176, 'mean_word_length': 4.565604977249513, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 11, '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}}
Local Electronic Structure of Defects in Superconductors

The electronic structure near defects (such as impurities) in superconductors is explored using a new, fully self-consistent technique. This technique exploits the short-range nature of the impurity potential and the induced change in the superconducting order parameter to calculate features in the electronic structure down to the atomic scale with unprecedented spectral resolution. Magnetic and non-magnetic static impurity potentials are considered, as well as local alterations in the pairing interaction. Extensions to strong-coupling superconductors and superconductors with anisotropic order parameters are formulated.


I. INTRODUCTION
Low-temperature superconductors almost always have a high concentration of non-magnetic impurities. Even in the dirty limit, however, where the mean free path is shorter than the coherence length, superconductivity endures 1 . This phenomenon can be understood by generalizing BCS pairing 2 to be between degenerate Kramers partnered states in a time-reversal invariant system 3 . Magnetic impurities, which do break time-reversal invariance, have more profound effects on the superconductor in dilute concentrations than non-magnetic impurities, lowering the critical temperature T c 4-6 and producing localized states within the gap 7-13 which at sufficient concentrations hybridize to produce gaplessness 14 . In the course of the investigation of the effects of impurities on superconductivity during the last four decades, the primary emphasis has been the influence of impurities on bulk properties. These effects have been treated within a strong-coupling formalism (e.g. Ref. 15), but only very recently self-consistently and beyond the Born approximation 16 . The above work was primarily concerned with bulk or impurity-averaged characteristics and ignored the spatial structure of electronic properties near to the impurity.
Among the first local properties calculated in the vicinity of an impurity in a superconductor were the structures of screening clouds around a charged impurity 17,18 and a magnetic impurity 18,19 in a superconductor (characterized by exponentially-decaying Friedel-like 20 oscillations). The oscillation of the order parameter around a magnetic impurity was first evaluated without selfconsistency [21][22][23] . A self-consistent calculation of the order parameter at the impurity and very far away for weak impurity potentials was reported by Schlottmann 24 for magnetic impurities and Fetter 25 for non-magnetic impurities.
Interest in local properties near impurities in superconductors has been revived by the capability of scanning tunneling microscopy (STM) to perform localized spectroscopic measurements. The spatially-dependent differential conductivity around a non-magnetic impurity at the surface of a superconductor has been considered theoretically 26 for both isotropic and anisotropic order parameters. This differential conductivity is proportional to the local density of states (LDOS) around the impurity. In Ref. 26 the impurity was modeled as a point defect, and spatial oscillations in the LDOS at various voltages were calculated. These oscillations are the superconductor's analog of oscillations in the LDOS created by step edges and defects on noble metal surfaces 27,28 . The superconductor's LDOS oscillations would allow one to measure the anisotropy of the superconductor's order parameter. The conductance oscillations for voltages just above a gap minimum or maximum are strongly pronounced in the real-space directions corresponding to the momenta of the gap minimum and maximum. Calculations followed which considered sharp energy features in the scattering process, such as resonant states 29,30 . Nevertheless, an important assumption of these calculations has remained unchecked in detail, that the electronic distortions induced by the impurity are local, including the deformation of the order parameter. Self-consistent calculations using the Bogoliubov-de Gennes (BdG) equations 31 followed for two-dimensional systems [32][33][34][35] , but have been hampered by finite-size effects.
A magnetic impurity differs from a non-magnetic impurity in that a localized state exists around it [7][8][9][10][11][12][13]36 . The first calculations of the LDOS of the localized and continuum states around a magnetic impurity were performed recently both with a simplified analytic model and numerically via the new self-consistent technique 37 described in this paper. Calculations of the LDOS of the localized state with angular momentum quantum number ℓ = 0 were performed with a slightly different analytic model and numerically in two dimensions via the self-consistent BdG equations 38 . These calculations were motivated by preliminary experimental results around a Gd atom on a niobium surface 39 .
Another type of defect is a local change in the pairing interaction. The resulting spatially-dependent order parameter can then distort the LDOS. An order parameter suppression can even localize states, as in a vortex core 40,41 .
The local electronic properties of all of these defects can be calculated self-consistently from the Gor'kov equation 42 without further approximation with the new technique introduced in Ref. 37. Our technique for calculating the electronic structure around a defect in a superconductor is related to the Koster-Slater inversion techniques for determining the local electronic structure of impurities in metals 43,44 . Since its original application to localized vibrational modes 45 , this algorithm has been applied to numerous problems including deep levels in semiconductors 46 and impurity states in magnets 47 . The new Koster-Slater technique separates space around the defect into two regions: the near field and the far field. The far field is a region distant enough from the defect that the potential is insignificant and the order parameter has relaxed back to its homogeneous value. The near field is the region close to the defect where the potential is finite or the order parameter is distorted. In essence, the Gor'kov equation that determines the Green's functions of the inhomogeneous superconductor is inverted in the real space region of the near field. This paper describes the technique in detail, expands on an analytic model introduced in Ref. 37 and reports several calculations of the properties of the defects described above.
Previous attempts to calculate the structure of the near field have often used the BdG equations 31 . Other theories of inhomogeneities, such as Ginzburg-Landau theory 48 or the Eilenberger equations 49 , treat the spatial degrees of freedom as coarse-grained over the superconductor's coherence length. Coarse-grained approximations are not appropriate for considering electronic structure on the atomic scale near a defect. The BdG equations are generalized Schrödinger equations for the electron and hole wavefunctions of a quasiparticle, and are valid for a superconductor with an arbitrarily-varying order parameter, only constrained by the validity of BCS theory.
Unfortunately, these equations have significant practical difficulties as well. Despite qualitative success modeling STM measurements of a single vortex in superconducting NbSe 2 50 , calculations of the electronic structure 51-55 using the BdG equations are hampered by the difference in energy scales between the Fermi energy and the order parameter. Since the BdG equations are solved numerically for a finite system, the difficulty of the calculation is determined by the necessary spectral resolution. The key energy scale which must be resolved is the superconducting gap. Thus, the numerical difficulty increases as the ratio of the Fermi energy to the gap becomes large. Hence the band structure assumed for the superconductor must be somewhat unrealistic (for Refs. 54, 55 the Fermi wavelength was approximately 100Å, which is inappropriately large 56 ). This limitation extends to calculations of the interaction between a vor-tex and an impurity 57 , the characteristics of the vortex lattice 58,59 , and work on a non-magnetic impurity [32][33][34][35] and a magnetic impurity 38 in a two-dimensional s-wave or d-wave superconductor. In contrast, the computational requirements of the Koster-Slater technique are determined by the range of the impurity potential, rather than the necessary spectral precision.
In the Section II of this paper we first describe the BdG formalism for local defect potentials and then compare with the new, Koster-Slater formalism. Section III describes an analytic model, based on a delta-function potential, which reproduces some of the quantitative behavior of the numerical results. Section IV discusses the results of the numerical calculations for magnetic impurities, non-magnetic impurities, impurities incorporating both magnetic and non-magnetic potentials, and inhomogeneities in the pairing interaction. A heuristic picture of the electronic structure near these impurities will be presented here, and the calculations will be compared with the analytic model of Section III. Section V generalizes the formalism of Section II to the case of strong-coupling and anisotropically-paired superconductors.

A. Bogoliubov-de Gennes Equations
To place our new formalism in context we will contrast it with the BdG equations, which are Schrödinger-like equations for the electron and hole components of the quasiparticle wavefunctions u(x) and v(x) respectively. These are, for a free electron band structure with mass m, the positive-energy (E) solutions to (1) Here σV S is a position-dependent, spin-dependent potential, such as one originating from an impurity with a classical spin. V 0 is a position-dependent non-magnetic potential and ∆(x) is the inhomogeneous order parameter. ∆(x) can be chosen real since the defect potential is real 60 . The quantization direction of the electronic spins in the superconductor (σ = ±1/2) is chosen parallel to the classical spin. A classical spin has no quantum dynamics, and cannot flip the quasiparticle spin. Hence spin is a good quantum number for the quasiparticles and only two coupled equations (Eqs. (1)) are required. The combinations σV S ± V 0 have physical significance: σV S + V 0 is the potential felt by an electron of spin σ, while σV S − V 0 is the potential felt by a hole of spin σ. The spatially-dependent order parameter is determined self-consistently: where n labels the states for each spin σ, T is the temperature, and k B is Boltzmann's constant. γ(x, E nσ ) is the effective electron-electron interaction potential, which is For a spherically-symmetric defect the wavefunctions are eigenstates of the angular momentum with quantum numbers ℓ and m. Typically the defect is placed in a sphere of radius R with appropriate boundary conditions. The value of R is determined by the spectral resolution necessary for accurately evaluating Eq. (2) and the spectral width of features measurable by (for example) the STM. The typical complications resulting from approximating an infinite system by a finite-size system apply, such as discrete states above the energy gap and the heavy investment of computer time required for large values of R. For example, in the calculations for the vortex in NbSe 2 54,55 , ǫ F /∆ o = 32 was the largest ratio of the Fermi energy to the homogeneous order parameter considered. This value is unrealistic, and is a result of inappropriately fitting the coherence length and upper critical field of NbSe 2 to a free-electron model. A more realistic band structure 56 has a bandwidth to energy gap ratio at least an order of magnitude greater, and a highly anisotropic Fermi velocity. A system where the free-electron model is appropriate, such as niobium, has an ǫ F /∆ o = 705 and is numerically inaccessible for general potentials V S and V 0 .

B. Koster-Slater Inversion Formalism
We now introduce a self-consistent method which works within a sphere whose radius is determined by the range of the defect's potential and utilizes the continuum spectrum of the homogeneous superconductor. In essence, we invert the Gor'kov equation in real space. The Gor'kov equation 42 for a defect in a superconductor can be written in the Nambu formalism 61 as: where the inhomogeneous retarded Green's function, .
The elements of this matrix are The explicit subscripts ↑ and ↓ do not refer to the spin of the excitation in the superconductor but rather to the spin band of the normal state used to construct the excitation. The key concept is that the spin-up band contains both up electrons and down holes just as the spin-down band contains both down electrons and up holes. The convention employed here is standard in the theory of semiconductors where a spin-up electron excited above the Fermi energy leaves a spin-down hole below the Fermi energy. This is convenient for magnetic potentials since if spin-up electrons are attracted to a magnetic impurity spin-down holes should be repelled by the impurity. In the presence of a single classical impurity spin, the quasiparticle spin is a good quantum number despite electronhole mixing. Our convention determines the composition of a spin-up quasiparticle to be part spin-up electron and part spin-up hole, rather than part spin-up electron and part spin-down hole.
For ω > 0, G(x, x ′ ; ω) describes spin-up excitations, involving the mixing of electrons in the spin-up band with holes in the spin-down band. For ω < 0, G(x, x ′ ; ω) describes spin-down excitations, involving the mixing of electrons in the spin-down band with holes in the spin-up band. Since spin is a good quantum number it is not necessary to use a 4 × 4 formalism, such as that of Ref. 62. The notation here has been simplified and improved relative to Ref. 37.
The homogeneous Green's function g is independent of σ, so .
The potential where δ∆(x) = ∆(x) − ∆ o and V S , V 0 , and ∆(x) have similar meaning above as in the BdG equations. The factor of one-half from the electron spin has been incorporated into the potential V S . To be concrete, and without loss of generality, the spin direction attracted by the potential will be called spin up, and the spin direction repelled will be called spin down 63 . The self-consistency equation for the order parameter is where n(ω) is the Fermi occupation function at temperature T . Self-consistent spin-dependent and charge-dependent potentials can also be constructed for V(x) using the calculated spatial structure of the spin s(x) and charge ρ(x) around the defect, No calculations are reported for such potentials 31 in this paper.
For a BCS superconductor with an isotropic gap in a parabolic band, for ω much smaller than the Fermi energy, where r = |x − x ′ | and N o is the density of states for each spin at the Fermi level. The coherence length, ξ =hv F /π∆ o , where v F is the Fermi velocity. These expressions are valid for ω above and below ∆ o so long as the imaginary parts of both f and g are multiplied by sgnω.
One strength of this formalism is its reliance on the short-range nature of the defect's potential. Solution of Eq. (4) requires inverting the frequency-dependent realspace matrix The structure of M(x, x ′ ; ω) allows for a precise description of the difference between the near field and the far field. We require that the defect's potential V(x) is zero for |x| > R. The space |x| > R belongs to the far field, whereas the space |x| ≤ R belongs to the near field. We can then separate any real-space matrix A symbolically into four pieces: where n and f label the near-field region and far-field region, respectively. The particular example of M is blocktriangular: It is clear from Eq. (18) that the computational effort in inverting M, and thus finding the inhomogeneous electronic structure, is entirely determined by the complexity of inverting M n→n . We also obtain the useful result that the electronic structure in the far field is easily determined once the electronic structure in the near field is known.
The local density of states in the near field is merely while in the far field it is In the limit that the size of the near-field approaches a point, Eq. (20) is The inhomogeneous part of the right-hand side of Eq. (21) has the following form if the off-diagonal elements of V eff (ω) can be ignored (this is appropriate when self-consistency is not important): Hence the spatially-dependent features reported in Ref. 26 are qualitatively retained in our more complete model. For a magnetic potential in the Born approximation there would be no such signal, since so the distortions in the LDOS of the two spin directions are equal and opposite. The local density of states under ideal conditions is directly measured by a scanning tunneling microscope. The differential conductivity measured at a point x, voltage V and temperature T can be related to the local density of states at the tip location as follows: Here e is the charge of the electron and n ST M (ω) is the Fermi function in the STM tip The local density of states is proportional to the imaginary part of the retarded Green's function fully dressed by the interaction of the electronic system with the impurity. Certain features of the above equations simplify their numerical implementation. Each angular momentum channel constitutes an independent block-diagonal submatrix in M(ω). Since the bare Green's functions in Eqs. (15) have analytic expansions in spherical harmonics, M(ω) can be calculated quickly. These expansions are detailed in the Appendix.
The value of R is governed by the longest-range potential. In this paper that will be δ∆(x). We model V S and V 0 with Gaussians of radius k −1 F and evaluate Eq. (12) at Figure 1. While oscillating with wavelength ∼ πk −1 F , δ∆(x) falls off to a negligible potential within 10k −1 F . A typical radial grid of 100 points provides a numerically robust solution. The self-consistent solution depends on the value of which for niobium is 3.6 × 10 −5 (Ref. 64). This is the single dimensionless parameter required to parametrize the free-electron model of a superconductor.
In contrast to the Koster-Slater technique described here, which exploits a physical distinction between the near field and the far field to accelerate the numerical calculation, the BdG equations treat the near-field and the far-field on the same level. A large R is desired for decent spectral resolution, but the possible size of R is limited by numerical constraints. Hence a numerical implementation of the BdG equations, by comparison, is typically slower and substantially less accurate than the Koster-Slater technique.

A. Magnetic and Non-Magnetic Point Potentials
Approximating the local potential by a delta function leads to a simple expression for Eq. (4), In principle M(ω) can be found from the Green's functions in Eq. (15), however there is a divergence in the real part of g(r; ω) as r → 0. This divergence is coped with in Ref. 8 by discarding the divergent piece. This approximation is essentially an assumption of strict particle-hole symmetry (not merely linearizing ǫ(k) around ǫ F ). We now discuss the effects of this approximation on the local properties of the system in the normal state. The heuristics will be simpler for the normal state, but the conclusions also apply to the superconducting state.

B. Particle-Hole Symmetry in the Normal State
In order to focus on the spin-dependent potential, the non-magnetic potential will be set to zero. The normal state properties can be obtained from Eq. (15) for ∆ o = 0. That yields a Green's function appropriate for an outgoing wave: where ω is considered close to the Fermi surface so the change in momentum due to ω = 0 is negligible. The inhomogeneous local density of states for a delta-function potential using this Green's function is unphysical, since it diverges at r = 0. Now we examine the local density of states when the divergent real part of g is ignored.
This expression yields the pathological behavior that the local density of states near a spin-dependent potential is exactly the same for spin-up electrons as for spin-down electrons. A spatial response in the LDOS to a spindependent potential that is identical for up and down spin electrons only occurs for certain band structures with ǫ F at special energies (e.g. half-filled tight-binding models). The local density of states for v s = πN o |V S | = 0.1 is plotted in Figure 2 as a function of distance from the potential. Although it is somewhat distorted from its homogeneous value, it does not show the spin-dependent asymmetry of a more realistic potential. An exact calculation for a Gaussian of range k −1 F is also shown in Figure 2 for comparison.
A more realistic approach to coping with the divergence in Eq. (29) without yielding the pathological result of Eq. (30) is to average the real part of g over a range given by the assumed range of the potential. That yields a finite value for the local density of states at the potential, but does not control the behavior for small r.
To perform that task, we consider a "muffin-tin" Green's function. This function has the form where α is the average of the real part of g over the range of the potential, and R o is chosen so that the spatiallyintegrated spectral weight of the Green's function is unchanged (as required by probability conservation). We show in Figure 2 the local density of states calculated with this Green's function. In particular, the asymmetry in the response of spin-up electrons and spin-down electrons to a spin-dependent delta-function potential is governed by this phenomenological parameter α. The agreement with the exact solution is good at the origin and far from the impurity. The muffin-tin Green's function is discontinuous, unfortunately, but yields a better approximation of the response of the system than the particlehole symmetric approximation (discarding the real part of the Green's function). It is also possible to generate an asymmetry between spin-up and spin-down electrons by adding a non- Spectral weight in the normal state around a magnetic impurity with vs = 0.1 as a function of the distance from the impurity r. The solid lines are exact solutions for spin-up (attracted by the impurity) and spin-down (repelled by the impurity) electrons for a Gaussian potential of range k −1 F in a metal with a free-electron dispersion relation. The dashed line, which is the same for both spin-up and spin-down electrons, is the calculated spectral weight for the particle-hole symmetric delta-function potential model. The dot-dashed line is the result for a delta-function potential calculated using the muffin-tin Green's functions. The muffin-tin Green's functions fix the pathological result of the symmetric model that the spectral weight is the same for spin-up and spin-down electrons. The muffin-tin parameter, α = 0.704, is determined by the range of the Gaussian potential and is therefore not a free parameter. magnetic potential v 0 with the v s to parametrize the impurity, but still maintaining a particle-hole symmetric band structure. Since the response of a particle-hole symmetric system to a potential does not depend on the sign of that potential, the v 0 is required to distinguish between electrons and holes. An electron feels a potential v s + v 0 , whereas a hole feels a potential v s − v 0 . This additional non-magnetic potential only breaks particle-hole symmetry locally (within the range of the potential), whereas for a realistic band structure the particle-hole symmetry is broken everywhere in the solid.
For a Gaussian potential with range a, For the Gaussian potentials numerically calculated in this paper, a = k −1 F , so α = 0.704.

C. Self-consistency within the Analytic Model
As seen from Figure 1, the distortion of the order parameter is short-ranged around an impurity. We may then consider the effect of the order-parameter distortion on the electronic structure to be parametrized by an effective delta-function potential δ∆δ(x) similarly motivated to the delta-function potentials for the magnetic and nonmagnetic potentials. The potential V(x) is changed in the following way: The relative effect of the δ∆ compared to the other two potentials is likely to be small for the potentials considered in this paper.

D. Energies of localized states in the Superconductor
The energies of the localized states of angular momentum ℓ correspond to the positive energies ω ℓ = |Ω|, where detM n→n (Ω) = 0, (34) and the solution is traced to the ℓ-channel block of M (see Appendix). For the analytic model, M n→n (ω) is the matrix shown in Eq. (28), where g(0, 0; ω) in the superconducting state is constructed similarly to that of the normal state, The anomalous Green's function is given by Eq. (15), since it does not have a divergence problem as r → 0, This analytic model only has localized states in the ℓ = 0 angular momentum channel, as expected for a deltafunction potential. Those energies are where Eq. (37) reduces to a result obtained by Shiba 8 when v 0 = α = δ 0 = 0, a result obtained by Rusinov 9 when α = δ 0 = 0, a result obtained by Salkola, Balatsky and Schrieffer 38 when α = 0, and a result obtained by us 37 when v 0 = δ 0 = 0. One solution of Eq. (37) is a spin-up quasiparticle and the other is a spin-down quasiparticle. There may be only one real solution to Eq. (37); then only one ℓ = 0 localized state exists around the impurity. This occurs for large v s . When v s = 0 the localized states are due to order parameter suppression, and the energies of the two spin states are degenerate. This follows from timereversal symmetry in the absence of a magnetic potential. For small v s the two energies are split by an amount

E. Spectral Weight Asymmetry in the Analytic Model
A spin-up quasiparticle consists of amplitudes for a spin-up electron (electron in a spin-up state), and a spinup hole (electron missing from a spin-down state). Therefore the spectral weight of a spin-up localized state will be divided between an electron-like pole in the spin-up band at ω = ω o (with weight A ↑ (r; ω o )) and a hole-like pole in the spin-down band at ω = −ω o (with weight A ↓ (r; −ω o )). These two types of excitation are independently resolvable by a scanning tunneling microscope since at positive sample voltage relative to the tip, the STM places electrons in the sample, whereas at negative sample voltage the STM places holes in the sample. We define the energy of the pole in the spin-up band to be ω ↑ and in the spin-down band to be ω ↓ . Even though ω o is always positive, ω ↑ can be positive or negative, and ω ↑ = −ω ↓ .
The spatial structure of the spectral weights of the spin-up band and spin-down band components of the localized state are given by The above expression is set up for use as a muffin-tin Green's function. The construction of such a Green's function requires that the spectral weight integrates to one. However, if R o → 0, the above expression can be integrated over all space, and still yields one-half for the spin-up band and one-half for the spin-down band to order ∆ o /ǫ F . Hence the quasiparticle is one-half electron and one-half hole. This provides additional confidence in the above expressions. The frequency-integrated weight at the defect of the two types of excitation can be calculated within the analytic model, and is (for δ 0 = 0) For a spin-down quasiparticle there is an electron-like pole in the spin-down band and a hole-like pole in the spin-up band and the relative weight is still given by the expression above. This expression for v 0 = 0 was reported in Ref. 37. In the following section these results are compared with the numerical calculations of properties in the superconducting state and in the normal state.

A. Non-Magnetic Impurity
Even strong non-magnetic impurities at moderate concentrations will not suppress the critical temperature of a superconductor 3 . Nevertheless, it was recognized early on 25 that the local order parameter may be affected. In Ref. 25 the effect of the non-magnetic impurity was calculated in the far field by modeling the impurity potential with a phase shift. The phase shifts were evaluated for two models: a spherical square-well potential and a deltashell potential. Self-consistency was ignored by only focusing on regions far from the impurity where the change in the order parameter is small compared to its homogeneous value. The order parameter change due to the impurity was found to oscillate with the Fermi wavelength, and decay as r −2 . Figure 3 shows the spectral weight A(r; ω) at several frequencies above the energy gap near a strong nonmagnetic impurity with a Gaussian potential of range k −1 F and strength v 0 = 7/8, calculated self-consistently for a superconductor with ξk F = 449 (niobium). The spectral weights are suppressed to approximately 30% of their homogeneous value at the center of the potential. Only continuum states are shown since no localized states Comparison of the spectral weight in the superconducting state at several frequencies to the spectral weight of the normal state and the spatial dependence of the order parameter. All of these curves are indistinguishable, but they have been displaced for better visibility. The spectral weights are normalized to their value at large distances from the impurity r, and the order parameter is normalized to the homogeneous order parameter ∆o. This figure indicates that normal-state properties drive the spatially-dependent features of the superconductor's spectrum.
were found for this potential. The curves showing the spectral weight have been displaced from each other so that they may be distinguished. Also shown displaced in Fig. 3 is the spectral weight in the normal state, normalized to the spectral weight in the homogeneous metal. All the quantities plotted in Fig. 3 are identical to the accuracy of the calculation. Fig. 3 is an illustration of a relationship between the spectral weight in the normal state and the spectral weight in the superconducting state, where A n (r) is the spectral weight in the inhomogeneous normal state for energies near the Fermi surface and A sc (ω) is the homogeneous superconductor's spectral weight as a function of frequency. 2N o is the normal state's spectral weight far from the impurity. This expression is valid for small r and small ω, the regime of interest for STM on a superconductor. For ω of order ∆ o , Eq. (42) is valid for r < ξ. We further illustrate the relationship of Eq. (42) in Figure 4(a), which shows the LDOS for this non-magnetic impurity as a function of voltage and position calculated from Eq. (24) with T = 2 15 ∆ o /k B . This temperature corresponds to 2K for niobium. There is no change in the energy gap due to this non-magnetic potential. Figure 4(a) shows that it is merely the local amplitude of the spectral weight which is reduced -this would manifest itself in a locally reduced oscillator strength for an optical transition, or the reduction in the tunneling current for an STM, which is directly proportional to the LDOS shown in Figure 4(a). Figure 4(b) is an identical calculation for a shorter coherence length, ξ = 10k −1 F . There appears to be little difference, although here a localized state exists near the continuum. Figures 3 and 4(ab) show for the non-magnetic impurity that the normal-state electronic structure determines the spatial dependence of the superconductor's A(r; ω) for all frequencies including near the energy gap. The potential strength of the impurity is orders of magnitude greater than ∆ o , and thus locally mixes in states far from the Fermi surface in the homogeneous metal. These states are required to construct probability densities which are suppressed by 70% near the impurity. The spatial structure of the spectral weight in the normal state is essentially identical to that seen in Fig. 3 over an energy range around the Fermi surface which is orders of magnitude greater than ∆ o .
Once the normal-state band structure has been distorted by the presence of the non-magnetic impurity, superconductivity is a small perturbation within a narrower frequency range. The formation of the gap in the single-particle excitation spectrum in the superconducting state is characterized by the "mixing" of electron and hole amplitudes to form quasiparticles near the gap edge. These quasiparticles, therefore, are constructed from single-particle eigenstates of the metal which have already been strongly distorted by the impurity potential.
Equation (42) has important implications for spectroscopy on a superconductor, for one of the procedures for normalizing spectra taken at different lateral positions on a superconducting surface is to assume that the LDOS at a particular voltage much larger than ∆ o is the same. This is an attempt to correct for possible changes in the tip-surface distance upon moving the tip laterally. A small change in the tip-surface distance can have a strong effect on the tunneling current. Unfortunately this procedure will prevent an experiment from seeing changes in the LDOS due to a non-magnetic impurity, including the conductance oscillations described in Ref. 26.
We now discuss the properties of the order parameter. ∆(r) is self-consistently determined, and is shown in Fig. 3 for small r to be identical in spatial structure to the normal-state spectral weight, Since a non-magnetic potential repulsive to electrons attracts holes, and ∆(r) depends equally on electron and hole amplitudes, one might expect a non-magnetic potential to have little effect on the spatial dependence of the order parameter. However the allowable maximum spectral density of holes depends on the spectral density of the electron band where the holes reside, so if most electrons are excluded from the site, holes will be effectively excluded as well. To emphasize this point we note that the scaled anomalous spectral weight ImF (r, r; ω) is identical to the scaled A(r; ω) for the frequencies shown in Fig. 3 (and for all relevant frequencies for the selfconsistency equation Eq. (12)). Since ImF (r, r; ω) is proportional to the product of electron and hole amplitudes, and A(r; ω) is proportional to the electron amplitude squared, the spatial structure of the electron and hole spectral weights must be similar. They are, since the normal-state spectral weight A n (r; ω) is roughly frequency independent around the Fermi energy over an energy range much greater than ∆ o .
We now comment on the lack of localized states near the non-magnetic impurity for ξk F = 449 and the small binding energy for the quasiparticle for ξk F = 10. The suppression of the order parameter near the impurity may be considered to form an attractive off-diagonal potential which may bind quasiparticles. Localized states created by order-parameter suppression would be doubly degenerate, due to the two possible spin states (see Eq. (37)). Since the quasiparticle is half hole and half electron, if the electron part is attracted and the hole part is repelled, one might expect the effects of such a non-magnetic potential on the quasiparticle to cancel. However the binding energy of the localized state is an order of magnitude smaller (ω o = (1 − 2 × 10 −4 )∆ o ) than that found in Section IV.D for a suppressed order parameter via pairing suppression. This may be explained by the well-known repulsive effect (quantum reflection) of a strong attractive potential on a quantum-mechanical particle. We find also that in the case of the magnetic impurity that a large enough non-magnetic potential of either sign will suppress the binding of a quasiparticle to the impurity. We note here that the ratio of the nonmagnetic potential to the off-diagonal potential (δ∆(r)) is much larger (v 0 /N o ∆ o ∼ 10 4 ) than the ratios of the non-magnetic potentials to the magnetic potentials considered below.

B. Magnetic Impurity
Recently, we have presented calculations of the LDOS (and thus the differential conductivity in an STM experiment) in the vicinity of a magnetic impurity 37 . These calculations indicated that the spatial structure of the electron amplitude of the localized state differed strongly from the hole amplitude of the localized state. A further result was that the spectrum should recover to the homogeneous spectrum within a few atomic spacings. Similarly motivated calculations of the LDOS due to the ℓ = 0 localized state have been presented since then 38 , although these calculations did not address the continuum LDOS. The two models used in Ref. 38 were (1) a δ-function model solved using particle-hole symmetric Green's functions, but not self-consistently, and (2) self-consistent calculations for a two-dimensional tightbinding s-wave superconductor within the BdG equations. The first method can only model the normal-state properties properly for a particle-hole symmetric band structure, such as at the Van Hove singularity in a twodimensional tight-binding band structure. The second method must contend with numerical finite-size effects, which make it difficult to calculate the continuum states. A result obtained from the first method which is only true for special band structures is that the spatial structure of the electron and hole components of the quasiparticle are the same. The authors of Ref. 38 did raise the possibility of an additional non-magnetic potential as a source of electron-hole amplitude asymmetry in the spatial structure of the localized state. We found 37 and will explore below that there is, for realistic band structures, electron-hole asymmetry without a non-magnetic potential.
We will begin with a discussion of the energies of the localized states around a magnetic impurity and the spin character of those states. Solutions to Eq. (34) can be evaluated numerically. Figure 5 shows the dependence of the energies of the localized state poles for the first two angular momentum channels on the strength of the magnetic potential. Results for a short coherence length (ξ = 10k −1 F ) are shown (solid line) as well as results for a long coherence length (ξ = 449k −1 F ). The localized quasiparticle state for small v s is the spin state attracted to the classical spin, which we will label up (↑) 63 . As the potential strength increases, the excitation energy of each angular momentum state decreases. At some critical value v * sℓ the localized state becomes a spin-down excitation, the energy changes abruptly, and then increases with increasing v s . This behavior can be extracted from the analytic model (Eq. (37)) as well.
Also shown in Figure 5 are the analytic results for the pole energies for α = 0 (Ref. 8) and α = 0.704. The muffin-tin model is no better than the particle-hole symmetric model in predicting the localized state energies. The muffin-tin model will prove more successful at predicting the asymmetry between the electron and hole amplitudes of the localized state.
The unimportance of self-consistency for determining the energies of the localized states can be gauged by the small difference between the short-coherence length result and the long-coherence length result. The most important feature it determines in Fig. 5 is the size of the discontinuity in the localized state energy at v * s0 . This discontinuity is due to a discontinuous change in the order parameter at this magnetic potential strength, a result pointed out in Ref. 37 which will be discussed more below.
A magnetic impurity affects each member of a pair oppositely. An up-spin electron is attracted to the impurity and a down-spin electron is repelled. The difference in energy between the two time-reversed electron states leads to a locally suppressed pairing. Therefore, the energy needed to break a pair and create a localized quasiparticle with angular momentum ℓ in the vicinity of the magnetic impurity, when v s < v sℓ * , is reduced from 2∆ o to ∆ o + ω ℓ . One member of the broken pair is a delocalized spin-down quasiparticle at the gap edge (with energy ∆ o ). The other member of the broken pair is a localized spin-up quasiparticle with energy ω ℓ . For v s > v * sℓ there is a spin-up quasiparticle with angular momentum ℓ in the ground state of the superconductor 65 . An essential point about the ground state of a superconductor containing classical magnetic impurities is that at T = 0 when v s < v sℓ * (for all ℓ) the ground state is composed entirely of paired electrons. When v s > v sℓ * for any ℓ then the T = 0 ground state contains localized quasiparticles as well as pairs and there are new low-energy excitations including the re-formation of a pair as well as the excitation of a localized quasiparticle to a higher energy localized or continuum state. Evident from Eq. (5) is that the electron pole involves entirely single-particle states within the spin-up band for v s < v * sℓ , but involves entirely states within the spindown band for v s > v * sℓ . The hole-like pole always involves single-particle states in the opposite spin band from the electron-like pole. The source of the quasiparticle amplitude for the various poles is indicated in Fig. 5. Figures 6(a), 6(b), and 6(c) show LDOS results for v s = 0.52, 0.875, and 1.75 respectively. They show the state split off from the continuum, with a larger electronlike amplitude than hole-like amplitude (6(a)), and then lower in energy with an increased electron/hole asymmetry (6(b)). Finally the larger peak becomes hole-like (6(c)). The asymmetry between the electron-like and hole-like peaks becomes more pronounced as v s increases. We note that the larger peak is always associated with the spin-up band, whereas the other is associated with the spin-down band. Despite the apparent differences in peak size, the spatially integrated electron spectral weight of the quasiparticle is equal to the spatially inte- grated hole spectral weight. The localized quasiparticle is always half electron and half hole for all potentials examined here. For v s < v * s0 the spin-up band amplitude is electron-like and the spin-down band amplitude is holelike. At v * s0 (1.32 for niobium), the spin-up component becomes hole-like and the spin-down component becomes electron-like, as required by the change in the spin of the excitation. Figure 6(d) shows the LDOS for a markedly different coherence length, ξ = 10k −1 F , and v s = 0.875. It is evaluated for the same values of ∆ o /k B T as Figs. 6(abc) and looks almost identical to Fig. 6(b). Since the localized state is broadened by temperature through Eq. (24), this is a manifestation of the proportionality of the spectral weight to N o ∆ o (Eq. (40)). Figure 7 shows the spectral weight at the origin for the ℓ = 0 state and for the ℓ = 1 state at its first maximum for v s = 0.875 as a function of the inverse of the coherence length, which is proportional to N o ∆ o . Figure 8 shows the spectral weight for ξ = 449k −1 F as a function of v s for the spin-up and spin-down poles at the origin for the ℓ = 0 state and at the first maximum for the ℓ = 1 state. It is clear that a nonmagnetic potential is not necessary to obtain an electronhole asymmetry.
In Figure 9 we show the asymmetry at the impurity as a function of v s for two values of ξ -a long coherence length appropriate for niobium, and a short coherence length. From Figure 7 it should be evident that the asymmetry is not sensitive to ξ. It is, however, predicted extremely well by the normal-state spin-up and spin-down band spectral weight asymmetry at the impurity (also shown in Figure 9). We can therefore conclude that as with the non-magnetic impurity (Figure 4), the spatial structure of the spectral weight is a normal-state property. We further show in Figure 10 for v s = 0.875 the ℓ = 0, spin-up band and spin-down band projections of the normal-state spectral weight to compare with the localized state spin-up band and spin-down band spectral weights for two values of ξk F . The normal-state and long-coherence length calculation are practically indistinguishable. The insets show r 2 A(r), which removes the rapid power-law decay of the state. The localized states for all angular momenta ℓ will decay as the power law r −2 . For the short-coherence length calculation the  9. Ratio of the spectral weight in the spin-up band and in the spin-down band at the impurity (r = 0) as a function of magnetic potential strength. This is plotted for the normal state ℓ = 0 projected Green's functions (solid line) as well as for the localized states for niobium (ξkF = 449, long dashed line), for ξkF = 10 (dot-dashed line), and for the muffin-tin model. The muffin-tin model is only successful for vs < 0.5, but that is due to a breakdown in describing the normal state. The normal-state electronic structure is a good predictor of the superconductor's electronic structure for the entire range of vs.
effect of an exponential envelope is also visible. In the analytic result the exponential envelope should have a range R = πξ/2 1 − (ω o /∆ o ) 2 , which corrects to better than 1% the discrepancy in Figure 10. The power-law The inset shows the spectral weight multiplied by r 2 to remove the algebraic decay. The normal-state and long-coherence length results are practically indistinguishable. The deviation shown in the short-coherence length superconductor's spectral weight is fit to within 1% by the exponential decay factor described in the text. Hence the spatial structure of the spectral weight of the superconductor's localized state is well-predicted by the normal state spectral weight.
fall-off and exponential envelope can be seen directly from Eq. (15) and Eq. (40).
We can summarize these comments with a general equation, similar in concept to that for the non-magnetic impurity, Eq. (42). That is, for a localized quasiparticle state of spin σ ′ , the spectral weight of a localized state with angular momentum ℓ would be where B is a normalization factor so that the spectral weight of the state integrates to one, and A nσ (r, ℓ) is the angular-momentum ℓ projection of the normal-state spectral weight in the spin σ band. We note that for small r there is an approximate relationship between the superconducting state's spectral weight and the normal state's spectral weight in each spin band, where ∆ o ≪ E ≪ ǫ F . This, in connection with Eq. (44), implies a dependence on the normal-state structure of the continuum spectral weight around the magnetic impurity. We now return to the structure of ∆(x). This quantity, which is not directly observable, formed the focus of several investigations of the local structure around a magnetic impurity. The oscillation of the order parameter around a magnetic impurity was first evaluated without self-consistency [21][22][23] . A self-consistent calculation of the order parameter at the impurity and very far away for weak impurity potentials was done by Schlottmann 24 .
As shown in Figure 1, for large values of v s , ∆(x = 0) < 0. Sign changes in ∆, as seen in pair tunneling, have been suggested for magnetic impurities in the barriers of Josephson junctions [66][67][68] . The sign change in ∆(0) occurs (at T = 0) precisely at v * s0 . Due to the spin and frequency symmetries of Eqs. (4)(5)(6)(7)(8)(9)(10)(11), the anomalous spectral weight ImF (r, r, ω) associated with the spin-up pole is always equal and opposite to the anomalous spectral weight associated with the spin-down pole. As the pole in the spin-up band goes from electron-like (ω > 0) to hole-like (ω < 0) and the pole in the spin-down band goes from hole-like to electron-like the contribution to ∆(0) changes sign abruptly at T = 0. The ∆(r) resulting from several values of v s and two values of the coherence length are shown in Figure 11(ab). The discontinuity at v * s0 is more pronounced for shorter coherence lengths since the localized state's spectral weight is more concentrated at the impurity (Eq. (44)). ∆(0) as a function of v s is shown in Figure 12 for two values of the coherence length. At T > 0 the transition would be smoothed somewhat.
The behavior of ∆(0) as a function of v s comes from the introduction at v * s0 of a quasiparticle into the ground state of the system. The spin-up quasiparticle localized near the impurity in the ground state suppresses the local order parameter. For time-reversal invariant potentials one cannot make ∆(r) negative by inserting a single quasiparticle, since the suppression from one quasiparticle is cancelled by the lack of suppression from its unexcited Kramers doublet partner. For a spin-dependent potential, however, the anomalous spectral weight near the impurity may be almost entirely contributed by the single low-energy localized state. When a quasiparticle is present in the ground state, the ground state has spin 1 up 65,38 and a negative ∆(r) 37 . Exciting the low-energy state for v s > v * s0 removes the spin-up quasiparticle, and therefore increases ∆(0), whereas excitation of quasiparticles typically reduces ∆(x) (which is the case for v s < v * s0 ). Also, exciting the low-energy state reduces the induced spin of the superconductor at the impurity.

C. Combined Magnetic and Non-Magnetic Potentials
We now discuss the addition of a non-magnetic potential to the magnetic potential. It had been suggested 38 that introducing a v 0 with a v s would provide electronhole asymmetry. We find that it does change the asymmetry, which we show in Figure 13 for a particular v s , but that once again this is a normal-state property. The relationship between the normal-state spectral weights and the superconducting-state spectral weights of Eq. (44) still holds. Introducing v 0 also alters the localized-state energies (see Eq. (37)), which we show in Figure 14 for v s = 0.875 and ξk F = 449. The presence of a nonmagnetic potential may affect the value of v * s0 38 . We show in Figure 15 a partial diagram of the ground state as a function of the parameters v s and v 0 for ξk F = 449. We note that the boundary between the two ground states is not shifted much from ξk F = 449 to ξk F = 10, hence the condensation energy is not very significant in determining this boundary. This fact, and the observation that the energy of the localized state (in the absence of selfconsistency) relative to the homogeneous order parameter does not depend on any parameter of the superconducting state (Eq. (37)) suggests that an understanding of the transition at v * s0 may be assisted by normal-state concepts. We hope to address this issue in a future In a similar result to that seen in Fig. 9, the normal-state spectral weights are good predictors of the superconducting state's spectral weight. We note that the curves are not symmetrical around v0 = 0, which results from a realistic band structure without particle-hole symmetry.

publication.
We now mention an example designed to showcase the drawbacks of Ginzburg-Landau theory. We present in Figure 16 the order parameter for a mixed magnetic and non-magnetic impurity (v s = 7/8, v 0 = −7/4, and ξk F = 449). The order parameter is everywhere larger than in the homogeneous superconductor, however, the presence of a localized state within the gap indicates that superconductivity has been weakened around the impurity. Since Ginzburg-Landau theory focuses on the order parameter and ignores the quasiparticle structure, a Ginzburg-Landau perspective would incorrectly predict an enhancement of superconductivity in the region.
We conclude by attempting to provide some guidance for attempts to extract impurity potentials from STM measurements. Figure 17 shows normal-state dI/dV 's for various potentials. These curves should also represent    16. The change in order parameter around an attractive non-magnetic potential combined with a magnetic potential for ξkF = 449. The order parameter is larger at the impurity than in the homogeneous superconductor, yet there exists a localized state in the gap.
the frequency-averaged spectral weight measured in the superconducting state (see Eq. (45)) 69 .
The enhancement or suppression of spectral weight near the origin is particularly sensitive to v 0 . A measurement of this quantity, the energy of the localized state and the asymmetry of the electron and hole amplitudes at the impurity overconstrains v 0 and v s , given an assumption of the shape of the potential. To extract information about the potential's detailed shape would require a fitting procedure using the differential conductivity at various positions. If, for some reason, the spin-down amplitude were too small to measure, it may remain possible to constrain the potential strength using the frequencyintegrated spectral weight and the localized-state energy.
It seems appropriate to mention again the tendency to normalize spectra according to the LDOS measured at energies much larger than ∆ o . Since the normalstate LDOS near the impurity changes substantially in the presence of magnetic or non-magnetic impurities, an experiment performed using such a normalization procedure would yield impurity parameters of questionable validity.

D. Pairing Suppression
The pairing potential, γ(x) in Eq. (12), may also have spatial structure. When this parameter is changed it induces a change in the order parameter which produces an off-diagonal potential felt by the quasiparticles. We set v s = v 0 = 0 so that there is no magnetic or nonmagnetic potential to compete with the order parameter change. Figure 18 shows the order parameter around a short-range suppression, for two values of the coherence length. The order parameter is strongly suppressed and since γ(0) = 0, ∆(0) = 0. For long coherence lengths, this change in the order parameter has no effect on the local density of states (shown in Figure 19(a)). It is possible to localize quasiparticle states, however, at shorter coherence lengths. These FIG. 17. Spectral weights in the normal state as functions of the distance from the impurity for several combinations of magnetic and non-magnetic potentials (a) only a magnetic potential, and (b) a fixed magnetic potential vs = 0.875 and a varying non-magnetic potential. By making measurements around the impurity in the normal state (or integrating the superconductor's spectrum over a frequency much larger than ∆o), information about the structure of the impurity may be obtained.
can produce features in the local density of states which are visible. One such case is shown in Figure 19  density of states without significantly changing the energy gap, a pairing suppression has a very weak effect on both, especially in the long-coherence length limit.

V. STRONG-COUPLING AND ANISOTROPIC ORDER PARAMETERS
A few observations are in order concerning the extension of this formalism to systems where the homogeneous order parameter has important frequency or spatial structure. The Gor'kov equation (Eq. (4)) changes due to the more general form for the off-diagonal potential originating from the order parameter. Taking this opportunity to generalize, where The diagonal terms are general potentials, possibly frequency-dependent, effective on spin-up electrons (V e↑ ) and spin-up holes (V h↑ ) for ω > 0. Since this potential is diagonal in frequency, as is the Gor'kov equation, the frequency structure of the order parameter does not add any additional complication to numerically solving the Gor'kov equation. However, the addition of spatial structure to the pairing has added another integral over the volume to the Gor'kov equation, and thus dramatically increased the size of the matrix which needs to be inverted. Fortunately, the range of the order parameter in |x − x ′ | is truncated by the range of the pairing interaction. For the isotropic-gap superconductor we have considered for most of this paper, the effective pairing interaction is modeled by a delta function in space. It is possible to obtain anisotropic order parameters, including the d-wave order parameters possibly appropriate for high-temperature superconductors, merely by considering pairing with nearest neighbors on a square tightbinding lattice. This model has been implemented with the BdG formalism for non-magnetic impurities 32-35 and magnetic impurities 38 and vortices 70 in d-wave superconductors. When formulated on a lattice the addition of nearest-neighbor pairing multiplies the rank of the matrix M n→n (ω) by 1 + z, where z is the coordination number of the lattice.
The order parameter's frequency dependence complicates the self-consistency equation (Eq. (12)). It must now be solved for each frequency: ∆(r, r ′ ; ω)Z(r, r ′ ; ω) where K ± are kernals of the pairing interaction and are different for each mechanism of superconductivity. They can be determined from the homogeneous solution. U c is a Coulomb factor. Z, the quasiparticle weight, is solved for self-consistently. Incorporating these strong coupling effects allows a determination of the effect of the frequency-dependence of the pairing interaction on the electronic structure around a defect. While selecting the particular location where the STM tip is similar to selecting the momentum of the quasiparticles of interest, selecting the STM voltage indicates which order parameter frequency one wishes to probe.

VI. SUMMARY
The local electronic structure around an defect reflects both the properties of the defect and the medium it is embedded in. For a strong non-magnetic or magnetic impurity in a superconductor, the distortion of the normalstate properties by the strong impurity plays a vital role in the response of the superconducting medium. The LDOS for the inhomogeneous superconductor can be related to the LDOS for the inhomogeneous normal metal via equations like Eq. (42) and Eq. (44). In the case of a non-magnetic impurity with no localized states around it, the LDOS for the inhomogeneous superconductor is merely the normalized LDOS for the inhomogeneous normal metal multiplied by the homogeneous superconductor's density of states. This should suggest some caution for the method used to normalize STM spectra taken at different places on a superconductor's surface.
For the case of a localized state (such as around a magnetic impurity) with angular momentum ℓ, the LDOS for the state is the angular-momentum-ℓ-projected LDOS of the inhomogeneous normal metal multiplied by a decaying exponential whose range is determined by the energy of the localized state. Since the spin-up band LDOS in the normal state differs from the spin-down band LDOS in the normal state, the electron-like pole of the localized quasiparticle will have different spatial structure than the hole-like pole of the localized quasiparticle.
The self-consistent calculations described here have been performed with a new, powerful, Koster-Slater technique which allows the Gor'kov equation to be solved in principle exactly. Although we have only presented calculations for weak-coupling isotropic order parameters within a free-electron model, we have formulated the extension of the new technique to strong-coupling pairing and general band structures and order parameter symmetries.
where J ℓ , H (1) ℓ , and H (2) ℓ are standard Bessel functions, r < (r > ) is the smaller (larger) of r and r ′ , ω is in units of ∆ o and r and r ′ are in units of k −1 F . The Green's functions are in units of N o .
Thus the three-dimensional integral has been reduced to a one-dimensional radial integral. Since the numerical inversion procedure depends on inverting a matrix whose rank is proportional to the number of spatial points considered, this reduction to a one-dimensional integral dramatically increases the speed of this calculation over a calculation for a three-dimensional potential which is not spherically symmetric.