17 March 2000
H. Akai

Korringa–Kohn–Rostoker Method
1. Brief summary of the theory of scattering
1.1. Radial Schrödinger equation
The Schrödinger equation describing the motion of the particle scattered by a central,spherical potential in units of h̄/2m = 1 is−∇2 + V (r) ψ(r ) = Eψ(r ),where the Laplacian operator ∇2 in spherical coordinate is


L2
1 ∂
2
2 ∂
r
− 2,
∇ = 2
r ∂r
∂r
r
with
1 ∂
L =−
sin θ ∂θ



2

∂
sin θ
∂θ


−

1 ∂2
.
sin2 θ ∂ϕ2

(1.1)

(1.2)

(1.3)

A solution of eq.(1.1) is then expressed in terms of a solution of the radial Schrödinger
equation ul (r) as
ψlm (r ) =ul(r)Ylm (θ, ϕ),
r

(1.4)

where ul (r) satisﬁes
l(l + 1)d2 ul (r)+ E − V (r) −ul (r) = 0
dr2
r2

(1.5)

and Ylm (θ, ϕ) are the spherical harmonics deﬁned to be eigenfunctions of L2 and Lz .

1.2. Boundary condition and the phase shifts
Outside the potential range (r > r0 ), the radial Schrödinger equation reduced to that of
the free space and its solution may be any linear combinations of the two independent
solutions (one regular, the other irregular) of
l(l + 1)d2 wl (r)wl (r) = 0.
+ E−
dr2
r2
1

(1.6)

The two solutions P and Q of eq.(1.6) are related to the spherical Bessel (regular solution)
and Neumann (irregular solution) functions, jl(z) and nl (z), respectively, by√√
Erjl ( Er)

(1.7)

√
√
Ql (r; E) = − Ernl ( Er).

(1.8)

Pl (r; E) =and

The negative sign is attached to the deﬁnition of Ql for the later convenience. At large
distances, these functions behave as

and

√
lπ
Pl (r; E) −→ sin( Er − )
2

(1.9)

√
lπ
Ql (r; E) −→ cos( Er − ).
2

(1.10)

The asymptotic form of the solution of eq.(1.4) at large distances is therefore
√
√
lπ
lπ
rr
ul (r) −→0 A cos( Er − ) + B sin( Er − )
2
2

