Energy of step formation on metal surfaces from the stabilized-jellium model

ARTICLE IN PRESS

Physica B 369 (2005) 187–195 www.elsevier.com/locate/physb

Energy of step formation on metal surfaces from the stabilized-jellium model M. Farjam Faculty of Physics and Nuclear Sciences, Amirkabir University of Technology, P.O. Box 15875-4413, Tehran, Iran Received 15 May 2005; received in revised form 10 August 2005; accepted 12 August 2005

Abstract We describe stepped metal surfaces by the stabilized-jellium model, and calculate the step formation energy within the local density approximation of the density-functional theory. Using fourth-order gradient expansion for the kinetic energy, we determine the electron density variationally from a parameterized trial density. We compare our numerical results for a monatomic step on the Al(1 1 1) surface with other theoretical calculations and estimates. r 2005 Elsevier B.V. All rights reserved. PACS: 68.35.Dv; 68.35.Md; 71.15.Mb Keywords: Step formation energy; Metal surfaces; Stabilized-jellium model; Density-functional theory

1. Introduction Recently, regularly stepped vicinal surfaces have received considerable attention, partly due to their relevance to nanoscience [1–4]. In these efforts, theoretical electronic structure calculations can map the energetics of stepped surfaces [5], including step formation energy and step–step interaction, providing useful data for understanding the role of steps in surface processes. Various approaches exist, including the tight-binding Tel.: +98 21 6400603, fax: +98 21 6495519.

E-mail addresses: [email protected], [email protected].

model [6], effective medium theory [7], ab initio pseudopotential methods [8], and the jellium model. The jellium or the uniform background model in conjunction with density-functional theory was one of the early models to be used to calculate the surface properties of simple sp-bonded metals [9,10]. It affords a simple and nonempirical description of a wide number of metals, and with some pseudopotential corrections it gives reasonable values for the surface energy. The stabilizedjellium model incorporates a structureless pseudopotential into the jellium model [11], making it a more realistic model of metals, while retaining the simplicity and universality of the jellium model.

0921-4526/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2005.08.014

ARTICLE IN PRESS M. Farjam / Physica B 369 (2005) 187–195

188

The treatment of flat surfaces by the jellium model, using density-functional theory within local density approximation, both by the self-consistent Kohn–Sham method [10,12,13] and variational techniques [9,11,14], is extensive and well-established. Addition of adatoms, steps, or other defects to the jellium surface, makes the self-consistent equations much more difficult to solve. Thus, although Kohn–Sham calculations for jellium stepped surfaces have been reported [15], numerical values for the step formation energy were not given. In another set of calculations [16–18], the step formation energy was derived from a more general theory of curved surfaces. We present a variational calculation of the energy of an isolated step on a metal surface using the stabilized-jellium model. We have formulated our model similarly to the ‘simple analytic model’ of Perdew [11,14] for flat surface properties. Although the resulting model cannot be solved analytically, its numerical solution requires little computational time. We present a description of the model in Section 2, sample numerical results for Al(1 1 1) in Section 3, and concluding remarks in Section 4. We gather detailed expressions for numerical calculations of the step energy in Appendix A, and for the electrostatic potential in Appendix B.

2. The model 2.1. Total energy In the stabilized-jellium model, the total energy is a functional of the electron density nðrÞ, and can be written as [11,16] E½n; nþ  ¼ T s ½n þ E es ½n; nþ  þ E xc ½n þ E ps ½n; nþ . (1) Here, T s is the noninteracting kinetic energy given, to fourth-order in gradient expansion, by T s ½n ¼ T 0 ½n þ T 2 ½n þ T 4 ½n,

(2)

where 3 T 0 ½n ¼ ð3p2 Þ2=3 10

Z

3

d rn

5=3

,

(3)

1 T 2 ½n ¼ 72

Z

d3 r

jrnj2 , n

(4)

Z ð3p2 Þ 2=3 d3 rn1=3 T 4 ½n ¼ 540 (   ) 2    r2 n 9 r2 n rn2 1 rn4

