Position-dependent exact-exchange energy for slabs and semi-infinite jellium

The position-dependent exact-exchange energy per particle $\varepsilon_x(z)$ (defined as the interaction between a given electron at $z$ and its exact-exchange hole) at metal surfaces is investigated, by using either jellium slabs or the semi-infinite (SI) jellium model. For jellium slabs, we prove analytically and numerically that in the vacuum region far away from the surface $\varepsilon_{x}^{\text{Slab}}(z \to \infty) \to - e^{2}/2z$, {\it independent} of the bulk electron density, which is exactly half the corresponding exact-exchange potential $V_{x}(z \to \infty) \to - e^2/z$ [Phys. Rev. Lett. {\bf 97}, 026802 (2006)] of density-functional theory, as occurs in the case of finite systems. The fitting of $\varepsilon_{x}^{\text{Slab}}(z)$ to a physically motivated image-like expression is feasible, but the resulting location of the image plane shows strong finite-size oscillations every time a slab discrete energy level becomes occupied. For a semi-infinite jellium, the asymptotic behavior of $\varepsilon_{x}^{\text{SI}}(z)$ is somehow different. As in the case of jellium slabs $\varepsilon_{x}^{\text{SI}}(z \to \infty)$ has an image-like behavior of the form $\propto - e^2/z$, but now with a density-dependent coefficient that in general differs from the slab universal coefficient 1/2. Our numerical estimates for this coefficient agree with two previous analytical estimates for the same. For an arbitrary finite thickness of a jellium slab, we find that the asymptotic limits of $\varepsilon_{x}^{\text{Slab}}(z)$ and $\varepsilon_{x}^{\text{SI}}(z)$ only coincide in the low-density limit ($r_s \to \infty$), where the density-dependent coefficient of the semi-infinite jellium approaches the slab {\it universal} coefficient 1/2.


