Asymptotics of the metal-surface Kohn-Sham exact exchange potential revisited

The asymptotics of the Kohn-Sham (KS) exchange-correlation ($xc$) potential $V_{xc}(z)$ of density-functional theory (DFT) at metal surfaces remains an important and open field of research. In their seminal DFT investigation of the electronic structure of metal surfaces, Lang and Kohn LK70 pointed out that far outside the surface $V_{xc}(z)$ should behave, such as the classical image potential $-\;e^{2}/4z$, with $z$ being the distance from the metal-vacuum surface located at $z=0$. This has been widely accepted; but more than 50 yr later there is no rigorous proof of its validity yet. As a first step towards this goal, here we split $V_{xc}(z)$ into its exchange ($x$) and correlation ($c$) contributions, i.e., $V_{xc}(z)=V_{x}(z)+V_{c}(z)$, and we focus on a rigorous evaluation of $V_{x}(z)$ in the framework of the so-called optimized-effective-potential (OEP) formalism Grabo00 ; KK08 . These rigorous $x$-only results should then serve as a basis for the determination of the still unknown KS correlation potential. Besides, and from a more general perspective, exact results are always most welcome in DFT, since a success strategy for the systematic building of progressively more predictive $xc$ energy functionals is based on the satisfaction of as many exact features as possible P05 ; PRSB14 ; PSMD16 .

Along these lines, the asympotics of $V_{x}(z)$ have been analyzed in the case of (i) metal slabs –with a discrete energy spectrum HPR06 ; HPP08 ; HCPP09 ; HPP10 -, and (ii) an extended semi-infinite (SI) metal – with a continuous energy spectrum HCPP09 ; HPP10 –. In the vacuum region of a metal slab of width $d$, one finds the following universal asymptotics for the KS exchange potential HPR06 :

$$V_{x}^{\text{Slab}}(z/d\gg 1)\rightarrow-\frac{e^{2}}{z}\;$$

(1)

and the exchange energy per particle HCPP09 :

$$\varepsilon_{x}^{\text{Slab}}(z/d\gg 1)\rightarrow-\frac{e^{2}}{2z}\;,$$

(2)

which do not depend on the metal electron density. This exact result is in sharp contrast with the much faster exponential decay displayed by $V_{x}(z)$ in all local or semi-local approximations, such as local density approximation (LDA), generalized gradient approximation (GGA), and meta-GGA ed11 , which precludes an accurate evaluation of the electronic structure of metal surfaces NP11 . Equations. (1) and (2), which were obtained for jellium slabs, have been validated recently for real slabs in a series of works by Engel Engel14a ; Engel14b ; Engel18a ; Engel18b ; Engel18c . In these publications, the jellium approximation was replaced by the real atomic configuration of the slabs by using a plane-wave pseudopotential approach. After averaging over the plane parallel to the surface, the universal Eqs. (1) and (2) were found again as in the case of jellium slabs. These rigorous slab asymptotics are very much in line with the corresponding asymptotics for finite systems, such as atoms, molecules, and metal clusters, where $V_{x}(r)\sim-\;e^{2}/r$ and $\varepsilon_{x}(r)\sim-\;e^{2}/(2r)$ far within the vacuum. In the case of finite systems, correlation is known not to contribute to the leading asymptotics, so far enough in the vacuum one can also write $V_{xc}(r)\sim-\;e^{2}/r$ and $\varepsilon_{xc}(r)\sim-\;e^{2}/(2r)$. AvB85 ; RDGG04 ; Kraisler20 It should be noted here that all these rigorous asymptotics are the result of the fact that far enough from the center of the finite system or slab the electron density is dominated by the highest occupied KS orbital. This is not so for the SI metal, as in this case the energy spectrum is continuous and one cannot, therefore, isolate one single KS orbital at the Fermi level. Hence, it is important, in the case of the more subtle SI metal, to combine analytical derivations with fully numerical calculations. This is precisely what we report here, by focusing on the SI metal, as the slab exchange results [Eqs. (1) and (2)] are known to be final and free from any controversy.

In this paper, we revisit the analytical calculations reported in Ref. [HPP10, ], which we combine here with fully numerical OEP calculations to reach the conclusion that the KS exchange potential far outside a SI metal surface scales indeed as

$$V_{x}(z/\lambda_{F}\gg 1)\to{e^{2}\over z}\,\ln{(az)}\;,$$

(3)