þ   . n 8 n n 3 n ð5Þ The electrostatic energy E es is given by Z 1 d3 rfð½n; nþ ; rÞ½nðrÞ nþ ðrÞ, E es ½n; nþ  ¼ 2 where fð½n; nþ ; rÞ ¼

Z

d3 r0 ½nðr0 Þ nþ ðr0 Þ=jr0 rj

(6)

(7)

is the electrostatic potential. The exchange-correlation energy E xc is given, in the local-density approximation, by Z E xc ½n ¼ d3 rnðrÞexc ðnðrÞÞ, (8) where exc ðnÞ ¼ 3kF =4p þ ec ðnÞ.

(9)

(We use atomic units, _ ¼ e ¼ m ¼ 1:) Here, ec ðnÞ is the correlation energy per electron in a uniform gas of density n, for which we use the Vosko et al. expression [19]. The positive background nþ equals a constant, n¯ , inside a sharp surface, and vanishes outside. For a neutral metal, Z Z 3 d rnþ ðrÞ ¼ d3 rnðrÞ. (10) 2

So far the description is identical to that of the jellium model. The energy, E ps , is the stabilizedjellium, or structureless pseudopotential, contribution, and can be written as Z C d3 rnþ ðrÞ½nðrÞ nþ ðrÞ, E ps ½n; nþ  ¼ (11) n¯ where C ¼ k2F =5 þ kF =4p þ ðrs =3Þdec =drs .

(12)

We calculate an approximation to E½n; nþ  by variation of parameters of a suitable form of the electron density.

ARTICLE IN PRESS M. Farjam / Physica B 369 (2005) 187–195

2.2. Electron density Our system consists of two flat terraces with their outward normal in the positive z-direction separated by a single ledge (monatomic step) extending along the y-axis, as illustrated in Fig. 1. For the flat surface, the electron density profile is written as [14] nðzÞ ¼ n¯ f ðgks zÞ,

(13)

