First-principles study of the Na(110) surface

481

Surface Science 209 (1989) 481-491 Nosh-Holl~d, Amsterdam

FIRST-PRINCIPLES T. RODACH,

STUDY OF THE Na(ll0)

SURFACE

K.-P. BOHNEN

Kernforschungszentrum Karlsruhe, Institut D-7500 Karlsruhe, Fed. Rep. of Germany

ftir Nukleare FestkGrperphysik,

P.0.

Box 3640,

and K.-M. HO Ames Laboratory, US Department Ames, IA 5OQ11, USA

of Energy and Department

of Physics,

Iowa State

University,

Received 5 September 1988; accepted for publication 25 October 1988

Using first-principles total-energy calculations the geometry of the Na(ll0) surface has been determined. We found a small contraction of the first interlayer spacing by 1.6% f 0.5%. The electronic band structure of Na(ll0) shows no significant difference to the projected bulk band structure especially there are no occupied electronic surface states. Use of the HeEmann-Feynman theorem allowed the calculation of forces and thus the determination of force constants. We found an increase of up to 60% for the surface force constants perpendicular to the surface in comparison to the corresponding bulk values. This has considerable influence on the surface phonon frequencies.

1. Iff~uetion The determination of structure and dynamics of surfaces is the basic question in surface science. The experimental techniques, such as LEED, ion scattering, EELS or He scattering [I], are well developed. First-principles total-energy calculations have been used successfully for the study of bulk properties. Applying this technique to the theoretical treatment of surfaces requires a considerable amount of computational effort. But the new supercomputers now available made it possible to apply these methods also to surfaces. In this paper we present a study of the Na(llO) surface using first-principles self-consistent total energy calculations. The alkali metals are often claimed to be the best understood elements on the basis of nearly free electron gas models. Nevertheless, there are open questions. Especially for the Na(ll0) surface the interpretation of the photoemission experiment by Jensen and ~39-6028/89/$03.50 Q Elsevier Science Publishers B.V. {North-Holland Physics ~blis~ng Division)

482

T. Rodach et al. / First-principles

study of Na(Il0)