where $\lambda_{F}$ is the Fermi wavelength and $a$ represents a coefficient that depends on the average electron density of the metal. This is in contrast with the analytical OEP derivations reported in Ref. [Qian12, ] where it is wrongly concluded that the KS exchange potential $V_{x}(z)$ coincides far away from the surface with the exchange energy per particle $\varepsilon_{x}(z)$, which is known to be of the image-like form $\varepsilon_{x}(z)\sim-\;e^{2}R_{x}/(2z)$, with $R_{x}$ being a dimensionless metal-dependent constant to be defined below.

The paper is organized as follows: In Sec. II, our SI jellium model is introduced, together with the basics of the OEP formalism, and a few known asymptotics are revisited. Fully numerical calculations of $V_{x}(z)$ and its various contributions are given in Sec. III, with a particular emphasis on the difficult large-$z$ limit. In Sec. IV, our main results are summarized.

We consider a SI jellium model of a metal surface, where the discrete character of the positive ions inside the metal is replaced by a SI uniform positive neutralizing background (the jellium). This positive jellium density is given by

$$n_{+}(z)=\overline{n}\,\,\theta(-z),$$

(4)

which describes a sharp jellium $(z<0)$ - vacuum $(z>0)$ interface at $z=0$. $\overline{n}$ is a constant with the dimensions of a three-dimensional (3D) density that through the global neutrality condition fixes the electron density. The model is invariant under translations in the $x-y$ plane (of normalization area $A$), whose immediate consequence is that the 3D KS orbitals of DFT can be factorized as

$$\varphi_{{\bf k}_{\parallel},k}({\bf r})=\frac{e^{i{{\bf k}_{\parallel}\cdot}\bm{\rho}}}{\sqrt{A}}\,\frac{\xi_{k}(z)}{\sqrt{L}},$$

(5)

where $\bm{\rho}$ and ${\bf k}_{\parallel}$ are the in-plane coordinate and wave vector, respectively, and $L$ represents a normalization length along the $z$ direction. The limit $L\rightarrow\infty$ will be taken below, both for the analytical derivations and for the numerical calculations. The spin-degenerate orbitals $\xi_{k}(z)$ are the normalized solutions of the effective one-dimensional KS equation:

$$\widehat{h}_{\text{KS}}^{k}(z)\xi_{k}(z)=\left[-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial z^{2}}+\overline{V}_{\text{H}}(z)+V_{\text{xc}}(z)-\varepsilon_{k}\right]\xi_{k}(z)=0,$$

(6)

where $\varepsilon_{k}$ are the KS energies, and $\overline{V}_{\text{H}}(z)$ is the electrostatic effective Hartree potential note1

	$\displaystyle\overline{V}_{\text{H}}(z)$	$\displaystyle:=$	$\displaystyle V_{\text{ext}}(z)+V_{\text{H}}(z)\;,$		(7)
		$\displaystyle=$	$\displaystyle 2\pi e^{2}\int_{-L/2}^{L/2}dz^{\prime}\left[n(z^{\prime})-n_{+}(z^{\prime})\right]\int_{0}^{\infty}\frac{{\rho}^{\prime}d{\rho}^{\prime}}{\sqrt{(\rho-\rho^{\prime})^{2}+(z-z^{\prime})^{2}}}\;,$
		$\displaystyle=$	$\displaystyle-2\pi e^{2}\int_{-L/2}^{L/2}dz^{\prime}\left|z-z^{\prime}\right|\left[n(z^{\prime})-n_{+}(z^{\prime})\right]\;,$

with $n(z)$ being the ground-state electron density

$$n(z)=2\sum_{{\bf k}_{\parallel},k}^{occ}|\varphi_{{\bf k}_{\parallel},k}({\bf r})|^{2}=\frac{1}{4\pi^{2}}\int_{-k_{F}}^{k_{F}}(k_{F}^{2}-k^{2})|\xi_{k}(z)|^{2}dk\;.$$

(8)

The Fermi wave number $k_{F}$ is determined by imposing the ”bulk” neutrality condition that the integral of $n(z)$ along $z$ ($|z|\leq L/2$) be equal to $\overline{n}L$; thus, $k_{F}=(3\pi^{2}\overline{n})^{1/3}$. The Fermi wavelength is simply $\lambda_{F}=2\pi/k_{F}$. Note that (i) in passing from the second to the third line in Eq.(7) the bulk neutrality condition has been used to nullify a divergent contribution arising from the $\rho^{\prime}$ integration and (ii) in Eq.(8) the discrete sums over ${\bf k}_{\parallel}$ and $k$ have been converted into the usual integrals using the conversion factors $A/(2\pi)^{2}$ and $L/(2\pi)$, respectively.

