Metal Surface Energy: Persistent Cancellation of Short-Range Correlation Effects beyond the Random-Phase Approximation

The role that non-local short-range correlation plays at metal surfaces is investigated by analyzing the correlation surface energy into contributions from dynamical density fluctuations of various two-dimensional wave vectors. Although short-range correlation is known to yield considerable correction to the ground-state energy of both uniform and non-uniform systems, short-range correlation effects on intermediate and short-wavelength contributions to the surface formation energy are found to compensate one another. As a result, our calculated surface energies, which are based on a non-local exchange-correlation kernel that provides accurate total energies of a uniform electron gas, are found to be very close to those obtained in the random-phase approximation and support the conclusion that the error introduced by the local-density approximation is small.


I. INTRODUCTION
The widely-used Kohn-Sham formulation of densityfunctional theory (DFT) 1 requires approximations to the exchange-correlation (xc) energy E xc [n(r)]. The simplest approximation to this functional is the so-called localdensity approximation (LDA), where E xc [n(r)] is given at each point by the xc energy of a uniform electron gas at the local density. This approximation was found to be remarkably accurate in some rather inhomogeneous situations, 2 and its widespread use in condensed-matter physics led to the early success of DFT.
Hence, it is important that the LDA be tested against benchmark systems, such as the jellium surface, and that new functionals be developed. Nevertheless, more than 30 years after Lang and Kohn reported the first selfconsistent LDA calculation of the jellium surface energy, 3 the question of the impact of non-local xc effects on the surface energy and their interplay with the strong charge inhomogeneity at the surface has remained a puzzle. 4,5 The simple LDA and more advanced density functionals such as generalized gradient approximations (GGA's) 6 and meta-GGA's 7 all predict the same jellium xc surface energy within a few percent, but show no such agreement with the available wave-function based methods: Fermi hypernetted chain (FHNC) 8 and diffusion Monte Carlo (DMC); 9 see Table I of Ref. 10.
An alternative formally exact way to find the xc energy of an arbitrary inhomogeneous system is provided by the adiabatic connection formula and the fluctuationdissipation theorem. 11 Within this approach, the exchange energy is fully determined from the exact Kohn-Sham (KS) orbitals and the correlation energy is obtained in terms of the xc kernel f xc . 12 In the randomphase approximation (RPA), f xc is taken to be zero. Full RPA or corrected-RPA calculations are now feasible not only for bulk jellium 13 but also for jellium surfaces 14 and molecules. 15,16,17 In this paper, we take a non-local xc kernel that provides accurate ground-state energies of a uniform electron gas, 13 and evaluate the jellium surface energy through the use of the adiabatic connection formula. We analyze the correlation surface energy into contributions from dynamic density fluctuations of various two-dimensional wave vectors, and find that short-range xc effects on intermediate and short-wavelength fluctuations nearly compensate. Hence, while RPA is known to be a poor approximation for the total correlation energy, our calculations show that it is a surprisingly good approximation for those changes in the correlation energy that arise in surface formation. This is in contrast with FHNC and DMC slab calculations, 8,9 which predict surface energies that are significantly higher than those obtained either in the LDA 3 or in a fully non-local RPA. 14 The FHNC variational equations have been shown to provide, in the homogeneous limit, reasonable agreement with the known properties of a uniform electron gas, 18 and DMC calculations are often regarded as essentially exact. 19 Furthermore, one may expect that when applied to non-uniform systems these wave-function-based approaches will lead to results whose accuracy is comparable to the high accuracy obtained for uniform systems. Nevertheless, we show that surface formation energies obtained from slab calculations either by a linear fit in the slab thickness 8 or as differences between slab energies and an independently determined bulk energy 9 may result in substantial imprecision, and conclude that wavefunction-based estimates need to be reconsidered.

II. THEORETICAL FRAMEWORK
We consider a jellium slab normal to the z axis, which is translationally invariant in the surface plane. Subtract-ing from the slab energy the corresponding energy of a uniform electron gas and using the adiabatic connection formula, one obtains the xc surface energy q being a wave vector parallel to the surface.k F is the Fermi momentum and γ xc q represents the surface energy associated with the good quantum number q: being the Fourier transform of the bare Coulomb interaction. n xc q,λ (z, z ′ ) andn xc q,λ (|z − z ′ |) represent Fourier components of the xc-hole density of a fictitious jellium slab at coupling strength λe 2 and the corresponding xc-hole density of a uniform electron gas of densityn =k 3 F /3π 2 , respectively. In the LDA, the xc surface energy is obtained by simply replacing n xc q,λ (z, z ′ ) by the xc-hole density of a uniform electron gas of density n(z). A parametrization of the uniform-gas xc-hole density has been given by Perdew and Wang, 20 which yields the DMC ground-state energy of a uniform electron gas. 21 According to the fluctuation-dissipation theorem, where χ q,λ (z, z ′ ; ω) is the interacting density-response function. Time-dependent DFT (TDDFT) shows that this function obeys the Dyson-type equation 22,23 where χ 0 q (z, z ′ ; ω) is the density-response function of noninteracting KS electrons 24 and f xc q,λ [n](z, z ′ ; ω) involves the functional derivative of the KS xc potential at coupling constant λ.
Using the coordinate-scaling relation for the λdependence of the xc kernel derived in Ref. 13, we find is the xc kernel at λ = 1. In order to derive an approximation for this quantity, we assume that the density variation [n(z) − n(z ′ )] is small within the short range of f xc q [n](z, z ′ ; ω) and write 13 wheref xc q (n; |z − z ′ |; ω) is the Fourier transform of the xc kernelf xc (n; k, ω) of a uniform electron gas of density n. Here k = (q, k z ) represents a three-dimensional wave vector.

III. RESULTS AND DISCUSSION
We have carried out simplified surface-energy calculations withf xc q (n; |z − z ′ |; ω) replaced byf xc (n; k = q, ω)δ(z − z ′ ) [thus assuming that the dynamic density fluctuation is slowly varying in the direction perpendicular to the surface] and using the parametrization of Richardson and Ashcroft forf xc (n; k, ω), 26 and have found that neglect of the frequency dependence of the xc kernel does not introduce significant errors. We have also carried out adiabatic LDA (ALDA) surfaceenergy calculations withf xc q (n; |z − z ′ |; ω) replaced bȳ f xc (n; k = 0, ω = 0)δ(z − z ′ ) [thus assuming that the dynamic density fluctuation is slowly varying in all directions], and have found that the spacial range of the xc kernel cannot be neglected. These conclusions also apply to the uniform electron gas. 13 Hence, we neglect the frequency dependence of the xc kernel and exploit the accurate DMC calculations reported in Ref. 27 for the static xc kernel of a uniform electron gas. A parametrization of this data satisfying the known small-and large-wavelength asymptotic behavior has been reported, 28 which allows us to writē where C, B, g, α, and β are dimensionless functions of the electron density (see Ref. 28), n = k 3 F /3π 2 , andz = z −z ′ . The finite q → 0 limit of Eq. (7) will be dominated by the q → 0 divergence of v q (|z−z ′ |), making RPA exact in this limit. In the large-q limit, where short-wavelength excitations tend to be insensitive to the electron-density inhomogeneity, introduction of Eq. (7) into Eq. (6) is expected to yield an xc kernel that is essentially exact.
If the interacting density-response function χ q,λ (z, z ′ ; ω) entering Eq. (3) is replaced by χ 0 q (z, z ′ ; ω), Eqs. (1) and (2) yield the exact exchange surface energy, as obtained in Ref. 14. Here we focus our attention on the correlation surface energy, which for comparison we also calculate in LDA by replacing n c q,λ in Eq. (2) by the uniform-gas correlation-hole density at the local density n(z). Fig. 1 shows the wave-vector analysis γ c q of both LDA and non-LDA correlation surface energies of a jellium slab of thickness a = 7.17r s and r s = 2. 29 First of all, we focus on our LDA calculations, which have been carried out either by using the uniform-gas correlation-hole density withf xc q (n; |z − z ′ |) = 0 (RPA-based LDA) or the xc kernel of Eq. (7) (Corradini-based LDA), or by using the Perdew-Wang (PW) parametrization of Ref. 20. We observe that in the long-wavelength limit (q → 0) both RPA and Corradini calculations coincide with the PW parametrization. At shorter wavelengths, the Corradini scheme predicts a substantial correction to the RPA and accurately reproduces the PW wave-vector analysis of the correlation energy. Hence, armed with some confidence in the accuracy of our choice of the xc kernel, we apply it to the more realistic non-local scheme described above. Fig. 1 shows that our non-LDA beyond-RPA ('best' non-local) calculation coincides in the q → 0 limit 30 with the non-LDA RPA, which is exact in this limit. In the large-q limit, local and non-local calculations coincide, and our 'best' non-local calculation accurately reproduces the PW-based LDA (thin solid line), which is expected to be esentially exact in this limit.
Consequently, our 'best' non-local calculation provides both the exact small-q limit, where LDA fails badly, and the exact large-q limit, where RPA is wrong. The LDA largely underestimates our non-local correlation surface I: Non-local xc (σ xc ) and total (σ) surface energies and their local (LDA) counterparts, as obtained from Eqs. (1) and (2) with the non-local xc kernel of Eq. (7). Also shown are non-local RPA xc surface energies (σ xc RPA ) and PW-LDA total surface energies (σPW−LDA). Tiny differences between these RPA xc surface energies and those reported before 14 are entirely due to differences in the parametrization of the xc potential. Units are erg/cm 2 . energy, but Fig. 1 shows that the difference between RPA and beyond-RPA γ c q (the short-range part of the correlation surface energy) is fairly insensitive to whether the LDA is used or not. This supports the assumption made in Ref. 31 that the short-range part of the correlation energy can be treated within LDA or GGA. However, shortrange xc effects on intermediate and short-wavelength contributions to the surface energy tend to compensate, and this cancellation happens to be even more complete than expected from LDA or GGA. Thus, our non-local scheme yields surface energies which are still closer to RPA than is the RPA + of Ref. 31.
To extract the surface energy of a semi-infinite medium, we have considered three different values of the slab thickness: the threshold width at which the n = 5 subband for the z motion is completely occupied and the two widths at which the n = 5 and n = 6 subbands are half occupied, and have followed the extrapolation procedure of Ref. 14. In Table I we show our extrapolated local (LDA) and non-local surface energies, as obtained from Eqs. (1) and (2) either withf xc q (n; |z − z ′ |) = 0 (RPA) or with the xc kernel of Eq. (7). These calculations indicate that the introduction of a plausible non-local xc kernel yields short-range corrections to RPA surface energies that are negligible. For comparison, also shown in Table I are PW-LDA surface energies, as obtained either from Eqs. (1) and (2) with the PW uniform-gas xc-hole density 20 or from the PW parametrization of the uniform-gas xc energy. 25 Corradini and PW-based LDA surface energies (σ LDA and σ PW−LDA ) are found to be very close to each other, and we expect our Corradinibased non-local surface energies (σ) to be close to the exact jellium surface energy, as well.
We close this paper with an analysis of the available wave-function-based surface-energy calculations. Krotscheck et al. 39 considered slabs of four different thickness a, and obtained both a bulk energy per particle ε ∞ and a surface energy σ from the FHNC slab energy per particle ε as a function of the particle number per unit areana, by a linear fit ε(na) = ε ∞ + 2σ/na that becomes exact in the limit of infinite thickness. These authors showed that the extrapolated bulk energies (ε ∞ ) and those obtained from a separate FHNC bulk calculation (ε) agree within about 1%, and claimed that this comparison lent credibility to their numerical treatment. However, these small differences in the bulk calculation yield an uncertainty in the surface energy ∆σ = (ε − ε ∞ )na/2, which for r s = 2.07 and 4.96 can be as large as 280 and 11 erg/cm 2 , respectively. Moreover, since for the smallest/largest width under study and due to oscillatory quantum-size effects the quantity ε(na) is larger/smaller than expected for a semi-infinite medium, the extrapolated bulk and surface energies are found to be too negative and too large, respectively.
Li et al. 32 calculated the fixed-node DMC surface energy of a jellium slab with r s = 2.07 and found σ = −465 erg/cm 2 , which is ∼ 40 erg/cm 2 larger than the RPA value. They also performed LDA calculations with either the Wigner or the Ceperley-Alder form for the uniform-gas xc energy, and found LDA surface energies that are also about 40 erg/cm 2 larger than the corresponding LDA surface energies of a semi-infinite jellium, which suggests that finite-size corrections might bring the DMC surface energy into close agreement with RPA. These fixed-node DMC calculations were extended by Acioli and Ceperley to study jellium slabs at five different densities, 9 but these authors extracted the surface energy from release-node bulk energies. Both Li et al. 32 and Acioli and Ceperley 9 claimed that the releasenode correction of the uniform electron gas at r s = 2.07 is 0.0023 eV/electron, and argued that this correction would only yield a small error in the surface energy. However, the unpublished uniform-gas fixed-node energy reported and used in Ref. 32 is actually 0.0123 eV/electron higher than its release-node counterpart; hence, by combining fixed-node slab and release-node bulk energies Acioli and Ceperley produced for r s = 2.07 a surface energy that is too large by 138 erg/cm 2 . Furthermore, had these authors used fixed-node bulk energies (see, e.g., Ref. 33), they would have obtained surface energies that are close to RPA.

IV. SUMMARY AND CONCLUSIONS
We have investigated the role that short-range correlation plays at metal surfaces, on the basis of a wave-vector analysis of the correlation surface energy. Our non-local calculations, which are found to provide the exact small and large-q limits, indicate that a persistent cancellation of short-range correlation effects yields surface energies that are in excellent agreement with RPA, and support the conclusion that the error introduced by the LDA is small. Although this conclusion seems to be in contrast with available wave-function based calculations, we have shown that a careful analysis of these data might bring them into close agreement with RPA. This is consistent with recent work, where jellium surface energies extracted from DMC calculations for jellium spheres are also found to be close to RPA. 34 We have found that the RPA xc surface energy displays an error cancellation between short-and intermediaterange correlations. A different (and less complete) error cancellation between long-and intermediate-range xc effects explains 10 why the LDA works for the surface energy. The GGA corrects only the intermediate-range contributions, 10 and so gives surface energies slightly lower and less accurate than those of LDA. Usually GGA works better than LDA, but not for the jellium surface energy where long-range effects are especially important. In the present work, as in four others, 10,31,34,35 we have found that the jellium xc surface energy is only a few percent higher than it is in LDA. These closely-agreeing methods include two different short-range corrections to RPA (present work using a non-local xc kernel and Ref. 31 using an additive GGA correction), a long-range correction to GGA (Ref. 10), extraction of a surface energy from DMC energies for jellium spheres (Ref. 34), and a meta-GGA density functional (Ref. 7). The corresponding correction to the LDA or GGA surface energy has been transferred 36 successfully to the prediction of vacancy formation energies 37 and works of adhesion. 38 We have almost reached a solution of the surfaceenergy puzzle, but one piece still does not fit. Krotscheck and Kohn 39 examined a "collective RPA" which approximates our full RPA, and also used several xc kernels to correct for short-range effects. When they used an isotropic xc kernel derived from the uniform gas, in the spirit of our Eq. (6), they found surface energies very close to RPA, as we do. When they used FHNC, 8 corresponding to an anisotropic (but to the eye not very different) kernel constructed explicitly for the jellium surface, they found a large positive correction to the RPA surface energy, amounting at r s = 4 to as much as 35% of the RPA σ xc or 60% of the RPA total σ. The story of the jellium surface energy cannot reach an end until their work is reconciled with the other work on this subject. The DMC energies of jellium slabs should also be re-considered. 40 Measured surface energies of real metals have been compared with calculations in Refs. 41,42,43. However, experimental and calculational uncertainties and differences between jellium and real-metal surfaces seem to preclude a solution to the surface-energy puzzle from these comparisons.