Plummer [2] has been controversial. Overhauser [3] claimed it as evidence for a charge density wave in sodium in contradiction to other authors [4]. Another explanation for the results of Jensen and Plummer could be the existence of occupied electronic surface states. For a closed-packed surface like the (110) surface in the bee structure, a lot of properties can be well described by jellium calculations [5-71. On the other hand, a successful description of, e.g. electronic surface states, lattice relaxation or reconstruction and surface phonons requires a full three-dimensional treatment. We present here self-consistent total energy calculations in the framework of density functional theory. For the Na(ll0) surface we studied multilayer relaxation, work function, surface energy, electronic surface band structure, surface force constants and surface phonons. The outline of the paper is as follows: A brief description of the method is followed by a discussion of the results for the equilibrium geometry. Section 4 is devoted to the discussion of surface phonons. A brief summary concludes this paper.

2. Calculation Self-consistent total energy calculations are performed within the local-density-functional formalism [8]. Details of our method can be found in refs. [9,10]. A norm-conserving [ll] pseudopotential is used. For electronic exchange and correlation the interpolation formula of Hedin-Lundqvist [12] is applied [13]. Plane waves with kinetic energy up to E, = 8.5 Ry are used in the expansion of the electronic wave function and plane waves with energy up to E, = 12.0 Ry are included with second-order perturbation. In addition to the total energy, forces exerted on each atomic layer are calculated using the Hellmann-Feynman (HF) theorem [14]. For the surface calculations the periodic slab geometry is used. The use of HF forces minimizes the number of trial geometries needed to determine the equilibrium geometry, especially when a number of layers are relaxed simultaneously. Having determined the equilibrium geometry an extension of the “frozen phonon” method is used to determine surface phonon frequencies. Details can be found in refs. [9,10]. This method allows us to calculate all intra- and interplanar force constant matrices at certain high symmetry points in the surface Brillouin zone (SBZ). Using these force constants phonon modes can be calculated. So far the results for the phonon modes are restricted to high symmetry points. Extension to arbitrary points in the SBZ requires the knowledge of atomic force constants. These can be obtained from the results at high symmetry points by a proper fitting procedure.

T. Rodach et al. / First-principles

483

study of Na(1 IO)

3. Equilibrium geometry We first calculated some bulk properties for sodium. This allowed us to check the reliability of our pseudopotential and at the same time gave some reference quantities for comparison with the surface results. Our bulk results are summarized in table 1 and compared with the experimental results and the calculated results of Dacorogna and Cohen [15]. Our results are in very close agreement with those of Dacorogna and Cohen. Besides the static properties, frequencies of selected phonon modes have also been determined. The phonon frequencies given in table 1 are calculated with the frozen phonon method [16]. The low frequency modes are less accurate than the higher modes due to the fact that they result from a near cancellation of two big contributions, the ionic and the electronic one. Nevertheless there is satisfactory agreement with the measured ones. Having determined the bulk quantities successfully the surface structure was studied. We employed the supercell method. Detailed tests are carried out to study the dependence of our results on the number of slab layers and the number of grid points sampled in the irreducible part of the SBZ. For the calculation slabs up to 15 atomic layers have been used together with a vacuum region up to 9 layers. The irreducible part of the SBZ is sampled by a uniform grid containing up to 64 k-points. Partial occupancy of states near the Table 1 (a) Calculated bulk properties for sodium in comparison with experimental results Calculated Lattice constant (A)

4.01

Bulk modulus (kbar) Cohesive energy (Ry)

90 - 0.0846

Experiment 4.225

[21]

73.6 1221 - 0.0833 [23]

Frequency of the zone boundary phonon in the (110) direction in THZ: Longitudinal 4.06 3.82 Transversal 2.29 2.56 0.98 0.93

[20]

(b) Detailed comparison of our calculated results for the lattice constant a, (A) and the bulk modulus B,, (kbar) with the calculations of Dacorogna and Cohen [15]. Different approximations for exchange and correlation are used. The correction for the zero-point motion (ZPMC) is described in their paper Our results

Hedin-Lundqvist, without ZPMC Hedin-Lundqvist, with ZPMC Wigner, without ZPMC Wigner, with ZPMC

Dacorogna and Cohen

a0

Bo

a0

Bo

4.01 4.03 4.05 4.07

90 85 85 80

3.99 4.01 4.04 4.06

105 100 96 91

484

T. Rodach et al. / First-principles

study

of Na(ll0)

Table 2 Convergence of relaxations (in percent interlayer distance, “-” inward, “+” outward) on Na(ll0) as a function of slab thickness and size of the sampling grid. The equilibrium geometry has been estimated using bulk force constants. A Gaussian smearing width of 0.25 eV (for the numbers in brackets of 0.05 eV) has been used for Fermi surface weighting Slab thickness

Number of k-points 8

16

36

64 - 1.9 -0.4 0.1 0.3

9 atomic layers, 5 vacuum layers

Ad,, Ad,, Ad,, Ad,,

-1.5 -1.5 -2.2 -0.5

-0.8 1.4 1.0 1.4

-1.9(-1.8) 0.0 (0.9) 0.7 (1.6) 0.5 (0.7)

11 atomic layers, 9 vacuum layers

Ad,, Ad,, Ad,, A d,,

-1.3 -1.5 -2.1 0.1

-0.9 1.4 1.1 1.0

- 2.1 0.0 0.7 0.3

13 atomic layers, 7 vacuum layers

Ad,, Ad,, A d,, %s

-1.7 -1.9 - 2.4 - 0.5

-0.9 1.4 1.0 1.0

-2.1 -0.1 0.6 0.2

(- 1.9) (-0.5) (0.0) (0.1)

Fermi level is taken into account by a Gaussian smearing scheme which broadens each energy level and calculates the Fermi energy and fractional occupancy for each state from the resultant density-of-states [17]. A Gaussian smearing width of 0.25 eV or 0.05 eV has been used for Fermi-surface weighting. The results summarized in table 2 were obtained starting from the ideal geometry and deducing the equilibrium via a force constant matrix. As a first estimate the bulk values have been used for these force constants. From the results of table 2 it is obvious that a slab of 9 atomic layers and 5 vacuum layers gives converged results. For the sampling grid a number of 36 k-points and a smearing width of 0.25 eV is sufficient. These are the parameters which have been used in all the following calculations. Our results for the force-free situation are summarized in table 3. These differ from the results given in table 2 due to changes in force constants at the surface. We found an inward relaxation of the first layer by 1.6% f 0.5% interlayer spacing, for the second layer 0.0% f 0.58, and for the third and fourth layer an outward relaxation of 0.6% & 0.8% 0.4% & 0.8%, respectively. Thus only the relaxation of the first layer differs from zero by more than the error bar. Our result of 1.6% is somewhat bigger than the fairly old LEED result of 0.33% f 0.33% by Andersson et al. [IS]. For the work function we obtained a value of 3.1 eV. The work function is very insensitive to the relaxation, the sampling grid and the number of layers (once the slab is thick enough). Our result is somewhat bigger than the experimental value of 2.9 eV

T. Rodach et al. / First-principles Table 3 Na(ll0):

Comparison of calculated results with experiment Calculated

Surface energy: (eV/surface atom) (erg/cm’) Work function (eV) Relaxation (in %interlayer distance, “ - ” inward, “ + ‘* outward): A 12 A 23 A 34 A 4.5

48.5

study of Na(i IO)

0.174

0.185 1191

245

261

3.1

2.9

-1.6+0.5 0.0 + 0.5 0.6

Experiment

*

(191 [18]

- 0.33 f 0.33

1181

0.8

0.4 * 0.8

[lS]. Lang and Kohn [S] found with a jellium calculation 3.1 eV, Monnier et al. [6] with a generalization of the Lang-Kohn method 3.13 eV. For the surface energy we found 245 erg/cm2 in close agreement with the experimental value of 261 erg/cm* [19]. Again for comparison, the Lang-Kohn result [5] is 230 erg/cm*, Rose and Dobson [7] found 250 erg/cm*, using a simple ansatz for the inhomogeneous screened linear electron gas susceptibility. Fig. 2a shows our results for the electronic band structure of the Na(ll0) surface along symmetry directions of the SBZ as given in fig. 1. For comparison fig. 2b shows the electronic band structure of Na bulk, projected on the SBZ. There are no additional occupied electronic surface states. Surface states occur only at c and Y, but lie about 2 eV above the Fermi energy. The electronic band structure of the Na(ll0) surface has been investigated with normal-e~ssion angle-resolved photoe~ssion by Jensen and Plummer [Z].

Fig. 1. Irreducible part of the SBZ.

486

T. Rodach et al. / First-principles

study

ofNa(l IO)

a 3 2 1 0 -1

-2 -3 -4

b

E

3 2 1 0 -1 -2

-3 -4 I i=

f

5

-

c

-

Y

a

I

i=

Fig. 2. (a) Electronic band structure of the Na(ll0) surface. (b) Electronic band structure of Na bulk, projected on the SBZ (the Fermi energy is set equal to zero).

The interpretation of their data was controversial: Overhauser [3] took their results as an indication that the band structure of sodium is distorted by a charge density wave. On the other hand, Shung and Mahan [4] explain the data with a more detailed model for photoemission by a nearly free electron band structure. Our results support the interpretation based on a nearly free electron metal.

4. Surface phonons As already described in section 2 the phonon spectrum (bulk and surface modes) is calculated at high symmetry points in the SBZ. This is done by calculating the force constant matrices coupling the layers in the slab. For the calculation of the phonon modes we used a slab of 50 layers. Calculations with a 100 layer film showed no significant differences to the 50 layer results. The force constant matrices coupling these layers were determined in two ways. Microscopically we calculated the inter- and intraplanar force constant matrices coupling the outermost layers. The coupling of the inner layers was described

T. Rodach et al. / First-principles Table 4 Bulk atomic force constants (dyn/cm)

study of Na(lI0)

481

as obtained by a fit to the data of Woods et al. [20]

Neighbour

Radial (F)

Tangential (C)

1 2

3830 588.5

- 181.5 86.9

by force constant matrices obtained from bulk atomic force constants. It is well established that the interatomic forces in sodium are very short ranged [20]. The bulk atomic force constants were determined by fitting inelastic neutron scattering data [20] by an axially symmetric Born-von-K&-man model including forces up to second nearest neighbours. The measured phonon frequencies could be reproduced with an rms error of 0.055 THz. In table 4 we give these atomic force constants expressed in terms of radial (F) and tangential (G) force constants. We have checked that these atomic force constants are consistent with our microscopically determined interplanar force constants for the bulk. The microscopic force constants are determined via the calculation of force differences. For the surface the two geometries used for these calculations are always the relaxed geometry as specified in table 3 and a situation where the atoms in the outermost layer are distorted according to the wave vector under consideration. Such calculations were done for the r-, S- and Y-point in the SBZ. The interplanar force constant matrices are given in table 5 for the interlayer coupling in the first layer (index 11) and the coupling of the first to the second layer (index 21). The coupling of the first to the third and fourth layers can be neglected, because in these cases we calculated force constants typically two orders of magnitude smaller than those given in table 5. At the Y-point we only calculated the Z&component. Such force constant calculations require a considerable amount of computing time and the data given in table 5 are sufficient for fitting the atomic surface force constants. We also show in table 5 the interplanar force constant matrices calculated with the bulk values of F and G given in table 4. While parallel to the surface the differences between the surface and bulk values of the force constants are small, perpendicular to the surface there is an increase of up to 60% for the surface force constants. So parallel to the surface the situation is very bulk-like, which is not astonishing for a closed-packed surface. But although the equilibrium geometry is very close to the ideal one, there is a remarkable increase in the force constants perpendicular to the surface. This shows that the assumption underlying a number of model calculations for surface phonons, bulk-like geometry means also bulk-like force constants for a surface, is not always correct. So far we have discussed only the results for high symmetry points in the irreducible part of the SBZ. This is due to the fact that our microscopic

488

T. Rodach et al. / First-principles

study of Na(l IO)

Table 5 Interplanar force constants for the surface (microscopically calculated) and the bulk (calculated with the F’s and G’s from table 4) at F, S and 7 in dyn/cm. The directions of the axes are X: (liO), Y: (OOl), Z: (110). At y only the ZZ-component was calculated Coupling

Surface

Bulk

464

0

0

312

0

0

0

2216 0

0 6926

0

2485 0

( 2662 0 )

- 464

0

0

- 312

0

0

0

-2216 0

0

- 2485 0

13030

8906

10632

7565

0

i 8906 0

11428 0

7565 0

9463 0

0 i 5283

0

0

- 502

3782 0

0 1 - 3782

11 r

21

11

(

0 ) - 6926

i

0 0 H 8332

0 ) - 5662

S 0

0

-444

4756 0

0 ) - 4756

