A full-potential approach to the relativistic single-site Green's function

0.77 MB·

by  Xianglin Liu

A full-potential approach to the relativistic
single-site Green’s function ∗
6
1 Xianglin Liu
0
2 Department of physics, Carnegie Mellon University

r Yang Wang
a
M
Pittsburgh supercomputing center, Carnegie Mellon University
Markus Eisenbach
5

] Center for Computational Sciences, Oak Ridge National Laboratory
i
c G. Malcolm Stocks
s
-
l Materials science and technology division, Oak Ridge National Laboratory
r
t
m
March 8, 2016
.
t
a
m
-
d Abstract
n
Onemajorpurposeofstudyingthesingle-sitescatteringproblemis
o
c toobtainthescatteringmatricesanddifferentialequationsolutionsin-
[ dispensable to multiple scattering theory (MST) calculations. On the

2 other hand, the single-site scattering itself is also appealing because
v it reveals the physical environment experienced by electrons around
6
the scattering center. In this paper we demonstrate a new formalism
7
3 to calculate the relativistic full-potential single-site Green’s function.
4
Weimplementthismethodtocalculatethesingle-sitedensityofstates
0
and electron charge densities. The code is rigorously tested and with
.
1
the help of Krein’s theorem, the relativistic effects and full potential
0
6 effects in group V elements and noble metals are thoroughly investi-
1 gated.
:
v
i ∗This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-
X
AC05-00OR22725 with the U.S. Department of Energy. The United States Government
r retains and the publisher, by accepting the article for publication, acknowledges that the
a
UnitedStatesGovernmentretainsanon-exclusive,paid-up,irrevocable,world-widelicense
to publish or reproduce the published form of this manuscript, or allow others to do so,
for United States Government purposes. The Department of Energy will provide public
access to these results of federally sponsored research in accordance with the DOE Public
Access Plan(http://energy.gov/downloads/doe-public-access-plan).
1
1 Introduction
The multiple-scattering theory (MST) is a DFT based ab-initio method that
iswidelyappliedtothecalculationoftheelectronicstructureofmetals, alloys
and impurities. One crucial step in MST calculation is to solve the single-site
scattering problem. Combined with the position information of the atoms,
solutions of the single-site scattering problem can be used to construct the
Green’s function of the whole system, from which most physical quantities
can be extracted. Because of the essential role played by the single-site
scattering in MST, there have been a constant effort to improve it in the
last few decades [1]. The earliest MST solves the Schr¨odinger’s equation
and employs the muffin-tin potential approximation, i.e. the potential is
spherically symmetric within a bounding sphere and is constant outside.
Then the generalization to full potential (FP)[2] and relativistic cases [3, 4,
5, 6], and eventually combination of the two, i.e. relativistic full potential
(RFP), were proposed and implemented using various schemes [7, 8, 9, 10].
AmongthemonewidelyusedRFPMSTcodeisdevelopedbyHuhneetal[11]
to either solve the coupled integral equations iteratively using Born’s series
expansion or to directly solve the coupled differential equations. In a recent
paper[12], the coupled differential equations are also solved by matching the
regular solutions at the boundary of the Wigner–Seitz cell.
In this paper we present a alternative formalism to tackle the RFP single-
site scattering problem by directly solving the Dirac differential equation.
This method is a relativistic generalization of the non-relativistic FP MST
method[13]. Compared to other RFP methods, this new formalism has the
feature that the differential equation, the t matrix and the single-site Green’s
function are all expressed in terms of the r dependent sine and cosine matri-
ces. We would like to point out that the sine and cosine matrices technique
to obtain the solutions for relativistic scattering theory was first proposed
by X. D. Wang et al. [8], based on the phase integral technique[14], and
later implemented by S. B. Kellen and A. J. Freeman [10]. The major dif-
ference between our method and the Kellen & Freeman’s approach is in the
calculation of observables. They find the observables by searching for energy
eigenvalues and eigenfunctions of the KKR secular equation, while in our
case all physical quantities are expressed in terms of the Green’s functions,
without the compute-intensive eigenvalue searching and the wavefunction
orthonormalization procedures.
As a test of our code, the single-site DOS are calculated using both the
Green’s function method and the Krein’s theorem method and two results
are compared. To investigate the relativistic effects and the full potential
effects, we study the density of states of noble metals and group V elements.
2
Finally, the charge density of tantalum is calculated and some interesting
relativistic features are discussed.
2 Theory
In the following discussion, we use the atomic Rydberg units; i.e., (cid:126) = 1,m =
1/2, c ≈ 274, etc. To construct the Green’s function, we first need to solve
the single-particle Dirac equation:
{−cαi∇+βmc2 +V(r)}ψ(E,r) = Wψ(E,r), (1)
where α, β are the Dirac matrices, W and E are the relativistic energy and
relativistic kinetic energy, respectively. W and E are related to momentum
p by
(cid:112)
W = (mc2)2 +(pc)2, (2)
E = W −mc2. (3)
The solution ψ(E,r) is a Dirac spinor. Within the framework of spin po-
larized relativistic density functional theory [15], the effective potential is
written as
(cid:18) (cid:19)
v(r)+σ ·B(r) 0
V(r) = , (4)
0 v(r)−σ ·B(r)
where σ are the Pauli matrices and B are effective magnetic fields. In the
following discussion we will focus on nonmagnetic systems and discard the B
terms. For r > R , where R is the radius of the circumscribed sphere of the
c c
Wigner-Seitz cell, the effective potential V(r) vanishes, and the solutions of
the corresponding free-space Dirac equation are well known. The right-hand
solutions are
(cid:18)W +mc2(cid:19)1/2(cid:18) j (pr)χ (ˆr) (cid:19)
J (E,r) = l Λ , (5)
Λ c2 Wip+cmSκc2j¯l(pr)χΛ¯(ˆr)
(cid:18)W +mc2(cid:19)1/2(cid:18) n (pr)χ (ˆr) (cid:19)
N (E,r) = l Λ , (6)
Λ c2 Wip+cmSκc2n¯l(pr)χΛ¯(ˆr)
where j (pr) and n (pr) are the usual spherical Bessel functions of the first
l l
kind and the second kind, with angular momentum index l. Λ stands for the
pair of relativistic angular-momentum indices (κ,µ) and S is the sign of κ
κ
¯
index. Λ = (−κ,µ) and
(cid:40)
l+1 if κ < 0
¯
l = . (7)
l−1 if κ > 0
3
χ and χ are the spin-angular functions
Λ Λ¯
(cid:88) 1
χ (ˆr) = C(l,j, |µ−m ,m )Y (ˆr)φ , (8)
Λ 2 s s l,µ−ms ms
ms=±1
(cid:88) 1
χ (ˆr) = C(¯l,j, |µ−m ,m )Y (ˆr)φ , (9)
Λ¯ 2 s s ¯l,µ−ms ms
ms=±1
where C(l,j, 1|µ − m ,m ) are the Clebsch-Gordan coefficients, Y are
2 s s l,µ−ms
spherical harmonics and φ are Pauli spinors
ms
(cid:18) (cid:19) (cid:18) (cid:19)
1 0
φ = , φ = . (10)
1/2 0 −1/2 1
In addition to the right-side solutions, we also need the left-hand solutions
to construct the Green’s function (see Appendix A and B). In free-space, the
left-hand solutions are given by
(cid:18)W +mc2(cid:19)1/2(cid:18) −ipcS (cid:19)
J+(E,r) = j (pr)χ†(ˆr), κ j (pr)χ†(ˆr) , (11)
Λ c2 l Λ W +mc2 ¯l Λ¯
(cid:18)W +mc2(cid:19)1/2(cid:18) −ipcS (cid:19)
N+(E,r) = n (pr)χ†(ˆr), κ n (pr)χ†(ˆr) . (12)
Λ c2 l Λ W +mc2 ¯l Λ¯
For r < R , the Dirac equation solutions can only be obtained numeri-
c
cally. By using the phase integral technique[14] , the solutions are written in
terms of the free-space solutions as (see Appendix A)
(cid:88)
ψ (E,r) = {S (E,r)N (E,r)−C (E,r)J (E,r)}, (13)
Λ Λ(cid:48)Λ Λ(cid:48) Λ(cid:48)Λ Λ(cid:48)
Λ(cid:48)
where the r dependent cosine matrix C (E,r) and sine matrix S (E,r)
Λ(cid:48)Λ Λ(cid:48)Λ
are defined as
(cid:90)
C (E,r) = p d3r(cid:48)N+(E,r(cid:48))V(r(cid:48))ψ (E,r(cid:48))−δ , (14)
Λ(cid:48)Λ Λ(cid:48) Λ ΛΛ(cid:48)
0<r(cid:48)<r
(cid:90)
S (E,r) = p d3r(cid:48)J+(E,r(cid:48))V(r(cid:48))ψ (E,r(cid:48)). (15)
Λ(cid:48)Λ Λ(cid:48) Λ
0<r(cid:48)<r
Note that this expression is also valid for r > R , because S (E,r) and
c Λ(cid:48)Λ
C (E,r) will become constants outside the circumscribed sphere. To dis-
Λ(cid:48)Λ
tinguish these constants from S (E,r) and C (E,r) we denote them as
Λ(cid:48)Λ Λ(cid:48)Λ
S (E) and C (E). Equations (13), (14) and (15) form a set of coupled
Λ(cid:48)Λ Λ(cid:48)Λ
4
integral equations. C (E,r) and S (E,r) can be obtained by solving the
ΛΛ(cid:48) ΛΛ(cid:48)
corresponding differential equations
(cid:90)
d
S (E,r) = r2 dˆr pJ+(E,r)V(r)ψ (E,r), (16)
dr Λ(cid:48)Λ Λ(cid:48) Λ
(cid:90)
d
C (E,r) = r2 dˆr pN+(E,r)V(r)ψ (E,r). (17)
dr Λ(cid:48)Λ Λ(cid:48) Λ
Note that the integral is only upon the angular part. For regular solutions,
the boundary conditions are
ψ (E,r) = J (E,r). (18)
Λ r→0 Λ
or equivalently
S (E,0) = 0, (19)
Λ(cid:48)Λ
C (E,0) = −δ . (20)
Λ(cid:48)Λ ΛΛ(cid:48)
To see the boundary conditions are well defined, we note the effective electric
potential has a spherically symmetric 1/r behavior at the origin. As a result
the integrals on the right side of equation (14) and (15) vanish at the origin
because spherical bessel functions with different l indices can not be coupled
by spherical potential.
Next, we proceed to discuss technical details on solving the coupled dif-
ferential equations. The explicit expressions of the differential equations (16)
and (17) are given by
d (cid:88) (cid:88)
S (E,r) = a (E,r)S (E,r)− b (E,r)C (E,r), (21)
Λ(cid:48)Λ Λ(cid:48)(cid:48)Λ(cid:48) Λ(cid:48)(cid:48)Λ Λ(cid:48)(cid:48)Λ(cid:48) Λ(cid:48)(cid:48)Λ
dr
Λ(cid:48)(cid:48) Λ(cid:48)(cid:48)
where
(cid:90)
a (E,r) = r2 dˆr pJ+(E,r)V(r)N (E,r), (22)
Λ(cid:48)(cid:48)Λ(cid:48) Λ(cid:48) Λ(cid:48)(cid:48)
(cid:90)
b (E,r) = r2 dˆr pJ+(E,r)V(r)J (E,r), (23)
Λ(cid:48)(cid:48)Λ(cid:48) Λ(cid:48) Λ(cid:48)(cid:48)
and
d (cid:88) (cid:88)
C (E,r) = c (E,r)S (E,r)− d (E,r)C (E,r),
Λ(cid:48)Λ Λ(cid:48)(cid:48)Λ(cid:48) Λ(cid:48)(cid:48)Λ Λ(cid:48)(cid:48)Λ(cid:48) Λ(cid:48)(cid:48)Λ
dr
Λ(cid:48)(cid:48) Λ(cid:48)(cid:48)
(24)
5
where
(cid:90)
c (E,r) = r2 dˆr pN+(E,r)V(r)N (E,r), (25)
Λ(cid:48)(cid:48)Λ(cid:48) Λ(cid:48) Λ(cid:48)(cid:48)
(cid:90)
d (E,r) = r2 dˆr pN+(E,r)V(r)J (E,r). (26)
Λ(cid:48)(cid:48)Λ(cid:48) Λ(cid:48) Λ(cid:48)(cid:48)
For full potential calculation, the effective potential V(r) is expanded in
terms of the spherical harmonics
(cid:88)
V(r) = V (r)Y (ˆr). (27)
Lv Lv
Lv
The angular integral in equation (22),(23),(25) and (26) can be done analyt-
ically and written in terms of the Gaunt coefficients
(cid:90)
CL(cid:48)(cid:48) = dˆr Y (ˆr)Y∗(ˆr)Y (ˆr). (28)
L,L(cid:48) L L(cid:48) L(cid:48)(cid:48)
In practice, the differential equations are solved on exponential radial grid
r = r exp(x), using the fourth-order Bashforth-Adams-Moulton predictor-
0
corrector method.
After we obtain the solutions of the Dirac equation, we can use them to
construct the single-site Green’s function, which has the following expression
(see Appendix B)
(cid:88) (cid:88)
G(E,r,r(cid:48)) = Z (E,r)t (E)Z+(E,r(cid:48))− Z (E,r)J+(E,r(cid:48)). (29)
Λ ΛΛ(cid:48) Λ(cid:48) Λ Λ
ΛΛ(cid:48) Λ
whereJ+(E,r)arethesolutionsofequation(1)withtheboundarycondition
Λ(cid:48)
that J+(E,r) = J+(E,r) outside the Wigner-Seitz cell. Z (E,r) is another
Λ(cid:48) Λ(cid:48) Λ
set of regular solutions defined as
(cid:88)
Z (E,r) = p ψ (E,r)S−1 (E). (30)
Λ Λ(cid:48) Λ(cid:48)Λ
Λ(cid:48)
The t matrix is given by equation (70) in Appendix B and the explicit ex-
pression is
1 (cid:88)
t (E) = − S (E)(C (E)−i S (E))−1. (31)
ΛΛ(cid:48) ΛΛ(cid:48)(cid:48) Λ(cid:48)(cid:48)Λ(cid:48) Λ(cid:48)(cid:48)Λ(cid:48)
p
Λ(cid:48)(cid:48)
FromtheGreen’sfunction,it’sstraightforwardtocalculatephysicalquan-
tities of the system. Here we focus on the charge density and the density of
6
states. The first one is necessary to calculate the new potential in a SCF cy-
cle and the latter one is needed to determine the Fermi energy of the system.
The charge density is given by
1 (cid:90) EF
ρ(r) = − Im Tr G(E,r,r)dE, (32)
π
0
where the Fermi energy E is chosen to give the correct number of total
F
electrons. The density of states is given by
(cid:90)
1
n(E) = − Im Tr G(E,r,r)dr, (33)
π
Ω
where Ω denotes the wigner seitz cell. For real energy, the imaginary part of
the ZJ+ term in equation (29) vanishes, so we can instead write the density
of states n(E) and the charge density ρ(r) as
(cid:90)
1 (cid:88)
n(E) = − Im Tr Z (E,r)t (E)Z+(E,r)dr, (34)
π Λ ΛΛ(cid:48) Λ(cid:48)
Ω ΛΛ(cid:48)
1 (cid:90) EF (cid:88)
ρ(r) = − Im Tr Z (E,r)t (E)Z+(E,r)dE. (35)
π Λ ΛΛ(cid:48) Λ(cid:48)
0 ΛΛ(cid:48)
Thebenefitofthisexpressionisthatitdoesnotcontaintheirregularsolutions
J (r), which are difficult to evaluate precisely near the origin. Even though
Λ
the focus of this paper is on single-site Green’s function , at the end of
this section we would like to make a few comments about the MST Green’s
function, which has the expression [16]
(cid:88) (cid:88)
G(E,r,r(cid:48)) = Z (E,r)τ (E)Z+(E,r(cid:48))− Z (E,r)J+(E,r(cid:48)). (36)
Λ ΛΛ(cid:48) Λ(cid:48) Λ Λ
ΛΛ(cid:48) Λ
Compare it with equation (29) we see the only difference is the replacement
of the t matrix by the scattering-path matrix τ [17], which describes the
ΛΛ(cid:48)
scattering of the electron along all paths in the solid. The scattering-path
matrix can be constructed from the t matrix and the structure constants
G , which describe the propagation of free electrons between two atoms.
ΛΛ(cid:48)
Therefore, it’s straightforward to implement our formalism to the MST cal-
culations.
3 Krein’s Theorem
For single-site scattering, the density of states can be calculated from either
theGreen’sfunctionortheKrein’stheorem[18][19]. Acomparisonofthetwo
7
results will be a good test of our code. First we focus on the Krein’s theorem
method. Details of the relation between the Krein’s theorem and the Green’s
function in the non-relativistic case has been derived in a previous paper [20]
and most of them can also be applied to the relativistic case. Therefore,
in this section we will just establish notation and quote known results. By
applying the Krein’s theorem to the scattering theory, it has been proved
that the integrated density of states (IDOS) is given by
N (E) = −2ξ(E) = N(E)−N (E)+n , (37)
K 0 c
where N(E) is the single-site density of states, N (E) is the free electron
0
integrated density of states, n is the total number of core electrons, which is
c
irrelevant here since we are only interested in valence electrons. ξ(E) is the
Krein’s spectral shift function [18] related to the standard unitary S-matrix
S(E) by [21]
e−i2πξ(E) = detS(E). (38)
The S-matrix is obtained from the t matrix using the relation
S (E) = δ (E)−2ipt (E). (39)
ΛΛ(cid:48) ΛΛ(cid:48) ΛΛ(cid:48)
The Krein DOS n (E) is obtained by taking the derivative of N ,
K K
dN (E)
K
n (E) = . (40)
K
dE
Note that n (E) is the DOS for the entire space. To compare the DOS
K
obtained from the Green’s function with n (E), we express the DOS inside
K
and outside the region bounded by a sphere of radius R in terms of Green’s
c
function as
(cid:90)
1
n (E) = − Im Tr(G(E,r,r)−G (E,r,r))dr, (41)
in 0
π
0<r<Rc
(cid:90)
1
n (E) = − Im Tr(G(E,r,r)−G (E,r,r))dr, (42)
out 0
π
r>Rc
and we should have n (E) = n (E) + n (E). The G(E,r,r) term in n
K in out in
is evaluated numerically and the G (E,r,r) term is evaluated analytically
0
using the following expression of the free space Green’s function
1 eip(r−r(cid:48))
Tr G (E,r,r) = lim Tr− (α·p+βmc2 +W) (43)
0 r(cid:48)→r c2 4π(r−r(cid:48))
ipW
=− . (44)
πc2
8
Therefore, the contribution to the DOS from the inside integral is
1 (cid:90) Rc (cid:88) 4pW
n (E) = − Im Tr Z (E,r)t (E)Z+(E,r)dr− R 3. (45)
in π Λ ΛΛ(cid:48) Λ(cid:48) 3c2π c
0 ΛΛ(cid:48)
For r > R , using another expression of the Green’s function
c
(cid:88)
G(E,r,r) = G (E,r,r)−p2 H (E,r)t (E)H+(E,r), (46)
0 Λ ΛΛ(cid:48) Λ(cid:48)
ΛΛ(cid:48)
where H (E,r) = J (E,r)+iN (E,r), the outside contribution to the DOS
Λ Λ Λ
is
2p2W (cid:90) ∞(cid:88)
n (E) = Im h (pr)t h (pr)r2dr. (47)
out πc2 l ll l
Rc l
This double Hankel function integral also occurs in the non-relativistic case
[20] and can be done analytically using [22]
(cid:90) ∞ R3
(cid:0) (cid:1)
h (pr)2r2dr = h (pr)h (pr)−h (pr)2 . (48)
l l−1 l+1 l
2
R
To compare the Krein’s theorem method and the Green’s function method,
the DOS of copper and vanadium are calculated and shown in figure 1.
The input potentials are obtained from non-relativistic self-consistent full-
potential KKR calculations. In both calculations 256 energy points are used.
In general the relative error is within the order of 10−4, which indicates an
excellent agreement between the two methods. Most of the errors are from
the small number of energy points and the relative primitive way to calculate
n from the numerical derivative of N .
K K
4 Singel-site DOS
The full potential effects and the relativistic effects can be directly observed
in the single-site DOS. As the first example, the DOS of noble metals, i.e.,
copper, silver and gold are calculated using both relativistic and nonrela-
tivistic Green’s function methods. The results are shown in figure 2. In all
calculations the angular momentum expansion cut-offs l are set to be 4
max
and the expansion cut-offs for potentials are set to be 2×l to satisfy the
max
angular momentum triangular relation. Both the relativistic effects and the
fullpotentialeffectsarewelldemonstratedbythethreepeaksintherelativis-
tic DOS plots of silver. The large energy differences between the leftmost
9
30
150
25
Cu V
S S
O O 20
SiteD 100 nin!nout SiteD 15 nin!nout
" "
Single 50 nK Single 10 nK
5
0 0
0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2
E Ry E Ry
Figure 1: (Color onl!in"e) Comparison of the DOS from! " Green’s function
method and the Krein’s theorem method. The blue solid lines show the
Krein DOS n . The dashed lines show n +n calculated using the Green’s
K in out
function. Because of the good agreement the two lines actually overlap.
peak and the right two peaks are due to relativistic effects, mostly spin-orbit
coupling. The smaller difference between the right two peaks is due to full
potential effect. The full potential effects of copper and gold, however, do
not cause any visible splitting of DOS because their peaks on the right are
relatively broad and merge into one single peak. A better demonstration
of the full potential effects and the relativistic effects is the splitting of the
IDOS components of the d electrons. According to Krein’s theorem, the
IDOS components are given by the generalized phase shifts divided by 2π,
where the generalized phase shifts are obtained by diagonalizing the S-matrix
[23]. As an example we calculated the Krein IDOS components of copper
in different cases and the results are shown in figure 3. For non-relativistic
muffin-tin calculation, all the d electrons degenerate. For relativistic muffin-
tin calculation, even though the input potential is spherically symmetric, the
IDOS still splits into two parts due to spin-orbit coupling. The introduc-
tion of asymmetric potential leads to further splitting in the relativistic full
potential calculation. Because of the cubic structure symmetry, the split-
ting is not complete and the degeneracies still exist. Similar calculations
are also performed for silver and gold and the Krein IDOS for relativistic
full-potential calculation are shown in figure 4, where the splittings due to
spin-orbit coupling are more significant. The width of the spitting provides
an estimate of the strength of the spin-orbit coupling or full potential effect.
For example, for copper, silver and gold the splittings caused by spin-orbit
coupling are 0.020 Ry, 0.042 Ry and 0.13 Ry, respectively. This agrees with
our expectation that the magnitude of spin-orbit coupling tend to increase
for heavier elements.
All the noble elements have face centered cubic structure, therefore the
10