$V_{xc}(z)$ is the so-called $xc$ potential, which is obtained as the functional derivative of the $xc$-energy functional $E_{xc}[n(z)]$; for our effective one-dimensional semi-infinite system : fd

$$V_{xc}(z)=\frac{1}{A}\;\frac{\delta E_{xc}}{\delta n(z)}\;.$$

(9)

Within the OEP framework, $V_{xc}(z)$ may be expressed as follows:

$$V_{xc}(z)=V_{xc}^{\text{KLI}}(z)+V_{xc}^{\text{Shift}}(z)\;,$$

(10)

where $V_{xc}^{\text{KLI}}(z)$ represents the so-called Krieger-Li-Iafrate (KLI) contribution KLI ,

$$V_{xc}^{\text{KLI}}(z)=\int_{0}^{k_{F}}\frac{|\xi_{k}|^{2}}{2\pi^{2}n(z)}\left[u_{xc}^{k}(z)+\overline{\Delta V}_{xc}^{k}\right]\widetilde{dk}\;,$$

(11)

with $\widetilde{dk}=(k_{F}^{2}-k^{2})dk$, $\Delta V_{xc}^{k}(z)=V_{xc}(z)-u_{xc}^{k}(z)$, and $u_{xc}^{k}(z)=[4\pi/A(k_{F}^{2}-k^{2})\xi_{k}^{*}(z)]\delta E_{xc}/\delta\xi_{k}(z)$; the $u_{xc}^{k}(z)$’s are usually referred to as orbital-dependent $xc$ potentials. Besides HPP10 ,

$$V_{xc}^{\text{Shift}}(z)=-\int_{0}^{k_{F}}\frac{\left[(k_{F}^{2}-k^{2})\Psi_{k}(z)\xi_{k}(z)+\Psi^{\prime}_{k}(z)\xi^{\prime}_{k}(z)\right]}{2\pi^{2}n(z)}\widetilde{dk}\;,$$

(12)

with primes denoting derivatives with respect to the coordinate $z$, and $\Psi_{k}(z)$ are the so-called orbital shifts; they are defined as the solutions of the following differential equation note4 :

$$\widehat{h}_{\text{KS}}^{k}(z)\Psi_{k}(z)=-\left[\Delta V_{xc}^{k}(z)-\overline{\Delta V}_{xc}^{k}\right]\xi_{k}(z)\;,$$

(13)

or, equivalently, by

$$\Psi_{k}(z)=\int_{-L/2}^{L/2}\xi_{k}(z^{\prime})\Delta V_{xc}^{k}(z^{\prime})G_{k}(z^{\prime},z)~{}dz^{\prime}\;,$$

(14)

with $G_{k}(z,z^{\prime})$ being the KS Green function note5 :

$$G_{k}(z,z^{\prime})=\frac{1}{\pi}P\int_{0}^{\infty}\frac{\xi_{k^{\prime}}^{*}(z)\xi_{k^{\prime}}(z^{\prime})}{(\varepsilon_{k}-\varepsilon_{k^{\prime}})}~{}dk^{\prime}.$$

(15)

Equations (6) - (13) form a closed set of equations, whose self-consistent solutions must be obtained numerically, once the key $xc$-energy functional $E_{xc}[n]$ is provided. Importantly, this set of equations determines $V_{xc}(z)$ up to an additive constant that should be fixed by imposing some suitable boundary condition; this fact is a consequence of the metal surface being represented as a closed system RHP15 . Since $E_{xc}[n]=E_{x}[n]+E_{c}[n]$, which in turn implies that $V_{xc}(z)=V_{x}(z)+V_{c}(z)$, we will now focus on the exchange contribution $V_{x}(z)$ to the KS $xc$ potential of Eq. (10), leaving the more complicated correlation contribution for further studies. The main reason for this procedure is that since $E_{x}[n]$ is known exactly in terms of the KS orbitals $\xi_{k}(z)$, it is feasible to obtain the exact $V_{x}(z)$ as the self-consistent solution of the exchange-only version of Eqs. (6) - (13).

In the case of a SI system,

$$E_{x}[n]=A\int_{-L/2}^{L/2}dz\,n(z)\,\varepsilon_{x}(z)\;,$$

(16)

where $\varepsilon_{x}(z)$ is the position-dependent exchange energy per particle at plane $z$. This quantity, which has been studied in detail in Ref. [HCPP09, ], is given by

