Exact-exchange Kohn-Sham potential, surface energy, and work function of jellium slabs

Exact-exchange self-consistent calculations of the Kohn-Sham potential, surface energy, and work function of jellium slabs are reported in the framework of the Optimized Effective Potential (OEP) scheme of Density Functional Theory. In the vacuum side of the jellium surface and at a distance $z$ that is larger than the slab thickness, the exchange-only Kohn-Sham potential is found to be image-like ($\sim -e^2/z$) but with a coefficient that differs from that of the classical image potential $V_{im}(z)=-e^2/4z$. The three OEP contributions to the surface energy (kinetic, electrostatic, and exchange) are found to oscillate as a function of the slab thickness, as occurs in the case of the corresponding calculations based on the use of single-particle orbitals and energies obtained in the Local Density Approximation (LDA). The OEP work function presents large quantum size effects that are absent in the LDA and which reflect the intrinsic derivative discontinuity of the exact Kohn-Sham potential.


I. INTRODUCTION
The analysis of the electronic structure of metal surfaces poses a big theoretical challenge; a suitable calculational tool is needed for large, interacting, and strongly inhomogeneous many-electron systems. More than 30 years since its first application by Lang and Kohn 1,2 to the surface problem, little doubt exists that one method of choice for fulfilling this goal is density functional theory ͑DFT͒. 3,4 DFT aims to a microscopic understanding of atoms, molecules, clusters, surfaces, and bulk solids starting from the fundamental laws of quantum mechanics. In the Kohn-Sham 5 ͑KS͒ implementation of DFT, the complicated many-body problem is mapped to an effective single-particle problem, with particles subjected to an effective single-particle potential ͑the KS potential͒. Although this mapping is exact, it gives no clue as to how to calculate in practice the so-called exchange-correlation ͑xc͒ contribution to the KS potential. Lang and Kohn 1,2 solved this problem by using the local-density approximation ͑LDA͒ for the surface problem. In LDA, the xc potential at each point is taken to be that of a homogeneous interacting electron gas with the local density. Since then, many authors have calculated the electronic properties of metal surfaces by using either the LDA ͑Ref. 6͒ or further elaborations that incorporate nonlocal ingredients to the unknown xc functional. 7,8 Other schemes of the computational electronicstructure tool kit available for the investigation of solid surfaces are the Fermi hypernetted-chain ͑FHNC͒ method, 9,10 the GW approximation, 11 quantum Monte Carlo ͑QMC͒, [12][13][14] and the inhomogeneous Singwi-Tosi-Land-Sjölander ͑ISTLS͒ approach. 15 In the framework of the optimized effective potential ͑OEP͒ scheme of DFT, 16,17 which had been first used in the context of atomic physics, 18 correlation is ignored altogether and the exact-exchange KS potential is obtained. Several advantages are associated with the use of the exact-exchange energy functional of DFT: ͑i͒ it corrects the self-interaction problem inherent in approximate treatments of the exchange energy 19 ͑this problem is particularly acute for localized systems such as atoms and molecules, although it is not relevant for extended systems like bulk solids and solid surfaces͒; ͑ii͒ it yields great improvements in the study of the KS eigenvalue spectrum, 20 semiconductor band structures and excitations, 21 and nonlinear optical properties; 22 ͑iii͒ it yields the correct asymptotics; 23 ͑iv͒ it reproduces the derivative discontinuity which should be present in the KS exchange potential each time the number of particles crosses through an integer value; [24][25][26][27][28] and ͑v͒ it yields the correct twodimensional ͑2D͒ exchange energy per particle in the case of a quasi-2D electron gas. 29 It is the aim of this paper to provide benchmark exact-exchange OEP calculations for jellium slabs, with the expectation that more accurate DFT schemes that include correlation be developed by starting from a well founded exchange analysis and tested once reduced to their exchange-only ͑x-only͒ counterparts.
The rest of the paper is organized as follows. We give in Sec. II the general theoretical background which will be used in the following sections. Section III is devoted to a discussion of the asymptotic behavior of the exact-exchange KS potential of jellium slabs; in Secs. IV and V we give the results that we have obtained for the OEP surface energy and work-function, respectively, and in Sec. VI we present the conclusions.