I. INTRODUCTION
The jellium model of a metal surface, introduced by Bardeen in 1936, 1 is the simplest model which reproduces qualitatively, and sometimes quantitatively, the physical properties of real metal surfaces. 2 While in his work Bardeen applied an approximated Hartree-Fock (HF) theory for the study of the electronic structure, since the seminal work of Lang and Kohn (LK) 3 the standard theoretical tool applied to the study of the electronic structure of metal surfaces has been Density-Functional Theory (DFT). 4  This approach has been highly successful, and routinely yields good results for global surface properties such as work functions, surface energies, crystal-structure relaxation and reconstruction, etc. 5 At a more basic level, however, some problems still remain to be solved, concerning for instance the asymptotic behavior of the exchange-correlation (xc) potential of the widely used Kohn-Sham (KS) approach to DFT. In the LDA, this potential decays exponentially when evaluated in the vacuum region, instead of the expected image-like ∝ − e 2 /z behavior. 6 This qualitative failure of the LDA xc potential translates to a similar failure of the positiondependent xc energy per particle, ε xc (r), which is defined through 7 E xc [n] = n(r) ε xc (r) dr, with E xc [n] being the xc-energy contribution to the universal energy functional of DFT, and n(r) representing the electron density. Three aspects of Eq. (1) are worth emphasizing: (i) it represents the basic expression for the LDA, in which the exact ε xc (r) of an arbitrary inhomogeneous electron system is replaced at each point r by that of a homogeneous electron gas at the local density n(r), and for a plethora of generalizations of the LDA, 8 (ii) since E xc [n] can be split as the sum of exchange (E x [n]) and correlation (E c [n]) contributions, one can write ε xc (r) = ε x (r) + ε c (r), and (iii) the position-dependent xc energy per particle ε xc (r) entering Eq. (1) is not unique. One can always add to ε xc (r) an arbitrary function ε(r) with the condition that weighted by the electron density n(r) integrates to zero. Here we have chosen ε xc (r) to represent the interaction between a given electron at r and its xc hole.
The goal of this work is to provide exact analytical and numerical calculations of the exact-exchange ε x (r) for jellium slabs and the semi-infinite (SI) jellium. In particular, we analyze the asymptotic behavior of the exact ε x (r) in the vacuum region far away from the surface, and we find that there is a qualitative difference between ε Slab x (z → ∞) and ε SI x (z → ∞): both exhibit an image-like behavior of the form −a e 2 /z (a > 0), but with a coefficient a that while in the case of jellium slabs is universal and equal to 1/2 in the case of a semi-infinite jellium depends on the density of the bulk material and only approaches 1/2 in the low-density limit. The results reported here should help to settle the still controversial issue of the asymptotic behavior of the position-dependent xc energy per particle and KS xc potential at metal surfaces. 7,9,10,11 Besides, being the results presented here exact at the exchange level, they should also serve as a benchmark against to which DFT xc calculations could be confronted, and hopefully improved, once reduced to their exchange-only version. In this context, very recently Luo et al. 12 have used a HF scheme to report self-consistent calculations of the surface energy and work function of jellium slabs.
The rest of the paper is organized as follows: In Section II we present the general theoretical background for both jellium slabs and the semi-infinite jellium, and we derive exact analytical expressions for the position-dependent exchange energy per particle in the vacuum region far away from the surface. Numerical calculations that are valid at all positions, from the bulk region to the vacuum, are reported in Section III. Section IV is devoted to the conclusions. In the case of jellium slabs, with the discrete character of the positive ions inside the slab being replaced by a uniform distribution of positive charge (the jellium background), the positive jellium density is which describes a slab of width d, number density n, 13 and jellium edges at z = −d and z = 0. θ(z) represents the Heaviside step function: θ(z) = 1 if z > 0 and θ(z) = 0 if z < 0.
The size of the slab is infinite in the x − y plane. The jellium slab is taken to be invariant under translations in the x − y plane, so the KS eigenfunctions can be rigorously factorized as follows where ρ and k are the in-plane coordinate and wave-vector, respectively, and A represents a normalization area in the x−y plane. ξ i (z) are normalized spin-degenerate eigenfunctions for electrons in slab discrete levels (SDL) i (i = 1, 2, ...) with energies ε i ; they are the solutions of the effective one-dimensional KS equation with m e being the electron mass. It is important to remark here that the factorization of the 3D wave-function as proposed in Eq.
(3) is only valid for the case of a local potential, as is the case of the KS implementation of DFT. On the other hand, in the HF approximation the nonlocality of the Fock potential introduces a coupling between k and i quantum numbers 12 . As a consequence, HF numerical calculations are more time-consuming than the ones presented here, either LDA, KLI 14 , or OEP. 15 The local KS potential V KS (z) entering Eq. (4) is the sum of two distinct contributions: where V H (z) is the classical (electrostatic) Hartree potential, given by 16 Here, n Slab (z) is the electron number density 17 where k i F = 2m e (ε F − ε i )/ , and ε F = ε F (n, d) is the Fermi energy or chemical potential, which in turn is determined from the neutrality condition for the whole system by the condition occ.
i (k i F ) 2 = 2π d n. V xc (z) is the nonclassical xc potential, which is obtained as the functional derivative of the xc-energy functional E xc [n(z)]: 18 Both for the slab and semi-infinite [d → ∞ limit of Eq.
(2))] geometries, and as a consequence of the translational symmetry in the x − y plane, Eq. (1) simplifies to where ε xc (z) is the position-dependent xc energy per particle at plane z. In the case of jellium slabs, the exchange-only contribution to E xc [n] (which is originated in the Pauli exchange hole with all other correlation effects excluded) is known to be given by the following where with J 1 (x) being the cylindrical Bessel function of first order. 20 Comparison of Eq. (9) and (10) yields the following expression for the exchange-only contribution to ε Slab xc (z): ε Slab which can be interpreted as the energy due to the interaction of an electron at z and its exchange-only Pauli hole. In order to demonstrate that the exact-exchange energy per particle ε slab x (z) of Eq. (12) represents indeed the interaction between an electron at z and its exact-exchange hole, we appeal to the following expression for the exchange-hole for our slab geometry 21 which represents the density of the exchange hole at point r + R (observational point) due to the presence of an electron located at r. Owing to the translational symmetry in the x − y plane, without loss of generality we have choose r = (0, z), r + R = (ρ, z + Z). Using Eq. (13), and defining z ′ = z + Z, Eq. (12) may be rewritten as which justifies the physical interpretation of ε Slab x (z) as the interaction energy of an electron located at z and the "charge distribution" given by h x (z; ρ, z ′ ). Equations (13) and (14) can also be used as a sort of alternative definition of the ε Slab x (z) investigated in this work, as they solve the non-uniqueness of ε Slab x (z) which results from its definition through the exchange-only version of Eq. (9). 22,23