$$\varepsilon_{x}(z)=-\frac{e^{2}}{2\pi^{2}n(z)}\int\limits_{-k_{F}}^{k_{F}}dk\int\limits_{-k_{F}}^{k_{F}}dk^{{}^{\prime}}(k_{F}^{2}-k^{2})^{1/2}(k_{F}^{2}-k^{{}^{\prime}2})^{1/2}\int\limits_{-L/2}^{L/2}dz^{\prime}\gamma_{k}(z,z^{\prime})\gamma_{k^{{}^{\prime}}}(z^{\prime},z)F_{k,k^{{}^{\prime}}}(z,z^{\prime}),$$

(17)

where $\gamma_{k}(z,z^{\prime})=\xi_{k}(z)^{*}\xi_{k}(z^{\prime})$ and $F_{k,k^{{}^{\prime}}}(z,z^{\prime})$ is the Kohn-Mattsson function KM98 , whose explicit expression is

$$F_{k,k^{\prime}}(z,z^{\prime})=\frac{1}{4\pi}\int_{0}^{\infty}\frac{d\rho}{\rho}\frac{J_{1}\left[\rho(k_{F}^{2}-k^{2})^{1/2}\right]J_{1}\left[\rho(k_{F}^{2}-{k^{\prime}}^{2})^{1/2}\right]}{[\rho^{2}+(z-z^{\prime})^{2}]^{1/2}}\;,$$

(18)

with $J_{1}(x)$ being the cylindrical Bessel function of the first kind. In most DFT standard approximations, $\varepsilon_{x}(z)$, which is a functional of the electron density $n(z)$, is obtained from the knowledge of the exchange energy per particle of a uniform electron gas at $z$ and nearby, thus leading to local (LDA) or semi-local (GGA) approximations of $\varepsilon_{x}(z)$ and $V_{x}(z)$ that decay exponentially far into the vacuum side of the surface and fail badly to describe the actual asymptotics of these quantities which are strongly non-local in nature.

The exact KS exchange potential is obtained [see Eq. (9)] as the functional derivative of the exact exchange energy of Eq. (16). As $\varepsilon_{x}(z)$ of Eq. (17) depends explicitly on the KS orbitals but only implicitly on the electron density $n(z)$, we follow the OEP method, which allows us to deal with orbital-dependent energy functionals. Indeed, the OEP method allows us to obtain the exact KS exchange potential and, therefore, its actual asymptotics far into the vacuum side of the surface.

In order to make contact with previous discussions, we write

$$V_{x}(z)=V_{x}^{\text{KLI}}(z)+V_{x}^{\text{Shift}}(z)\;,$$

(19)

and we split the exchange part of the KLI $xc$ potential of Eq. (11) as follows:

$$V_{x}^{\text{KLI}}(z)=V_{x}^{\text{S}}(z)+V_{x}^{\Delta}(z)\;,$$

(20)

where $V_{x}^{\text{S}}(z)=2\,\epsilon_{x}(z)$ is the so-called Slater potential involving the orbital-dependent exchange potentials $u_{x}^{k}(z)$, whereas $V_{x}^{\Delta}(z)$ represents the contribution from $\overline{\Delta V}_{x}^{\,k}$. The corresponding asymptotics are SS97 ; Nastos ; QS05 as follows:

$$V_{x}^{\text{S}}(z/\lambda_{F}\gg 1)\rightarrow-\;\left(\frac{\pi+2\alpha\ln{\alpha}}{\pi(1+\alpha^{2})}\right)\frac{e^{2}}{z}=-\;{R_{x}}\frac{e^{2}}{z}\;,$$

(21)

where $\alpha=a_{0}k_{F}/\sqrt{2W}$, $a_{0}$ and $W$ being the Bohr radius and the metal-vacuum work function —in atomic units—, and HPP10

$$V_{x}^{\Delta}(z/\lambda_{F}\gg 1)\rightarrow\frac{e^{2}}{2\pi\alpha z}[\ln(\alpha k_{F}z)+C]\;,$$

(22)

where $C\sim 0.96351$. Equation (22) is obtained by replacing $\Delta V_{x}^{k}(z)$, entering the calculation of $\overline{\Delta V}_{x}^{k}$, by its bulk value [$\Delta V_{x}^{k}(z)\simeq-\;e^{2}k_{F}/\pi-u_{x}^{k}(z/\lambda_{F}\rightarrow-\infty)$], which is fully justified, as for a SI metal the mean value is not sensitive to the actual form of the KS orbitals near the metal surface and far into the vacuum. An explicit analytical expression for $u_{x}^{k}(z/\lambda_{F}\rightarrow-\infty)$ is obtained through a ${\bf k}_{\parallel}$ Fourier transform of the orbital-dependent exchange potential $u_{x}$ of a 3D electron gas HPP10 ; bardeen . The all-important boundary condition commented above is satisfied by imposing this choice for the bulk value of the exchange potential in the OEP numerical procedure. By doing so, the related additive constant to $V_{x}(z)$ is fixed in such a way that $V_{x}(z/\lambda_{F}\gg 1)\rightarrow 0$. Below, we show how these analytical asymptotics [Eqs. (21) and (22)] are accurately confirmed by our numerical OEP calculations.