*

0

0

11

i 0

* 0

21

* 0 0

0 * 0

21

(444 0

0 1 6549

i 502 0

i

20 255

0

0

11730 0

1038 0 0

0 - 2138 0

0 0 1 4209

Y 0 0 - 6265

0 0 -4310

calculations can give directly only interplanar force constant matrix elements. Extension of our results to arbitrary points in the SBZ requires the knowledge of atomic force constants. Since the interplanar force constants are well defined linear combinations of the atomic force constants, we can obtain the latter from the interplanar one. With a fitting procedure we obtained the values of the F’s and G’s given in table 6 under the assumption of an axially symmetric Born-von-KQrman model including forces up to second nearest

Table 6 Atomic surface force constants (dyn/cm) as obtained by a fit from the microscopically calculated interplanar surface force constants from table 5 Coupling

Neighbour

Radial (F)

11

1 2

4407 865

Tangential (G) -47 309

21

1 2

4911 559

-34 - 299

489

T. Rodachet al. / FirsI-principles study of Na(i IO)

a

b

t

5

sv

a

F

Fig. 3. (a) Catculated surface-phonon dispersion curves along symmetry directions of the SBZ. The shaded regions indicate the projection of the bulk-phonon dispersion curves onto the SBZ. Different shadings indicate phonons with different symmetries. (b) Surface-phonon dispersion curves obtained in the truncated bulk model.

neighbours. Having determined in this way atomic surface force constants the full phonon dispersion curves could be calculated. In fig. 3a we show our results for the relaxed surface using surface force constants. For comparison we present in fig. 3b the results for the truncated bulk, i.e., no relaxation and use of bulk force constants.

490

T. Rodach et al. / First-principles

study

of Na(ll0)

Due to the increased surface force constants there is an increase in the frequency of the surface modes by about 0.1 to 0.2 THz in comparison to the values of the truncated bulk. Furthermore, there appear a few new surface modes with frequencies just above the region of bulk modes in the projected phonon spectrum. These can be also traced back to the increased surface constants. Unfortunately, there are no experimental results available for comparison.

5. Summary We have used first-principles self-consistent total-energy calculations to investigate the surface structure and surface phonon spectrum of the Na(ll0) surface. Our results for work function and surface energy are close to that obtained from jellium calculations. Our investigations of electronic band structure, surface force constants and phonons go beyond the capabilities of the jellium model. We found only a small contraction of the first interlayer spacing by 1.6% + 0.5% and a bulk-like electronic band structure. These results confirm the general observation that relaxation effects are small for close-packed surfaces. But on the other hand, even for a nearly unrelaxed geometry, there is an increase of up to 60% for the surface force constants perpendicular to the surface in comparison to the bulk values. This has considerable influence on the surface phonon frequencies.

Acknowledgements Ames Laboratory is operated for the US Department of Energy by Iowa State University under contract NO. W-7405-Eng-82. Part of this work is supported by the Director for Energy Research, Office of Basic Energy Sciences (including a grant of computer time on the Cray computers at the Lawrence Livermore Laboratory), and by NATO Grant No. RG(86/0516).

References [l] M.A. van Hove and S.Y. Tong, The Structure of Surfaces (Springer, New York, 1985); H. Ibach and T.S. Rahman, in: Chemistry and Physics of Solid Surfaces, Eds. R. Vanselow and R. Howe (Springer, New York, 1984); J.P. Toennies, J. Vacuum Sci. Technol. A 2 (1984) 1055. (21 E. Jensen and E.W. Plummer, Phys. Rev. Letters 55 (1985) 1912. [3] A.W. Overhauser, Phys. Rev. Letters 55 (1985) 1916. [4] K.W.-K. Shung and G.D. Mahan, Phys. Rev. Letters 57 (1986) 1076; Comment by A.W. Overhauser and response by K.W.-K. Shung and G.D. Mahan, Phys. Rev. Letters 58 (1987) 959;

T. Rodach et al. / First-principles

[5] [6] [7] [8] (91 [lo] [ll] [12] [13]

[14] [15] [16] [17] [18] [19] [20] [21] [22] [23]

study of Na(l10)

491

R. Taylor and A.H. MacDonald, Phys. Rev. Letters 57 (1986) 1639. N.D. Lang and W. Kohn, Phys. Rev. B 1 (1970) 4555; B 3 (1971) 1215. R. Monnier, J.P. Perdew, D.C. Langreth and J.W. Wilkins, Phys. Rev. B 18 (1978) 656. J.H. Rose and J.F. Dobson, Solid State Commun. 37 (1981) 91. W. Kohn and P. Vashishta, General Density Functional Theory, in: Theory of the Inhomogeneous Electron Gas, Eds. S. Lundqvist and N.H. March (Plenum, New York, 1983). K.-P. Bohnen and K.-M. Ho, Surface Sci. 207 (1988) 105. K.-P. Bohnen and K.-M. Ho, Phys. Rev. B 32 (1985) 3446; Phys. Rev. Letters 56 (1986) 934; Phys. Rev. B, submitted. D.R. Hamann, M. Schhuter and C. Chiang, Phys. Rev. Letters 43 (1979) 1494. L. Hedin and B.I. Lundqvist, J. Phys. C 4 (1971) 2064. Unlike previous calculations in ref. [15], corrections for the nonlinear exchange and correlation interaction between the core and the valence charge densities, as discussed by S.G. Louie, S. Froyen and M.L. Cohen, Phys. Rev. B 26 (1982) 1738, are not included in our calculations. The close agreement between the present results and the results of ref. [15] suggests that such corrections have only small effects on the bulk properties of Na. H. Hellmann, Einftihrung in die Quantenchemie (Deuticke, Leipzig, 1937); R.P. Feynman, Phys. Rev. 56 (1939) 340. M.M. Dacorogna and M.L. Cohen, Phys. Rev. B 34 (1986) 4996. See K.-M. Ho, C.L. Fu and B.N. Harmon, Phys. Rev. B 29 (1984) 1575, and references therein. C.L. Fu and K.-M. Ho, Phys. Rev. B 28 (1983) 5480. S. Andersson, J.B. Pendry and P.M. Echenique, Surface Sci. 65 (1977) 539. W.R. Tyson and W.A. Miller, Surface Sci. 62 (1977) 267. A.D.B. Woods, B.N. Brockhouse, R.H. March and A.T. Stewart, Phys. Rev. 128 (1962) 1112. C.S. Barrett, Acta Cryst. 9 (1956) 671. M.S. Anderson and C.A. Swenson, Phys. Rev. B 28 (1983) 5395. T.M. DiVincenzo, R. Kalkan and L.A. Girifalco, Phys. Rev. B 9 (1974) 3180.