√
lπ
+ ηl ).
= A2 + B 2 sin( Er −2
Here the so-called phase shift ηl is introduced.

(1.11)The regularity of ul at the origin

determines the ratio B/A and the asymptotic boundary condition at r  r0 determines√
the normalization A2 + B 2 . In the theory of scattering, the proper asymptotic boundary
condition for the wave function is
rr0

√
i Ez

ψ(r) −→ e√ei Erf(θ, ϕ).
+r(1.12)

Correspondingly, the condition on the normalization of the radial wave function is√lπrrul (r) −→0 eiηl sin( Er −+ ηl )
2(1.13)

or equivalentlyul (r) = eiηl [cos ηl Pl (r; E) + sin ηl Ql (r; E)] ,2

r > r0 .

(1.14)

1.3. t-matrix
A more familiar expression of ul (r) with the same normalization as implied by eq.(1.14)
may beul (r) = Pl (r; E) −√Etl (E)R+l (r; E),r > r0 ,(1.15)(1)

where R+l = Ql + iPl is related to the spherical outgoing Hankel function hl (z) =jl (z) + inl (z) by

√
√
(1)
R+
(r;
E)
=
i
Erh
(
Er)
l
l

(1.16)and tl is the so-called atomic t-matrix deﬁned by1(1.17)
tl (E) = − √ sin ηl eiηl .
E
It is also possible to express ul (r) in terms of the k-matrix (corresponding to eq.(1.14)) or−s-matrix (if R+l and Rl = Ql − iPl are used instead of Pl and Ql ).

Practically, the phase shift ηl is calculated in the following way. First, given a potential
V (r) and an energy parameter E, we integrate the radial Schrödinger equation numericallyfrom the origin with the regular boundary condition to obtaine a solution ul (r), which
satisﬁesP  (r) + tan ηl Ql (r)ul(r)= l.
(1.18)ul(r)Pl (r) + tan ηl Ql (r)at r > r0 irrespective of the normalization, Now, equation (1.18) is readily solved withrespect to ηl .
1.4. Integral expression for the t-matrix
Another expression of the phase shifts, or equivalently the t-matrices, is obtained fromthe integral form of the Schrödinger equation. The expression is never used for actual
calculations but is often useful in formal discussions. First we deﬁne a normalization
constant Λl by(2l + 1)!! ul (r)√.
(1.19)
r→0( Er)l+1
Then from the Wronskian relation between ul and wl , which is a direct consequence of theΛl = lim

Schrödinger eqs.(1.5) and (1.6),

r2
r2




W [wl (r), ul (r)] = {wl(r)ul (r) − wl (r)ul (r)} 
r1

r1

r2

=

dr wl (r)V (r)ul (r),r13

(1.20)

where wl is either Pl or Ql , it follows that∞tl (E) =E −1 0 dr Pl (r; E)V (r)ul (r)∞Λl + E −1/2 0 dr R+l (r; E)V (r)ul (r).

(1.21)

If ul is normalized as eq.(1.13), the denominator of eq.(1.21) simply gives a unity. Similarly,
k-matrix (or s-matrix) is obtained if R+
l (or Pl ) is suitably replaced by a linear combination−of Pl and Ql (or R+l and Rl ).

2. Radial Green’s functions
2.1. Integral form of the Schrödinger equation
From the Wronskian relation eq.(1.20), it follows thatr

W [Pl (r), ul (r)] =dr Pl (r)V (r)ul (r)

(2.1)

0

and

∞

C − W [Ql (r), ul (r)] =r

dr Ql (r)V (r)ul (r),(2.2)

where the constant C depends on the normalization of ul . For the normalization impliedby eq.(1.13), C satisﬁes
C = W [Ql (r), ul (r)]r>R = W [Ql , Pl ] eiηl cos ηl√= E eiηl cos ηl ,where the relation W [j(z), n(z)] = 1/z 2 is used.

(2.3)

In the following, we adopt thisnormalization for ul . Now we multiply eq.(2.1) by Ql , multiply eq.(2.2) by Pl and add
them, obtaining
1ul (r) = eiηl cos ηl Pl (r) − √E∞dr Ql (r> )Pl (r< )V (r )ul (r ),

(2.4)

0

where r< and r> denote the smaller and the greater of r and r . Equation (2.4) is reduced
to the expression given by eq.(1.14) when r > r0 is satisﬁed.

2.2. Green’s function
The quantityg(r, r ) = −Pl (r< )Ql (r> )
√
E
4

(2.5)

appearing in eq.(2.4) is a radial Green function.
The expression which is reduced to eq.(1.15) for r > r0 is also obtained by replacing
Ql with R+l − iPl :1ul(r) = Pl (r) − √
E

∞
0



dr R+l (r> )Pl (r< )V (r )ul (r )(2.6)

Here the relation eq.(1.21), with the denominator of r.h.s. being a unity, is used. The
quantityg + (r, r ) = −(r> )Pl (r< )R+
√ l
E

(2.7)

in eq.(2.6) again is a radial Green function.

2.3. Boundary condition
The diﬀerence between two Green functions, g(r, r ) (standing wave Green function) andg + (r, r ) (outgoing wave Green function), becomes clear if we consider the integral equationof the form∞ul (r) = Pl (r) +dr g + (r, r )V (r )ul (r ).

(2.8)

0

The solution will satisfy eq.(1.14). On the other hand, if we use g, instead of g + , ineq.(2.8), the solution will satisfy ul (r) = Pl (r) + tan ηl Ql (r) for r > r0 . This means thatthe diﬀerence between g and g + is only in the normalization of the solution or, in otherwords, in the asymptotic boundary conditions. Note that the analytic properties of g and
g + nevertheless are quite diﬀerent. For example, g + is analytic in the upper half complex
plane of the physical energy-sheet while g is not.

2.4. Diﬀerential equation for the Green function
We introduce a Green function which satisﬁes an inhomogeneous diﬀerential equation asis the case of usual Green functions for the linear diﬀerential equations;l(l + 1)d2 gl (r, r )gl (r, r ) = δ(r − r ).
+ E−
dr2
r2

(2.9)

The discontinuity arising at r = r in the ﬁrst derivative Pl (r< )Ql (r> ) yields the δ-function
on the r.h.s. of eq.(2.9).
5

2.5. Eigenfunction expansion
Seemingly a quite diﬀerent expression of the radial Green function is obtained byintroducing an eigenfunction expansion of the Green function. A set of eigenfunctions{Pl (r; En )} are deﬁned by the radial Schrödinger equation eq.(1.6) with a boundary
condition wl (R) = 0 at r = R > r0 . In the limit R  1, the energy eigenvalues En ’sof Pl (r; En ) satisfy (see eq.(1.9))


lπ
.
2

(2.10)Gn (r )Pl (r; En )(2.11)Pl∗ (r ; En )Pl (r; En ),E − En(2.12)

En R = nπ +Expanding g(r, r ) asg(r, r ; E) =n

and inserting it into eq.(2.9), we obtaing(r, r ; E) =n

where we assume that Pl is normalized as
R|Pl (r; En )|2 dr = 1.

(2.13)

02.6. Radial Green function in general cases
Though all the arguments so far have regarded the free space as the reference system, thewhole story goes as well for any other choice of the reference system. Thus if we take thesystem described by the Schrödinger equation with a potential term, eq.(1.5), as a reference
system, we will obtain the Green function describing the propagation of the electron waves
in the system with a central, spherical potential. The expression for such Green functions
clearly has the same form as that for the free space Green function. In the later sections wewill refer to such general cases by the same formulae, simply replacing Pl (r; E), Ql (r; E)
and R±l (r; E) by the corresponding radial wave functions obtained in each case.

63. Free space Green’s function
3.1. Real space representation
The free space Green function is deﬁned by the diﬀerential equation(∇2 + E)g(r , r  ; E) = δ(r − r  ),(3.1)

whose momentum space representation(−k 2 + E)g(k ; E) = 1(3.2)

readily gives the solutiong(r , r  ; E) =d3 k exp [ik · (r − r  )].
(2π)3E − k2(3.3)

Equation (3.3) is nothing but the eigenfunction expansion of the Green function. The
integration on the r.h.s is easily performed as∞π2π1exp(ik|r − r  | cos θ)
2
k
dk
sin
θdθ
dϕ
(2π)3 0E − k200∞1
sin(k|r − r  |)=kdk
4π 2 |r − r  | −∞E − k2√exp(i E|r − r  |)
=−
.
(3.4)
4π|r − r  |(r.h.s. of eq.(3.3)) =In the last expression the limit E + iδ(δ → 0), corresponding to the outgoing-wave

boundary condition (retarded Green function), is assumed.

3.2. Partial wave expansion
Partial wave expansion of the free space Green function is obtained from the similar
expansion of the plane waves:il jl (kr)YL (k̂ )YL (r̂ ) = 4πexp(ik · r ) = 4π
LilL

Pl (r; E)
YL (k̂ )YL (r̂ )
kr

(3.5)

(notice: Pl in eq.(3.5) are not the Legendre polynomials.) Here YL are the real sphericalharmonics with the usual notation (use cos mφ and sin mφ instead of exp(imφ) for real7

harmonics) and L represents (l, m). By substituting the exponential term in eq.(3.3) andperforming the angular part integration, we obtain∞2g(r , r ; E) =
πrr


L

0

Pl(r; E)YL (r̂ ) Pl (r ; E)YL (r̂  )dk
E − k2

YL (r̂ )GL (r, r ; E)YL (r̂  )

=

(3.6)

L

By expanding the radial Green function g(r, r ; E) by the radial eigenfunctions
PL (r; En ) satisfying the boundary condition PL (R; En ) = 0, where R is some ﬁxed large
value (R  1, see Sec. 2), and using the relation
−→
n

R
π

∞

R

dk

and

0

0

PL2 (r; k 2 )dr −→

R
,
2

(3.7)

where k = {nπ + (lπ/2)}/R, both being valid for R  1, we ﬁnd
GL (r, r ; E) =g(r, r ; E)
.
rr

(3.8)

Finally, using the expression for g(r, r; E) given in Sec. 2 we obtain the partial wave
expansion of g(r , r  ; E):√g(r , r ; E) = − E
YL (r̂ )√= −i EL

Pl (r< ; E)
√
Er<



R+
l (r> ; E)
√
Er>



YL (r̂  )JL (r< ; E)HL (r> ; E),

(3.9)

√
JL (r; E) = jl ( Er)YL (r̂ )(3.10)√(1)
HL (r; E) = hl ( Er)YL (r̂ ).

(3.11)

L

where we deﬁneand3.3. General cases
Similarly the Green function of the system with a central, spherical potential may be
expressed (see the last paragraph of Sec. 2) as
 +


√
Rl (r> ; E)
Pl (r< ; E)√√
YL (r̂  )YL (r̂ )G(r , r ; E) = − E
Er<
Er>
L
√
JL (r< ; E)HL (r> ; E),= −i EL

8(3.12)

where Pl and R+
l = Ql + iPl are the regular and the irregular solution of the radial
Schrödinger equation, eq.(1.5), and the deﬁnition of JL and HL are understood; Pl andQl (hence Pl and R+l ) satisfy the Wronskian relation
W [Pl , Ql ] =√E.

(3.13)

4. Multiple scattering by muﬃn-tin potentials
4.1. Integral equation for the Green function
Suppose that the system is composed of an assembly of non-overlapping spherical potential(muﬃn-tin potential) centered on each lattice site of a crystal.
The crystal Green function G(r , r  ; E) can be calculated from the free space Greenfunction g(r , r  ; E) by
G(r , r  ) = g(r , r  ) +dr  g(r , r  )V (r  )G(r  , r  ).

(4.1)

Hereafter, the explicit notation of the energy dependence of the Green function, etc. isomitted for brevity.

4.2. Cell-centered expansion of the Green function
In eq.(4.1) the potential V (r ) denotes the assembly of muﬃn-tin potentials. Practically,eq.(4.1) can not be used for the calculation of G because of the convergence problems.
Instead we introduce the cell-centered representation of eq.(4.1). First we consider thecell-centered representation of g;g(r + R m , r  + R n ) = g(r , r  + R n − R m )√JL (r )HL (r  + R n − R m )=− E( for m = n). (4.2)

L

Using the addition theorem for H,
H(x + x  ) = 4π

il−l +l CLL L JL (x < )HL (x > ),L L

9(4.3)

where CLL L is the Gaunt number
CLL L =dr̂ YL (r̂ )YL (r̂ )YL (r̂ )and x > and x < is the bigger and the smaller (in the absolute value) of x and x  , we obtaing(r + R m , r  + R n ) =mnJL (r )gLL
 JL (r )( for m = n),

(4.4)

LL

where

√
mn
gLL
E
 = 4πiil−l +l CLL L HL (R m − R n ).

(4.5)

L
mm
Deﬁning gLL
 = 0, cell-centered representation of g may be expressed asg(r + R m , r  + R n ) = δmn g(r , r  ) +mnJL (r )gLL
 JL (r ).

(4.6)

LL

Similarly the crystal Green function has the cell-centered representationG(r + R m , r  + R n ) = δmn Gms (r , r ) +nJLm (r )Gmn
LL JL (r ).

(4.7)

LL
m
Here Gm
at m-ths is the Green function for the system of a single muﬃn-tin potential V

site given by√Gms (r , r ) = −i E

JLm (r< ; E)Hm
L (r> ; E),

(4.8)

L

which is exactly the same Green function as deﬁned by eq.(3.12).

4.3. Dyson-type equation for multiple scattering
In eq.(4.7) the second term represents the contribution of the multiple scattering (or “backscattering”) to the Green function and giving rise to the band structure. The relationmnbetween Gmn
LL and gLL is obtained from the cell-centered representation of eq.(4.1). By

direct insertion of eqs.(4.6) and (4.7) into eq.(4.2) and using the integral equations whichmmn
are satisﬁed by JLm and Gm
s for the single muﬃn-tin potential V , it is proven that GLLmn
is related to gLL
 by



mn
Gmn
LL = gLL +





mm m
mn
gLL
 tL GL L .
L m

10(4.9)

For periodic systems eq.(4.9) is solved by Fourier transform as−1

GLL (k ) =LgLL (k ) [1 − tg(k )]L L .

(4.10)

The crystal Green function is then given byG(r + R m , r  + R n ) = δmn Gms (r , r ) +JLm (r )LL×dkexp(−ik · (R m − R n ))
τ τ

−1

L

gLL (k ) [1 − tg(k )]L L JLn (r  ),(4.11)

where τ is the Brillouin zone volume.

5. Band structure calculations
5.1. Traditional KKR method
Standard KKR band structure calculation is performed by ﬁnding the zero of thedeterminant of the matrix appearing in eq.(4.11) (called KKR matrix) for each k ;det |δLL − tL (E)gLL (k ; E)| = 0.

(5.1)

This procedure gives the energy eigenvalues as a function of k and hence determines anenergy dispersion relation E(k ). The reason why the energy eigenvalues are obtainedin this way is the following: The zeros of the KKR matrix gives the poles of the Greenfunction (it actually is a one-to-one correspondence); each pole however corresponds tothe eigenstate of the Hamiltonian, as understood from the eigenfunction expansion of theGreen function.
In actual calculations, a numerical diﬃculty arises when the E(k ) falls near the freestate branches or the atomic energy levels. Both the free states and the atomic statescause poles in the determinant (poles of the structural Green function gLL (k ; E) or thoseof the t-matrix), changing the sign of the determinant. This can be a big trouble (however
diﬃcult to imagine!) in a root-ﬁnding procedure. For this reason, the KKR-method issometimes referred to as an unstable method. Though such a statement may not be quite
faire, it certainly is true that careless use of KKR sometimes leads to funny conclusions.
115.2. An alternative
A more direct way to obtain the electronic structure may be to calculate the crystal Greenfunction following eq.(4.11). Such a method is especially attractive because it enable us
to apply the KKR technique to various other classes of problem, e.g., impurity problems,disordered alloys, surfaces, interfaces, LEED states, and other spectroscopic properties.
The k-space integration appearing in eq.(4.11), however, causes a tough problem: Since
the imaginary part of the integrand is the δ-function, it is very hard to perform theintegration literally (actually impossible!). A conceivable way to overcome this problem is
to work with complex energies. For the complex energy the integrand, as a function of k ,becomes quite dull and the integration can be performed accurately.
Such a trick is possible since, ﬁrst, the back scattering term is the diﬀerence between
two Green functions, both of them being analytic, and hence analytic in the upper complex
half plane and, secondly, all we need in the self-consistent band structure calculation is not
the Green function but its energy integral. Owing to these facts, we can safely deform theenergy integration path toward the upper complex plane and still can hold all necessary
information. This also increases the accuracy of the energy integration because the Greenfunction is quite structureless oﬀ the real axis as a function of E.
Contrary to the traditional method mentioned in the previous subsection, thisalternative way is free from the spurious poles. Moreover numerical steps needed forthis procedure increase, roughly speaking, as N 2 , where N is the number of atoms perunit cell, which is compared with N 3 needed for the traditional method.
Another big advantage of this method is that the procedure is quite eﬃciently
vectorized in the sense of the vector processing of the supercomputers. This point, thoughnothing much physics, is important in actual calculations.

5.3. Total-energy calculation
Once we obtain the Green function we can calculate, in principle, all the ground stateelectronic properties in the framework of the muﬃn-tin potential model. For instance, thespin dependent local charge distributions are calculated from the imaginary part of energy
integral of the site diagonal Green function:
ρ(r ) = −

1

π

EF
−∞

12

dEG(r , r ; E),(5.2)

where G(r , r ; E) is given by eq.(4.11) with m = n = 0. In order to calculate the totalenergy of the system we exploit the muﬃn-tin charge model where the spherical symmetriccharge distribution inside the muﬃn-tin sphere and the uniform charge distribution forinterstitial region is assumed. In this case we can express the total energy per atom as asum of various separate contributions:
Etot = Eband − Epot + Estat + Exc(5.3)

To calculated the band energy, we introduce a function Z(E), whose imaginary part gives
the total number of states (per atom) below a real energy E:
Z(E) =dkτ τ

ln det | − E + (k + g )2 | + ln det |1 − tg(k ; E)| exp(iηL (E)) − ln det ΛL (E) (5.4)

Here ΛL (E) is deﬁned by eq(1.19), giving a normalization for the phase shift ηL . Sucha normalization of the phase shift makes the analytic continuation of Z(E) to the uppercomplex half plain possible; the procedure removes all singularities which exp(iηL (E)) will
have in the upper half plain at the resonances in the single site scattering.
The band energy is now calculated by integrating Z(E) with respect to E up to the
Fermi level:1
Eband = −  Z(EF )EF −
π



EF

Z(E)dE(5.5)−∞

We subtract from the band energy the potential energyr0

Epot =

r2 drρ(r)V 0 (r),

(5.6)

0

which is the expectation value of the single site potential at 0-th site V 0 , to obtain the
kinetic energy of the non-interacting system in the sense of the local density functional
theory. The electrostatic energy Eelectrostatic (par atom, in atomic unit) is obtained as
r0

Estat = 32π 2r

dr rρ(r)0r0dr r2 ρ(r ) − 8π Q

0rρ(r)dr0

1 Q2out
− C
,
2
a

(5.7)

where a is the lattice constant and the constant C, depending on the crystal structure,is given as, e.g., C =4.8320664, 4.085521 and 3.1166857 for fcc, bcc and sc lattice,13

respectively; Q is the nuclear charge and Qout is the number of electrons per atom
accumulated at the interstitial region, i.e.,r0

Qout = Q −r2 drρ(r).

(5.8)

0

The exchange-correlation energy (per atom) is given by the LSD approximation as
r0

Exc =

r2 drρ/xc (ρ) + Qout /xc (ρout ),

(5.9)

0

where /xc is the exchange-correlation energy suitable for the homogeneous electron gas.
All the core contributions should be included in the above expression.

5.4. Pressure
For pressure we use a mixture of the pressure cell-boundary relation of Liberman–Pettifor
and the usual expression which relates the pressure to a sum of twice the kinetic energy
and the potential energy with an additional exchange-correlation contribution. We apply
the former inside the muﬃn-tin sphere and use the latter for the interstitial region:EF

3P Ω =
L−∞dE ρL (r0 , E)r0 (E − V (r0 )) r02 + (DL − l)(DL + l + 1)
EF

+ Eband −

r0

dE E
L

−∞

r2 dr ρL (r, E)

0





− 4πr03 ρ(r0 ) /xc ρ(r0 ) − µxc ρ(r0 )−1 Q2out
C− 3Qout {/xc (ρout ) − µxc (ρout )}2
a(5.10)

Here, ρL (r, E) is the partial density of states inside the muﬃn-tin sphere, i.e., the Lth component of the negative imaginary part of the Green function divided by π (seeeqs.(4.11) and (5.2)), DL = r0 JL (r0 , E)/JL (r0 , E) is related to the logarithmic derivativeof the radial wave functions at r = r0 and µxc is the exchange correlation potential.
An advantage of this mixture is that in this way we can avoid the numerical problem
caused by a big cancellation between kinetic and potential energies without introducing
the atomic sphere approximation which the Liberman–Pettifor approach (usually used inLMTO) had to use for a surface integral. Furthermore, we can use this method both for
non-relativistic and relativistic cases. The usual treatment, on the otherhand, fails for a14

relativistic case since the sum of twice the kinetic energy and the potential energy nevercorresponds to the energy cost for the volume dilation in this case. The reason why thepresent method works even for the relativistic treatment is that we calculate the pressure
at the muﬃn-tin sphere from the momentum ﬂow across the sphere; the procedure is valid
both for relativistic and non-relativistic cases. For the interstitial region the relativistic
eﬀects are not important any more, which allows us to apply the second method. A
key point is that we need not use the atomic sphere approximation, in contrast to the
Liberman–Pettifor approach, for the second method and thus can avoid the additional
approximation.

6. Application of the KKR method
6.1. Impurity problems
Suppose an impurity atom located, say, at the origin is embedded in an otherwise perfect
host crystal. Now let us assume that the potentials at the neighboring host sites are notaﬀected by the presence of the impurity atom (single site approximation). In this case, all
we have to know is the site diagonal Green function at the origin G(r , r ), where G is given
by eq.(4.7). This problem is immediately solved if the back scattering term of eq.(4.7) forthe perfect host crystal G̃00
LL (i.e. m = n = 0) is known:

−1
00
00
1
−
(t
=
−
t
)
G̃
G̃
G00


impurity
host
LL
LL

L L

L

,

(6.1)

where
G̃00
LL =

dk
G̃LL (k ) =
τ τ

dk
τ τ

−1

L

gLL (k ) [1 − thost g(k )]L L(6.2)

may be calculated for the host crystal once and while. The above treatment might not
be enough when the perturbation caused by the impurity extends beyond the single site.
More sophisticated treatment which can deal with not only the impurity site but also thesurrounding host sites simultaneously is required in these cases; the extension of eq.(6.1)
is rather straightforward.

6.2. KKR-coherent potential approximation
For disordered systems, the coherent potential approximation (CPA) is one of the mosteﬃcient way to solved the problem. What are discussed about the disordered systems are15

averaged properties. In the present context, it is the conﬁguration average of the Greenfunction. In CPA, averaged Green function is expressed in terms of the eﬀective medium; itis the assembly of the ﬁctitious atoms whose atomic potential is speciﬁed by the coherent
t-matrix t̃(E). Once the coherent t-matrix is known, any ground state properties arecalculated within the single site approximation (actually it is proven to be the best single
site treatment).
To determine the coherent t-matrix, we consider the following impurity problem: Let
the disordered system be A1−x Bx . Suppose we put either A or B atom in the eﬀective
medium as an impurity. The impurity problems in these cases are solved if t̃ is known (seethe preceding subsection); the site diagonal Green function is then given by
A (or B)
=
GLL

G̃00
LL
L


−1
00
,
1 − (tA (or B) − t̃) G̃
 
L L

where
G̃00
LL =

dk
G̃LL (k ) =
τ τ

dk
τ τ

L


−1
gLL (k ) 1 − t̃g(k ) L L .

(6.3)

(6.4)

Now the CPA prescription determining t̃ is
B00(1 − x)GA
LL + xGLL = G̃LL .

(6.5)

In other words, the t̃ should be determined in such a way that the scattering disappears
if the single site scattering by A and by B impurity located at the origin in the eﬀective
medium is conﬁgurationally averaged.
All the average quantities are easily calculated if t̃ is obtained by solving the above
self-consistent equation. As for the problems associated to the numerical treatment of the
CPA equation and some eﬃcient techniques to get the procedure run, see elsewhere.

16