Our numerical strategy to solve the SI metal-surface OEP equations starts with the definition of the following three regions: far-left bulk region ($-L/2<z<z_{1}$, with $L\to\infty$), central region (a few $\lambda_{F}$’s to the left and to the right of the jellium edge), and far-right vacuum region ($z_{N}<z<L/2$). In the bulk region, the KS eigenfunctions are taken to be of the form $\xi_{k}(z)=\sin(kz-\delta_{k})$, where $\delta_{k}$ are phase shifts, thus fixing an overall normalization constant LK70 . By doing that, the system becomes effectively infinite along the negative-$z$ direction. 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 the KS potential $\overline{V}_{\text{H}}(z)+V_{xc}(z)\sim 0$ for $z\geq z_{N}$, the KS eigenfunctions can be approximated as $\xi_{k}(z)=s\,e^{-k^{*}z}$, where $s$ is a constant and $k^{*}=(-2m\varepsilon_{k}/\hbar^{2})^{1/2}$. At the mesh points, the KS orbitals $\xi_{k}(z)$ are calculated by using the Numerov integration procedure liebsch . By doing that, the system becomes effectively infinite also along the positive-$z$ direction. Regarding the orbital shifts, they are calculated from Eq. (14), with the KS Green function being computed by using the procedure described in Ref. [1liebsch, ]. For our systems we have found it convenient to place $z_{1}$ at a distance of 5 $\lambda_{F}$ to the left of the jellium edge ($z_{1}=-5\lambda_{F}$), $z_{N}$ at a distance of 25 $\lambda_{F}$ to the right of the jellium edge ($z_{N}=25\lambda_{F}$) and a mesh of $1200$ points. The approximation $\xi_{k}(z)=s\,e^{-k^{*}z}$ quoted above is used only for $z>z_{N}$, but we have checked that it is already valid for smaller values of $z$. See Ref. [HCPP09, ] for more numerical details.

As an example of our numerical OEP calculations, in Fig. 1 we display $V_{x}(z)$ and its three contributions: $V_{x}^{\text{S}}(z)$, $V_{x}^{\Delta}(z)$, and $V_{x}^{\rm Shift}(z)$, for a high-density metal with $r_{s}=2.07$ and a low-density metal with $r_{s}=6.00$. note2 All contributions have been scaled with the factor $(a_{0}k_{F})^{-1}$, as the scaled $V_{x}(z)$ yields then an ”universal” (material independent) value in the bulk: $V_{x}(z/\lambda_{F}\to-\infty)/(a_{0}k_{F})=-1/\pi$. It is interesting to note the presence of marked Friedel-like oscillations in the case of $r_{s}=6.00$, which are known to be due to the poor screening capability of the electron gas at this low electron density. Within this context, the ”shoulder” exhibited, for $r_{s}=6.00$ but not for $r_{s}=2.07$, by $V_{x}(z)$ right after the interface on the vacuum side of the surface may be considered as a last enhanced oscillation that becomes indeed a shoulder.

First of all, we proceed with the OEP numerical study of the two contributions entering the KLI exchange potential of Eq. (20) as this serves as a strong test for our OEP calculations and allows us to estimate the regime where $V_{x}(z)$ reaches its asymptotic behavior.

Figure 2 shows the result of our OEP calculations of $V_{x}^{\text{S}}(z)$ and $V_{x}^{\Delta}(z)$ for $r_{s}=2.07$ and $r_{s}=6.00$, together with the corresponding analytical asymptotics of Eqs. (21) and (22). In the case of $V_{x}^{\Delta}(z)$, the agreement between our fully numerical calculations and the analytical asymptotics is excellent for both electron densities at distances from the surface that are beyond 6 $\lambda_{F}$. In the case of $V_{x}^{\text{S}}(z)$ and for high electron densities ($r_{s}=2.07$), however, the analytical asymptotics deviate, even at distances from the surface that are beyond 6 $\lambda_{F}$, from the actual numerical calculations. A detailed numerical study shows that this slight deviation is a consequence of the extension of the KS orbitals close to the Fermi energy beyond the interior of the metal. More explicitly, to obtain the analytical asymptotics of Eq. (21), it is assumed that the main contribution in Eq. (17) comes from inside the metal and, accordingly, the vacuum region in the $z^{\prime}$ integral is neglected SS97 ; Nastos ; QS05 . This approximation starts to fail at high electron densities ($r_{s}=2.07$) where the orbitals extend considerably into the vacuum, as in this case the neglected contribution to the integral of Eq. (17) starts to play a role. It is well known that at high electron densities the electron spill out of the surface is more pronounced, as seen in Fig. 3, where the electron-density spill out is represented for both $r_{s}=2.07$ and $r_{s}=6.00$. Indeed, when the numerical contribution from the vacuum region near the surface is added (for $r_{s}=2.07)$ to Eq. (21), we find the dashed curve in the lower panel of Fig. 2, which nicely agrees with the full numerical calculation for all distances beyond 6 $\lambda_{F}$.

