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

C. M. Horowitz  Instituto de Investigaciones Fisicoquímicas Teóricas y
Aplicadas, (INIFTA), UNLP, CCT La Plata-CONICET, Sucursal 4, Casilla de Correo
16, (1900) La Plata. Argentina  C. R. Proetto  Centro Atómico Bariloche and
Instituto Balseiro, 8400, S. C. de Bariloche, Río Negro, Argentina  J. M.
Pitarke  CIC nanoGUNE BRTA, Tolosa Hiribidea 76, E-20018 Donostia, Basque
Country, Spain  Materia Kondentsatuaren Fisika Saila, Centro Física Materiales
CSIC-UPV/EHU, and DIPC, 644 Posta Kutxatila, E-48080 Bilbo, Basque Country,
Spain

###### Abstract

The asymptotics of the Kohn-Sham (KS) exact exchange potential $V_{x}(z)$ of a
jelliumlike semi-infinite metal is investigated, in the framework of the
optimized-effective-potential formalism of density-functional theory. Our
numerical calculations clearly show that deep into the vacuum side of the
surface $V_{x}(z)\propto e^{2}\ln(az)/z$, with $a$ being a system-dependent
constant, thus confirming the analytical calculations reported in Phys. Rev. B
81, 121106(R) (2010). A criticism of this work published in Phys. Rev. B 85,
115124 (2012) is shown to be incorrect. Our rigorous exchange-only results
provide strong constraints both for the building of approximate exchange
functionals and for the determination of the still unknown KS correlation
potential.

##  I Introduction

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.

##  II Model and basics of the OEP method

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.

Figure 1: OEP self-consistent calculations of
$\tilde{V}_{x}(z)=V_{x}(z)/(a_{0}k_{F})$ and its three contributions
$\tilde{V}_{x}^{\text{S}}(z)$, $\tilde{V}_{x}^{\Delta}(z)$, and
$\tilde{V}_{x}^{\text{Shift}}(z)$, for $r_{s}=2.07$ (solid lines) and
$r_{s}=6.00$ (dotted lines). The jellium-vacuum interface is at $z=0$. The
horizontal dashed lines on the vertical axis represent the corresponding
universal bulk limits (in Hartree atomic units):
$V_{x}(z/\lambda_{F}\to-\infty)/(a_{0}k_{F})=-1/\pi\simeq-0.318$,
$V_{x}^{\text{S}}(z/\lambda_{F}\to-\infty)/(a_{0}k_{F})=-3/2\pi\simeq-0.477$
bardeen , $V_{x}^{\Delta}(z/\lambda_{F}\to-\infty)/(a_{0}k_{F})=1/2\pi\simeq
0.159$ GLLB95 , and $V_{x}^{\text{Shift}}(z/\lambda_{F}\to-\infty)=0$ note3 .

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.

##  III Numerical asymptotics

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: Upper panel: OEP self-consistent calculations of $V_{x}^{\Delta}(z)$
(in Hartree atomic units) for $r_{s}=2.07$ and $r_{s}=6.00$ (solid lines), and
their analytic asymptotics of Eq. (22) (dotted lines). Lower panel: OEP self-
consistent calculations of $V_{x}^{\text{S}}(z)$ (in Hartree atomic units) for
$r_{s}=2.07$ and $r_{s}=6.00$ (solid lines), and their analytic asymptotics of
Eq. (21) (dotted lines). The dashed curve for $r_{s}=2.07$ represents the sum
of the analytic asymptotics of Eq. (21) and the contribution from the
penetration of the vacuum orbitals (see the main text) that is neglected in
the derivation of Eq. (21).

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}$.

Figure 3:  Electron density $n(z)$ for $r_{s}=2.07$ (solid line) and
$r_{s}=6.00$ (dotted line). The jellium edge is located at $z=0$. In both
cases, the electron density is normalized with the bulk density
$\overline{n}$. By using $\lambda_{F}$ as the length unit, the comparison
(Friedel oscillations and spill out) between both electron densities becomes
feasible.

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.

