The single-site Green’s function and Krein’s theorem

An important step in electronic structure calculations using multiple-scattering theory is obtaining the density of states for the central site from the Green’s function for that site. We have found that the Krein’s spectral displacement function for the central site contributes significantly to the understanding of these calculations. We argue that these insights can lead to improvements in the robustness of MST electronic structure codes without negatively impacting their performance.


Introduction
During the academic year of 1968-69, Balazs Gyorffy held a research position at the University of Sheffield. Malcolm Stocks was finishing his PhD dissertation, and Sam Faulkner had taken a year away from his position as head of a Theory Group at the Oak Ridge National Laboratory (ORNL) to visit Sheffield as a Senior Fulbright Scholar. The three developed a mutual scientific rapport, and Malcolm and Balazs came to ORNL the following summer. This was the first of many visits that Balazs made to ORNL, and Malcolm stayed there. Sam had been introduced to multiple-scattering theory (MST) and alloy theory by his Professor, Jan Korringa. Balazs and Malcolm became interested in MST, and quickly produced important papers [1,2] in the rapidly burgeoning field of calculating the electronic states of disordered solids using the coherent potential approximation (CPA). The three scholars who began collaborating in 1969 have had many graduate students, post-doctoral students, and collaborators who have helped develop the MST, CPA, and related fields to their Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. present state. An author of the present paper, Yang Wang, is among that group. It is with great sadness that the authors note the passing of our friend and colleague, Balazs Gyorffy.
The proposal to use multiple-scattering theory (MST) to study the quantum states of electrons moving in a local one-electron potential was put forward by Korringa [3] as a technique for doing band theory calculations on periodic solids. Kohn and Rostoker published another derivation of the technique [4], and it is in use as the Korringa-Kohn-Rostoker (KKR) method [5] at many institutions in Europe, Japan, and the USA. MST also underpins the locally self-consistent multiple-scattering (LSMS) method that has been developed by some of the authors of this work [6]. The LSMS method is a real space MST code that can perform ab initio studies of extended defects and complex magnetic structures in O [N] time for models comprising a number of atoms (N) in the range 10 4 -10 5 using current leadership-class massively parallel computers (of [7]).
The local potentials V (r) used in MST are obtained from density functional theory [8]. They are written as the sum of non-overlapping potentials associated with each atom (1.1) Each of these potentials is considered a scattering center for the electron. The MST equations couple the incoming and outgoing waves from each center to form a continuous eigenfunction. The division of V (r) into non-overlapping v i (r) is achieved by constructing planes that bisect the lines connecting the nuclei to obtain space-filling Voronoi polyhedra i . The potential v i (r) is defined to be zero for r / ∈ i . Originally the v i (r) were approximated by muffin-tin (MT) potentials, which are spherically symmetric within spheres inscribed within the i , and zero outside of them. MST works quite well for v i (r) that are not spherically symmetric and whose domain of definition fills the interior of i . Such calculations are typically referred to as being full-potential (FP).
The following section contains a brief review of MST and a demonstration of the key role played by the singlesite Green's function G n (E, r, r ) in MST calculations. In section 3, G n (E, r, r ) is analyzed with the help of one of the most powerful tools of formal scattering theory, Krein's theorem. Krein's spectral shift function ξ(E) is used to obtain the single-site density of states and integrated density of states associated with a potential v n (r) that can be compared with those predicted from G n (E, r, r ). The other contributions from the full MST Green's function G(E, r, r ) are treated adequately in other publications and we do not reiterate that discussion. In section 4, the algebra in section 3 is illustrated by calculations using G n (E, r, r ) and ξ(E) with simple model single-site potentials. In section 5, similar calculations are done using non-spherical single-site potentials v n (r) obtained from self-consistent full-potential LSMS calculations on crystalline aluminum, copper, and molybdenum.