II. OEP APPROACH
Our calculations are restricted to a jellium-slab model of metal surfaces, where the discrete character of the positive ions inside the metal is replaced by a uniform distribution of positive charge ͑the jellium͒. The positive jellium density is defined as which describes a slab of width d, number density n, 30 and jellium edges at z =−d and z =0. ͑x͒ represents the Heaviside step function: ͑x͒ =1 if x Ͼ 0 and ͑x͒ =0 if x Ͻ 0. A schematic view of our jellium slab is given in Fig. 1. Besides, and for convenience for the numerical calculations, infinite barriers are located far from the jellium edges, well inside the left and right evanescent vacuum regions. We have checked that these infinite barriers are located far enough for all the numerical calculations presented here to be independent of their precise location. 31 The jellium-slab model 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 represents a normalization area. i ͑z͒ are the normalized spin-degenerate eigenfunctions for electrons in slab discrete levels ͑SDLs͒ i ͑i =1,2,...͒ with energy i . They are the solutions of the effective one-dimensional KS equation, with m e as the bare electron mass.
The KS potential V KS entering Eq. ͑3͒ is the sum of two distinct contributions, where V H ͑z͒ is the classical ͑electrostatic͒ Hartree potential, given by 32 Here, n͑z͒ is the electron number density, 33 where k F i = ͱ 2m e ͑ − i ͒ / ប, and = ͑n , d͒ is the chemical potential, which in turn is determined from the neutrality condition for the whole system by the condition ͚ i occ. ͑k F i ͒ 2 =2dn. V xc ͑z͒ is the nonclassical xc potential, which is obtained as the functional derivative of the so-called xc energy functional E xc ͓n͑z͔͒, 34 Applications of DFT typically proceed from explicit density-dependent forms of E xc , as obtained using a variety of local or semilocal approximations. However, in the last few years increasing attention has been devoted to orbitaldependent forms of E xc : E xc = E xc ͓͕ i ͖ , ͕ i ͖͔, which are only implicit functionals of the electron density n͑z͒. In this case, one resorts to the OEP method 16 or, equivalently, uses repeatedly the chain rule for functional derivatives to obtain the following expression for the xc potential of Eq. ͑7͒ ͑Ref. 35͒: where G i KS ͑zЈ , z͒ is Green's function of noninteracting KS electrons. In the calculation of KS ͑z , zЈ͒, the chain rule for functional derivatives has been used; now we are considering the density itself as a functional of the occupied SDL. In obtaining Eq. ͑13͒ from Eq. ͑12͒, we have used Eq. ͑11͒ and also that which follows from Eq. ͑6͒. Introducing Eqs. ͑11͒ and ͑13͒ into the central Eq. ͑10͒, we obtain the final and compact version of the OEP integral equation for V xc ͑z͒, ͚ i occ. and Here, ⌬V xc i ͑z͒ = V xc ͑z͒ − u xc i ͑z͒, where u xc i ͑z͒ are SDLdependent xc potentials of the form 36 The magnitudes ⌿ i ͑z͒ are called the "shifts," as they can be physically interpreted as the first-order corrections of the KS eigenfunctions i ͑z͒ under the perturbation ⌬V xc i ͑z͒. These shifts also provide a useful and practical tool for the numerical solution of the OEP equation. 37,38 From Eq. ͑15͒, we find the orthogonality constraint between the KS eigenfunctions and the shifts: ͐ i ͑z͒ ‫ء‬ ⌿ i ͑z͒dz = 0. It is also immediate that the shifts are invariant under the replacement V xc ͑z͒ → V xc ͑z͒ + ␣, with ␣ being an arbitrary constant. This means that the above set of equations determines V xc ͑z͒ up to an additive constant, which should be fixed by imposing a suitable boundary condition. Moreover, the shifts ⌿ i ͑z͒ are easily found to satisfy the following inhomogeneous differential equation: 20 Here, mean values are defined as Equations ͑3͒-͑5͒ and ͑14͒, which determine the local V xc ͑z͒ corresponding to a given SDL-dependent E xc , form a closed system of equations ͑the OEP equations͒, which should be solved in a self-consistent way. In order to accomplish some contact with other useful versions of the OEP equations for the present problem, a few additional steps are required. First of all, we write which is easily obtained from Eq. ͑3͒. Second, we multiply the left-hand side of Eq. ͑16͒ by i ͑z͒ ‫ء‬ to obtain Then, we start from the self-evident identity, we eliminate the factor ͓⌬V xc i ͑z͒ − ⌬V xc i ͔͉ i ͑z͉͒ 2 by using Eq. ͑17͒, and we obtain