Single occupied slab discrete level
In order to obtain the asymptotic behavior of ε Slab x (z) at z → ∞, we first restrict our analysis to the case where there is one single occupied SDL. 24 In this special case, in which i = j = 1, Eq. (12) yields [see also Eqs. (7) and (11)]: or, equivalently (see Appendix): with I 1 and L 1 being the modified Bessel and Struve functions, respectively. 20 We note that Eq. (16) is valid for all z, both inside and outside the jellium slab. Also, the cancellation of n Slab (z) which occurs in passing from Eq. (12) to Eq. (15) allows for the numerical calculation of ε Slab x,1 (z) for arbitrarily large values of z. a. Asymptotic behavior. For the slab geometry, it is permissible (and rigorous) to take the asymptotic limit k 1 F |z − z ′ | ≃ z k 1 F ≫ 1, although the integral over z ′ runs from −∞ to +∞. This is due to the fact that for a given z the main contribution to the integral in Eq. (16) comes from values of z ′ inside the slab (−d z ′ 0), as ξ 1 (z) decays exponentially one or two λ F 's from each jellium edge.
At this point, we define F (x) ≡ 1−I 1 (2x)/x+L 1 (2x)/x and use the asymptotic expansions of I 1 and L 1 in the limit x ≫ 1. We obtain 20 Hence, as z k 1 F ≫ 1, Eq. (16) yields the following asymptotic behavior: where The r s dependence of z 1 (d, r s ) and  (16), we refer to the Appendix.

General situation
For the general situation where more than one SDL is occupied, we obtain the asymptotic limit of Eq. (12) by using the fact that for z → ∞ (i) the electron density is dominated by the slowest decaying KS orbital, which corresponds to the highest occupied SDL (i = m), and (ii) the numerator of Eq. (12) is dominated by the term i = j = m, since all ξ i (z) with i = m decay exponentially two or three λ F 's from each jellium edge. Hence, in the vacuum region far away from the surface we find: and or, equivalently [see Eq. (11)]: Finally, following the same procedure as in the case of a single occupied SDL, we find: 25 This result represents a straightforward generalization of the result presented above [Eq. (18)] for the case of one single occupied SDL. At this point, it is interesting to note that the leading contribution to ε Slab x (z → ∞) → − e 2 /2z can be easily obtained directly from Eqs. (15) or (21), by approximating the argument inside the square root by z (in the large z limit), and using the normalization of the KS orbitals ξ i (z) and the identity By considering a slab of thickness sufficiently large to make the energy spectrum continuous, Solamatin and Sahni 10 reached the conclusion that far away from the slab the so-called Slater potential V S (z) [which is twice the exchange energy per particle: V S (z) = 2ε x (z)] decays as −e 2 /z 2 , in contrast with the asymptotic structure dictated by Eq. (22). This result is, however, not correct due to the fact that for a finite jellium slab (no matter how thick it is) the slab intrinsic discrete spectrum [corresponding to the eigenvalues ε i entering Eq. (4)] can never be replaced by a continuous one.
Equation (22) leads us to the conclusion that in the vacuum region of a finite jellium slab and at distances from the surface that are large compared to 1/k m F (which is typically larger than the slab thickness d), ε Slab x (z → ∞) → − e 2 /2z, which is exactly half the corresponding KS exact-exchange potential V x (z → ∞) → − e 2 /z. 26 Hence, as in the case of finite systems, 27 the Slater potential V S (z) of jellium slabs [or, equivalently, twice the exchange-energy per particle ε x (z)] embodies the asymptotics of the KS exchange potential In contrast, Solamatin and Sahni 9,10 concluded that in the case of a semi-infinite jellium only half the Slater potential embodies the asymptotics of the KS exchange potential, i.e., is still something remaining to be clarified on this issue. Work along these lines is now in progress. 28