Green's function
The MST Green's function was derived some time ago [9], and is applicable to potentials that have no shape approximation. When r and r are in the nth cell and r > r , it is written where L stands for the pair of angular-momentum indices , m. The Z n L (E, r) are obtained by integrating outward from the origin at the nucleus of atom n, and are regular there. The J n L (E, r) satisfy a boundary condition away from the nucleus and are obtained by integrating inward. There is no reason for them to be regular at the origin, and they typically approach infinity there 4 . The scattering-path matrices τ nn L L (E) describe the propagation of the electron around all possible paths in the system, which is infinite in the case of standard KKR band theory and confined to a finite cluster around each atom in the system in the case of the LSMS method. The dot notation in the above expression indicates that the complex conjugates of any spherical harmonics are taken, but the remainder of the function is left unchanged.
Quantities that are used to describe the electronic structure of condensed matter can be obtained from the Green's function. The ones that of interest in the present context are the density of states n(E) = − 2 π Im G(E, r, r) dr, (2.2) and the integrated density of states where E B is the lowest eigen-energy in the conduction band. The straightforward application of equation (2.1) to calculate these parameters requires an amount of computational time that is almost prohibitive because G(E, r, r ) varies rapidly with energy. The resolution of this problem follows from the observation that, when viewed as a function of complex energy, G(z, r, r ), is holomorphic everywhere except for poles at the bound states and a cut on the part of the real axis starting at E B and extending to plus infinity. Because G(z, r, r ) is a much smoother function of z than of E, quantities that require an integral over energy are evaluated by integrating over a contour in the upper half-plane that is restricted only by the conditions that it starts at E B and ends at the Fermi energy E F [10]. The Green's function is usually broken up in two parts where G n (E, r, r ) is the single-site Green's function corresponding to the nth site. It is useful to define the t-matrix in terms of sine and cosine matrices where α = √ E. It has been shown [11] that c n and s n can best be calculated with the help of an ancillary function φ n L (E, r) that is the solution of and satisfies the boundary condition lim r →0 φ n L (E, r) = Y L (r) j l (αr ). (2.8) In this equation j l (αr ) is the spherical Bessel function that is regular at the origin and Y L (r) is the spherical harmonic with the angles (ϑ, ϕ) defined by the direction of the vector r. The subscript on φ n L (E, r) refers only to the boundary condition equation (2.8). It can be expanded in spherical harmonics in the usual way It is obvious that for r > R n c , where R n c is the radius of the sphere that circumscribes n , φ n L L (E, r ) can be written as a linear combination of the j l (αr ) and the spherical Bessel functions of the second kind, n l (αr ), (2.10) By comparing this function with a solution of equation (2.7) in the scattered wave form, it can be shown that the coefficients are the elements of the sine and cosine matrices. It is by no means obvious, but it has been shown by exhaustive computational and algebraic studies [12] that equation (2.10) holds true for r < R n c as long as r is on the surface or outside of the polyhedron n . Formulas for the sine and cosine matrices are The functions Z n L (E, r) that appear in equation (2.1) are linear combinations of the φ n L (E, r) so they are also regular at the origin. The J n L (E, r), however, are solutions of that take on the asymptotic values when r > R n c . These functions are typically infinite at r = 0. The second term in equation (2.4) describes the scattering from every site except the central one. In many ways this is the easiest part to calculate because contour integrals may always be used in dealing with it. Since the techniques for doing this have been discussed extensively, we will not consider it further. The major problem is caused by the second term in G n (E, r, r ). The inward integration required in the evaluation of J n L (E, r) is very difficult to converge near r = 0, particularly for complex E. A method for dealing with this problem that is used extensively is to stop calculating the term within a sphere surrounding the origin and to extrapolate by some means within that sphere [13]. It has been shown [14] that the correct charge density for the conduction electrons has undulations that reflect the orthogonality of the wavefunctions to the core functions and are appreciable at radii as small as the radius of the 1s electrons. Misrepresenting these undulations can lead to significant uncontrolled errors in the charge density, total energy, Hellmann-Feynman forces, and other physical quantities.
From equations (2.6) and (2.13), the central-site Green's function can be written One of the authors, Wang, derived a useful lemma that the sine and cosine matrices satisfỹ where the tilde indicates the transpose of the matrix. In a previous paper [14], Wang's lemma was used to separate the first term into real and imaginary parts, where =c n· c n +s n· s n .