͑19͒
Finally, we proceed with the elimination from Eq. ͑19͒ of the term proportional to ‫ץ‬ 2 ⌿ i ͑z͒ / ‫ץ‬z 2 , the subsequent elimination of ‫ץ‬ 2 i ͑z͒ ‫ء‬ / ‫ץ‬z 2 proceeds via the KS equations, and as a result of all these manipulations we obtain the following expression for the DFT xc potential: 39 with primes denoting derivatives with respect to the z coordinate. It is important to note that Eqs. ͑14͒ and ͑20͒ are just two different, but fully equivalent, ways to obtain the OEP xc potential for the present problem. If the shifts ⌿ i ͑z͒ are ͑arbitrarily͒ forced to be identically equal to zero, the only term that survives is V xc,1 ͑z͒. This is exactly the Krieger-Li-Iafrate 24 ͑KLI͒ approximation, which brings the identification V xc,1 ͑z͒ϵV xc KLI ͑z͒. 40 As before, Eqs. ͑3͒-͑5͒ and ͑20͒ form a closed set of equations, which should be solved self-consistently.
Both exchange and correlation have been included so far. Unless stated otherwise, we will now focus on the x-only case, where E xc , V xc ͑z͒, and u xc ͑z͒ are replaced by E x , V x ͑z͒, and u x ͑z͒, respectively. We have achieved the self-consistent numerical solution of the x-only version of the OEP equations by two different methods: ͑i͒ direct calculation of the shifts of Eq. ͑15͒, by solving Eq. ͑16͒, 37 and ͑ii͒ direct solution of the OEP integral equation for V x ͑z͒, as given by the x-only version of Eq. ͑14͒. 38 Both methods yield results that agree within numerical accuracy, although the first approach is found to be computationally more efficient than the second. Both methods face numerical instabilities beyond a critical coordinate z in the vacuum region.
Finally, we note that the exact-exchange energy of a jellium slab is given by the following expression: where u x i ͑z͒ represent the SDL-dependent exchange potentials, being the "universal" ͑that is, independent of V KS ͒ function introduced by Kohn and Mattsson 41 and J 1 ͑x͒ being the firstorder cylindrical Bessel function. 42

III. ASYMPTOTICS OF THE EXACT-EXCHANGE KS POTENTIAL
The long-range behavior of V xc ͑z͒ in the vacuum region is an important and open issue in DFT studies of metal surfaces. 43 The aim of this section is to present a detailed derivation of the analytical asymptotic limit of V x ͑z͒ reported in Ref. 39 for a slab geometry. First of all, we note that by making the choice that V KS ͑z → ϱ͒ → 0, Eq. ͑3͒ leads us to the conclusion that i ͑z → ϱ͒ → e −z ͱ −2m e i /ប for all occupied i ͑disregarding a factor involving powers of z͒. We also remark the following points: ͑i͒ due to the exponential decay of V H ͑z → ϱ͒, the assumption V KS ͑z → ϱ͒ → 0 implies that V x ͑z → ϱ͒ → 0; ͑ii͒ for this choice of the zero of energy, one finds i Ͻ 0 for all occupied states; ͑iii͒ the slowest decaying of all the occupied SDL corresponds to i = m, where m is the highest occupied SDL. Now we look at the asymptotic behavior of the shifts ⌿ i ͑z͒. Turning to the x-only version of Eq. ͑16͒, we focus on the asymptotic behavior of the three terms on the right-hand side ͑rhs͒ of this equation, 20 follows from an inspection of Eq. ͑22͒ in the limit z → ϱ; in this limit, the sum over j is exponentially dominated by the term j = m, and the result of Eq. ͑26͒ follows at once. Hence, for i m Eq. ͑24͒ yields i.e., ⌿ i ͑z → ϱ͒ → e −z␤ m . For i = m, all three terms in the rhs of Eq. ͑24͒ decay equally ͑to exponential accuracy͒, and further analysis is necessary. Equation ͑14͒ can be rewritten as follows,

͑29͒
and by studying its asymptotic limit, it is clear that its rhs can be approximated by the term i = m −1 ͑with exponential ac-curacy͒. Given that both m ͑z → ϱ͒ and ⌿ m−1 ͑z → ϱ͒ decay as e −z␤ m−1 , it follows that ⌿ m ͑z → ϱ͒ decays as m−1 ͑z → ϱ͒, that is, ⌿ m ͑z → ϱ͒ → e −z␤ m−1 . Armed with these results, the asymptotic limit of V x ͑z͒ is immediate from Eq. ͑20͒: V x,2 ͑z → ϱ͒ tends exponentially to zero, while The leading contribution to u x m ͑z → ϱ͒ is easily obtained from Eq. ͑22͒ by considering once again that in this regime the sum over j is exponentially dominated by the term j = m. For this case, the integral over the coordinate t can be evaluated analytically, yielding