Now we focus on the last -and most subtle- exact-exchange contribution $V_{x}^{\text{Shift}}(z)$ entering Eq. (19) [beyond the KLI contribution of Eq. (20)], which we plot by a solid line in Figs. 4 and 5, together with the second exact-exchange contribution $V_{x}^{\Delta}(z)$ entering Eq. (20), well outside the surface up to many Fermi wavelengths, after the scaling $\tilde{V}^{i}(x)=(2\pi\,x/a_{0}k_{F})\,V^{i}(x)$ with $x=\alpha k_{F}z$ . With this scaling, the so-called Slater potential $\tilde{V}_{x}^{\text{S}}(x)$ (also plotted in Figs. 4 and 5) yields simply [see Eq. (21)] a material-dependent constant value at a few Fermi wavelengths outside the surface.

Unlike the scaled $\tilde{V}_{x}^{\text{S}}(x)$ Slater potential, $\tilde{V}_{x}^{\Delta}(x)$ and $\tilde{V}_{x}^{\text{Shift}}(x)$ do not yield a constant value in the asymptotic region outside the surface. Instead, they are found to be well fitted, in the asymptotic region, to the following fitting equation:

$$(e^{2}/a_{0})[a^{i}+b^{i}\ln(x)]\;,$$

(23)

with $a^{i}$ and $b^{i}$ being dimensionless fitting parameters ($i=\Delta,\text{Shift}$). In the case of the second contribution entering Eq. (20) ($i=\Delta$), one finds $a^{\Delta}=0.94\;(\pm\;0.03)$ and $b^{\Delta}=1.02\;(\pm\;0.02)$ for $r_{s}=2.07$, and $a^{\Delta}=0.96\;(\pm\;0.03)$ and $b^{\Delta}=0.99\;(\pm\;0.02)$ for $r_{s}=6.00$, which perfectly agrees, within error bars, with the universal coefficients $a_{\text{e}xact}^{\Delta}=0.96351$ and $b_{\text{e}xact}^{\Delta}=1$ entering Eq. (22). In the case of the last contribution to Eq. (19) ($i=\text{Shift}$), one finds $a^{\text{Shift}}=1.70\;(\pm\;0.03)$ and $b^{\text{Shift}}=-0.77\;(\pm\;0.02)$ for $r_{s}=2.07$, and $a^{\text{Shift}}=0.34\;(\pm\;0.05)$ and $b^{\text{Shift}}=-0.70\;(\pm\;0.03)$ for $r_{s}=6.00$. Figs. 4 and 5 clearly show that the numerical fit is essentially not distinguishable from the actual OEP calculations for $x>100$ and $x>70$ for $r_{s}=2.07$ and $r_{s}=6.00$, respectively. For the electron densities under study ($r_{s}=2.07$ and $r_{s}=6.00$), the full coefficient ($b^{\Delta}+b^{\text{Shift}}$) is positive, so we conclude that for these densities the leading asymptotic contribution to $V_{x}(z)$ is indeed of the form given by Eq. (3), as anticipated in Ref. [HPP10, ] and in contrast with the main result of Ref. [Qian12, ].

In Ref. [Qian12, ], by taking only the contribution to the KS exchange potential from the KS eigenfunction at $k=k_{F}$ and after a number of approximations, Qian concluded that in the case of a SI metal surface the asymptotics of $V_{x}(z)$ are embodied by half the Slater potential $V_{x}^{\text{S}}(z)$, i.e., $V_{x}^{\text{Q}}(z/\lambda_{F}\gg 1)=V_{x}^{\text{S}}(z/\lambda_{F}\gg 1)/2$ (see Eq. (1) of Ref. [Qian12, ]). Qian also reproduced Eq. (22) above (originally reported in Ref. [HPP10, ]); but he stated (without proof, see Eq. (88) of Ref. [Qian12, ]) that