(2.19)
As in all MST calculations that use angular-momentum expansions, the formulas in this section are made tractable by ignoring terms with l values greater than some l max .
It was pointed out [14] that this function could be integrated along the real axis to get, for example, (2.20) because the poles on the real axis that occur in the last term are exactly canceled by those in the second term. This approach was not pursued because it appeared that contour integrations based on the equivalence of the poles in the last term with those in the second would be a faster method to eliminate the problems caused by the function J n L (E, r). The method works and the calculations are fast, but a multiplicity of poles appear as l max was increased. Many of those poles are fictitious, but they slow the calculation to the point that it now appears that it is better integrate along the real axis.

Krein's theorem
A very general theorem that has found many uses in scattering theory, statistical mechanics, and mathematics was derived by Krein [15]. He proved that, for any two self-adjoint operators H and H 0 whose difference V = H − H 0 is trace class, there exists a real function of real variables ξ(E) such that where φ(X ) is any sufficiently smooth function. For applications to scattering theory, φ(H ) = (z − H ) −1 and therefore where H and E are a Hermitian operator and a real scalar and z is a complex number. It follows that For this case, Birman and Krein [16] proved that Krein's spectral shift function ξ(E) is given by is demonstrated in the papers quoted, and in the massive literature on the subject that followed. In simple terms, the definition includes the conditions that ξ(E) is a continuous function of E and approaches zero as E approaches plus or minus infinity. The spectral shift function defined in this way is consistent with the known properties of the S-matrix. From this definition and Cauchy's integral theorem, it follows that Recall that the potentials used in this paper are such that v n (r) is defined to be zero for r outside the finite domain n . It follows that the energy spectrum for a system with one such potential embedded in a vacuum is discrete for E < 0, being the set of bound state energies E ν . The system has a continuous spectrum for E ≥ 0.
The trace operation indicated on the left side of equation (3.2) can be carried out in many ways. Taking it to be the integral over all space, leads to and hence where d ν is the degeneracy of the bound state, n n (E) is the density of states, and n 0 (E) is the free-electron density of states. From the discussion in the preceding paragraph, n n (E) and n 0 (E) are zero for E < 0. The 2 is included to account for the spin of the electrons. Integrating equation (3.6) from minus infinity to E leads to where N K (E) is designated the Krein integrated density of states (IDOS).
Integrating equation (3.6) from minus infinity to zero leads to a generalization of Levinson's theorem 2 ν d ν = n c = −2ξ(0) = N K (0), (3.8) where n c is the total number of core electrons. This equation, which holds for potentials that are not spherically symmetric, can be checked for consistency with the ordinary Levinson's theorem by applying it to potentials that are. The S-matrix for such potentials is diagonal with elements repeated 2l + 1 times. Considering only the states that correspond to a given l leads to, the standard version of the theorem. It is obvious that the advantage of Krein's formulation is its generality. The trace operation in equations (3.1) and (3.2) is independent of the orthonormal basis used, as long as it is complete. The versions familiar in quantum theory are spatial integrations, as in equation (2.1), or integrations over momenta. The requirement that the potential is trace class is essentially equivalent to the condition that the trace on the left side of equation (3.2) exists, with no shape approximation.
As with the Green's function formulas in the preceding section, the ones based on Krein's theorem are made tractable by ignoring contributions from terms with l > l max .
Krein's theorem, as stated in equation (3.2), holds for any energy. Calculations in the range −∞ < E < 0 are very difficult for the specific kind of potentials being treated in this paper. For this reason, a numerical test of the generalized Levinson theorem, equation (3.8), for these potentials is carried out by calculating ξ(E) for 0 ≤ E < E max , where E max is sufficiently large that ξ(E) has become a constant. Then ξ(E max ) is set equal to zero, and integers are added to ξ(E) over various energy intervals is such a way as to make it a continuous function. Following the function in this way from E max to zero is a way of counting the bound states of the potential. The result can be checked by finding the bound states with the numerical techniques that have historically been used in atomic physics. This calculation is done for a simple model potential in the next section.