͑31͒
where I 1 and L 1 are the modified Bessel and Struve functions, respectively. 42 Noting now that in this regime k F m ͉z − zЈ͉Ӎk F m z ӷ 1, 44 it is permissible to expand the integrand of Eq. ͑31͒ as follows: Using the normalization of the orbitals m ͑z͒, we obtain 45 Since the exchange potential V x ͑z͒ has been chosen to vanish at large distances from the surface into the vacuum ͓V x ͑ϱ͒ =0͔, Eq. ͑30͒ leads us to the important constraint, which fixes the undetermined constant in V x ͑z͒ discussed above. All numerical results presented here have been obtained by using this constraint. From Eqs. ͑30͒, ͑33͒, and ͑34͒, we conclude that which is the main result of this section. At this point, we emphasize that the asymptotics dictated by Eq. ͑35͒ hold only at z coordinates that are larger than 1 / k F m . As k F m is of the order of 1 / d ͑or smaller, depending on the actual value of d͒, Eq. ͑35͒ shows that the x-only KS potential happens to be four times larger than the classical image potential ͓V im ͑z͒ =−e 2 / 4z͔ only at a distance z that is considerably larger than the slab thickness. Furthermore, the arguments leading to Eqs. ͑26͒ and ͑31͒ are only valid for a discrete slab spectrum, such that there is a finite-energy-gap between m and the remaining occupied energy levels i ͑i Ͻ m͒. An extension of the present OEP framework to treat the case of a semi-infinite jellium surface 46 is now in progress. 47 Finally, we note that under the condition V KS ͑ϱ͒ = 0 Eq. ͑35͒ for the asymptotics of V x ͑z͒ remains valid when correlation is included in the evaluation of the shifts ⌿ i ͑z͒. The point here is that the shifts are separable in their exchange and correlation components, and they also satisfy separated differential equations ͓like Eq. ͑24͒ for exchange͔. Once exchange and correlation contributions are split, the analysis of the asymptotic behavior of V x ͑z͒ follows the same lines as above, and the asymptotic limit of Eq. ͑35͒ remains the same.