$$V_{x}^{\text{Shift,Q}}(z/\lambda_{F}\gg 1)=-\frac{1}{2}\;V_{x}^{\text{S}}(z/\lambda_{F}\gg 1)-V_{x}^{\Delta}(z/\lambda_{F}\gg 1)\;,$$

(24)

in such a way that [see Eqs. (19) and (20) above]

$$V_{x}^{\text{Q}}(z/\lambda_{F}\gg 1)=\frac{1}{2}\;V_{x}^{\text{S}}(z/\lambda_{F}\gg 1)\;.$$

(25)

Figs. 6 and 7 clearly show that the actual $V_{x}^{\text{Shift}}(z/\lambda_{F}\gg 1)$ and Qian’s version $V_{x}^{\text{Shift,Q}}(z/\lambda_{F}\gg 1)$ given by Eq. (24) do not coincide, which means that Eqs.  (1) and (88) of Ref. [Qian12, ] are simply wrong. In these figures, we plot (for $r_{s}=2.07$ and $r_{s}=6.00$) our numerical OEP calculation of the actual $V_{x}^{\text{Shift}}(z)$ together with our numerical OEP calculation of Qian’s quantity $V_{x}^{\text{Shift,Q}}(z)=-V_{x}^{\text{S}}(z)/2-V_{x}^{\Delta}(z)$. From these plots, it is clear that $V_{x}^{\text{Shift}}(z)$ and $V_{x}^{\text{Shift,Q}}(z)$ are two different functions in the asymptotic region. We emphasize here that $V_{x}^{\text{Shift,Q}}(z)$ of Eq. (24) has been plotted by taking the numerically exact $V_{x}^{\text{S}}(z)$ and $V_{x}^{\Delta}(z)$ functions, as obtained with our OEP code for all possible values of $z$, from deep into the bulk ($z/\lambda_{F}\rightarrow-\infty$) to far into the vacuum ($z/\lambda_{F}\rightarrow\infty$). Also plotted in these figures (by a dotted line) is the corresponding $V_{x}^{\text{Shift,Q}}(z/\lambda_{F}\gg 1)$ that one obtains by introducing Eqs. (21) and (22) into Eq. (24). The asymptotic limit of $V_{x}^{\text{Shift,Q}}(z)$ is nicely reached at $z/\lambda_{F}\sim 2$; but it does not coincide with the actual $V_{x}^{\text{Shift}}(z/\lambda_{F}\gg 1)$

Figs. 8 and 9 show additional evidence of the fact that Eqs. (1) and (88) of Ref. [Qian12, ] are not correct. According to Eq. (24) above (Eq. (88) of Ref. [Qian12, ]):

$$V_{x}^{\Delta}(z/\lambda_{F}\gg 1)+V_{x}^{\text{Shift,Q}}(z/\lambda_{F}\gg 1)\rightarrow-V_{x}^{\text{S}}(z/\lambda_{F}\gg 1)/2=(R_{x}/2)e^{2}/z.$$

(26)

Hence, if the actual $V_{x}^{\text{Shift}}(z)$ were to coincide far outside the surface with $V_{x}^{\text{Shift,Q}}(z)$, as claimed by Qian in Ref. [Qian12, ]), the scaled quantity $[V_{x}^{\Delta}(z)+V_{x}^{\text{Shift}}(z)](z/a_{0}R_{x})$ [see left-hand side of Eq. (26)] should approach the constant value 1/2 (in atomic units) as $z/\lambda_{F}\gg 1$. Figs. 8 and 9 clearly show that this is not the case. Indeed, the scaled quantity $[V_{x}^{\Delta}(z)+V_{x}^{\text{Shift}}(z)](z/a_{0}R_{x})$ exhibits a minimum at $z/\lambda_{F}\sim 9$ and $z/\lambda_{F}\sim 12$, for $r_{s}=2.07$ and $r_{s}=6.00$, respectively. As the scaled $V_{x}^{\Delta}(z)$ has a positive slope in the displayed range of $z$, while the opposite is true for the scaled $V_{x}^{\text{Shift}}(z)$, the existence of a minimum implies that the large-$z$ behavior of $V_{x}^{\Delta}(z)$ overcomes the large-$z$ contribution of $V_{x}^{\text{Shift}}(z)$, and in fact provides the asymptotic limit for $V_{x}(z)$, as stated in Ref. [HPP10, ].

