- Email: [email protected]
Non-local pseudopotentials in the multi-slice method of calculating electron wavefunctions in crystals
ultramicroscopy ELSEVIER
Ultramicroscopy 55 (1994) 221-227
Non-local pseudopotentials in the multi-slice method of calculating electron wavefunctions in crystals Kyozaburo Kambe, Catherine Stampfl Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, D-14195 Berlin (Dahlem), Germany Received 3 December 1993 Dedicated to Gunter Lehmpfuhl on the occasion of his 65th birthday
Abstract
A method for including non-local potentials in the multi-slice formalism for calculating electron wavefunctions in crystals is presented. With the assumption that non-local potentials appear only in atomic core regions as correction terms having a so-called separable form, it is shown that these corrections can be included as an additional procedure carried out after each passage through a core region. Possible ways of deriving pseudopotentials of the present type for various fields of electron diffraction and electron microscopy are briefly discussed.
1. Introduction
After the multi-slice method was introduced by Cowley and Moodie [1], it has been widely used for theoretical interpretation of electron-diffraction patterns and electron micrographs. The range of application has been expanded to the low-energy regime, and thus to low-energy electron diffraction [2] and current image diffraction [3]. The method has been used successfully to calculate convergent b e a m [4] and usual [5] reflection high-energy electron diffraction intensities, and was shown recently by Stampfl et al. [6] to be well applicable to calculations of usual bulk band structures and surface electronic states, and consequently to, for example, angle-resolved ultraviolet photoemission spectroscopy [7]. The applicability of the method to these last cases relies on the use of pseudopotentials, which enables us to
choose relatively small numbers of beams and slices per unit lattice period. Pseudopotentials have been implemented in band theory to replace true atomic potentials for the purpose of keeping calculation time manageable without changing the obtained results. Recent developments in the theory of pseudopotentials [8], however, tend to have non-local forms (see below) which cannot be used directly in the usual multi-slice method. In order to utilize fully the advantage of the method, in particular for surface problems, it is necessary to modify the method to incorporate non-local potentials. In the present p a p e r we report the method to include non-local corrections which has been developed for a special type of non-local pseudopotential [9-12] which is currently in use in self-consistent band-structure calculations. Hence, the electron is assumed to lie in the energy range of
0304-3991/94/$07.00 © 1994 Elsevier Science B.V. All rights reserved SSDI 0 3 0 4 - 3 9 9 1 ( 9 4 ) 0 0 0 5 6 - S
222
K. Kambe, C. Stampfl / Ultramicroscopy 55 (1994) 221-227
valence bands. It is believed, however, that similar types of pseudopotentials can be constructed appropriately for the use in low- and high-energy electron-diffraction and -microscopy problems when a reduction of calculation time is desired. We discuss these possibilities in Section 7.
2. Form of the pseudopotential
After application of the method of Kleinman and Bylander [9] (see also Refs. [10,11]) we may choose the pseudopotential to give, in the Schr6dinger equation 2m V2~+ 7(E-
V)~=0,
(1)
the potential term in the form
V(r)&(r) = vl°cal(r)&(r) + EUilm(r-ri) ilm
xfWam(r'-ri)O(r' )
3. Two-dimensional expansion of the Schr6dinger equation
For obtaining wavefunctions for crystal slabs or crystals with surfaces we calculate at first those solutions of the Schr6dinger equation (1) for which the values of energy E and the component of the wavevector parallel to the surface kll are given. It is convenient to expand 0 ( r ) in a two-dimensional Fourier series parallel to the surface
O(r)
= El]Ira(Z)
(4)
where r t is the tangential component of r to the surface, z is taken normal to the surface and is positively increasing into the crystal, m represents the two-dimensional Miller indices (h, k), and kmt is the tangential component of the wavevector k to the surface given by (5)
kmt = kot + bmt ,
d3r '.
(2)
VL°Cat(r) is the main part of the potential and is local; the second term is the correction having a "separable" non-local form. r i represents the position of the ith atom center. Expressions for the functions g i l m ( r - r i) and W i l m ( r ' - - r i) are given in Section 7. For brevity we may write
where k0t = kll , and brat is a two-dimensional reciprocal lattice vector. Substitution of Eq. (4) into Eq. (1) gives d2
dz2Om(Z)
+ k 2 m z t ~ m ( z ) --
EVm_.(Z)O.(Z) n
- ~_. ~UmL(Z)fbLWL~(Z')~O~(Z ') n
E Uilm(r - ri)Wilm(r'- ri) = ~-, Ul_(r)Wl_(r'), ilm
exp(ikmt "rt),
m
dz' = O,
~al•
L
(6)
L
(3) L representing ilm. The solution may be obtained at first by using V ~°ca' for V, and the non-local correction term may be introduced for modifying the solution. It is to be noted that the non-local correction term in Eq. (2) is represented in a so-called "degenerate" (Courant-Hilbert [13]) or "finite" (MorseFeshbach [14]) form, that is, a sum of a finite number of terms of a product of functions of r and r' only. In the theory of integral equations it is known that an integral equation having a kernel in a degenerate form can be solved relatively easily. We are going to derive these types of integral equations by applying the method of variation of constants.
where k2z --
2m h-----TE-]kmt[ 2,
(7)
and Um( Z ) -- _2m _
1 f vl°cal(r) e x p ( - i b m t h 2 A "A
"rt)
d2rt,
(8) 1 fAUL(r ) WL"(Z')
e x p ( - i b m t "rt) d2rt,
(9)
2m 1
h 2 A fA WL(r')
e x p ( i b n t "r't) d2r't'
(10)
K. Kambe, C. Stampfl / Uhramicroscopy 55 (1994) 221-227
where A indicates the area of the two-dimensional unit cell. The ranges of integrals (a L, b L) in Eq. (6) are taken to include the ranges of non-vanishing values of WLn(Z'), and may depend on L, particularly on i through the position vector r i. In Eq. (6) the summation over i includes only the atoms in a two-dimensional unit cell. It is convenient to reduce Eq. (6) to a system of first-order differential equations by putting
223
/2z~0. We see that the expression (13) can be interpreted as a linear combination of these solutions
xO(z)
=
E 6°(Zo) 2N
+
%(z))" E x°_~(Zo)(~'x%(z)
v=N+l
d
dz~bm(Z)
= Xm(Z),
(14)
(11)
and writing Eq. (6) as d
dzX,,( z) + k2mz$m(z)
5. Variation of constants - 2 U m _ n ( Z)~ln( Z ) n
- ~n Y'~UmL(Z)f~:
dz'=O. (12)
We assume that it is sufficient to take into account only N equations of Eq. (6), that is, 2 N equations of the system (11), (12).
Using the solutions (O°,(z), x O v ( g ) ) T, which are assumed here to be already obtained, we are going to solve Eqs. (11), (12) by using the method of variation of constants. We put
( )2N
= E c.(z)[ ~'%(z)
Xm(Z)
.=l
lxOv(Z)
)
(15)
and iook for the solutions for given initial values Substituting Eq. (15) into Eqs. (11), (12) and solving for (d/dz)c,(z), defining the determinants
c,(Zo). 4. Solutions for a local potential
At first a set of solutions ~O°(z), X°(z) is derived for the system of Eqs. (11), (12) without the non-local correction term. We adopt here the set of solutions which are implicitly used in the so-called omega-matrix method [6]. In this method an arbitrary solution can be written, in matrix notation, as
x~(z)
Z°lx°(zo) '
(13)
where (O°(Zo),)(m0(Z0)) T (T represents the transpose) are the values of the solution at z = z0, assumed to be known, and g2Zz0 is a 2 N × 2 N matrix having elements which are functions of z and z o. The set of 2 N solutions, (qt°,(z), X°m,(Z))T (v = 1 - 2 N ) is represented by the vth column of
A(z)
= det a Zo z
(16)
and A~m(z) as the co-factor of the element in A(z), we obtain - -
d c.(z) = Y'.F.L(z)J dz L
"bl
~
~ ~GL.,(z
r
X°,(z)
)c.,(z') dz',
aL v'
(17) where
U aCm(z)
F.L(Z)= ~_,--UmL(Z)
(18)
,.=1 a ( z )
and N
oL~,( z') = E W,~n(z ' )~0..,(z'). 0 n=l
(19)
K. Kambe, C. Stampfl / Ultramicroscopy 55 (1994) 221-227
224
Integration of Eq. (17) gives cv(z) =cv(zo)
6. Correction matrix
+ ~_~H,,L(Z, Zo)
L
X~bL~,GL,,,(Z')C,~,(Z ') dz',
(20)
where
In the case of non-local potentials we evaluate a matrix Z{., corresponding to ~ zZ(I in Eq. (13), defined by
Xm( Z )
HvL( Z, Zo) = fz,IFvL( Z ) dz.
(21)
It is noted that the equivalence of Eq. (20) to integral equations with a degenerate kernel can be seen by putting gc~(z) = c~(z) - c~(z o) in Eq. (20), giving a set of Fredholm's integral equations of the second kind,
=
z, Zo) EJuc L
Just as for J2Zo z the/zth column of Z z~o represents the/zth solution (~b.w(z), Xm~(Z)) T which has the same initial conditions as for (q~O (z) ' X m) u (0Z ) T (V being replaced by p.). We apply to each of the solutions (0mu(z), Xmu(Z))T the method of variation of constants described above, putting
,( z0)
v'
+
Xm.(Z)
Zo) L
tX%(Z)
where bl
SLy, = £LGL~,( Z) dz.
(23)
By applying the standard method of solution [13,14] for these types of integral equations we arrive at
c~(z) =c~(z0) + ~ K d ( z ,
Zo)C~,(Zo),
(24)
1/
where D~L' .
K,,,,,(z, Zo) = ~_, ~_,HvL(Z, Zo)---D--JL,~,, L
(25)
L'
where D is the determinant 1 --111
--I12
"'"
--I21
1 -- 122
...
=
,
(26)
where bl
'
the initial condition being given by c~.(z0) = 6v..
D
(2S)
z,,t Xm( ZO) .
ILC: ~ £LGL~(Z)H~c(z, z0) dz
(27)
and D~L, is the co-factor of the L'Lth element of D.
(30)
Introducing the correction matrix C~o having c~,(z) as its elements, we may write Eq. (29) as Z[,, = n~,CZ , (31)
and Eq. (30) as C ZO ~'' = E,
(32)
E being the unit matrix. We assume that the non-local correction given in Eq. (2) has a non-vanishing value only in small regions inside the atomic cores. Correspondingly, the functions U,,L(Z), WLn(Z') appearing in Eq. (6) may have non-vanishing values only in small ranges of z = (a L, bL). Outside these ranges of z the solution may well be represented by using omega-matrices. We consider here only the simplest case of one atomic layer, extension to more layers can be easily done. The correction appears only in one range (a, b) including all ranges (aL, b L) for that layer, and lying inside a certain range of z, let us say one lattice period (0, c). The solution is derived successively in the regions (0, a), (a, b), (b, c). Obviously we have Z~ = o c 7 b o a -- o c obtQbc')a ~Ub~a~O
-- ~t, b a ~ a V a~,() ,
(33)
K. Kambe, C. Stampfl / Ultramicroscopy 55 (1994) 221-227
using Eq. (31). Matrix G b has the elements cry(b), given from Eq. (24) by
c . . ( b ) = 6~. + K~u(b, a).
(34)
The evaluation of Gab is most conveniently carried out in parallel to the evaluation of ~ b (by using discrete steps of z in the multi-slice scheme). That is, the non-local corrections are included as an additional procedure during and after passage through a core region. The integrals of Eqs. (21), (23) and (27) may be evaluated by accumulating the integrands in the course of evaluating J2~. The expressions ACm(Z)/A(z) in Eq. (18) and D ~ c / D in Eq. (25) may be evaluated by Gauss' algorithm. As can be seen clearly from Eq. (33) the effect of the non-local correction is represented as the insertion of Gab at the starting point z = a of the range (a, b) in which the deviation of the solutions takes place. Phenomenologically the effect is that a rearrangement of the initial values occurs at z = a, otherwise the calculation follows the omega-matrix procedure.
7. Discussion
The present formulation gives exact solutions of the problem, except for the use of a finite number of beams and finite steps of z in the multi-slice method. Various approximative methods may be conceived. For example, Eqs. (11) and (12), or (20) may be solved by simple iteration, beginning with the solutions (qJ°~(z), X°v(z)) T, or c ~ ( z ) = c,(zo), respectively. Alternatively, using the discrete values of z used for the multi-slice method, we may convert Eqs. (11), (12) or (17) into algebraic equations, replacing derivatives by differences and integrals by sums, and solve them by means of numerical programs developed for large matrices. In these methods the non-local term essentially needs not be separable, but a separable form may require much less calculation time. As mentioned already, the present formulation has been developed for electrons in the energy range of valence bands. The non-local pseudopo-. tentials [8-12] are considered only for atomic
225
cores. The electron under consideration sees the pseudopotential of atomic cores with the non-local correction term and the potential of other valence electrons. The latter emerges from the results of self-consistent calculations and is usually regarded as local. The correction term of Kleinman and Bylander [9] has essentially the form
UL( F -- ri ) = CilXil( r ) Yim( O, ~9),
(35)
WL(r -- ri) = X ~ ( r ) Y f f ( O, q~),
(36)
where r, 0, and , are the polar coordinates of r - ri, Cil is a constant, X a ( r ) is an l dependent radial part, and Ylm(O, ~O) are spherical harmonics. Usually l only up to l = 3 is necessary [8]. Obviously the method is immediately applicable to electron-diffraction and -microscopy problems in the energy range of band theory, e.g. for the initial states of photoemission electron microscopy [15], and for the final states of scanning tunnelling microscopy [16]. The advantage lies in the possibility of taking over the results of selfconsistent calculations using non-local pseudopotentials without changing the model potentials. In the energy ranges of L E E D and H E E D the form of the effective potential becomes different, so that the data derived for valence electrons cannot be used. Furthermore, due to the increase of the necessary number of lm's with increasing energy, the expansion in spherical harmonics becomes more and more ineffective. A possible way around this may be to construct a new pseudopotential. It has been pointed out by B16chl [10] and Vanderbilt [11] that the separable non-local forms of pseudopotentials can be derived in quite general cases. Apparently it may be regarded as one of the generally applicable methods of deriving pseudopotentials at all. We may try to apply their method, in the context of the multi-slice method, for example, using the coordinates rt(= r sin 0), q~, z in place of r,O,~ around one atom in the following way. At first we need an exact solution of the scattering of an electron at a given energy by a single atom. This may be done, for example, by taking a sufficiently large number of partial waves. It is assumed here that the solution exists. Assuming
K. Kambe, C. Stampfl/ Ultramicroscopy55 (1994) 221-227
226
that the electron is incoming in the positive z-direction, we separate from the wavefunction O(r) the rapidly oscillating factor exp(ik0~z), with a properly chosen reference wavevector k0z, into the form
~b(r) = w ( r ) exp(ik0~z ).
(37)
A small range z = (a, b) is chosen around the atomic core and divided into very thin slices z =(zs, z i + Azi). The average values of the atomic potential Va(r) is constructed inside the slices
1
zi+azi 2m vi(rt, q ~ ) = ~ z / f z , h2 Va(rt, q~,z) d z
(38)
and similarly wi(rt, q~) is derived from w(r). The partial derivatives of w(r) in z may be replaced by differences, for example, Ow
0z
1 - --½(ws+,-
W,_l),
(az )
632W
D{w,} - viw i = 0,
(40)
(41)
where D is the differential-difference operator acting on w i r t ~t2 -Or,
- -
Ort
(Az,) 2
(Wi+l--2Wiq-Wi t)
8.= ffx? ir,
d r t dq~.
(45)
Then the non-local correction term is given by
uNLwi = E x k ( r t , @)( B-1)kj kj
× f fx7( r;, q~')wi(r;,
q~')r; dr~ d~', (46)
ik0~ (42) where (43)
(47)
We see that the potential c i in Eq. (41) is replaced by the pseudopotential v/es = t i,local _[_ I!NL in Eq. (47), and that ~bi's are solutions of Eq. (47). For going back to continuous z we define t:l°cal(r), x(r), and u(r) in such a way that their z-averages in the sense of Eq. (38) above are v] °cal, Xi, and S, kXk(B l)kiA Z i, respectively. Then the integral
l)kjX*(r ~, q£)
= £l'u(rt, q~, z ' ) x * ( r ~ , qY, z") dz"
q- AZ--~(Wi+I--Wi_,) -}-(~:-k2z)Wi,
2m e = ~-E.
and evaluate the matrix
E X k ( r t , q~)( g kj
1
q-
(44)
D{ wi} - v]°Calwi - uNLwi = 0.
( Azi)z (wi+l -- 2wi + wi_t).
=
Xi( rt, q~) = D{ q~i} - b']°cal~i
leading to
The Schr6dinger equation is then written as
D{wi}
(The smoothening should occur with respect to rt, q: and i.) Generalizing slightly the method of B16chl [10] and Vanderbilt [11] to be applied to Eq. (41), where wi's are coupled and t,]°cal depends on i, we construct
(39)
1
OZz
local pseudopotential U)°cal (independent of Oi)-
(48)
may be satisfactorily approximated by a summation formula taking a much smaller number of discrete points z~' than z i above. It follows then from Eq. (47) that V2O q- EI~ -- ul°cal(r)O - exp(ik0~z ) ~ u ( r t ,
q~, z';)Az'/
l From the existing exact solution of the scattering problem the solution wi can be derived. We replace the solution by a smoothened pseudo-wavefunction ~i, and the potential v~ by a smoothened
× f x * ( r ~ , q£, z'l' ) e x p ( - i k o ~ z ' ) ~ b ( r ' ) d3r '. (49)
K. Kambe, C. Stampfl / Ultramicroscopy 55 (1994) 221-227
A f t e r adding up the potential terms over atoms, this e q u a t i o n leads to an equation having the form of Eq. (1) with Eq. (2). (Two q u a n t u m n u m b e r s replacing lm in Eq. (2) are included, for simplicity, in the indices i and l above.) T h e local pseudopotential U]°cal and the pseudo-wavefunction ~bi may be chosen in such a way that the summation in Eq. (2) contains only a limited n u m b e r of terms. Having b r o u g h t the Schr6dinger equation into the form of Eqs. (1), (2), with an effective non-local pseudopotential suitable for energies in the ranges of L E E D and H E E D , the m e t h o d of solution using the multi-slice scheme follows exactly as outlined above and summarized in Eq. (33). T h e p r o c e d u r e given above may a p p e a r to be rather complicated, but it is r e g a r d e d only as an example to d e m o n s t r a t e that the flexibility allowed in construction of the present types of pseudopotentials opens new possibilities of introducing them into problems in the L E E D and H E E D energy ranges. Finally, in the special cases of normal incidence in H E E D and higher energy ranges (highvoltage electron microscopy, electron channeling) [17] it is known that only the potential averaged in the z-direction plays a role. It is also known that b o u n d states may a p p e a r in two-dimensional core regions. These states may be treated separately by using the tight-binding m e t h o d [18], and used for deriving two-dimensional O P W - t y p e pseudopotentials (see, e.g., Ref. [9]) in the field of the z-averaged potential.
8. Conclusions It is shown that the multi-slice m e t h o d can be modified to include pseudopotentials in the separable non-local form which is an important imp r o v e m e n t over the f o r m e r r e q u i r e m e n t that the potential be local. In the framework of the S2-matrix m e t h o d only a correction matrix is necessary to be introduced for the ranges of z in which non-local correction terms appear. In view of this
227
extension of the applicability of the multi-slice m e t h o d it is d e m o n s t r a t e d that, by the help of new developments in the theory of p s e u d o p o t e n tials, the present types of pseudopotentials may be constructed for problems in the usual energy ranges of electron diffraction and -microscopy in o r d e r to save the computational time avoiding the necessity of using the so-called muffin-tin model. A m o r e detailed description of the theory will be reported in a later p a p e r with numerical exampies.
References [1] J.M. Cowley and A.F. Moodie, Acta Cryst. 10 (1957) 609. [2] D.F. Lynch and A.F. Moodie, Surf. Sci. 32 (1972) 422. [3] A.E. Smith and D.F. Lynch, Surf. Sci. 184 (1986) 189. [4] A.E. Smith, J. Vac. Sci. Technol. 5 (1987) 1262. [5] T.C. Zhao, H.C. Poon and S.Y. Tong, Phys. Rev. B 38 (1988) 1172; P.A. Maksym and J.L. Beeby, Surf. Sci. 110 (1981) 423; G. Meyer-Ehmsen, Surf. Sci. 219 (1989) 177; A. Ichimiya, Jpn. J. Appl. Phys. 22 (1983) 176. [6] C. Stampfl, K. Kambe, J.D. Riley and D.F. Lynch, J. Phys. (Cond. Matter) 4 (1992) 8461. [7] C. Stampfl, K. Kambe, J.D. Riley and D.F. Lynch, J. Phys. (Cond. Matter) 5 (1993) 8211. [8] G.B. Bachelet, D.R. Hamann and M. Schliiter, Phys. Rev. B 26 (1982) 4199. [9] L. Kleinman and D.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. [10] P.E. Bl6chl, Phys. Rev. B 41 (1990) 5414. [11] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892. [12] X. Gonze, R. Stumpf and M. Scheffler, Phys. Rev. B 44 (1991) 8503. [13] R. Courant and D. Hilbert, Methods of Mathematical Physics, Vol. I (Interscience, New York, 1953) p. 114. [14] P. Morse and H. Feshbach, Methods of Theoretical Physics, Vol. I (McGraw-Hill, New York, 1953) p. 956. [15] O.H. Griffith and W. Engel, Eds., Emission Microscopy and Related Techniques, Special Issue, Ultramicroscopy 36 (1991) 1-274. [16] R. Wiesendanger and H.-J. Giintherodt, Eds., Scanning Tunneling Microscopy III, Vol. 29 of Springer Series in Surface Sciences, (Springer, Berlin, 1993). [17] M.V. Berry, B.F. Buxton and A.M. Ozorio de Almeida, Rad. Eft. 20 (1973) 1. [18] K. Kambe, G. Lehmpfuhl and F. Fujimoto, Z. Naturforsch. 29a (1974) 1034.
Ultramicroscopy 55 (1994) 221-227
Non-local pseudopotentials in the multi-slice method of calculating electron wavefunctions in crystals Kyozaburo Kambe, Catherine Stampfl Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, D-14195 Berlin (Dahlem), Germany Received 3 December 1993 Dedicated to Gunter Lehmpfuhl on the occasion of his 65th birthday
Abstract
A method for including non-local potentials in the multi-slice formalism for calculating electron wavefunctions in crystals is presented. With the assumption that non-local potentials appear only in atomic core regions as correction terms having a so-called separable form, it is shown that these corrections can be included as an additional procedure carried out after each passage through a core region. Possible ways of deriving pseudopotentials of the present type for various fields of electron diffraction and electron microscopy are briefly discussed.
1. Introduction
After the multi-slice method was introduced by Cowley and Moodie [1], it has been widely used for theoretical interpretation of electron-diffraction patterns and electron micrographs. The range of application has been expanded to the low-energy regime, and thus to low-energy electron diffraction [2] and current image diffraction [3]. The method has been used successfully to calculate convergent b e a m [4] and usual [5] reflection high-energy electron diffraction intensities, and was shown recently by Stampfl et al. [6] to be well applicable to calculations of usual bulk band structures and surface electronic states, and consequently to, for example, angle-resolved ultraviolet photoemission spectroscopy [7]. The applicability of the method to these last cases relies on the use of pseudopotentials, which enables us to
choose relatively small numbers of beams and slices per unit lattice period. Pseudopotentials have been implemented in band theory to replace true atomic potentials for the purpose of keeping calculation time manageable without changing the obtained results. Recent developments in the theory of pseudopotentials [8], however, tend to have non-local forms (see below) which cannot be used directly in the usual multi-slice method. In order to utilize fully the advantage of the method, in particular for surface problems, it is necessary to modify the method to incorporate non-local potentials. In the present p a p e r we report the method to include non-local corrections which has been developed for a special type of non-local pseudopotential [9-12] which is currently in use in self-consistent band-structure calculations. Hence, the electron is assumed to lie in the energy range of
0304-3991/94/$07.00 © 1994 Elsevier Science B.V. All rights reserved SSDI 0 3 0 4 - 3 9 9 1 ( 9 4 ) 0 0 0 5 6 - S
222
K. Kambe, C. Stampfl / Ultramicroscopy 55 (1994) 221-227
valence bands. It is believed, however, that similar types of pseudopotentials can be constructed appropriately for the use in low- and high-energy electron-diffraction and -microscopy problems when a reduction of calculation time is desired. We discuss these possibilities in Section 7.
2. Form of the pseudopotential
After application of the method of Kleinman and Bylander [9] (see also Refs. [10,11]) we may choose the pseudopotential to give, in the Schr6dinger equation 2m V2~+ 7(E-
V)~=0,
(1)
the potential term in the form
V(r)&(r) = vl°cal(r)&(r) + EUilm(r-ri) ilm
xfWam(r'-ri)O(r' )
3. Two-dimensional expansion of the Schr6dinger equation
For obtaining wavefunctions for crystal slabs or crystals with surfaces we calculate at first those solutions of the Schr6dinger equation (1) for which the values of energy E and the component of the wavevector parallel to the surface kll are given. It is convenient to expand 0 ( r ) in a two-dimensional Fourier series parallel to the surface
O(r)
= El]Ira(Z)
(4)
where r t is the tangential component of r to the surface, z is taken normal to the surface and is positively increasing into the crystal, m represents the two-dimensional Miller indices (h, k), and kmt is the tangential component of the wavevector k to the surface given by (5)
kmt = kot + bmt ,
d3r '.
(2)
VL°Cat(r) is the main part of the potential and is local; the second term is the correction having a "separable" non-local form. r i represents the position of the ith atom center. Expressions for the functions g i l m ( r - r i) and W i l m ( r ' - - r i) are given in Section 7. For brevity we may write
where k0t = kll , and brat is a two-dimensional reciprocal lattice vector. Substitution of Eq. (4) into Eq. (1) gives d2
dz2Om(Z)
+ k 2 m z t ~ m ( z ) --
EVm_.(Z)O.(Z) n
- ~_. ~UmL(Z)fbLWL~(Z')~O~(Z ') n
E Uilm(r - ri)Wilm(r'- ri) = ~-, Ul_(r)Wl_(r'), ilm
exp(ikmt "rt),
m
dz' = O,
~al•
L
(6)
L
(3) L representing ilm. The solution may be obtained at first by using V ~°ca' for V, and the non-local correction term may be introduced for modifying the solution. It is to be noted that the non-local correction term in Eq. (2) is represented in a so-called "degenerate" (Courant-Hilbert [13]) or "finite" (MorseFeshbach [14]) form, that is, a sum of a finite number of terms of a product of functions of r and r' only. In the theory of integral equations it is known that an integral equation having a kernel in a degenerate form can be solved relatively easily. We are going to derive these types of integral equations by applying the method of variation of constants.
where k2z --
2m h-----TE-]kmt[ 2,
(7)
and Um( Z ) -- _2m _
1 f vl°cal(r) e x p ( - i b m t h 2 A "A
"rt)
d2rt,
(8) 1 fAUL(r ) WL"(Z')
e x p ( - i b m t "rt) d2rt,
(9)
2m 1
h 2 A fA WL(r')
e x p ( i b n t "r't) d2r't'
(10)
K. Kambe, C. Stampfl / Uhramicroscopy 55 (1994) 221-227
where A indicates the area of the two-dimensional unit cell. The ranges of integrals (a L, b L) in Eq. (6) are taken to include the ranges of non-vanishing values of WLn(Z'), and may depend on L, particularly on i through the position vector r i. In Eq. (6) the summation over i includes only the atoms in a two-dimensional unit cell. It is convenient to reduce Eq. (6) to a system of first-order differential equations by putting
223
/2z~0. We see that the expression (13) can be interpreted as a linear combination of these solutions
xO(z)
=
E 6°(Zo) 2N
+
%(z))" E x°_~(Zo)(~'x%(z)
v=N+l
d
dz~bm(Z)
= Xm(Z),
(14)
(11)
and writing Eq. (6) as d
dzX,,( z) + k2mz$m(z)
5. Variation of constants - 2 U m _ n ( Z)~ln( Z ) n
- ~n Y'~UmL(Z)f~:
dz'=O. (12)
We assume that it is sufficient to take into account only N equations of Eq. (6), that is, 2 N equations of the system (11), (12).
Using the solutions (O°,(z), x O v ( g ) ) T, which are assumed here to be already obtained, we are going to solve Eqs. (11), (12) by using the method of variation of constants. We put
( )2N
= E c.(z)[ ~'%(z)
Xm(Z)
.=l
lxOv(Z)
)
(15)
and iook for the solutions for given initial values Substituting Eq. (15) into Eqs. (11), (12) and solving for (d/dz)c,(z), defining the determinants
c,(Zo). 4. Solutions for a local potential
At first a set of solutions ~O°(z), X°(z) is derived for the system of Eqs. (11), (12) without the non-local correction term. We adopt here the set of solutions which are implicitly used in the so-called omega-matrix method [6]. In this method an arbitrary solution can be written, in matrix notation, as
x~(z)
Z°lx°(zo) '
(13)
where (O°(Zo),)(m0(Z0)) T (T represents the transpose) are the values of the solution at z = z0, assumed to be known, and g2Zz0 is a 2 N × 2 N matrix having elements which are functions of z and z o. The set of 2 N solutions, (qt°,(z), X°m,(Z))T (v = 1 - 2 N ) is represented by the vth column of
A(z)
= det a Zo z
(16)
and A~m(z) as the co-factor of the element in A(z), we obtain - -
d c.(z) = Y'.F.L(z)J dz L
"bl
~
~ ~GL.,(z
r
X°,(z)
)c.,(z') dz',
aL v'
(17) where
U aCm(z)
F.L(Z)= ~_,--UmL(Z)
(18)
,.=1 a ( z )
and N
oL~,( z') = E W,~n(z ' )~0..,(z'). 0 n=l
(19)
K. Kambe, C. Stampfl / Ultramicroscopy 55 (1994) 221-227
224
Integration of Eq. (17) gives cv(z) =cv(zo)
6. Correction matrix
+ ~_~H,,L(Z, Zo)
L
X~bL~,GL,,,(Z')C,~,(Z ') dz',
(20)
where
In the case of non-local potentials we evaluate a matrix Z{., corresponding to ~ zZ(I in Eq. (13), defined by
Xm( Z )
HvL( Z, Zo) = fz,IFvL( Z ) dz.
(21)
It is noted that the equivalence of Eq. (20) to integral equations with a degenerate kernel can be seen by putting gc~(z) = c~(z) - c~(z o) in Eq. (20), giving a set of Fredholm's integral equations of the second kind,
=
z, Zo) EJuc L
Just as for J2Zo z the/zth column of Z z~o represents the/zth solution (~b.w(z), Xm~(Z)) T which has the same initial conditions as for (q~O (z) ' X m) u (0Z ) T (V being replaced by p.). We apply to each of the solutions (0mu(z), Xmu(Z))T the method of variation of constants described above, putting
,( z0)
v'
+
Xm.(Z)
Zo) L
tX%(Z)
where bl
SLy, = £LGL~,( Z) dz.
(23)
By applying the standard method of solution [13,14] for these types of integral equations we arrive at
c~(z) =c~(z0) + ~ K d ( z ,
Zo)C~,(Zo),
(24)
1/
where D~L' .
K,,,,,(z, Zo) = ~_, ~_,HvL(Z, Zo)---D--JL,~,, L
(25)
L'
where D is the determinant 1 --111
--I12
"'"
--I21
1 -- 122
...
=
,
(26)
where bl
'
the initial condition being given by c~.(z0) = 6v..
D
(2S)
z,,t Xm( ZO) .
ILC: ~ £LGL~(Z)H~c(z, z0) dz
(27)
and D~L, is the co-factor of the L'Lth element of D.
(30)
Introducing the correction matrix C~o having c~,(z) as its elements, we may write Eq. (29) as Z[,, = n~,CZ , (31)
and Eq. (30) as C ZO ~'' = E,
(32)
E being the unit matrix. We assume that the non-local correction given in Eq. (2) has a non-vanishing value only in small regions inside the atomic cores. Correspondingly, the functions U,,L(Z), WLn(Z') appearing in Eq. (6) may have non-vanishing values only in small ranges of z = (a L, bL). Outside these ranges of z the solution may well be represented by using omega-matrices. We consider here only the simplest case of one atomic layer, extension to more layers can be easily done. The correction appears only in one range (a, b) including all ranges (aL, b L) for that layer, and lying inside a certain range of z, let us say one lattice period (0, c). The solution is derived successively in the regions (0, a), (a, b), (b, c). Obviously we have Z~ = o c 7 b o a -- o c obtQbc')a ~Ub~a~O
-- ~t, b a ~ a V a~,() ,
(33)
K. Kambe, C. Stampfl / Ultramicroscopy 55 (1994) 221-227
using Eq. (31). Matrix G b has the elements cry(b), given from Eq. (24) by
c . . ( b ) = 6~. + K~u(b, a).
(34)
The evaluation of Gab is most conveniently carried out in parallel to the evaluation of ~ b (by using discrete steps of z in the multi-slice scheme). That is, the non-local corrections are included as an additional procedure during and after passage through a core region. The integrals of Eqs. (21), (23) and (27) may be evaluated by accumulating the integrands in the course of evaluating J2~. The expressions ACm(Z)/A(z) in Eq. (18) and D ~ c / D in Eq. (25) may be evaluated by Gauss' algorithm. As can be seen clearly from Eq. (33) the effect of the non-local correction is represented as the insertion of Gab at the starting point z = a of the range (a, b) in which the deviation of the solutions takes place. Phenomenologically the effect is that a rearrangement of the initial values occurs at z = a, otherwise the calculation follows the omega-matrix procedure.
7. Discussion
The present formulation gives exact solutions of the problem, except for the use of a finite number of beams and finite steps of z in the multi-slice method. Various approximative methods may be conceived. For example, Eqs. (11) and (12), or (20) may be solved by simple iteration, beginning with the solutions (qJ°~(z), X°v(z)) T, or c ~ ( z ) = c,(zo), respectively. Alternatively, using the discrete values of z used for the multi-slice method, we may convert Eqs. (11), (12) or (17) into algebraic equations, replacing derivatives by differences and integrals by sums, and solve them by means of numerical programs developed for large matrices. In these methods the non-local term essentially needs not be separable, but a separable form may require much less calculation time. As mentioned already, the present formulation has been developed for electrons in the energy range of valence bands. The non-local pseudopo-. tentials [8-12] are considered only for atomic
225
cores. The electron under consideration sees the pseudopotential of atomic cores with the non-local correction term and the potential of other valence electrons. The latter emerges from the results of self-consistent calculations and is usually regarded as local. The correction term of Kleinman and Bylander [9] has essentially the form
UL( F -- ri ) = CilXil( r ) Yim( O, ~9),
(35)
WL(r -- ri) = X ~ ( r ) Y f f ( O, q~),
(36)
where r, 0, and , are the polar coordinates of r - ri, Cil is a constant, X a ( r ) is an l dependent radial part, and Ylm(O, ~O) are spherical harmonics. Usually l only up to l = 3 is necessary [8]. Obviously the method is immediately applicable to electron-diffraction and -microscopy problems in the energy range of band theory, e.g. for the initial states of photoemission electron microscopy [15], and for the final states of scanning tunnelling microscopy [16]. The advantage lies in the possibility of taking over the results of selfconsistent calculations using non-local pseudopotentials without changing the model potentials. In the energy ranges of L E E D and H E E D the form of the effective potential becomes different, so that the data derived for valence electrons cannot be used. Furthermore, due to the increase of the necessary number of lm's with increasing energy, the expansion in spherical harmonics becomes more and more ineffective. A possible way around this may be to construct a new pseudopotential. It has been pointed out by B16chl [10] and Vanderbilt [11] that the separable non-local forms of pseudopotentials can be derived in quite general cases. Apparently it may be regarded as one of the generally applicable methods of deriving pseudopotentials at all. We may try to apply their method, in the context of the multi-slice method, for example, using the coordinates rt(= r sin 0), q~, z in place of r,O,~ around one atom in the following way. At first we need an exact solution of the scattering of an electron at a given energy by a single atom. This may be done, for example, by taking a sufficiently large number of partial waves. It is assumed here that the solution exists. Assuming
K. Kambe, C. Stampfl/ Ultramicroscopy55 (1994) 221-227
226
that the electron is incoming in the positive z-direction, we separate from the wavefunction O(r) the rapidly oscillating factor exp(ik0~z), with a properly chosen reference wavevector k0z, into the form
~b(r) = w ( r ) exp(ik0~z ).
(37)
A small range z = (a, b) is chosen around the atomic core and divided into very thin slices z =(zs, z i + Azi). The average values of the atomic potential Va(r) is constructed inside the slices
1
zi+azi 2m vi(rt, q ~ ) = ~ z / f z , h2 Va(rt, q~,z) d z
(38)
and similarly wi(rt, q~) is derived from w(r). The partial derivatives of w(r) in z may be replaced by differences, for example, Ow
0z
1 - --½(ws+,-
W,_l),
(az )
632W
D{w,} - viw i = 0,
(40)
(41)
where D is the differential-difference operator acting on w i r t ~t2 -Or,
- -
Ort
(Az,) 2
(Wi+l--2Wiq-Wi t)
8.= ffx? ir,
d r t dq~.
(45)
Then the non-local correction term is given by
uNLwi = E x k ( r t , @)( B-1)kj kj
× f fx7( r;, q~')wi(r;,
q~')r; dr~ d~', (46)
ik0~ (42) where (43)
(47)
We see that the potential c i in Eq. (41) is replaced by the pseudopotential v/es = t i,local _[_ I!NL in Eq. (47), and that ~bi's are solutions of Eq. (47). For going back to continuous z we define t:l°cal(r), x(r), and u(r) in such a way that their z-averages in the sense of Eq. (38) above are v] °cal, Xi, and S, kXk(B l)kiA Z i, respectively. Then the integral
l)kjX*(r ~, q£)
= £l'u(rt, q~, z ' ) x * ( r ~ , qY, z") dz"
q- AZ--~(Wi+I--Wi_,) -}-(~:-k2z)Wi,
2m e = ~-E.
and evaluate the matrix
E X k ( r t , q~)( g kj
1
q-
(44)
D{ wi} - v]°Calwi - uNLwi = 0.
( Azi)z (wi+l -- 2wi + wi_t).
=
Xi( rt, q~) = D{ q~i} - b']°cal~i
leading to
The Schr6dinger equation is then written as
D{wi}
(The smoothening should occur with respect to rt, q: and i.) Generalizing slightly the method of B16chl [10] and Vanderbilt [11] to be applied to Eq. (41), where wi's are coupled and t,]°cal depends on i, we construct
(39)
1
OZz
local pseudopotential U)°cal (independent of Oi)-
(48)
may be satisfactorily approximated by a summation formula taking a much smaller number of discrete points z~' than z i above. It follows then from Eq. (47) that V2O q- EI~ -- ul°cal(r)O - exp(ik0~z ) ~ u ( r t ,
q~, z';)Az'/
l From the existing exact solution of the scattering problem the solution wi can be derived. We replace the solution by a smoothened pseudo-wavefunction ~i, and the potential v~ by a smoothened
× f x * ( r ~ , q£, z'l' ) e x p ( - i k o ~ z ' ) ~ b ( r ' ) d3r '. (49)
K. Kambe, C. Stampfl / Ultramicroscopy 55 (1994) 221-227
A f t e r adding up the potential terms over atoms, this e q u a t i o n leads to an equation having the form of Eq. (1) with Eq. (2). (Two q u a n t u m n u m b e r s replacing lm in Eq. (2) are included, for simplicity, in the indices i and l above.) T h e local pseudopotential U]°cal and the pseudo-wavefunction ~bi may be chosen in such a way that the summation in Eq. (2) contains only a limited n u m b e r of terms. Having b r o u g h t the Schr6dinger equation into the form of Eqs. (1), (2), with an effective non-local pseudopotential suitable for energies in the ranges of L E E D and H E E D , the m e t h o d of solution using the multi-slice scheme follows exactly as outlined above and summarized in Eq. (33). T h e p r o c e d u r e given above may a p p e a r to be rather complicated, but it is r e g a r d e d only as an example to d e m o n s t r a t e that the flexibility allowed in construction of the present types of pseudopotentials opens new possibilities of introducing them into problems in the L E E D and H E E D energy ranges. Finally, in the special cases of normal incidence in H E E D and higher energy ranges (highvoltage electron microscopy, electron channeling) [17] it is known that only the potential averaged in the z-direction plays a role. It is also known that b o u n d states may a p p e a r in two-dimensional core regions. These states may be treated separately by using the tight-binding m e t h o d [18], and used for deriving two-dimensional O P W - t y p e pseudopotentials (see, e.g., Ref. [9]) in the field of the z-averaged potential.
8. Conclusions It is shown that the multi-slice m e t h o d can be modified to include pseudopotentials in the separable non-local form which is an important imp r o v e m e n t over the f o r m e r r e q u i r e m e n t that the potential be local. In the framework of the S2-matrix m e t h o d only a correction matrix is necessary to be introduced for the ranges of z in which non-local correction terms appear. In view of this
227
extension of the applicability of the multi-slice m e t h o d it is d e m o n s t r a t e d that, by the help of new developments in the theory of p s e u d o p o t e n tials, the present types of pseudopotentials may be constructed for problems in the usual energy ranges of electron diffraction and -microscopy in o r d e r to save the computational time avoiding the necessity of using the so-called muffin-tin model. A m o r e detailed description of the theory will be reported in a later p a p e r with numerical exampies.
References [1] J.M. Cowley and A.F. Moodie, Acta Cryst. 10 (1957) 609. [2] D.F. Lynch and A.F. Moodie, Surf. Sci. 32 (1972) 422. [3] A.E. Smith and D.F. Lynch, Surf. Sci. 184 (1986) 189. [4] A.E. Smith, J. Vac. Sci. Technol. 5 (1987) 1262. [5] T.C. Zhao, H.C. Poon and S.Y. Tong, Phys. Rev. B 38 (1988) 1172; P.A. Maksym and J.L. Beeby, Surf. Sci. 110 (1981) 423; G. Meyer-Ehmsen, Surf. Sci. 219 (1989) 177; A. Ichimiya, Jpn. J. Appl. Phys. 22 (1983) 176. [6] C. Stampfl, K. Kambe, J.D. Riley and D.F. Lynch, J. Phys. (Cond. Matter) 4 (1992) 8461. [7] C. Stampfl, K. Kambe, J.D. Riley and D.F. Lynch, J. Phys. (Cond. Matter) 5 (1993) 8211. [8] G.B. Bachelet, D.R. Hamann and M. Schliiter, Phys. Rev. B 26 (1982) 4199. [9] L. Kleinman and D.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. [10] P.E. Bl6chl, Phys. Rev. B 41 (1990) 5414. [11] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892. [12] X. Gonze, R. Stumpf and M. Scheffler, Phys. Rev. B 44 (1991) 8503. [13] R. Courant and D. Hilbert, Methods of Mathematical Physics, Vol. I (Interscience, New York, 1953) p. 114. [14] P. Morse and H. Feshbach, Methods of Theoretical Physics, Vol. I (McGraw-Hill, New York, 1953) p. 956. [15] O.H. Griffith and W. Engel, Eds., Emission Microscopy and Related Techniques, Special Issue, Ultramicroscopy 36 (1991) 1-274. [16] R. Wiesendanger and H.-J. Giintherodt, Eds., Scanning Tunneling Microscopy III, Vol. 29 of Springer Series in Surface Sciences, (Springer, Berlin, 1993). [17] M.V. Berry, B.F. Buxton and A.M. Ozorio de Almeida, Rad. Eft. 20 (1973) 1. [18] K. Kambe, G. Lehmpfuhl and F. Fujimoto, Z. Naturforsch. 29a (1974) 1034.