IV. SURFACE ENERGY
In this section, surface-energy calculations are presented, as obtained at the x-only level. The surface energy is the work required, per unit area of the new surface formed, to split the crystal in two along a plane. 1 For our slab geometry, where E͑d͒ is the total ground-state energy for each half of the slab after it is split ͑width d͒, and E͑2d͒ is the total ground-state energy of the unsplit slab ͑width 2d͒, both the split and unsplit systems with the same jellium density. Following the standard DFT energy-functional partitioning, the surface energy ͑without correlation contribution͒ can be written as the sum of three terms, 6 where K ͑d͒ is the noninteracting kinetic contribution to the surface energy, el ͑d͒ is the electrostatic surface energy due to all noncompensated positive and negative charge distributions in the slab, and x ͑d͒ is the exchange contribution to the surface energy. From elementary physical arguments, it follows that K ͑d͒ Ͻ 0, while el ͑d͒ and x ͑d͒ are both positive. 2 Also, the stability of the slab against spontaneous fragmentation is accomplished if ͑d͒ Ͼ 0. From Eqs. ͑36͒ and ͑37͒, one writes with l = K ,el,x, and and ͓see Eqs. ͑21͒ and ͑22͔͒ The dependence on the slab width d in Eqs. ͑39͒-͑41͒ enters through the self-consistent KS eigenvalues ͑ i ͒ and eigenfunctions ͓ i ͑z͔͒.
Alternatively, one can define the effective single-slab surface energies, 48 where E l unif ͑d͒ is the ground-state energy of a uniform slab of electron gas of size d and l = K ,el,x. 49 Note that Eq. ͑42͒ only reproduces the surface-energy definition of Eq. ͑36͒ as d → ϱ. However, for a correct extrapolation of finite-slab calculations to the infinite-width limit, 48 here we calculate numerically the three components of the surface energy from the single-slab ͓Eq. ͑43͔͒. We have checked that the differ-ences between surface energies obtained from Eqs. ͑36͒ and ͑42͒ are quite small even for the narrowest slabs studied and that both agree in the extrapolation toward the semi-infinite limit.
Being the ground-state density the basic ingredient of DFT, we found interesting to compare the differences between the different density profiles that we have obtained. We exhibit in Fig. 2 the self-consistent electron-density profiles that we have obtained within the x-only LDA and OEP schemes for r s = 2.07 and d =8 F . 50 It is expected that the amplitude of the difference between both densities diminishes as z approaches the slab center, where both n LDA ͑z → −d / 2͒ and n OEP ͑z → −d / 2͒ should approach n as d → ϱ. Figure 2 shows that there are noticeable differences between both densities: n LDA ͑z͒ extends further into the vacuum region than n OEP ͑z͒, which is a result of the LDA orbitals being more extended or "diffuse" than their OEP counterparts, and the amplitude of the Friedel oscillations near the surface is larger for n OEP ͑z͒ than for n LDA ͑z͒. We have found the same behavior for other values of r s . Figure 3 shows the results that we have obtained for the slab kinetic surface energy, as a function of the slab width d, for r s = 2.07. As in the case of the electron density, we have performed these calculations within the x-only LDA and OEP schemes. In the LDA, the kinetic surface energy K ͑d͒ ͑LDA͒ is obtained by introducing the x-only self-consistent LDA eigenfunctions i LDA ͑z͒ and eigenvalues i LDA into the formally exact Eq. ͑39͒. In the OEP, the kinetic surface energy K ͑d͒ ͑OEP͒ is obtained by using the same equation ͓Eq. ͑39͔͒ but with the LDA eigenfunctions and eigenvalues replaced by their x-only OEP counterparts i OEP ͑z͒ and i OEP . The strong oscillations in both K ͑d͒ ͑LDA͒ and K ͑d͒ ͑OEP͒ are the result of the sequential filling of empty slab discrete levels as d increases. Maxima in K ͑d͒ correspond to the onset for the filling of a new slab discrete level. For this particular case, and following the extrapolation procedure of Ref. 48, we have obtained the infinite-width extrapolated surface energies K ͑LDA͒ = −4832 erg/ cm 2 ͑as reported in Ref. 48͒ and K ͑OEP͒ = −4720 erg/ cm 2 . Figure 4 displays the results that we have obtained for the electrostatic contribution to the surface energy, as a function of the slab width d and for r s = 2.07, again within the x-only LDA and OEP schemes. The electrostatic surface energies el ͑d͒ ͑LDA͒ and el ͑d͒ ͑OEP͒ are obtained from Eq. ͑40͒ by using either the x-only LDA electron density n LDA ͑z͒ or the x-only OEP electron density n OEP ͑z͒, respectively. In this case, the onset for the filling of a new slab discrete level is always associated with a minimum. Following the extrapolation procedure of Ref. 48, we have obtained the infinitewidth surface energies indicated by arrows in Fig. 4: el ͑LDA͒ = 1172 erg/ cm 2 ͑as reported in Ref. 48͒ and el ͑OEP͒ = 1103 erg/ cm 2 . In Fig. 5, we show the results that we have obtained for the exact-exchange contribution to the slab surface energy, as a function of the slab width d and for r s = 2.07, again within the x-only LDA and OEP schemes. As in the case of the kinetic and electrostatic surface energies, exact-exchange surface energies x ͑d͒ ͑LDA͒ and x ͑d͒ ͑OEP͒ ͓both derived LDA and OEP self-consistent electron densities and its difference for d =8 F and r s = 2.07. Note that n LDA ͑z͒ is slightly more diffuse than n OEP ͑z͒, as n LDA ͑z͒ − n OEP ͑z͒ Ͼ 0 for z outside the jellium edge ͑in the vacuum͒.  from the formally exact Eq. ͑41͔͒ are obtained by using either the x-only self-consistent LDA eigenfunctions i LDA ͑z͒ and eigenvalues i LDA or their x-only OEP counterparts i OEP ͑z͒ and i OEP . For comparison, we have also calculated standard LDA-exchange surface energies, 1 where x unif ͑n͒ is the exchange energy per particle of a uniform electron gas of density n, x unif ͑n͒ =−3e 2 ͑3 2 n͒ 1/3 / ͑4͒, and n LDA ͑z͒ represents the x-only LDA electron density.
All x ͑d͒ ͑LDA͒, x ͑d͒ ͑OEP͒, and x LDA ͑d͒ exhibit the characteristic oscillatory behavior also shown by the other components of the surface energy. As in the case of the electrostatic surface energy, the onset for the filling of a new slab discrete level is associated with a minimum. Figure 5 shows that while the LDA ͓see Eq. ͑44͔͒ considerably overestimates the exchange surface energy, which is a known result, the exact-exchange surface energy is not very sensitive to the actual shape of the single-particle orbitals and energies, i.e., to whether LDA or OEP orbitals are used. Following the extrapolation procedure of Ref. 48, we have obtained the infinite-width surface energies indicated by arrows in Fig. 5: x LDA = 2767 erg/ cm 2 , x ͑LDA͒ = 2390 erg/ cm 2 ͑both as reported in Ref. 48͒, and x ͑OEP͒ = 2316 erg/ cm 2 .
We have also computed kinetic, electrostatic, and exchange surface energies for other values of the electrondensity parameter r s , and we have obtained the infinite-width extrapolated results shown in Table I. A comparison of the LDA and OEP calculations presented in Table I shows that ͑i͒ LDA orbitals being more delocalized than the more realistic OEP orbitals, surface energies that are based on the use of LDA orbitals are too large relative to those obtained with the use of OEP orbitals, and ͑ii͒ the sum of kinetic, electrostatic, and exchange surface energies are not very sensitive to whether LDA or OEP is used in the evaluation of the singleparticle KS eigenfunctions and eigenvalues.