Figure 4: OEP self-consistent calculations of
$\tilde{V}_{x}^{\Delta}(x)=(2\pi\,x/a_{0}k_{F})\,V_{x}^{\Delta}(x)$ and
$\tilde{V}_{x}^{\text{Shift}}(x)=(2\pi\,x/a_{0}k_{F})\,V_{x}^{\text{Shift}}(x)$
for $r_{s}=2.07$ (solid lines), with $x=\alpha k_{F}z$. The dotted lines
represent numerical fits from Eq. (23) with the fitting parameters given in
the text. OEP self-consistent calculations of
$\tilde{V}_{x}^{\text{S}}(x)=(2\pi\,x/a_{0}k_{F})\,V_{x}^{\text{S}}(x)$ and
$\tilde{V}_{x}(x)=(2\pi\,x/a_{0}k_{F})\,V_{x}(x)$ are given for comparison.
The arrow on the right represents the material-dependent asymptotic constant
value of $(2\pi\,x/a_{0}k_{F})\,V_{x}^{\text{S}}(x\gg 1)$, which for
$r_{s}=2.07$ is: $-2\pi\alpha
R_{x}\;[(e^{2}/a_{0})]\simeq-\;4.85\;[(e^{2}/a_{0})]$. Figure 5: As in Fig. 4,
but for $r_{s}=6.00$. The arrow on the right represents the asymptotic
constant value of $(2\pi\,x/a_{0}k_{F})\,V_{x}^{\text{S}}(x>>1)$, which for
$r_{s}=6.00$ is: $-2\pi\alpha
R_{x}\;[(e^{2}/a_{0})]\simeq-\;2.98\;[(e^{2}/a_{0})]$. In this case
($r_{s}=6.00$), our OEP calculations do not go as deep into the vacuum as in
the case of $r_{s}=2.07$ (Fig. 4), which is due to the fact that for
$r_{s}=6.00$ the electron density decays faster outside the metal (see Fig. 3)
and instabilities associated with the vanishing electron density appear,
therefore, earlier in this case.

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)$

Figure 6: OEP self-consistent calculations of $V_{x}^{\text{Shift}}(z)$ and
$V_{x}^{\text{Shift,Q}}(z)$, for $r_{s}=2.07$, the latter computed as in Eq.
(24), i.e.,
$V_{x}^{\text{Shift,Q}}(z)=-V_{x}^{\text{S}}(z)/2-V_{x}^{\Delta}(z)$. The
dotted line represents the analytical asymptotics of
$V_{x}^{\text{Shift,Q}}(z)$, as obtained by introducing Eqs. (21) and (22)
into Eq. (24). The inset displays an enlarged view of the far-vacuum region.
Figure 7: As in Fig. 6, but for $r_{s}=6.00$. Figure 8: OEP self-consistent
calculations of the scaled quantities
$\tilde{V}_{x}^{\Delta}(z)+\tilde{V}_{x}^{\text{Shift}}(z)$ =
$[V_{x}^{\Delta}(z)+V_{x}^{\text{Shift}}(z)](z/a_{0}R_{x})$ and
$-\;\tilde{V}_{x}^{\text{S}}(z)$ = $-\;V_{x}^{\text{S}}(z)(z/a_{0}R_{x})/2$,
for $r_{s}=2.07$ (solid lines). Separate OEP self-consistent calculations of
the scaled quantities $V_{x}^{\Delta}(z)(z/a_{0}R_{x})$ and
$V_{x}^{\text{S}}(z)(z/a_{0}R_{x})$ are also plotted for comparison. The
dotted line represents the analytical asymptotics of
$V_{x}^{\Delta}(z)(z/a_{0}R_{x})$, as obtained from Eq. (22). The inset
displays an enlarged view of the minimum in the curve
$[V_{x}^{\Delta}(z)+V_{x}^{\text{Shift}}(z)](z/a_{0}R_{x})$. Figure 9: As in
Fig. 8, but for $r_{s}=6.00$.

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, ].

##  IV Conclusions

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.

##  V Acknowledgments

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.