We have presented a detailed analysis of the KS exchange potential $V_{x}(z)$ at a jellium-like SI metal surface, by numerically solving the corresponding OEP equations of DFT. Our numerical calculations clearly show that, at least for the metal electron densities under study, deep into the vacuum side of the surface $V_{x}(z)\propto e^{2}\ln(az)/z$, for $z>>\lambda_{F}$ and with $a$ being a system-dependent constant, thus confirming the analytical calculations reported in Ref. [HPP10, ]. This result is in sharp contrast with the universal (metal independent) asymptotics for metal slabs, where $V_{x}^{\text{Slab}}(z)\sim-\;e^{2}/z$, for $z\gg d$. An important difference between these two cases is that whereas for metal slabs the asymptotics are dominated by the so-called Slater potential $V_{x}^{\text{S}}(z)$ entering Eq. (20), for the SI metal the asymptotics are dominated by the two remaining contributions $V_{x}^{\Delta}(z)$ and $V_{x}^{\text{Shift}}(z)$ entering Eqs.  (19) and (20).

As in the case of the exact-exchange energy density HCPP09 , we attribute the qualitatively different behavior of $V_{x}(z\to\infty)$ and $V_{x}^{\text{Slab}}(z\to\infty)$ to the fact that these asymptotes are approached in different ranges. Although in the case of the SI metal the asymptote is reached at distances $z$ from the surface that are large compared to the Fermi wavelength $\lambda_{F}$ (the only existing length scale in this model), for metal slabs the asymptote is reached at distances $z$ from the surface that are large compared to the slab thickness $d$. For thick slabs with $d>>\lambda_{F}$, $V_{x}^{\text{Slab}}(z)$ first coincides with $V_{x}(z)$ [dictated by Eq. (3) at $z>>\lambda_{F}$] in the vacuum region near the surface, but at distances from the surface that are considerably larger than $d$, $V_{x}^{\text{Slab}}(z)$ turns to the slab image-like behavior of the form of Eq. (1). For increasingly wide slabs, the SI regime extends to larger values of $z$ (in the intermediate range $\lambda_{F}<z\ll d$, if feasible) , and finally for $d\rightarrow\infty$, the slab regime is never reached.

Our results and conclusions are in contrast with the claim of Ref. [Qian12, ] that for the SI metal the asymptotics of the KS exchange potential $V_{x}(z)$ are embodied by half the Slater potential $V_{x}^{\text{S}}(z)$. Our OEP numerical calculations clearly show that this is not the case. Indeed, the SI-metal actual asymptotics are dominated by $V_{x}^{\Delta}(z)$ and $V_{x}^{\text{Shift}}(z)$, which after scaling and in the asymptotic region are found to be well fitted by Eq. (23). In the case of $V_{x}^{\Delta}(z)$, the coefficients $a^{\Delta}$ and $b^{\Delta}$ entering Eq. (23) are universal, i.e., do not depend on the electron density. In the case of $V_{x}^{\text{Shift}}(z)$, however, $a^{\text{Shift}}$ and $b^{\text{Shift}}$ they both depend on the electron density even after scaling. For the metal electron densities under study ($r_{s}=2.07$ and $r_{s}=6.00)$, there is a net logarithmic contribution to Eq. (23) and the leading asymptotic contribution to $V_{x}(z)$ happens to be of the form of Eq. (3). Based on our experience with the OEP formalism, it should be feasible to perform ab-initio calculations, such as the ones presented here but for a real SI metal-vacuum surface, relaxing the jellium approximation. It would then be interesting to check whether the $e^{2}\ln(az)/z$ asymptotic scaling of the exchange potential found in Ref.[HPP10, ] and confirmed here survives to such strong test, as it happens in the slab case.

Now that we have a clear understanding of the asymptotics of the KS exchange potential of DFT for both metal slabs and a SI metal, the next ambitious step is to investigate the contribution to the asymptotics coming from correlation. In the case of three-dimensional finite systems, correlation is known not to contribute to the leading asymptotics. This consideration is not applicable, in principle, neither to the slab which is finite along $z$ but extended in the perpendicular plane nor to the SI metal-vacuum case, which is extended in the three spatial directions. Whether in the case of a SI metal surface correlation brings, in the asymptotic region, $V_{xc}(z)$ to the classical image potential of the form $-e^{2}/4z$ we do not know yet. Work in this direction is now in progress.

C.M.H. wishes to acknowledge financial support received from CONICET of Argentina through PIP 2014-47029. C.R.P. wishes to acknowledge financial support received from CONICET and ANPCyT of Argentina through Grants No.PIP 2014-47029 and PICT 2016-1087.