V. WORK FUNCTION
The work function W is the minimum work that must be done to remove an electron from the metal at zero temperature. In the context of DFT, the rigorous expression for the work function for a slab of thickness d is 51 where is the chemical potential. We note that as we are considering an electron system that is infinite in the x-y plane, electronic relaxation effects after removal of one electron are infinitesimal. For a slab geometry, the work function becomes size dependent through the chemical potential ͑n , d͒. We are imposing the boundary condition V KS ͑ϱ͒ = 0; accordingly, W͑d͒ =− Ͼ 0. Besides, the only energy of the full KS spectrum which has a physical significance is precisely the energy of the highest occupied level, which can be identified with . 52 The work function for a slab with r s = 2.07 and d =4 F is shown schematically in Fig. 1. For this particular case, nine SDLs are occupied and is between the ninth and tenth SDLs. Now we focus on the slab-width dependence of the work function. Figure 6 shows the result of the x-only calculations that we have performed within LDA and OEP ͓W LDA and W OEP ͔ for r s = 2.07. The weakly oscillating x-only W LDA ͑d͒ is  equivalent to the slab-width dependent work-function reported by Schulte a long-time ago. 53 As discussed by Schulte, the oscillations in W LDA ͑d͒ are the result of a combination of the shift of the bottom of the slab potential well and an effective film thickness shift, both effects suffering from an abrupt change each time the number of occupied SDL changes by one. The important point here, however, is the much stronger oscillations found in our W OEP ͑d͒ calculations, whose explanation is provided now with some detail. First of all, we note that, strictly speaking, the OEP work function W OEP ͑d͒ exhibits discontinuities of large size each time a new SDL becomes infinitesimally occupied. The first discontinuity in Fig. 6 appears at the 1 SDL→ 2 SDL transition ͑for d Շ F / 2͒, the second discontinuity appears at the 2 SDL→ 3 SDL transition ͑for d Շ F ͒, and so on. In order of clarify the source of such a discontinuous behavior, we have plotted in Fig. 7 the OEP exchange potential for slightly increasing values of the slab width d, around the 6 SDL → 7 SDL transition. Each slab width d is characterized by a "filling factor" of the last occupied SDL, which is defined as follows: Hence, ␣ m → 0 + ͑implying → m + ͒ corresponds to an infinitesimally small filling of the last occupied SDL ͑i = m͒, while ␣ m → 1 − corresponds to the threshold of occupancy of the next SDL ͑i = m +1͒. The key point here is the dramatic change in V x ͑z͒ when passing from the slab thickness corresponding to ␣ 6 =1 − to the infinitesimally thicker slab corresponding to ␣ 7 =0 + ͑ϳ10 −5 ͒. The remaining curves have been obtained for slab widths corresponding to the seventh SDL being progressively occupied; as ␣ 7 increases from 0 + to 1 − , V x ͑z͒ approaches the form it had at ␣ 6 =1 − , both in depth and asymptotic behavior, the only difference being a lateral shift of V x ͑z͒ to the right that is simply due to the larger value of d.
Second, we note that the potential barrier that forms at the interface, right after the jellium edge on the vacuum side of the surface, exhibits both V x,1 ͑z͒ and V x,2 ͑z͒ contributions ͓see Eq. ͑20͔͒, so the KLI approximation ͓which sets V x,2 ͑z͒ϵ0͔ cannot be used for the analysis of the characteristic discontinuous behavior of the work function. In all cases in Fig. 7, V x ͑z → ϱ͒ → 0. While this is clearly seen in the figure for the curves corresponding to ␣ 6 =1 − and ␣ 7 =1 − ͓in which case k F m ϳ 1 / d; see the asymptotics of Eq. ͑35͔͒, it is not evident at all for the set of potentials with small occupancies of the last occupied level, i.e., ␣ 7 Ӷ 1. In this case, k F m Ӷ 1 / d and the asymptotic regime only takes place at z coordinates that go to infinity ͑as ␣ 7 → 0 + ͒ far beyond the z coordinates considered in Fig. 7. This is the situation for ␣ 7 Ӎ 10 −5 , 10 −4 , and 10 −3 . As a final remark on this figure, it is important to realize that in the bulk and near the interface the exchange potentials V x ͑z͒ corresponding to ␣ 6 =1 − and ␣ 7 =0 + are simply related through a single vertical ͑constant͒ shift. This property, which can be verified numerically from Fig. 7, may also be derived analytically ͑see below͒. Finally, we note that although we have restricted our discussion to the case of a particular SDL transition, the same happens at every highest occupied→ lowest unoccupied SDL transition.
With the aim of understanding how this discontinuous behavior of V x ͑z͒ versus the slab width explains the results of Fig. 6 for the work function W OEP ͑d͒, we show in Fig. 8 the slab OEP electronic structure just before occupation of the SDL 7 ͑left panel͒, that is, at the slab width corresponding to ␣ 6 → 1 − , and just after occupation of the SDL 7 ͑right panel͒, i.e., at the slab width corresponding to ␣ 7 → 0 + . We note that while the Hartree potential approaches zero outside the surface exponentially and remains essentially unaffected by the infinitesimal population of the SDL 7 ͑compare left and right panels of Fig. 8͒, the OEP exchange potential ͓and therefore V KS ͑z͒ as well͔ suffers the abrupt jump explained in Fig. 7 which induces in turn the corresponding abrupt jump in the Fermi level. The net result in going from the left to the right panels of Fig. 8 is that the work function W OEP ͑d͒ suffers an  abrupt ͑discontinuous͒ decrease, as the boundary condition V KS ͑ϱ͒ = 0 is rigorously valid in both cases. This discontinuous behavior of W OEP ͑d͒, shown schematically in Fig. 8, represents precisely the origin of the jumps that are visible in Fig. 6 at every threshold for SDL occupation. It is evident from Fig. 6 that the size of the discontinuity decreases as d increases.
Finally, we investigate the size of the discontinuities that are visible in Fig. 6. For this, we rewrite the central OEP equation ͓as given by Eq. ͑14͔͒ in the following way: where i ͑z , zЈ͒ = i ͑z͒ ‫ء‬ i ͑zЈ͒. In writing Eq. ͑47͒ the contribution of all the m − 1 occupied SDLs has been split from the contribution of the last occupied ͑m͒ SDL. The label m in V x ͑z ; m͒ and u x i ͑z ; m͒ has been introduced in order to emphasize that they are solutions of a system with m occupied SDLs.
Let us now define a distance Z, such that for z Ͼ Z the electron density is dominated by the contribution of the last occupied ͑m͒ SDL, which is the one with the slowest decay. Equation ͑6͒ clearly shows that Z → ϱ when k F m → 0, which is the case whenever ␣ m → 0 + , i.e., whenever the filling of the last occupied SDL is infinitesimally small. We consider the following trial solution of Eq. ͑47͒: for z Ͻ Z and k F m → 0, with C x ͑m͒ being a constant which depends on the last occupied SDL. Introducing this trial solution into Eq. ͑47͒, we obtain In the limit k F m → 0, the second-line contribution of Eq. ͑49͒ is arbitrarily small; also, the KS wave functions i ͑z͒ and eigenvalue differences ͑denominators͒ entering G i KS ͑z , zЈ͒ should be extremely similar for the slab width corresponding to m − 1 occupied levels and ␣ m−1 → 1 − and the slab width corresponding to m occupied levels and ␣ m → 0 + . Therefore, an inspection of Eq. ͑22͒ leads us, using similar arguments, to the conclusion that u x i ͑z ; m͒ → u x i ͑z ; m −1͒ for all i Ͻ m, z Ͻ Z, and k F m → 0. Under these conditions, the first line of Eq. ͑49͒ reverts to the OEP equation for a slab width corresponding to m − 1 occupied states, and the proposal of Eq. ͑48͒ is proved. Considering now that C x ͑m͒ = V x ͑z ; m͒ − V x ͑z ; m −1͒ and taking the expectation value at the last occupied state ͑m −1͒ of m − 1 system, we find Now, for the m − 1 system we can use the boundary condition V x m−1 ͑m −1͒ = ū x m−1 ͑m −1͒ and, once again, approximate which has the nice feature that both the exchange potential V x and the orbital-dependent exchange potential u x are referred to the m system. For the m system V x m ͑m͒ = ū x m ͑m͒, which does not prevent the constant C x ͑m͒ from being nonzero ͑as shown in Fig. 7͒ since the KS orbitals m−1 ͑z͒ and m ͑z͒ are different. As the slab width increases, m also increases and the difference between m−1 ͑z͒ and m ͑z͒ decreases, thereby leading to the expectation that C x ͑m͒ → 0 as d → ϱ. This is explicitly shown in Fig. 9. While this analysis explains why C x ͑m͒ 0 for any finite m, it does not give a hint about its sign; Fig. 9 shows, however, that C x ͑m͒ is positive for all m. This positive jump in V x ͑z͒ is exchange driven: at each threshold width for the occupation of a new level, a barrier appears against the occupancy of an empty SDL. This is due to the fact that intra-SDL exchange is stronger than inter-SDL exchange. As a consequence, the slab gains exchange energy by restricting new SDL occupancies. On the other hand, correlation induces in general a negative jump in V c ͑z͒, so the net jump in V xc ͑z͒ depends on the relative weight of exchange and correlation for each particular system. 54 Finally, we have observed numerically that the average of the OEP work functions for slab widths corresponding to ␣ m−1 → 1 − and ␣ m → 0 + remains the same ͑within error bars͒ for all the m values that we have considered. Hence, we have taken the infinite-width extrapolated work function to be simply that average. Table II exhibits the infinite-width x-only LDA and OEP work functions that we have obtained in this way for various values of the electron-density parameter r s . OEP work functions are slightly and systematically smaller than their LDA counterparts.