Calculations on a model potential
The best way to illustrate the algebra in the preceding sections is to do some calculations on a spherical square well potential. The potential is v(r) = − for |r| < R and v(r) = 0 otherwise. The advantage of this model is that all integrals that appear in the calculations can be done analytically, and any hypothesis can be checked numerically in a matter of minutes. The radius R of the wells we use is chosen to be 3 because that is approximately the radius of atomic spheres in a metal in units of bohr radii. Because of this, the dimensionless energies in the model calculations have the same scale as the Rydbergs  is called a partial density of states because the integral is only over the sphere. This is equivalent to a partial trace. It follows that is a partial IDOS. The free particle Green's function is removed in equation (4.2) so that the result can be compared with the Krein IDOS. The potential with = 1.4 models a system with 8 s and p states, but is not close to binding d states. The potential with where the volume ∞ − is all of the space outside the domain in which the potential is nonzero. The corresponding partial IDOS is then  of Z L (E, r) in [11] is used,  In these equations, h l (αr ) is the Hankel function that describes outgoing waves. The contribution to the DOS from the outside integral is therefore For a spherical potential, this becomes This advantage in writing n out (E) in this form is that the definite integral in equation (4.9) has been worked out and published [17] ∞ R h l (αr ) 2 (4.11) The partial IDOS N out (E) obtained from these equations for the three values of are shown in figure 3. A visual  comparison of figures 2 and 3 confirms that as it should be. A numerical differentiation of N K (E) gives the Krein density of states This is plotted in figure 4  The last three graphs are truncated because the Krein DOS goes to plus or minus infinity at E equals zero. The reason for this can be seen from equation (4.1). It is well known that the (4.14) The Krein IDOS will be dominated by the l = 0 phase shift, which has the limit This behavior is illustrated in figure 7. It follows that the Krein DOS approaches infinity like The partial DOS n in (E) approaches zero smoothly, so it follows that n out (E) = n K (E) for small energies. The equality in equation (4.12) and the equality n K (E) = n in (E) + n out (E),  The two orders of magnitude difference in the errors for l max equal to 16 and 8 may be significant in some applications. Finally, the numerical methods used in the present calculations are fairly primitive. It is assumed that they will be improved over time so that, for example, the number of energy points can be reduced and the results will be less dependent on the choice of E min .