where g is a variational parameter of order unity and ks ¼ ð4kF =pÞ1=2 is the inverse of the Thomas– Fermi screening length. The form function that was used to calculate the surface energy in the stabilized-jellium model is [11] ( 1 þ B0 ez þ B1 eCz ; zo0; f ðzÞ ¼ (14) B2 e Dz ; z40; where B0 ¼ 0:621, B1 ¼ 0:085758, B2 ¼ 0:464758, C ¼ 2:98893 and D ¼ 0:784656. This form satisfies the charge neutrality constraint (10), and certain other requirements. Let us construct a parameterized electron density for the stepped surface. To do this, we consider the contours of constant electron density, which can be described by Zðz; xÞ ¼ C i ,

(15)

where fC i g is a set of arbitrary constant labels for the contours. As fC i g is not unique, any number of functions may satisfy Eq. (15), including the electron density itself. A particularly useful solution for our purpose is one that asymptotically approaches the z-coordinate relative to the terraces, or more precisely, Zðz; xÞ ! z  h=2

for x ! 1.

(16)

z h

x

Fig. 1. Step of height h. Origin of coordinate system is at the mark on the ledge.

189

The function Z that satisfies both Eqs. (15) and (16) is a natural coordinate for describing the electron density of our system, allowing us to write, nðz; xÞ ¼ n¯ f ðgks ZÞ,

(17)

where f is given by Eq. (14), i.e., the same function that was used for the flat surface. (We use n to denote both flat surface and stepped surface electron densities, however, the meaning should be clear from the context.) Our task is now reduced to seeking a suitable function for Zðz; xÞ. We expect the contours of constant electron density to be in the shape of rounded steps, and we ‘guess’ this function to be (1 x xo0; 2 ðe 1Þ; F ðxÞ ¼ 1 (18) x x40: 2 ð1 e Þ; This function tends to 12 as x ! 1, so a function Zðz; xÞ corresponding to our ledge can be written as Z ¼ z þ hF ½aðx bZÞ.

(19)

Here, a and b are the variational parameters used to optimize Z, and they serve to vary the steepness of the shape of the contour, and the placement of the contours in the zx-plane, respectively. To make the meaning of Eq. (19) clear, we can regard z as a function of x with Z being a constant; this describes a rounded step of height h, which is translated by Z along the z-axis and by bZ along the x-axis. The shape of the contours can be seen in Fig. 2(a), which is a contour plot of the electron density we obtained for Al(1 1 1). A similar idea, without b, i.e., b ¼ 0, was already used by Thompson and Huntington [20], in a study of adatom binding near surface ledges. Setting b ¼ 0 makes the computation of Z analytic, but intuitively we expect b to be closer to 1 rather than 0, and our numerical computations confirm this. With Z present on both sides of Eq. (19), however, we have to solve for it by numerical root finding methods [21]. It may not be obvious that the electron density thus defined satisfies the charge neutrality constraint described by Eq. (10). In fact, a mathematical property of the form of Z that we have chosen guarantees this. To show this property consider an

ARTICLE IN PRESS M. Farjam / Physica B 369 (2005) 187–195

190

ELECTRON DENSITY
2(z,x)

ELECTRON DENSITY ne(z,x)

5

ksz

ksz

5

0

0

-5

-5

-10

-5

(a)

0 ksx

5

10

-10

(b)

-5

0 ksx

5

10

Fig. 2. (a) Electron density nðz; xÞ defined by Eq. (17) for a ¼ 2:24 and b ¼ 0:79. The plot also illustrates the contours of constant Zðz; xÞ, Eq. (19). (b) Charge density r2 , the difference between nðz; xÞ and the ‘perfect’ electron density of each terrace. It shows a dipolar character localized around the ledge, with negative values on the left side of the ledge and positive values on the right. Its largest magnitude is about half of the bulk density.

algebraic function of the density, gðnÞ. Examples of interest here are gðnÞ ¼ n, gðnÞ ¼ n5=3 , and gðnÞ ¼ n4=3 , corresponding to the density itself, to the zeroth-order kinetic energy, and the exchange energy, respectively. We can write Z L Z 1 dx dzfg½nðz; xÞ nþ ðz; xÞg L 1   Z 1 Z L qz dx dZ fg½nðZÞ n¯ yð ZÞg, ¼ qZ x L 1 ð20Þ where L is a large number representing boundaries of the system in the x-direction. Here, we could replace nþ ðz; xÞ by n¯ yð ZÞ because they yield the same integral due to the symmetric form of F ðxÞ, Eq. (18). Noting that the integrand vanishes for Z ! 1, and with the Jacobian given by   qz ¼ 1 þ abhF 0 ½aðx bZÞ (21) qZ x and using Z 1 aF 0 ½aðx bZÞ dx ¼ 1,

(22)

1

we find Z L Z dx

1

dzfg½nðz; xÞ nþ ðz; xÞg L 1 Z 1 ¼ ð2L þ bhÞ dzfg½nðzÞ n¯ yð zÞg, 1

ð23Þ

where n¯ yð zÞ is the positive charge density of the flat surface. This establishes a simple relation between certain integrals of the stepped surface and the corresponding ones of the flat surface. For a charge neutral flat system the right-hand side of Eq. (23) vanishes, thus the derived stepped density also satisfies charge neutrality. For nonzero integrals, the term proportional to bh is the contribution due to the step formation. The relation expressed in Eq. (23), applies to T 0 , E x and E c , and can serve to check the numerical calculation of these terms. The terms T 2 and T 4 involve the gradient and Laplacian of the density, and E es is related to the density by the Poisson equation, and there is no such simple relation between their contribution to step formation energy and the flat surface terms. Neither does the above property apply to E ps , because the integral involved does not extend to the entire space. Although we can calculate the gradient and Laplacian of the density by numerical methods, it is more efficient computationally to use the following formulae to calculate them: jrnj2 ¼ ðgks Þ2 jrZj2 ðf 0 Þ2 ,

(24)

r2 n ¼ gks ðr2 ZÞf 0 þ ðgks Þ2 jrZj2 f 00 ,

(25)

ARTICLE IN PRESS M. Farjam / Physica B 369 (2005) 187–195

where the argument of f 0 and f 00 is gks Z. With Z given by Eq. (19), we have jrZj2 ¼

1 þ ðahF 0 Þ2 , ð1 abhF 0 Þ2

r2 Z ¼

(26)

a2 ð1 þ b2 ÞhF 00 , ð1 abhF 0 Þ3

(27)

where the argument of F 0 and F 00 is aðx bZÞ. Now, the evaluation of all energy terms, except for the electrostatic energy E es , is a matter of numerical integration. The electrostatic energy requires the calculation of the electrostatic potential first, which is related to the density by the Coulomb law, or, equivalently, by the Poisson equation. Thompson and Huntington [20] intro-

duced an efficient scheme for evaluating the electrostatic potential and energy appropriate for a stepped surface. We have used the same scheme, which we describe in outline here, giving expressions needed for calculations in Appendix B. In two dimensions Coulomb law takes the form ZZ fðrÞ ¼ 2 r dr dy ln jr r0 j½nðr0 Þ nþ ðr0 Þ. (28) First, we write the total charge density rðz; xÞ ¼ nðz; xÞ nþ ðz; xÞ as the sum rðz; xÞ ¼ r1 ðz; xÞ þ r2 ðz; xÞ, where r1 is ( r1 ðz; xÞ ¼

ELECTROSTATIC POTENTIAL φ1(Z,X)

ksz

ksz

xo0; x40:

(30)

5

0

0

-5

-5

-10

rs ðz h=2Þ; rs ðz þ h=2Þ;

(29)

ELECTROSTATIC POTENTIAL φ2(Z,X)

5

(a)

191

-5

0 ksx

5

10

-10

-5

0 k sx

(b)

5

10

ELECTROSTATIC POTENTIAL φ(Z,X)

ksz

5

0

-5

-10

(c)

-5

0 k sx

5

10

Fig. 3. (a) Electrostatic potential f1 , resulting from ‘perfect’ electron density for each terrace joined at the ledge. The potential is continuous, but its structure reflects the discontinuity of the charge at the ledge. (b) Potential f2 resulting from r2 , Fig. 2(b), showing the field due to a dipolar charge distribution. (c) The total potential f ¼ f1 þ f2 , showing expected contours in the shape of rounded steps.

ARTICLE IN PRESS 192

M. Farjam / Physica B 369 (2005) 187–195

Here, rs ðzÞ denotes the perfect total charge density of the flat surface. The difference r2 ¼ r r1 is a dipolar charge density distribution localized around the step (see Fig. 2(b)). By superposition law, the electrostatic potential is the sum of two contributions due to r1 and r2 , fðz; xÞ ¼ f1 ðz; xÞ þ f2 ðz; xÞ.

(31)

All three potentials from our numerical calculations are plotted in Fig. 3. Comparison with Fig. 2 shows that the potentials have the expected behaviour. It is convenient to write the electrostatic energy as Z Z 1 1 3 E es ¼ d r fr ¼ d3 r f1 r1 2 2 Z Z 1 þ d3 rf1 r2 þ ð32Þ d3 r f2 r2 . 2 The first term in Eq. (32) is constant with respect to the variational parameters and nearly cancels the electrostatic energy of the two terraces, while the second and third terms describe the interaction of the short-ranged r2 with r1 and with itself, respectively. The total energy of our system is a sum of bulk, surface and ledge energies, E ¼ EV þ sA þ lL,

(33)

where V, A and L are the volume of the system, the surface area of the two terraces and length of the ledge, respectively, and E, s and l denote bulk energy per unit volume, surface energy per unit area and ledge energy per unit length, respectively. We summarize appropriate expressions for calculating the various components of surface and ledge energies in Appendix A.

3. Results The stabilized-jellium model of the stepped surface is characterized by the metal density parameter, rs , defined by 4pr3s =3 ¼ 1=¯n, and the step height h. Here, we present the numerical results of applying this model to a monatomic step on Al(1 1 1), since we can compare it with other theoretical calculations [8,17,18]. In order to make the comparison

Fig. 4. Minimization of ledge energy.

meaningful, we use rs ¼ 2:04 corresponding to a lattice constant of a ¼ 7:52 bohr obtained in the references mentioned, although the experimental value of rs is 2.07. This p results in a monatomic ffiffiffi step height of h ¼ a= 3 ¼ 4:34 bohr, and a distancepffiffibetween atoms along the step line of ffi d ¼ a= 2 ¼ 5:32 bohr. We first minimize the surface energy of a flat surface with respect to the variational parameter g. Next, we minimize the ledge (step) energy with respect to the variational parameters a and b. A plot of the ledge energy in the neighbourhood of its minimum is shown in Fig. 4. The data is summarized in Table 1. Compared to 1.93 mhartree/bohr that we have obtained for the step energy, Stumpf and Scheffler [8], in detailed ab initio DFT-LDA calculations found 1.7 mhartree/bohr, and Ziesche et al. [18], from their theory of curved surfaces on stabilized-jellium, arrived at 1.8 mhartree/bohr. Considering that we have neglected atomistic detail, and that the variational calculation is only an upper bound to the exact energy of the model, the agreement with other theoretical calculations is surprisingly good. It is noteworthy to compare our value of surface energy, 0:526 mhartree=bohr2 or 819 erg=cm2 , with other theoretical calculations. Kohn–Sham calculations on the same stabilizedjellium model as ours have resulted in 925 erg=cm2 [13], while ab initio calculations of Ref. [8] have yielded 1120 erg=cm2 . Perdew et al. [17,18] have interpreted the step formation energy per unit length as a surface

ARTICLE IN PRESS M. Farjam / Physica B 369 (2005) 187–195 Table 1 Data for Al(1 1 1) stepped surface in atomic units System rs a h

2.04 7.52 4.34

Acknowledgements

Variational

Energy 103

g a b

s l

1.60 2.24 0.79

193

0.526 1.93

The author is grateful to Dr. Asgari for helpful suggestions and remarks, and to Prof. P. Ziesche, and Prof. A. Kiejna for valuable correspondence. Thanks are also due the reviewer for constructive criticisms.

Table 2 Term-by-term comparison of ledge and surface energies Index

li

103

hsi

103

0 2 4 es x c ps j sj

9.889 0.722 0.472 1.465 5.005 0.535 3.619 1.690 1.929

12.518 1.665 0.551 1.637 6.335 0.677 3.937 1.654 2.283

Index j refers to jellium and sj to stabilized-jellium.

energy corresponding to an additional surface area h plus two edge energies. It is illuminating to compare the various terms of the ledge energy with those of the surface energy times the step height. These are given in Table 2. A relation to notice from Table 2 is that for the terms i ¼ 0; x; c, we have li ¼ bhsi with b ¼ 0:79. This relationship is just what Eq. (23) established. We also see the well-known fact that the jellium model yields negative surface and ledge energies, for a highdensity metal such as aluminum, and therefore it is essential to use the stabilized-jellium model for this calculation.

4. Conclusions The stabilized-jellium model of a step on a metal surface is characterized by the parameters rs , the metal density parameter, and h, the step height. We proposed a form of trial function for the electron density, and presented calculations of the step formation energy. We gave numerical results for the Al(1 1 1) surface and found good agreement with other theoretical calculations.

Appendix A. Ledge energy expressions We write ledge and surface energies as sums of the various components as l ¼ l0 þ l2 þ l4 þ les þ lx þ lc þ lps ,

(A.1)

s ¼ s0 þ s2 þ s4 þ ses þ sx þ sc þ sps

(A.2)

and we give integral expressions for each component. We will use 1=ks as the unit of length, and introduce scaled densities g ¼ n=¯n, gþ ¼ nþ =¯n, and the various R ¼ r=¯n. (The bulk density is n¯ ¼ k3F =3p2 ). We denote by a0 ; a2 ; . . . integrals related to surface energy, and by b0 ; b2 ; . . . those related to the ledge energy. The numbers a are given by Z 1 1 a0 ¼ dz½f 5=3 ðzÞ yð zÞ, ðA:3Þ 20 1 Z 1 1 ðf 0 Þ2 , ðA:4Þ a2 ¼ dz 108p 1 f a4 ¼

Z 1 2 dzf 1=3 405p2 1 "    0 2  4 # 2 f 00 9 f 00 f 1 f0

þ , ðA:5Þ 8 f 3 f f f

Z 1 1 dzjðzÞ½f ðzÞ yð zÞ, 144p 1 Z 1 1 ax ¼ dz½f 4=3 ðzÞ yð zÞ, 8p 1 Z 1 1 ac ¼ dz½f ðzÞEc ð¯nf ðzÞÞ yð zÞEc ð¯nÞ, 6 1 Z 0 1 aps ¼ dz½f ðzÞ 1. 6 1 aes ¼

ðA:6Þ ðA:7Þ ðA:8Þ ðA:9Þ

ARTICLE IN PRESS M. Farjam / Physica B 369 (2005) 187–195

194

Here, j is the scaled electrostatic potential Z 1 Z 1 jðzÞ ¼ 4p dz0 dz00 ½f ðz00 Þ yð z00 Þ. z0

z

bc ¼

1 6

Z

1

Z dx

1

dz½gðz; xÞEc ð¯ngðz; xÞÞ  gþ ðz; xÞEc ð¯nÞ ac =g , 1

1

ðA:21Þ

(A.10) The surface energy is " ! 9=2 kF aes ax ac þ Caps 1 þ s ¼ 3=2 3 þ a0 g g kF p k2F # a2 g a4 g 3 þ þ 2 , kF kF where C is given in Eq. (12). The numbers b are given by Z 1 Z 1 1 dx dz½g5=3 ðz; xÞ b0 ¼ 20 1 1  gþ ðz; xÞ a0 =g , Z 1  Z 1 1 jrgj2 b2 ¼ a2 g , dx dz 108p 1 g 1 Z 1 2 b4 ¼ dxI 4 , 405p2 1 Z

1

I4 ¼

dzg

1=3

"

1

r2 g g

 4 # 1 rg þ   a4 g 3 , 3 g

2

I ð3Þ es ¼

Z

Z

0

1 bx ¼ 8p

Z

ðA:12Þ

ðA:13Þ ðA:14Þ

ðA:15Þ

ðA:16Þ ðA:17Þ ðA:18Þ ðA:19Þ

0 1

Z dx

1

dz½g4=3 ðz; xÞ  gþ ðz; xÞ ax =g , 1

1 6

Z

Z

1

z0

dx 1

 dzRðz; xÞ aps =g ,

(A.22)

1

where in (A.22) z0 is the z-coordinate of the terraces, ( þh=2; xo0; (A.23) z0 ¼ h=2; x40: Here, F is the scaled electrostatic potential, which we treat in Appendix B. Subtraction of the a’s in the brackets is necessary to make the result of zintegrations vanish for x ! 1. With the notation we have used we can write the ledge energy in a similar form as the surface energy as " # k4F bx bc þ Cbps b2 b4 bes þ b0 l¼ þ þ þ . 2p kF kF k2F k2F (A.24)

Appendix B. Electrostatic potential

1

r drF2 ðr; yÞR2 ðr; yÞ,

dy

ðA:11Þ

   9 r2 g rg2 8 g g

Z 1   1 ð2Þ ð3Þ dx I ð1Þ bes ¼ es þ I es þ I es , 144p 1 Z 1 ¼ dz½F1 ðz; xÞR1 ðz; xÞ aes g3 , I ð1Þ es 1 Z 1 ð2Þ I es ¼ 2 dzF1 ðz; xÞR2 ðz; xÞ, 1 2p

bps ¼

1

ðA:20Þ

In this appendix, we give expressions used to calculate the scaled electrostatic potential F ¼ F1 þ F2 , and defined by ZZ Fðz; xÞ ¼ dz0 dx0 ln½ðz0 zÞ2 þ ðx0 xÞ2 Rðz0 ; x0 Þ.

ðB:1Þ

Because of the simple x-dependence of R1 , part of F1 integration can be carried out analytically. Thus, integration by parts over z0 followed by integration over x0 yields Z 1 1 0 z0 F1 ðz; xÞ ¼ dz 0 tan 1 2p 1 jz j x 0

0 ½Eðz þ z h=2Þ jz j Eðz0 þ z þ h=2Þ

ARTICLE IN PRESS M. Farjam / Physica B 369 (2005) 187–195

Z 1 1 0 z0 dz 0 ½Eðz0 þ z h=2Þ 4 1 jz j þ Eðz0 þ z þ h=2Þ,

þ

ðB:2Þ

where E is the scaled electric field defined by Z z Z 1 0 0 dz Rs ðz Þ ¼ 4p dz0 Rs ðz0 Þ. EðzÞ  4p z

References

The density R2 is localized near the step, so it is more convenient to calculate F2 in polar coordinates and make use of the expansion of the kernel [22], ln jr r0 j2 ¼ lnðr24 Þ þ 2

1 X m¼1

  1 ro m

cos mðy y0 Þ. m r4

ðB:4Þ

Substitution of (B.4) into (B.1), in polar coordinates, yields F2 ðr; yÞ ¼ Q0 ðrÞ þ 2

1 X 1 ½Pm ðrÞ sin my m m¼1

þ Qm ðrÞ cos my, Z

ðB:5Þ

1

r0 dr0 lnðr24 Þc0 ðr0 Þ, Q0 ðrÞ ¼ 0  m Z 1 0 0 ro r dr cm ðr0 Þ, Qm ðrÞ ¼ r4 0  m Z 1 0 0 ro Pm ðrÞ ¼ r dr sm ðr0 Þ. r4 0

ðB:6Þ ðB:7Þ ðB:8Þ

Here, the functions cm ðrÞ and sm ðrÞ are the Fourier cosine and sine transforms of R2 , Z 2p cm ðrÞ ¼ R2 ðr; yÞ cos my dy, ðB:9Þ Z

0 2p

sm ðrÞ ¼

R2 ðr; yÞ sin my dy, 0

where m ¼ 0; 1; 2; . . . and y is measured from the z-axis. These can be calculated efficiently by fast Fourier transform routines [21].

1

(B.3)

where

195

ðB:10Þ

[1] C. Barreteau, F. Raoufi, M.C. Desjonque`res, D. Spanjaard, J. Phys.: Condens. Matter 15 (2003) S3171. [2] T.S. Rahman, A. Kara, S. Durukanog˘lu, J. Phys.: Condens. Matter 15 (2003) S3197. [3] N. Ne´el, T. Maroutian, L. Douillard, H.-J. Ernst, J. Phys.: Condens. Matter 15 (2003) S3227. [4] A. Mugarza, J.E. Ortega, J. Phys.: Condens. Matter 15 (2003) S3281. [5] M. Scheffler, J. Neugebauer, R. Stumpf, J. Phys.: Condens. Matter 5 (1993) A91. [6] M.C. Desjonque`res, D. Spanjaard, Concepts in Surface Physics, second ed., Springer, Berlin, 1996. [7] P. Stoltze, J. Phys.: Condens. Matter 6 (1994) 9495. [8] R. Stumpf, M. Scheffler, Phys. Rev. B 53 (1996) 4958. [9] J.R. Smith, Phys. Rev. 181 (1969) 522. [10] N.D. Lang, W. Kohn, Phys. Rev. B 1 (1970) 4555. [11] J.P. Perdew, H.Q. Tran, E.D. Smith, Phys. Rev. B 42 (1990) 11627. [12] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [13] A. Kiejna, Phys. Rev. B 47 (1993) 7361. [14] J.P. Perdew, Phys. Rev. B 21 (1980) 869. [15] H. Ishida, A. Liebsch, Phys. Rev. B 46 (1992) 7153. [16] C. Fiolhais, J.P. Perdew, Phys. Rev. B 45 (1992) 6207. [17] J.P. Perdew, P. Ziesche, C. Fiolhais, Phys. Rev. B 47 (1993) 16460. [18] P. Ziesche, J.P. Perdew, C. Fiolhais, Phys. Rev. B 49 (1994) 7916. [19] S.H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 58 (1980) 1200. [20] M.D. Thompson, H.B. Huntington, Surf. Sci. 116 (1982) 522. [21] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C, second ed., Cambridge University Press, Cambridge, 2002. [22] J.D. Jackson, Classical Electrodynamics, third ed., Wiley, New York, 1998.