VI. CONCLUSIONS
We have reported benchmark exact-exchange selfconsistent calculations of the KS potential, surface energy, and work function of jellium slabs in the framework of the OEP scheme. Special emphasis has been put into the asymptotical behavior of the exact-exchange KS potential far into the vacuum and the large quantum size effects that are present in the slab-width dependence of the surface energy and work function.
We have performed a detailed analysis of the asymptotics of the exact-exchange KS potential far into the vacuum, 39 showing that at a distance z that is larger than the slab thickness the exact-exchange potential takes an imagelike form, V x ͑z → ϱ͒ → −e 2 / z, but with a coefficient that differs from that of the classical image potential V im ͑z͒ =−e 2 / 4z. Although this result has been obtained in the x-only approximation, it is also true in the presence of correlation due to the separability of the basic OEP equations in their basic exchange and correlation components.
The OEP kinetic, electrostatic, and exchange contributions to the surface energy of jellium slabs have been obtained as a function of the slab width d and for a set of electron densities characterized by the parameter r s . We have shown that these components of the surface energy are all oscillating functions of d, with the oscillating period being Ϸ F / 2. By a suitable extrapolation procedure, we have found the values of the different components of the surface energy of a semi-infinite jellium. We have compared our OEP surface energies with those obtained from the same formally exact expressions ͓see Eqs. ͑39͒-͑41͔͒ but using single-particle LDA wave functions and energies; we have found small differences between these OEP and LDA surface energies, which appear as a consequence of the LDA orbitals being slightly more delocalized ͑diffuse͒ than their more realistic OEP counterparts.
Finally, we have performed x-only OEP calculations of the work function of jellium slabs, again as a function of the slab width d. We have found that the OEP work-function exhibits large quantum size effects that are absent in the LDA and which reflect the intrinsic derivative discontinuity of the exact KS potential. The amplitude of this discontinuity diminishes as the slab width increases and becomes arbitrarily small as d → ϱ, i.e., in the case of a semi-infinite system. This has been proved both analytically and numerically. We also note that although the precise value of the x-only OEP work functions reported here would change with the inclusion of correlation, the exact slab work function is expected to exhibit the large quantum size effects and discontinuities observed in the present work, barring possible accidental cancellations of exchange-driven and correlationdriven contributions to the total discontinuity. The presence of these large discontinuities in the x-only OEP slab work function ͑and presumably also in the actual work function that includes correlation͒ highlights the potential danger in which can be incurred by performing elaborated calculations for a restricted set of slab sizes without performing a suitable and reliable extrapolation toward the semi-infinite case. In summary, we expect that the benchmark exact-exchange OEP calculations reported here for jellium slabs will serve as motivation and as a starting point for the development of more realistic approximations for the exchange-correlation energy functional of jellium and real surfaces.