Calculations with real potentials
At this point, the way that Krein's spectral displacement function compares with densities of states and integrated densities of states obtained from Green's functions is understood. The crucial question is if the theory is applicable to full-potential multiple-scattering calculations, and if anything of value for such applications has been learned. The single-site potentials used in this section were extracted from self-consistent full-potential LSMS calculations on periodic crystals of aluminum, copper, and molybdenum. The lattice constants for the face-centered cubic crystals of Al and Cu were 7.65 and 6.76 bohr, while for body-centered cubic Mo it was 5.913 bohr. To obtain potentials suitable for the present purposes we used minimal settings of the usual l-truncation parameters, namely l max kkr = 2, l max phi = 2, l max rho = 4, l max pot = 4, and l max truc = 12, where l max kkr , l max phi , l max rho , l max pot , and l max truc are the angular-momentum expansion cut-offs for the KKR-matrix, wavefunction, charge density, single-site potential and σ (r)-function respectively (the latter is used when performing spatial integrations over a Voronoi polyhedron). The von Barth-Hedin exchangecorrelation functional was used in all of the calculations. Because Al, Cu, and Mo are all cubic crystals the output potentials contain both a spherical, l = 0, and non-spherical, l = 4, components. Of course, a fully converged full-potential calculation requires much larger values for these parameters. For studying the single-site Green's function and Krein's theorem we use l max kkr = 8 and l max truc = 12 together with the above l max pot = 4 single-site potential. When solving the single-site equations we further apply the σ (r)-function to the untruncated potential, terminating this expansion at l = 8. As a consequence, the actual potential used in solving the single-site equation contains l = 0, 4, 6 and 8 components. For performing the energy integral we tested a variety of linear and logarithmic grids with different lower limits for the integration within the range 10 −8 to 10 −3 Ryd (the upper limit being set at 1.5 Ryd) to test overall convergence of the various integrals. To test the small-energy asymptotic behavior of various partial quantities and the generalized phase we also used dense (∼1000-point) linear and logarithmic meshes with lower/upper limits of 10 −8 /10 −3 Ryd; the results of these exercises were consistent with what was already expected from previous studies of the square well (above) and muffin-tin potentials. The specific full-potential results shown below are for a logarithmic mesh with 1000-steps and a lower limit of 10 −3 Ryd, which are sufficient for the present purposes.
In FP calculations it is obviously necessary to use the general formula for the S-matrix when calculating the contribution to n(E) from inside the Voronoi polyhedron . The two formulas are mathematically equivalent, but equation (5.2) gets around the problem that thẽ c · c term in is badly behaved for large l max . The portion of n(E) obtained by integrating the Green's function over the space outside the Voronoi polyhedron, ∞ − , is called n out (E). To calculate this quantity, we consider two concentric spheres. The largest sphere that can be inscribed within has a radius R mt , while the smallest sphere that circumscribes has a radius R c . We introduce a step function σ (r) that is one inside the Voronoi polyhedron and zero outside. The technique for expanding this function has been well established [18]. Using these quantities we obtain The definite integral in equation (4.11) is used to evaluate the first term. Using these formulas, the various integrated densities of states for the Al, Cu, and Mo potentials were calculated at 1000 energy points on a logarithmic scale between E min = 10 −3 Ryd and E max = 1.5 Ryd. These data were used to generate figures 8-10. Careful calculations based on the data shown in the last three figures demonstrate that the sum in equation (4.12) holds as well for FP calculations on real solids as it does for model calculations. As explained above, the Mo IDOS is analogous to the IDOS for = 1.4, while those for Cu and Al are similar to the ones for = 2.1 and = 2.4.
The densities of states for the three elements are shown in figures 11-13. These data show that the Krein DOS is the sum of those calculated from the Green's function, as expected. In addition, n K (E) and n out (E) approach infinity as E approaches zero. There is not a simple proof of this like the one that led to equation (4.16) for spherical potentials, but an argument first introduced by Gyorffy [19] is helpful in this connection. Because the S-matrix is normal and unitary, it has eigenvalues λ n (E) and they fall on a unit circle, i.e. λ n (E) = e i2ϑ n (E) (5.5) with ϑ n (E) real numbers. Gyorffy coined the term generalized phase shifts for these functions. The generalized phase shifts  have been calculated for full-potential aluminum, copper, and molybdenum. At small energies they are dominated by an s-like phase shift that is proportional to the square root of E, just as in equation (4.15). The other phase shifts depend on energy like the ones in equation (4.14), except there are more than one for each l. This is interesting in itself, and it provides the explanation for the form of the Krein spectral displacement function of full potentials.

Conclusions
In MST calculations, the primary interest is in the density of states n in (E), the integrated density of states N in (E), and the charge density calculated from the Green's function inside the Voronoi polyhedron. These are difficult and time consuming calculations, especially when the density of states has sharp structure like copper. It is much easier to calculate Krein's spectral displacement function and the density of states from the Green's function outside the Voronoi polyhedron n out (E). At a minimum, these calculations can be used to check on the n in (E) and N in (E). The algebra and calculations in this paper demonstrate the meaning of formulations that are in common use in multiple-scattering theory. The original derivations of those formulations were not mathematically rigorous. It should be clear from the discussion that Krein's theorem is useful for non-spherical potentials because it makes it possible to generalize theorems that have only been proved for spherically symmetric potentials.
In a previous paper [20], one of the author's showed how Lloyd's formula for a cluster of scatterers can be derived from Krein's theorem. The usefulness of Lloyd's formula rests on the assumption that it gives the integrated density of states for one cell in the cluster. This might seem strange at first because the calculations in this paper indicate that Krein's SDF is not the same as the partial IDOS N in (E) obtained from the Green's function. This is not a problem because in [20] Krein's theorem is being applied to the entire cluster of scatterers, not just one. Model calculations show that the larger the scatterer, the closer the agreement is between Krein's IDOS and N in (E).