B. Two-dimensional electron gas
The exchange energy of a strict two-dimensional (2D) electron gas can be obtained from that of a jellium slab with a single occupied SDL [Eq. (10) with i = j = 1], by first performing the one-dimensional non-uniform scaling 29 and then taking the limit as λ → ∞. The scaling above preserves the total number of electrons. Noting that for a single occupied SDL the jellium-slab exchange energy takes the following form where we find: where N = A (k 1 F ) 2 /(2π) represents the total number of electrons. Previously, this scaling limit had been formulated in a different way, resulting in the much generous constraint that the exchange energy per particle in the 2D (λ → ∞) limit should be greater than − ∞. 29,30,31 It is interesting to note that the exchange energy functional as given by Eq. (24) is an explicit functional of the density, which is only possible in this single occupied SDL case, due to the simple (invertible) relation between density and wave-function, as given by Eq. (25). In the general, many SDL occupied case, Eq. (25) is replaced by Eq. (7), the direct inversion from wave-functions to density is not feasible anymore, and the exchange energy functional is an explicit functional of the KS orbitals, but an implicit functional of the density, as in Eq. (10).
We note at this point that the exchange energy of a strict 2D electron gas can also be obtained directly from Eq. (15) through the replacement the Dirac delta function, and taking z = 0: Either from Eq. (26) or (27), we find for the exchange energy per particle of the strict 2D homogeneous electron gas the well-known result ε 2D In the case of a semi-infinite jellium, a half-space filled with a uniform distribution of positive charge (the jellium background), the jellium density is with the jellium edge at z = 0 defining the surface of a metal. As in the case of jellium slabs, the semi-infinite jellium is invariant under translations in the x − y plane, so the KS eigenfunctions can be factorized as follows where ρ and k are the in-plane coordinate and wave-vector, respectively, and A (L) represents a normalization area (length). ξ kz (z) are spin-degenerate eigenfunctions for electrons with a continuous energy spectrum ε kz = V KS (−∞)+( k z ) 2 /2m e (k z is a continuum quantum number). They are the solutions of the effective one-dimensional KS equation The KS potential V KS (z) is given by Eq. (5), as in the case of a jellium slab but with the slab electron density of Eq. (7) being replaced by the SI electron density In the case of a semi-infinite jellium, the position-dependent exchange energy per particle at plane z is given by the following expression: where ϕ kz (z, z ′ ) = ξ kz (z) * ξ kz (z ′ ), and F kzk ′ z (z, z ′ ) is of the form of Eq. (11) but with k i F (k j F ) being replaced by [k 2 z ] 1/2 ). a. Asymptotic behavior. The derivation of the asymptotic limit of ε SI x (z) is more delicate than in the case of jellium slabs, due to the fact that in the present case we have a continuous energy spectra. Hence, the crucial argument that we have used to derive the asymptotic behavior for jellium slabs, concerning the fact that in that case the highest occupied SDL dominates in the vacuum region far away from the surface, is not so transparent when the spectrum is continuous. Stated in other words, contributions to the asymptotic position-dependent exchange energy of Eq. (32) come indeed from values of k z and k ′ z that approach k F , but not necessarily only from the highest occupied value, i.e., from The asymptotic behavior of Eq. (32) was analyzed by Nastos 11 by assuming that at z → ∞ the KS potential V KS (z) takes the image-like form V KS (z → ∞) → − α KS e 2 /z, with α KS positive, but otherwise arbitrary. One finds that in the vacuum region far away from the surface (z → ∞) the KS orbitals ξ kz can be expanded with respect to the KS orbital at k z = k F as follows 11,33 with α standing for the square root of the ratio between the Fermi energy and the work function W (α 2 = ε F /W ). By introducing Eqs. (33) and (34) into Eqs. (31) and (32), one finds 11,33 and Furthermore, the asymptote of Eq. (36) does not depend on the actual form of ξ k F (z → ∞).
This is due to a cancellation, when z is large, of the orbitals ξ k F (z → ∞) entering the numerator and denominator of Eq. (32). This is the reason why the leading term in the expansion of ε SI x (z → ∞) is independent of α KS . That means that the result remains valid also in the absence of this image-like contribution, i.e., even assuming that the KS potential V KS (z) entering Eq. (4) decays exponentially as z → ∞.