## References

  * (1) N. D. Lang and W. Kohn, Phys. Rev. B 1, 4555 (1970).
  * (2) T. Grabo, J. Kreibich, S. Kurth, and E. K. U. Gross, in Strong Coulomb Interactions in Electronic Structure Calculations: Beyond the Local Density Approximation, edited by V. I. Anisimov (Gordon and Breach, Amsterdam, 2000).
  * (3) S. Kümmel and L. Kronik, Rev. Mod. Phys. 80, 3 (2008).
  * (4) J. P. Perdew, J. Chem. Phys. 123, 062201 (2005).
  * (5) J. P. Perdew, A. Ruzsinszky, J. Sun, and K. Burke, J. Chem. Phys. 140, 18A533 (2014).
  * (6) J. P. Perdew, J. Sun, R. M. Martin, and B. Delley, Int. J. Quantum Chem. 116, 847 (2016).
  * (7) C. M. Horowitz, C. R. Proetto, and S. Rigamonti, Phys. Rev. Lett. 97, 026802 (2006).
  * (8) C. M. Horowitz, C. R. Proetto, and J. M. Pitarke, Phys. Rev. B 78, 085126 (2008).
  * (9) C. M. Horowitz, L. A. Constantin, C. R. Proetto, and J. M. Pitarke, Phys. Rev. B 80, 235101 (2009).
  * (10) C. M. Horowitz, C. R. Proetto, and J. M. Pitarke, Phys. Rev. B 81, 121106(R) (2010).
  * (11) E. Engel and R. M. Dreizler, in Density Functional Theory: An Advanced Course (Springer, Berlin, 2011).
  * (12) For a review on electronic surface calculations, see M. Nekovee and J. M Pitarke, Comput. Phys. Commun. 137, 123 (2001).
  * (13) E. Engel, J. Chem. Phys. 140, 18A505 (2014).
  * (14) E. Engel, Phys. Rev. B 89, 245105 (2014).
  * (15) E. Engel, Phys. Rev. B 97, 075102 (2018).
  * (16) E. Engel, Phys. Rev. B 97, 155112 (2018).
  * (17) E. Engel, Computation 6, 35 (2018).
  * (18) P. Rinke, K. Delaney, P. García-González, and R. W. Godby, Phys. Rev. A 70, 063201 (2004).
  * (19) C.-O. Almbladh and U. von Barth, Phys. Rev. B 31, 3231 (1985).
  * (20) E. Kraisler, Isr. J. Chem. 60, 805 (2020).
  * (21) Z. Qian, Phys. Rev. B 85, 115124 (2012).
  * (22) Note that this Hartree potential includes the contribution coming from the uniform positive background (proportional to $n_{+}$). Alternatively, this contribution could have been denoted separately as the ”external potential”.
  * (23) Due to the assumed translational invariance in the $x-y$ plane, functional derivatives in this work are conveniently defined as $\delta f:=\int\\{\delta f[g]/\delta g(z)\\}\delta g(z)dz$, where $\delta g(z)$ represents a uniform variation of the function g(r) in the plane ${\bf r}=z$. This is the origin of the factor $A^{-1}$ in Eq. (9).
  * (24) J. B. Krieger, Y. Li, and G. J. Iafrate, Phys. Rev. A 46, 5453 (1992).
  * (25) Mean values are defined as $\overline{O}^{\,i}=\int\xi_{i}^{*}(z)\,O\,\xi_{i}(z)\,dz$.
  * (26) The symbol “$P$” denotes the “principal value”.
  * (27) S. Rigamonti, C. M. Horowitz and C. R. Proetto, Phys. Rev. B 92, 235145 (2015).
  * (28) W. Kohn and A. E. Mattsson, Phys. Rev. Lett. 81, 3487 (1998).
  * (29) A. Solomatin and V. Sahni, Phys. Rev. B 56, 3655 (1997).
  * (30) F. Nastos, Ph.D. Thesis, Queen’s University, Kingston, Ontario, Canada, 2000.
  * (31) Z. Qian and V. Sahni, Int. J. Quantum Chem. 104, 929 (2005).
  * (32) J. Bardeen, Phys. Rev. 49, 653 (1936).
  * (33) A. Liebsch, in Electronic Excitations at Metal Surfaces (Springer, Boston, 1997).
  * (34) A. Liebsch, J. Phys. C 19, 5025 (1986).
  * (35) See, e.g., O. Gritsenko, R. van Leeuwen, E. van Lenthe, and E. J. Baerends, Phys. Rev. A 51, 1944 (1995).
  * (36) Deep in the bulk, the electronic density goes to its bulk value, $n(z\rightarrow-\infty)\rightarrow\overline{n}$. From Eq. (13), this implies that ${\Delta V}_{x}^{k}(z\rightarrow-\infty)\rightarrow\overline{\Delta V}_{x}^{k}$, which in turn implies that $\Psi_{k}(z\rightarrow-\infty)\rightarrow 0$.
  * (37) $r_{s}$ is a useful dimensionless parameter, defined as $r_{s}:=r_{0}/a_{0}$, with $r_{0}$ being the radius of a sphere containing on average one electron. It is related to $\overline{n}$ by $r_{s}=(3/4\pi a_{0}^{3}\overline{n})^{1/3}$ and to the Fermi wave number by $a_{0}k_{F}=(9\pi/4)^{1/3}/r_{s}$.