III. NUMERICAL RESULTS
All the numerical calculations presented below have been carried out by ignoring all correlation effects (beyond the Pauli exchange). By this we mean that the xc potential V xc (z) entering Eq. (4) has been replaced by the exchange-only contribution V x (z), disregarding V c (z), both for jellium slabs and for the semi-infinite jellium. In the case of jellium slabs, V x (z) and the corresponding KS orbitals of Eq. (4) can be obtained through the solution of the discrete version of the x-only optimized effective potential (OEP) method, as given for example by Eqs. (14) or (20)   The numerical methods that we have used to obtain exact-exchange (OEP) orbitals and exchange-only LDA orbitals suffer from instabilities in the vacuum region far from the surface. As in this region exchange-only LDA orbitals are stabler than their OEP counterparts and the results presented in this work do not depend significantly on whether exact-exchange (OEP) or exchange-only LDA orbitals are used in Eq. (12) (see Fig. 1), all the calculations presented below have been obtained with the use of exchange-only LDA orbitals. We emphasize, however, that this is not a crucial approximation, and that the magnitude of the error that it introduces is given by the almost indistinguishable difference between the full and dotted lines in Fig. 1.
The numerical self-consistent calculations presented in Fig. 1 and in  In the bulk region, the KS eigenfunctions are taken to be of the form ξ k (z) = sin(kz − γ k ), where γ k are phase shifts, and this fixes an overall normalization constant. In the central region, we define a mesh of N points between z 1 and z N (z 1 < z N ), the first point z 1 being chosen far enough from the jellium edge in the bulk so that the Friedel oscillations can be neglected, and the outer point z N being chosen to be far enough from the jellium edge into the vacuum so that the effective one-electron potential is negligibly small. Since V KS (z) ∼ 0 for z ≥ z N , the orbitals can be approximated as ξ k (z N ) = a e −k * z where a is a constant and k * = (−2m e ε kz / 2 ) 1/2 . The KS orbitals at the mesh points are calculated by using the Numerov integration procedure. 38 As in the vacuum region the orbitals follow exponential form, it is numerically most stable to integrate them inwards, so the Numerov integration procedure in this case is given by:  x (z) reaches the asymptote −e 2 /2z only at a distance from the jellium edge of several Fermi wavelengths (∼ 8λ F ). We remind here that point (i) above was our main assumption in the derivation of the slab asymptotic limit of Eq. (22). This assumption is fully justified after the numerical results shown in Fig. 3.
At this point, with a few algebraic manipulations, we rewrite the asymptote of Eq. (22) in the physically motivated image-like form where and z Slab x (d, r s , z), which represents the location of the so-called image plane, results in As z → ∞, z Slab jellium slabs. The same comparison has been done for other values of r s , and we have found that our full numerical calculation is always very close (as in Fig. 6) to the asymptote of Eq. (36).
Finally, we display in Fig. 7

IV. SUMMARY AND CONCLUSIONS
We have presented a detailed analysis of the position-dependent exchange energy per particle ε x (z) at jellium slabs and the semi-infinite jellium.
For jellium slabs, we have found that in the vacuum region far away from the surface ε Slab x (z → ∞) → − e 2 /(2z), independent of the bulk electron density. This is the equivalent to the well-known result ε x (r → ∞) → − e 2 /(2r), which holds in the case of localized finite systems like atoms and molecules. 4 The equivalence between these results is however not straightforward, since slabs have an extended character in the x − y plane, being "localized" only along the z-coordinate. In the vacuum side of the surface there is a region where ε Slab x (z) coincides with ε SI x (z) and this region increases as d increases. The fitting of our numerical calculations of ε Slab x (z) to a physically motivated image-like expression is feasible, but the resulting location of the image plane [z Slab For a semi-infinite jellium, we have found that our numerical calculations agree well with the analytical asymptote [see Eq. (36)] obtained in Refs. 9,10,11 and 33, which approaches the slab asymptote − e 2 /2z only in the extreme low-density limit (r s → ∞).
We attribute the qualitatively different behavior of ε Slab x (z → ∞) and ε SI x (z → ∞) to the fact that these asymptotes are approached in different ranges. While in the case of the semi-infinite jellium the asymptote is reached at distances z from the surface that are large compared to the Fermi wavelength (the only existing length scale in this model), for slabs the asymptote is reached at distances z from the surface that are large compared to 1/k m F (which is typically larger than the slab thickness d). For thick slabs with d >> λ F (λ F being the Fermi wavelength), ε Slab x (z) first coincides with ε SI x (z) [dictated by Eq. (36) at z >> λ F ] in the vacuum region near the surface (see Fig. 5), but at distances from the surface that are large compared to 1/k m F , ε Slab x (z) turns to the slab image-like behavior of the form of Eq. (22) [or, equivalently, Eq. (38) with α Slab x = 1/2]; in the limit as d → ∞ (i.e., when the jellium slab becomes semi-infinite), ε Slab x (z → ∞) coincides with ε SI x (z → ∞) everywhere. In the low-density limit, where λ F → ∞, the condition d >> λ F is never fulfilled and ε Slab x (z) reaches (at z >> 1/k m F ) one single asymptote: the slab image-like behavior of the form of Eq. (22) [or, equivalently, Eq. (38) with α Slab x = 1/2], which turns out to coincide with the semi-infinite-jellium asymptotic behavior dictated by Eq. (36).
Finally, we note that as ε xc (z) = ε x (z) + ε c (z), the same conclusion is expected to be valid for the exchange contribution to the position-dependent xc energy per particle. Recent developments concerning the asymptotic behavior of the correlation contribution to the KS exchange correlation potential V xc (z) of a semi-infinite jellium can be found in Ref. 33. Here we will shown in detail how to go through Eqs. (15)- (18) in the text. Starting from Eq. (15), one first notice that 39 with I 1 and L 1 being the modified Bessel and Struve functions, respectively. Substitution of Eq. (A1) in Eq. (15) yields at once Eq. (16). Now, in the asymptotic limit, 40 L 1 (x ≫ 1) → I 1 (x ≫ 1) − 2 π + 2 x 2 − ... , which inserted in the definition of the function F (x) yields Eq. (17). Substitution of this expansion for F (x) in Eq. (16) leads to the expression, Expanding the denominators of Eq. (A4) in the large z limit, the different contributions in Eq. (18) arise. For instance, the leading contribution −e 2 /(2z) comes from the first term on the r.h.s. in Eq. (A4), with |z − z ′ | approximated by z. It is interesting to note that one important feature of this result, as is its material and slab-size independence, is consequence (in this context) of the normalization of the SDL wave-functions. On more general grounds, and returning to the alternative definition of ε Slab x (z) given in Eq. (14), this is more physically understood as a consequence of that the integral of h x (z; ρ, z + Z) over all possible "observational" coordinates (ρ, Z) is exactly −1. 21 The next term in the expansion, proportional to β 1 and with decay z −2 , is obtained from the sub-leading contribution of the first term on the r.h.s. in Eq. (A4), together with the leading contribution from the second term. To the order explicitly displayed in Eq. (18), no contribution arises from the last (third) term in Eq. (A4), as the leading contribution coming from this term to ε Slab x,1 (z → ∞) is of the order z −4 . While this analysis has been performed for the single-occupied SDL case, it also applies to the general case where more than a SDL is occupied, as explained in Section II.A.2. * Permanent address: Centro Atómico Bariloche and Instituto Balseiro, 8400, S. C. de Bariloche, Río Negro, Argentina