# Orbital-free functional with sub-milliHartree errors for slabs

Pavel Okun  Department of Chemistry, University of California, Irvine, CA
92617, USA  Antonio C. Cancio  Department of Physics and Astronomy, Ball State
University, Muncie, IN 47304, USA  Kieron Burke  Departments of Physics and
Astronomy and of Chemistry, University of California, Irvine

###### Abstract

Using principles of asymptotic analysis, we derive the exact leading
corrections to the Thomas-Fermi kinetic energy approximation for Kohn-Sham
electrons for slabs. This asymptotic expansion approximation includes crucial
quantum oscillations missed by standard semilocal density functionals. Because
these account for the derivative discontinuity, chemical accuracy is achieved
at fourth-order. The implications for both orbital-free electronic structure
and exchange-correlation approximations are discussed.

††preprint: APS/123-QED

##  I Introduction

Almost all modern density functional theory (DFT) calculations use the Kohn-
Sham (KS) scheme [1], where only the exchange-correlation (XC) energy is
approximated as a density functional [2]. Accuracy comes at the cost of
solving self-consistently the KS equations for the orbitals. However, the KS
kinetic energy, $T_{s}$, is also a density functional [2, 3]. If it could be
approximated with sufficient accuracy, without incurring substantial
additional computational cost, one would bypass the KS equations [4], speeding
up every DFT calculation on the planet. The dream of orbital-free DFT lives on
[5].

Thomas-Fermi (TF) theory [6, 7], the original DFT, is orbital-free, but is too
crude for modern calculations. First, $T_{s}$ is typically far larger than the
XC energy, so a far smaller fractional error is required. Second, via the
Euler equation that follows from the variational principle [8], the functional
derivative of $T_{s}$ determines the density, which is typically very poor
[9].

Over the decades, many attempts have been made to construct orbital-free
density functional approximations [2, 4, 10, 11, 12], often aimed at a limited
set of circumstances [13]. These include the original gradient expansion from
slowly-varying densities [14, 15, 2, 16], generalizations of that expansion
[13, 17], functionals designed for weakly interacting subsystems such as water
molecules [18], two-point functionals for use in materials calculations [19,
20], and functionals for surfaces [21].

But Lieb and Simon (LS) showed long ago [22, 23] that, in a very specific
semiclassical limit, TF theory becomes relatively exact. One approach to the
LS limit is to take $\hbar\to 0$, keeping $\mu$, the chemical potential,
fixed. Expanding about this limit, functionals of the potential have been
developed for the total energy of one-dimensional (1D) problems [24, 25, 26,
27]. The dominant term is given by a TF calculation, while higher order terms
include not only the gradient expansion, but also oscillating corrections to
it, and depend on both the Maslov [28] indices and whether the system is
finite or extended. By resumming this asymptotic series for a linear half-
well, the total energy of 10 non-interacting fermions was found to 33 digits
[24].

Figure 1: Slab kinetic energy per unit area vs electron number per unit area,
with the exact (black), gradient expansion (gray), and asymptotic
approximation (red, see text). The TF kinetic energy has been subtracted from
all curves. Cusps, i.e., derivative discontinuities, occur when a new band
starts to be filled.

Here, we demonstrate the capabilities of a potential functional to yield
chemical accuracy [errors below 1 milliHartree (mH)] for a three-dimensional
slab geometry, with a potential that varies in only one direction. We find
excellent results when the leading correction is added to TF theory, and
chemical accuracy when the next two orders are included, providing a
systematic, parameter-free approach to orbital-free calculations. Figure 1
illustrates the deviation of the kinetic energy from its TF value, for a
Pöschl-Teller (PT) slab with well depth $10$. The standard density functional,
the gradient expansion approximation (GEA) yields the smooth gray curve, which
averages over the density of states. Our potential functional (red curve)
includes the quantum oscillations characteristic of systems with discrete
states, with errors almost too small to be visible here. These results are
characteristic of any system to which our functional applies, and it is
unlikely any existing orbital-free density functional can compete. Our
calculation (a) shows how the asymptotic expansion corrects standard semi-
local functionals to account for the derivative discontinuity [29, 30], (b)
connects the rich field of semiclassical analysis of eigenvalue problems,
including densities of states and the importance of Maslov indices [31, 32],
to the construction of approximate functionals for non-trivial 3D problems,
(c) explicitly connects the failures of semilocal functionals for stretched
bonds (symmetry breaking) to the divergence of the asymptotic expansion [33],
and (d) provides insight for creating functionals of chemical accuracy for
systems with more general potentials [34].

##  II Background

For any quantum problem, one can consider functionals of either the potential
or the density (thanks to Hohenberg-Kohn [35]). The KS scheme requires density
functionals for the XC energy, but orbital-free approximations can use either
variable for $T_{s}$. The derivation of the gradient expansion of DFT begins
with a potential functional approximation (PFA) for a slowly-varying gas [14],
and converts it, order-by-order, into a density functional approximation (DFA)
[2]. The lowest order is Thomas-Fermi theory, which is local in either the
density or the potential. For a slowly-varying gas, the procedure is extended
to include second-order gradients in either variable. The oscillating
contributions captured by the PFAs derived here depend on integrals over the
entire system, i.e., they have no known analog as density functionals.

Before discussing slab PFAs, our main topic in this paper, we will briefly
discuss PFAs in 1D. This will give us the opportunity to introduce our
notation and semiclassical jargon which may be unfamiliar to nonspecialists.
Moreover for non-interacting particles, one dimension is special, because the
semiclassical expansion of the eigenvalues can be performed to arbitrary order
(the celebrated WKB expansion [36]). The relation to DFT was first noted in
1956, when (non-interacting) TF theory was derived by summing WKB eigenvalues
[37]. Lastly we will use 1D quantities to construct our slab functionals.

Just as in Ref. [24], we define a quantum action as a function that yields the
eigenvalues $\epsilon_{j}$ via [36]

$s(\epsilon_{j})=j+\nu,\qquad j=0,1,\cdots,$ (1)

where $\nu$ is the Maslov index. In this paper we only consider smooth single
potential wells so $\nu=1/2$. Expanding $s$ in powers of $\hbar$ yields a
series expansion for $s$ in even gradients of the potential [36]:

$s=s^{(0)}+\Delta s^{(2)}+\Delta s^{(4)}+\cdots.$ (2)

We use atomic units ($\hbar=m_{e}=1$), so energies are in Hartrees and
distances are in Bohr radii. In these units the lowest order action (also
called the classical action) at energy $\mu$ is

$s^{(0)}(\mu)=\int dx\,\frac{p_{\scriptscriptstyle\rm F}(x)}{\pi},$ (3)

where $p_{\scriptscriptstyle\rm F}(x)=\sqrt{2[\mu-v(x)]}$ is the classical
momentum of a particle with energy $\mu$. All spatial integrals throughout
this paper are understood to run only between the classical turning points
satisfying $v(\pm x_{\scriptscriptstyle\rm F})=\mu$, as we consider only
symmetric single well potentials $v(x)$. Inserting Eq. (3) into Eq. (1) yields
the familiar WKB eigenvalues found in introductory textbooks on quantum
mechanics [38].

The higher-order corrections are well-defined, although less well-known. A
useful way of writing them was given in Ref. [39]. Define

$\displaystyle\begin{split}I(\mu)&=\int
dx\,\frac{v^{\prime\prime}(x)p_{\scriptscriptstyle\rm F}(x)}{8\pi},\\\
J(\mu)&=\int
dx\frac{7v^{\prime\prime}(x)^{2}-5v^{(4)}(x)p_{\scriptscriptstyle\rm
F}(x)^{2}}{\pi\,p_{\scriptscriptstyle\rm F}(x)},\\\ \end{split}$ (4)

where $v^{(j)}(x)=d^{j}v/dx^{j}$. $I(\mu)$ has a second-order gradient and
$J(\mu)$ has fourth order gradients. Then the next corrections can be written
as

$\Delta s^{(2)}(\mu)=-\frac{I^{\prime\prime}(\mu)}{3},\qquad\Delta
s^{(4)}(\mu)=\frac{J^{\prime\prime\prime}(\mu)}{5760},$ (5)

where the derivatives on $I$ and $J$ are with respect to $\mu$. This
particular way of writing the higher order action terms avoids dealing with
divergences at the turning points (compared to the forms in Ref. [36]).

The 1D number staircase $n(\mu)$ gives the number of eigenstates below $\mu$:

$n(\mu)=\sum_{j=0}^{\infty}\Theta(\mu-\epsilon_{j}),$ (6)

where $\Theta(x)$ is the Heaviside theta function. In terms of the quantum
action [25]

$n(\mu)=s(\mu)-\expectationvalue{s(\mu)},$ (7)

where $\expectationvalue{y}=y-{\lfloor}y+1/2{\rfloor}$ is periodic,
oscillating between $\pm 1/2$, and ${\lfloor}y{\rfloor}$ is the integer part
of $y$.

The exact 1D total energy can be written as a potential functional using the
1D number staircase:

$E(\mu)=\mu\,n(\mu)-\int_{0}^{\mu}d\epsilon\,n(\epsilon).$ (8)

Using the asymptotic expansion of $s$ in Eq. (2) we can construct the
asymptotic expansion of the 1D number staircase. Later we will use this
expansion of $n(\mu)$ to construct order by order approximations for our slab
potential functionals. We define the asymptotic expansion approximation (AEA)
to order $M$ (AEAM) as

$n^{(M)}(\mu)=s^{(M)}(\mu)-\expectationvalue{s^{(M)}(\mu)},\qquad{\rm AEAM}$
(9)

i.e., it includes terms up to $M$-th order in both the smooth and oscillating
pieces. Because the oscillating contribution is bounded in the LS limit, where
$j\to\infty$ in Eq. (1), it is typically of lower order than the smooth
contribution. Thus we define an alternative to AEA, denoted AEA’, as

$n^{(M)}(\mu)=s^{(M)}(\mu)-\expectationvalue{s^{(M-2)}(\mu)},\qquad{\rm
AEAM^{\prime}}$ (10)

which is easier to calculate and is correct to the same order. However, any
derivative with respect to $\mu$ may not be.

Figure 2: The eigenvalues of a slab with a Pöschl-Teller well of depth 6 as a
function of the parallel wave vector $K$.

We have shown in Eq. (8) how the 1D number staircase determines the total
energy. It also determines the kinetic and potential energy and the density.
If you have $n(\mu)$ you have all the information you need about a system.
This is not surprising since Eq. (6) shows that knowing $n$ is equivalent to
knowing all the eigenvalues. The structure of $n(\mu)$, in particular its
oscillatory term, is very interesting. If we ignore the oscillating term
entirely and approximate $n(\mu)$ with the lowest-order action, $n(\mu)\approx
s^{(0)}(\mu)$, we recover 1D TF theory. If we use a higher-order action,
$n(\mu)\approx s^{(M)}(\mu)$, we recover the 1D GEA to $M$-th order [15]. Thus
TF theory corresponds to the special case $M=0$, without the oscillating
contribution. Here we derive the GEA as a potential functional, but it is
always possible to invert it into the traditional density functional GEA.
However, the failure of the GEA, when applied to finite systems, is well known
[40]. In this case (single well 1D potentials), we can easily see why the GEA
fails: it neglects the oscillatory term $\expectationvalue{s(\mu)}$. By
including this term we derive PFAs which (unlike the GEA) yield the correct
asymptotic expressions. The resulting PFAs are much more accurate, and they
capture the derivative discontinuity. But converting PFAs into DFAs is much
more difficult when these oscillatory terms are included, and it may not even
be possible. The lesson is clear: the GEA fails to produce correct asymptotics
because it ignores this oscillatory piece. Modern generalized gradient
approximations (GGAs) that change the GEA coefficient [17] cannot compensate
for this missing oscillatory correction in general, so either their accuracy
or their range of applicability is limited. Restoration of the oscillations is
crucial on a path to a systematic, asymptotically correct, expansion of
functional approximations.

##  III Theory

Now we apply this technology to a three-dimensional non-interacting problem, a
slab with a potential profile $v(x)$, but uniform and infinite in the other
two directions. This geometry allows comparison with any existing orbital-free
kinetic energy functional [17]. By choosing a $v(x)$ that is the sum of two
potential wells (a dimer), separated by a distance $R$, we mimic certain
aspects of more realistic calculations, such as binding energies as a function
of bond lengths in diatomics [41]. We will only examine symmetric single well
potentials, $v(x)=v(-x)$, whose derivatives are finite up to at least the
fourth derivative. For simplicity we fix $v(0)=0$.

The eigenstates are two-dimensional free electron bands, with energies
$\epsilon_{j}(K)=\epsilon_{j}+K^{2}/2$, where $\epsilon_{j}$ is the $j$-th
eigenvalue in the 1D well, and $K$ is the parallel wavevector. Figure 2 shows
the simple shape of the PT slab bands, which are free-electron like in the two
directions perpendicular to $v(x)$ and begin at each of the eigenvalues of the
1D well. The system is a band metal. Summing over parallel directions in the
continuum limit, with double occupation, yields the number of electrons per
unit area from the 1D number staircase:

$N(\mu)=\int_{0}^{\mu}\frac{d\epsilon}{\pi}n(\epsilon),$ (11)

i.e., the density of states per unit area, $dN/d\mu$, is simply proportional
to the 1D number staircase. The slab energy (per unit area) is then

$E(\mu)=\mu\,N(\mu)-\int_{0}^{\mu}d\epsilon\,N(\epsilon).$ (12)

From now on all energy components are per unit area. Since $N(\mu)$ is
completely determined by $n(\mu)$ the slab $E(\mu)$ is also completely
determined by $n(\mu)$. Figure 1 shows its kinetic contribution. The density
of states, d$N/d\mu$, consists of constants (2D uniform gas), with steps up at
each $\epsilon_{j}$. We populate our slabs with non-interacting KS electrons.
Inversion of $N(\mu)$ then yields $E(N)$, either for the exact functional or
any approximation. The density is never calculated if not needed for other
purposes.

Inserting $s^{(0)}$ and dropping the oscillatory term in Eq. (7), and then
plugging this $n(\mu)$ into Eqs. (11) and (12) yields the familiar (non-
interacting) TF theory to lowest order:

$\displaystyle\begin{split}N^{\rm
TF}(\mu)&=\int\frac{dx}{3\pi^{2}}\,p_{\scriptscriptstyle\rm F}(x)^{3},\\\
E^{\rm
TF}(\mu)&=\int\frac{dx}{3\pi^{2}}\,\left[\frac{3}{10}p_{\scriptscriptstyle\rm
F}(x)^{2}+v(x)\right]p_{\scriptscriptstyle\rm F}(x)^{3}.\\\ \end{split}$ (13)

Both of these are potential functionals, and the first term in the energy is
the kinetic contribution. As we typically solve problems for a given $N$, one
inverts to find $\mu^{\rm TF}(N)$ and then $E^{\rm TF}(N)=E^{\rm TF}(\mu^{\rm
TF}(N))$.

To extract the density from such a PFA, we take the functional derivative of
the energy functional with respect to the potential while keeping the particle
number $N$ fixed:

$\rho^{\rm TF}(x)=\frac{\delta E^{\rm TF}}{\delta
v(x)}\bigg{|}_{N}=\frac{p_{\scriptscriptstyle\rm F}(x)^{3}}{3\pi^{2}}.$ (14)

We could also use this to eliminate $p_{\scriptscriptstyle\rm F}(x)$ from the
above expressions, thereby converting to the usual density functional
approximations. The same equations result from starting with the kinetic
energy as a density functional, and finding the density by minimizing the
total energy for fixed $N$.

A principal achievement of this work is the derivation of the leading post-TF
corrections. There is no correction to first order, but even orders appear to
have both smooth and oscillatory terms, while third order is purely
oscillatory. A full careful derivation, which we will give elsewhere, involves
subtle interchanges of orders of limits to avoid apparent divergences, but a
non-rigorous derivation is straightforward. The smooth contributions in Eq.
(7) from $s(\mu)$ can be easily integrated, as they are given by explicit
derivatives in Eq (5). The subtleties occur in the oscillating contributions.
These can be treated by replacing the integration variable $\epsilon$ in Eq.
(11) with the action itself. The change of variables allows the integrand to
be evaluated by elementary means (as in Eq. (16) of Ref. [25]), but only if
$s^{\prime}(\mu)$ is known as a function of $s(\mu)$. This smooth function can
be found, order-by-order, by expanding for large $\mu$. The resulting
expression for the slab staircase is

$\Delta N^{(2)}=-\frac{I^{\prime}(\mu)}{3\pi}+\frac{q}{2\tau},$ (15)

where

$\tau=\pi s^{(0)\prime}(\mu)=\int\frac{dx}{p_{\scriptscriptstyle\rm F}(x)},$
(16)

is the classical time to cross the well at energy $\mu$, and

$q=\frac{1}{12}-\expectationvalue{s}^{2},$ (17)

is a periodic function of the action, originating from an integral
proportional to the linear term in Eq. (16) of Ref. [25]. The smooth
contribution can be isolated by taking $\tau\rightarrow\infty$. These smooth
terms yield the standard gradient expansion when converted to density
functionals.

The total energy can then be deduced via the same procedure we used to find
$N(\mu)$, but applied to Eq. (12), yielding

$\Delta E^{(2)}=\mu\,\Delta N^{(2)}+\frac{I}{3\pi}.$ (18)

Given an expression for the total energy, one can extract the potential
contribution by taking a derivative with respect to the magnitude of $v(x)$.
We do this order-by-order in the 1D expressions, and find the slab
contribution via

$V=\int_{0}^{\mu}\frac{d\epsilon}{\pi}V^{\rm 1D}(\epsilon),$ (19)

which can be subtracted from the total energy, to yield

$\Delta T^{(2)}=-\frac{I}{6\pi}+\frac{\pi s^{(0)}q}{4\tau^{2}}.$ (20)

These expressions apply to any smooth symmetric $v(x)$ with one minimum and
two turning points. As we illustrate below, these are the exact asymptotic
corrections to TF theory for such problems. The red curve in Fig. 1 is the
AEA2 kinetic energy.

The process is repeated to find further higher-order terms, in which the
lengthy algebra and careful avoidance of apparent divergences, will be
reported elsewhere. Here we simply report the results:

$\displaystyle\begin{split}\Delta
N^{(3)}&=\frac{\pi\tau^{\prime}}{6\tau^{3}}h,\qquad\Delta E^{(3)}=\mu\Delta
N^{(3)}-\frac{\pi h}{6\tau^{2}},\\\ \Delta
T^{(3)}&=\frac{\pi^{2}s^{(0)}\tau^{\prime}}{4\tau^{4}}h,\\\ \end{split}$ (21)

where all primes are derivatives with respect to $\mu$ and

$h=\expectationvalue{s}\left(q+\frac{1}{6}\right).$ (22)

The fourth order corrections are

$\displaystyle\begin{split}\Delta
N^{(4)}&=\frac{J^{\prime\prime}}{5760\pi}+\frac{\pi
I^{\prime\prime\prime}}{6\tau^{2}}q-\frac{\pi^{2}w}{2\tau^{4}}\left[3\frac{\tau^{\prime
2}}{\tau}-\tau^{\prime\prime}\right],\\\ \Delta E^{(4)}&=\mu\,\Delta
N^{(4)}-\frac{J^{\prime}}{5760\pi}+\frac{3\pi^{2}\tau^{\prime}}{2\tau^{4}}w,\\\
\Delta T^{(4)}&=-\frac{J^{\prime}}{11520\pi}+\frac{\pi
q}{12\tau^{2}}\left[I^{\prime\prime}+\frac{2\pi
s^{(0)}I^{\prime\prime\prime}}{\tau}\right]+\\\
&\frac{\pi^{2}w}{4\tau^{4}}\left[3\tau^{\prime}+4\pi
s^{(0)}\frac{\tau^{\prime\prime}}{\tau}-15\pi
s^{(0)}\left(\frac{\tau^{\prime}}{\tau}\right)^{2}\right],\\\ \end{split}$
(23)

with $w=(7-240g)/2880$ and

$g=\expectationvalue{s}^{2}\left(q+\frac{5}{12}\right).$ (24)

We note that our separation of terms into smooth and oscillatory is based on
their origin in the derivation, rather than any process averaging over energy
levels.

##  IV Pöschl-Teller slabs

Our next step is to apply these formulas to a specific problem. Our
prototypical well is the Pöschl–Teller well, with $v(x)=D\,\tanh^{2}x$. Such
potentials are used in semiconductor physics [42, 43], and yield analytic
formulas. Here $s^{(0)}=\sqrt{2D}c$, $\Delta s^{(2)}=1/(8\sqrt{2D})$, and
$\Delta s^{(4)}=-1/(256\sqrt{2D^{3}})$, where $c=1-\sqrt{1-\mu/D}$. The TF PT
particle number and energies are

$\displaystyle\begin{split}N^{\rm
TF}&=\frac{\sqrt{2D^{3}}c^{2}}{\pi}\left(1-\frac{2}{3}c\right),\\\ E^{\rm
TF}&=\frac{\sqrt{2D^{5}}c^{3}}{\pi}\left(\frac{4}{3}-\frac{3c}{2}+\frac{2c^{2}}{5}\right),\\\
T^{\rm TF}&=\frac{3}{2}(\mu N^{\rm TF}-E^{\rm TF}).\\\ \end{split}$ (25)

The 2nd order contributions from Eqs. (15), (18), and (20) are

$\displaystyle\begin{split}\Delta
N^{(2)}&=\frac{\sqrt{2D}}{48\pi}[c\,(4-3c)-24\expectationvalue{s}^{2}(1-c)],\\\
\Delta
T^{(2)}&=-\frac{\sqrt{2D^{3}}}{192\pi}[c^{2}(4-3c)+96\expectationvalue{s}^{2}(1-c)^{2}]c,\\\
\Delta E^{(2)}&=\mu\left[\Delta
N^{(2)}+\frac{\sqrt{2D}}{96\pi}(4-6c+3c^{2})\right].\\\ \end{split}$ (26)

The third order corrections from Eq. (21) are

$\displaystyle\begin{split}\Delta N^{(3)}&=\frac{h}{6\pi},\\\ \Delta
E^{(3)}&=\mu\Delta N^{(3)}-\frac{Dh}{3\pi}(1-c)^{2},\\\ \Delta
T^{(3)}&=\frac{Dh}{2\pi}c(1-c).\\\ \end{split}$ (27)

The fourth order corrections from Eq. (23) are

$\displaystyle\begin{split}\Delta
N^{(4)}&=\frac{2-6c+3c^{2}}{768\pi\sqrt{2D}},\\\ \Delta E^{(4)}&=\mu\Delta
N^{(4)}+\frac{\sqrt{2D}}{2\pi}\beta,\\\ \Delta
T^{(4)}&=\frac{\sqrt{2D}}{4\pi}\left[\beta-\frac{(1-c)^{2}q+24cw}{8}\right],\\\
\end{split}$ (28)

where

$\beta=\frac{24-40c+80c^{2}-60c^{3}+15c^{4}}{7680}+3w(1-c).$ (29)

Figure 3 shows the number staircase (integrated density of states) for a given
PT slab. The function is rather smooth, making it difficult to see differences
between approximations. However, there are kinks in the exact curve whenever a
new band begins to be occupied. Both the TF and GEA2 curves have no such
kinks.

Figure 3: The PT slab number staircase with $D=6$ (black), and its TF (blue),
GEA2 (gray), AEA2’ (orange), and AEA2 (red) approximations. The lower panel
shows the errors in the upper panel.

Figure 3 also plots the errors (defined as approximate minus exact) in the
number staircase of the various approximations. The TF and GEA2 error curves
have kinks because they are smooth, while the exact curve is not. The orange
curve is AEA2’, which only accounts for the leading behavior of the
oscillating term, while the red curve is AEA2, which includes the next
contribution to the oscillating term. Their difference becomes negligible for
sufficiently large $\mu$ (both are asymptotically correct), but AEA2 clearly
has smaller errors for small $\mu$.

Figure 4: The, $D=6$, PT slab density of states $dN/d\mu$ (black), and its TF
(blue), GEA2 (gray), AEA2’ (orange), and AEA2 (red) approximations.

Figure 4 shows the density of states of a particular PT slab. This is just the
derivative of the number staircase given in Eq. (11) and shown in Fig. 3. Both
TF and GEA2 yield smooth approximations to it, and miss the discrete steps
(the origin of the infamous DFT derivative discontinuity [44]). Unlike how it
is treated in many semiclassical works [31], the smooth curve is not
synonymous with the TF contribution, as GEA2 makes a small but finite
correction. The asymptotic expansion approximation contains approximate steps,
with approximations to the plateau inbetween. The exact density of states
jumps discontinuously when $\mu=\epsilon_{j}$, where the $\epsilon_{j}$ are
the exact 1D eigenvalues. Using the definition of the saw-tooth function
$\expectationvalue{x}=x-{\lfloor}x+1/2{\rfloor}$, we can see that the AEA2’
approximation jumps whenever

$s^{(0)}(\mu)=j+\frac{1}{2},\qquad j=0,1,2,\cdots.$ (30)

This is just the lowest order WKB quantization rule for a single 1D well [38].
This means that AEA2’ jumps when $\mu=\epsilon^{(0)}_{j}$, the $j$th WKB
eigenvalue. Similar analysis shows that AEA2 jumps when

$s^{(2)}(\mu)=j+\frac{1}{2},\qquad j=0,1,2,\cdots,$ (31)

which is just the second order WKB quantization rule from Eq. (1). Thus AEA2
jumps when $\mu=\epsilon^{(2)}_{j}$, the second order WKB eigenvalue. AEA2 is
much more accurate, as the lower panel of Fig. 4 shows, because
$\epsilon^{(2)}_{j}$ is a better approximation to $\epsilon_{j}$ than
$\epsilon^{(0)}_{j}$. The inaccuracies in both second order AEAs vanish as
$\mu$ becomes large. Neither curve is quite flat, but AEA2 is flatter than
AEA2’.

Figure 5: Upper panel: Densities for $M=1$ (orange), 2 (blue), 4 (red), and
their TF limiting value (black). The areas under the curves are TF (0.0879), M
= 1 (0.0809), M = 2 (0.0854), M = 4 (0.0871). Lower panel: The deviation of
the exact scaled densities from the scaled TF density. The exact chemical
potential is $\mu$ and $\mu_{\rm{TF}}$ is its TF approximation.

$\begin{array}[]{|c|r|r|r|c|c|c|c|r|c|c|c|c|}\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Errors
(mH)}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Potential
Functionals}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Density
Functionals}\hfil\lx@intercol\vrule\lx@intercol\\\ \hline\cr
M&\lx@intercol\hfil D\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil
N\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil
T/N\hfil\lx@intercol\vrule\lx@intercol&\text{TF}&\text{GEA2}&\text{AEA2'}&\text{AEA2}&\lx@intercol\hfil\text{AEA4}\hfil\lx@intercol\vrule\lx@intercol&\text{TF}&\text{GEA2}&\text{MGE2}&\text{GEA4}\\\
\hline\cr 1&12.685&1.293&3.059&-87&-126&-29&-2.74&0.04097&-156&-41&-8&-2\\\
2&36.000&6.525&8.728&-85&-125&-14&-0.92&0.00487&-159&-35&1&-6\\\
3&70.971&18.318&17.233&-85&-125&-9&-0.46&0.00124&-162&-31&7&-7\\\
4&117.599&39.293&28.573&-84&-125&-7&-0.28&0.00045&-164&-28&11&-6\\\
5&175.883&72.075&42.748&-84&-125&-6&-0.18&0.00020&-165&-26&14&-6\\\
6&245.824&119.288&59.757&-84&-125&-5&-0.13&0.00010&-166&-25&16&-6\\\
7&327.422&183.555&79.602&-84&-125&-4&-0.10&0.00006&-167&-24&18&-6\\\
8&420.677&267.500&102.281&-84&-125&-3&-0.08&0.00003&-168&-23&19&-6\\\
9&525.589&373.746&127.795&-84&-125&-3&-0.06&0.00002&-168&-22&21&-5\\\
10&642.157&504.918&156.145&-84&-125&-3&-0.05&0.00001&-169&-21&22&-5\\\
\hline\cr\end{array}$

Table 1: Kinetic energy per particle for Pöschl-Teller slabs with
$D=2\epsilon_{M}$ and $\mu=D/2$. Here AEA$p$ is the $p$-th order asymptotic
expansion approximation, GEA$p$ is its analog without oscillating terms, while
MGE2 is a generalized gradient approximation (GGA) [17]. The rightmost columns
are density functionals evaluated on the exact density. Table 4 in the
appendix shows the analogous results for the total energy.

##  V Numerical results

In Table 1 we report results for a PT slab [45]. We approach the LS limit by
deepening the well while keeping $\mu/D=1/2$ fixed. Values of $D$ are chosen
where a new band begins to fill at $\mu=\epsilon_{j}$. Fig. 6 shows that our
approximations are poorest when $\mu$ is close to a 1D eigenvalue, but still
beat the competition. For each PFA, $N(\mu)$ is found to the given order and
inverted numerically to yield $E[\mu(N)]=E(N)$. Here $T/N$ is the exact
kinetic energy per particle in Hartrees, but errors are in milliHartrees. In
this paper errors are defined as approximate minus exact. The TF PFA error
starts at about 3%, and shrinks to less than 0.1%, consistent with the LS
theorem [22, 23]. TF is particularly good for our slab, due to the smoothness
of the density in the $x$-direction and uniformity in the other two (see Fig.
5). The (2nd-order) gradient expansion (GEA2) PFA worsens the result, but AEA2
has much smaller errors, especially for larger $D$. In fact, errors are less
than a mH for all but the shallowest well.

Figure 6: Kinetic energy errors for a shallow PT slab ($D=3$). Legend: TF
(black), GEA2 (gray), MGE2 (magenta), AEA2’ (orange), AEA2 (red), AEA4 (blue).
Dashed lines denote density functionals acting on the exact density.

$\begin{array}[]{|c|r|r|c|c|c|c|c|c|c|}\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Errors
(mH)}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil
e_{RE}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{$\mu$}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr M&\lx@intercol\hfil
e_{RE}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\mu\hfil\lx@intercol\vrule\lx@intercol&\text{TF}&\text{AEA2}&\text{AEA4}&\text{TF}&\text{AEA2'}&\text{AEA2}&\text{AEA4}\\\
\hline\cr 1&5.557&6.342&-63&-3&0.185&-242&-41&0.010&-0.0011911\\\
2&17.607&18.000&-172&-131&0.296&-239&-21&0.013&-0.0001006\\\
3&35.224&35.486&-204&-166&0.159&-238&-14&0.009&-0.0000218\\\
4&58.603&58.799&-217&-180&0.095&-237&-11&0.006&-0.0000072\\\
5&87.784&87.941&-224&-187&0.062&-237&-9&0.004&-0.0000030\\\
6&122.781&122.912&-227&-191&0.044&-237&-7&0.003&-0.0000015\\\
7&163.599&163.711&-230&-194&0.033&-237&-6&0.003&-0.0000008\\\
8&210.240&210.339&-231&-196&0.025&-237&-5&0.002&-0.0000005\\\
9&262.707&262.794&-232&-197&0.020&-237&-5&0.002&-0.0000003\\\
10&321.000&321.079&-233&-198&0.016&-237&-4&0.001&-0.0000002\\\
\hline\cr\end{array}$

Table 2: The removal energies, $e_{RE}$, and the chemical potentials for the
slabs of Table 1. The GEA2, AEA2’, and AEA2 total energy PFA’s, as functions
of $N$, are identical and thus yield the same $e_{RE}$ values.

The right-hand set of columns are DFAs evaluated on the exact density.
Unusually, TF does better self-consistently than on the exact density [46,
47], but GEA2 is improved. The MGE2 [17] (gradient expansion with coefficient
modified to recover asymptotics of atoms) is better for small wells, but
comparable to GEA2 in the asymptotic limit. GEA4 is better, but not nearly as
good as AEA2. AEA4 is the best approximation and always yields sub-
milliHartree errors, which are microHartree for $M>3$.

This table only tests one value of $\mu/D$, but similar trends hold for any
value. In Fig. 6, we plot errors for a shallow PT slab, $D=3$, as a function
of $N$, far from the LS limit. Our 2nd order approximations are typically
better than GEA2 and even MGE2 in this range, and all are beaten by fourth
order. Both this figure and Fig. 1 show that the AEA includes cusps in the
kinetic energy as a function of $N$, the infamous derivative discontinuity
that is the source of many errors in modern DFT calculations [29, 30]. Cusps
occur when $\mu=\epsilon_{j}$ both exactly and in approximations. The
resemblance with Fig. 9 of Ref. [48] (non-interacting electrons in a Coulomb
potential) is striking.

Figure 5 shows exact densities, from Table 1, and their TF approximation. We
scaled the densities so that the TF density is the same for all values of $M$.
As $M$ increases we approach the semiclassical limit and these densities
weakly approach their TF counterpart. Because the chemical potential (relative
to well-depth) is held fixed but the particle number is not (unlike in Fig. 2
of Ref. [49]), the normalization changes, but approaches that of TF in the LS
limit. The lower panel simply shows the differences from the TF curve, making
the weak approach to TF evident. Here, weak means that the integral over any
well-behaved function times the exact density approaches its TF counterpart
[50].

##  VI Energy differences

While it is important to demonstrate the validity of our expressions on total
energies, essentially all useful DFT calculations are of energy differences,
such as ionization potentials and binding energies of molecules (infinitesimal
differences determine bond lengths and lattice parameters). In the asymptotic
behavior of ionization potentials of neutral atoms [51], local exchange in a
KS calculation was found to match the exact result (Hartree-Fock), even
capturing variations across a row of the (very extended) non-relativistic
periodic table. Moreover, the average over such a row matched that of extended
interacting TF theory.

###  VI.1 Removal Energies

Figure 7: Upper panel: The exact energy (black) and its AEA2 (red) and AEA4
(blue) approximations (with the TF energy subtracted), for the $M=5$ PT slab
in Table 1. The middle panel shows a zoomed in view of the upper panel. The
gray line marks the particle number we examined in Table 1. Lower panel: The
exact chemical potential for a PT slab with $D=10$ (black), $dE^{\rm AEA2}/dN$
(red), and $\mu^{\rm AEA2}(N)$ (magenta). We have subtracted the TF value from
all curves.

In Table 2 we calculate electron removal energies (RE) in two distinct ways.
To simulate electron removal in a molecule, we introduce $e_{RE}$: the energy
per electron to remove $0.5$ electrons (each per unit area). In fact $\mu$ is
the true ionization potential of a metal, which $e_{RE}$ approaches in the LS
limit. The left-hand errors are those found from energy differences, in which
AEA2 is surprisingly poor. This is because AEA2 for $E(N)$ has no oscillations
(the cusps in Fig. 1 are in $T$, which does). On the other hand, $\mu^{\rm
AEA2}(N)$ does contain quantum oscillations, as shown in the lower panel,
which plots the exact chemical potential and two approximations to it. The
right-hand side of the table shows much better results from $\mu^{\rm
AEA2}(N)$, illustrating the subtleties of derivatives of the oscillating
terms. The derivatives of oscillations are typically one order larger than
those of the corresponding smooth terms.

Finally, Fig. 7 shows the exact and two approximate energy curves, each with
the TF curve subtracted, as a function of $N$. The black curve (exact) is
oscillating, as is the blue (AEA4) approximation. But the red curve (AEA2)
contains no oscillations, as the oscillations cancel out of the total energy
curve to 2nd order [which means that $E^{\rm
AEA2}(N)=E^{\rm{AEA2^{\prime}}}(N)$]. The middle panel in Fig. 7 zooms in on
the energy curve in the region of the particle number corresponding to
$\mu=D/2$, marked by the vertical line. Clearly the oscillations play a large
role in the energy change if you remove 1/2 an electron when $\mu=D/2$. The
red curve yields a very poor approximation to this energy difference. The
exact energy satisfies $E^{\prime}(N)=\mu$, but $dE^{\rm AEA2}/dN$ yields a
poor approximation to the chemical potential. Instead we derive $\mu^{\rm
AEA2}(N)$ by inverting $N^{\rm AEA2}(\mu)$. The blue AEA4 curve will clearly
yield almost exact answers.

###  VI.2 Binding energies

Our last test is the most stringent. We consider the sum of two identical PT
slabs (the PT dimer), each of well depth $D=3$, as a function of their
separation, $R$, to mimic the binding energy curve of a molecule. From 0 to
$R_{c}=2\,\text{asech}\sqrt{2/3}=1.31696$ the depth of each dimer is given by
$\mathcal{D}=6\,\sech^{2}(R/2)$ and there is a single minimum. Beyond $R_{c}$,
there is a double well and the formula for the total well depth becomes more
complicated (see Fig. 8). We keep $N$ fixed throughout. Table 3 shows the
result. For $R\leq R_{c}$, our results for energy differences are similar to,
but less accurate than, those for total energies in Table 1, with AEA4 still
yielding chemical accuracy. Figure 8 shows various potentials of the PT dimer
slabs as a function of their separation, given in units of the critical
separation ($R_{c}$) at which the second derivative of the midpoint potential
vanishes. Beyond this critical value, there are two wells, and our derivation
of the semiclassical asymptotic expansion no longer applies.

Figure 8: Several PT dimer potentials made from two PT slabs (each with $D=3$)
at various separations: $R/R_{c}=0$, 0.5, 1, 3 (black, blue, magenta, red),
where $R_{c}=2\,\text{asech}\sqrt{2/3}$ and $v(\pm\infty)=0$.

###  VI.3 Delocalization error

$\begin{array}[]{|c|c|c|c|c|c|r|r|c|c|c|}\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Errors
(mH)}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Potential
Functionals}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Density
Functionals}\hfil\lx@intercol\vrule\lx@intercol\\\ \hline\cr
R/R_{c}&T-2T_{A}&\text{TF}&\text{GEA2}&\text{AEA2'}&\text{AEA2}&\lx@intercol\hfil\text{AEA4}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{TF}\hfil\lx@intercol\vrule\lx@intercol&\text{GEA2}&\text{MGE2}&\text{GEA4}\\\
\hline\cr 0&0.792&27&28&32.2&6.1&-0.177&-23.8&-11&-7&-12\\\
0.1&0.784&28&30&32.1&6.1&-0.179&-23.0&-11&-7&-12\\\
0.2&0.760&31&34&31.7&6.2&-0.186&-20.6&-10&-7&-12\\\
0.25&0.743&33&37&31.4&6.3&-0.191&-18.8&-10&-7&-12\\\
0.3&0.722&35&41&31.0&6.2&-0.197&-16.8&-9&-7&-12\\\
0.4&0.672&38&48&30.0&6.1&-0.211&-11.8&-8&-7&-11\\\
0.5&0.613&41&55&28.5&5.7&-0.219&-6.0&-6&-7&-10\\\
0.6&0.547&43&60&26.6&5.0&-0.204&0.2&-5&-6&-9\\\
0.7&0.477&42&63&24.4&4.1&-0.141&6.5&-2&-5&-8\\\
0.75&0.441&41&64&23.3&3.7&-0.083&9.4&-1&-4&-7\\\
0.8&0.406&40&64&22.2&3.2&-0.001&12.3&0&-4&-6\\\
0.9&0.337&35&62&20.1&2.5&0.251&17.3&2&-3&-4\\\
1&0.271&30&59&18.4&2.1&0.644&21.3&3&-2&-3\\\
1.1&0.211&24&54&17.2&2.2&1.191&24.0&4&-2&-2\\\
1.2&0.157&18&48&16.8&2.8&1.876&25.3&4&-2&-1\\\
1.25&0.133&15&45&16.8&3.4&2.253&25.5&4&-2&-1\\\ \hline\cr\end{array}$

Table 3: The kinetic binding energies, $T-2T_{A}$, for a series of PT dimer
slabs with $3/\pi$ electrons per unit area. $T$ is the exact kinetic energy
and $T_{A}$ is the kinetic energy of an isolated PT slab ($D=3$). All
approximations are defined as in Table 1. We give more energy tables for these
dimers in the appendix.

While AEA4 yields chemical accuracy for any $R\leq R_{c}$, it fails to do so
beyond $R_{c}$. This may be because a different, more complicated asymptotic
expansion applies, one that accounts for the additional turning points below
$\mu$ (there are still only two turning points at $\mu$). The $R_{c}<R$ error
in AEA2 remains below that of $R=0$, although it is also growing slightly.
This is characteristic of an asymptotic expansion and presumably, for $R$
sufficiently large, AEA2 would fail. Eventually for $R$ very large, as the
transition to four turning points at $\mu$ is approached, even TF should fail.
This is highly analogous to the situation for XC [48], but beyond the scope of
the present work.

##  VII Conclusions

Our PFAs are the asymptotic expansion in $\hbar$ of the energy to a given
order, derived directly from the Schrödinger equation. Unlike any DFAs in use
today, they are the correct generalization of the gradient expansion (as a
PFA) to finite systems. Our results bring the proof-of-principle from earlier
1D studies closer to realistic orbital-free electronic structure calculations.
The leading terms in the semiclassical AEA achieve much greater accuracy than
modern orbital-free DFAs, and their quantum oscillations capture derivative
discontinuities. However, our results are limited to homogeneous slabs with
only two turning points.

For such homogeneous slab calculations, our PFAs can be tested immediately.
With any existing KS code and XC choice, find the self-consistent solution,
and then test our PFAs by feeding them the KS potential. Pseudopotentials are
smoother than full potentials, and so are likely to lead to higher accuracy
for a given order. To do this self-consistently, the XC potential requires the
ground-state density. There are highly accurate PFA’s for the 1D density [52,
53, 54], which can be generalized to slabs and then tested for sufficient
accuracy. Otherwise one can extract a PFA density via $\rho(x)=\delta E/\delta
v(x)$ for fixed $N$, a non-trivial derivation when oscillating terms are
included. Alternatively, the Lieb maximization principle,
$T[\rho]=\sup\\{E[v]-\int\rho\,v\\}$, where the supremum is over all
potentials [55, 56], converts a PFA into a DFA. This non-trivial derivation
for AEA2 has several caveats. Since the AEA is not variational, such a
maximization might not yield useful results (the bare AEA might need to be
generalized). Moreover, the resulting DFA might be less accurate than the
original PFA for finite systems.

Another interesting area is surfaces. As the interior region of a slab is
broadened, quantum oscillations will become less relevant, and the AEA should
reduce to the Airy-gas [57] approximation for $T_{s}$, and AEA2 should give
accurate surface kinetic energies. As for inhomogeneity in the perpendicular
directions, there are many approximations one could try (e.g., application to
surface averaged potentials, application in the direction of the local
gradient of the potential, and simple GGA-type generalizations of AEA
expressions) that could be constrained to recover our results in the uniform
limit. Another exciting area is to generalize the AEA to include complex
turning points when bonds are stretched [58]. This has not yet been done, even
in one-dimension. Real turning points determine the dominant form of the
semiclassical expansion, and complex turning points yield subdominant
corrections to this expansion.

Could these techniques be applied to the XC energy? Yes, but with difficulty.
Consider exchange. Our most accurate results come from fitting highly accurate
atomic data for atoms up to $Z=$ 120, because we lack analytic results [59].
Asymptotic expansions must be found for the double sum in the exchange energy,
and one also needs the semiclassical expansion for the KS density matrix with
real turning points [60]. A computational nirvana might be reached [24] if we
could overcome that one difficulty.

We thank the NSF (CHE-2154371) for funding. More rigorous derivations, exact
results for PT slabs, and limiting cases, will appear in future work.

## References

  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Dreizler and Gross [1990] R. M. Dreizler and E. K. Gross,  _Density Functional Theory: An Approach to the Quantum Many-Body Problem_, 1st ed. (Springer-Verlag Berlin Heidelberg, Berlin Heidelberg, 1990).
  * Parr and Yang [1994] R. G. Parr and W. Yang,  _Density-Functional Theory of Atoms and Molecules (International Series of Monographs on Chemistry)_ (Oxford University Press, USA, 1994).
  * Lignères and Carter [2005] V. L. Lignères and E. A. Carter, An introduction to orbital-free density functional theory, in  _Handbook of Materials Modeling: Methods_, edited by S. Yip (Springer Netherlands, Dordrecht, 2005) pp. 137–148.
  * Mi  _et al._ [2023] W. Mi, K. Luo, S. B. Trickey, and M. Pavanello, Orbital-free density functional theory: An attractive electronic structure method for large-scale first-principles simulations, Chemical Reviews 123, 12039 (2023).
  * Thomas [1927] L. H. Thomas, The calculation of atomic fields, Mathematical Proceedings of the Cambridge Philosophical Society 23, 542–548 (1927).
  * Fermi [1927] E. Fermi, Rend. Acc. Naz. Lincei 6 (1927).
  * Lieb [1983] E. H. Lieb, Density functionals for coulomb systems, International Journal of Quantum Chemistry 24, 243 (1983).
  * Kim  _et al._ [2013] M.-C. Kim, E. Sim, and K. Burke, Understanding and reducing errors in density functional calculations, Phys. Rev. Lett. 111, 073003 (2013).
  * Fabiano  _et al._ [2022] E. Fabiano, F. Sarcinella, L. A. Constantin, and F. Della Sala, Kinetic energy density functionals based on a generalized screened coulomb potential: Linear response and future perspectives, Computation 10, 10.3390/computation10020030 (2022).
  * Wesolowski and Wang [2013] T. A. Wesolowski and Y. A. Wang,  _Recent Progress in Orbital-free Density Functional Theory_ (World Scientific, 2013).
  * Witt  _et al._ [2018] W. C. Witt, B. G. del Rio, J. M. Dieterich, and E. A. Carter, Orbital-free density functional theory for materials research, Journal of Materials Research 33, 777 (2018).
  * Luo  _et al._ [2018] K. Luo, V. V. Karasiev, and S. B. Trickey, A simple generalized gradient approximation for the noninteracting kinetic energy density functional, Phys. Rev. B 98, 041111 (2018).
  * Kirzhnits [1957] D. Kirzhnits, Quantum corrections to the thomas-fermi equation, Sov. Phys. JETP 5, 64 (1957).
  * Šamaj and Percus [1999] L. Šamaj and J. K. Percus, Recursion representation of gradient expansion for free fermion ground state in one dimension, The Journal of Chemical Physics 111, 1809 (1999).
  * Yang [1986] W. Yang, Gradient correction in thomas-fermi theory, Phys. Rev. A 34, 4575 (1986).
  * Constantin  _et al._ [2011] L. A. Constantin, E. Fabiano, S. Laricchia, and F. Della Sala, Semiclassical neutral atom as a reference system in density functional theory, Phys. Rev. Lett. 106, 186406 (2011).
  * Lastra  _et al._ [2008] J. M. G. Lastra, J. W. Kaminski, and T. A. Wesolowski, Orbital-free effective embedding potential at nuclear cusps, The Journal of Chemical Physics 129, 074107 (2008).
  * Wang  _et al._ [1999] Y. A. Wang, N. Govind, and E. A. Carter, Orbital-free kinetic-energy density functionals with a density-dependent kernel, Phys. Rev. B 60, 16350 (1999).
  * Huang and Carter [2010] C. Huang and E. A. Carter, Nonlocal orbital-free kinetic energy density functional for semiconductors, Phys. Rev. B 81, 045206 (2010).
  * Shao  _et al._ [2021] X. Shao, W. Mi, and M. Pavanello, Revised huang-carter nonlocal kinetic energy functional for semiconductors and their surfaces, Phys. Rev. B 104, 045118 (2021).
  * Lieb and Simon [1973] E. H. Lieb and B. Simon, Thomas-fermi theory revisited, Phys. Rev. Lett. 31, 681 (1973).
  * Lieb and Simon [1977] E. H. Lieb and B. Simon, The thomas-fermi theory of atoms, molecules and solids, Advances in Mathematics 23, 22 (1977).
  * Berry and Burke [2020] M. V. Berry and K. Burke, Exact and approximate energy sums in potential wells, Journal of Physics A: Mathematical and Theoretical 53, 095203 (2020).
  * Burke [2020a] K. Burke, Leading correction to the local density approximation of the kinetic energy in one dimension, The Journal of Chemical Physics 152, 081102 (2020a).
  * Burke [2020b] K. Burke, Deriving approximate density functionals with asymptotics, Faraday Discuss. 10.1039/D0FD00057D (2020b).
  * Okun and Burke [2023] P. Okun and K. Burke, Semiclassics: The hidden theory behind the success of dft, in  _Density Functionals for Many-Particle Systems: Mathematical Theory and Physical Applications of Effective Equations_, edited by B.-G. Englert, M.-I. Trappe, and H. Siedentop (World Scientific Publishing Company, Singapore, 2023) pp. 179–249.
  * Maslov and Fedoriuk [2001] V. Maslov and V. Fedoriuk,  _Semi-Classical Approximation in Quantum Mechanics_ , Mathematical Physics and Applied Mathematics (Springer Netherlands, 2001).
  * Mori-Sánchez  _et al._ [2008] P. Mori-Sánchez, A. J. Cohen, and W. Yang, Localization and delocalization errors in density functional theory and implications for band-gap prediction, Phys. Rev. Lett. 100, 146401 (2008).
  * Kronik and Kümmel [2020] L. Kronik and S. Kümmel, Piecewise linearity, freedom from self-interaction, and a coulomb asymptotic potential: three related yet inequivalent properties of the exact density functional, Phys. Chem. Chem. Phys. 22, 16467 (2020).
  * Brack and Bhaduri [2003] M. Brack and R. Bhaduri,  _Semiclassical Physics_, 1st ed. (Routledge, Boca Raton, 2003).
  * Gutzwiller [1990] M. C. Gutzwiller,  _Chaos in Classical and Quantum Mechanics_ (Springer, New York, NY, 1990).
  * Cohen  _et al._ [2008] A. J. Cohen, P. Mori-Sánchez, and W. Yang, Insights into current limitations of density functional theory, Science 321, 792 (2008).
  * Cancio  _et al._ [2018] A. Cancio, G. P. Chen, B. T. Krull, and K. Burke, Fitting a round peg into a round hole: Asymptotically correcting the generalized gradient approximation for correlation, The Journal of Chemical Physics 149, 084116 (2018).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
  * Bender and Orszag [1999] C. M. Bender and S. A. Orszag,  _Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory_ (Springer, Verlag New York, 1999).
  * March and Plaskett [1956] N. H. March and J. S. Plaskett, The relation between the wentzel-kramers-brillouin and the thomas-fermi approximations, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 235, 419 (1956).
  * Griffiths [2005] D. J. Griffiths,  _Introduction to quantum mechanics_ (Pearson Prentice Hall, 2005).
  * Krieger  _et al._ [1967] J. B. Krieger, M. L. Lewis, and C. Rosenzweig, Use of the wkb method for obtaining energy eigenvalues, The Journal of Chemical Physics 47, 2942 (1967).
  * MA and BRUECKNER [1968] S.-k. MA and K. A. BRUECKNER, Correlation energy of an electron gas with a slowly varying high density, Phys. Rev. 165, 18 (1968).
  * Perdew  _et al._ [1995] J. P. Perdew, A. Savin, and K. Burke, Escaping the symmetry dilemma through a pair-density interpretation of spin-density functional theory, Phys. Rev. A 51, 4531 (1995).
  * Radovanović  _et al._ [2000] J. Radovanović, V. Milanović, Z. Ikonić, and D. Indjin, Intersubband absorption in pöschl–teller-like semiconductor quantum wells, Physics Letters A 269, 179 (2000).
  * Hartmann and Portnoi [2017] R. R. Hartmann and M. E. Portnoi, Two-dimensional dirac particles in a pöschl-teller waveguide, Scientific Reports 7, 11599 (2017).
  * Sagredo and Burke [2020] F. Sagredo and K. Burke, Confirmation of the pplb derivative discontinuity: Exact chemical potential at finite temperatures of a model system, Journal of Chemical Theory and Computation 16, 7225 (2020).
  * Nieto [1978] M. M. Nieto, Exact wave-function normalization constants for the ${B}_{0}tanhz-{U}_{0}{cosh}^{-2}z$ and pöschl-teller potentials, Phys. Rev. A 17, 1273 (1978).
  * Vuckovic  _et al._ [2019] S. Vuckovic, S. Song, J. Kozlowski, E. Sim, and K. Burke, Density functional analysis: The theory of density-corrected dft, Journal of Chemical Theory and Computation 15, 6636 (2019).
  * Sim  _et al._ [2022] E. Sim, S. Song, S. Vuckovic, and K. Burke, Improving results by improving densities: Density-corrected density functional theory, Journal of the American Chemical Society 144, 6625 (2022).
  * Burke  _et al._ [2016] K. Burke, A. Cancio, T. Gould, and S. Pittalis, Locality of correlation in density functional theory, The Journal of Chemical Physics 145, 054112 (2016).
  * Cangi  _et al._ [2010] A. Cangi, D. Lee, P. Elliott, and K. Burke, Leading corrections to local approximations, Phys. Rev. B 81, 235128 (2010).
  * Lieb [1976] E. H. Lieb, The stability of matter, Rev. Mod. Phys. 48, 553 (1976).
  * Constantin  _et al._ [2010] L. A. Constantin, J. C. Snyder, J. P. Perdew, and K. Burke, Communication: Ionization potentials in the limit of large atomic number, The Journal of Chemical Physics 133, 241103 (2010).
  * Ribeiro  _et al._ [2015] R. F. Ribeiro, D. Lee, A. Cangi, P. Elliott, and K. Burke, Corrections to thomas-fermi densities at turning points and beyond, Phys. Rev. Lett. 114, 050401 (2015).
  * Ribeiro and Burke [2017] R. F. Ribeiro and K. Burke, Leading corrections to local approximations. ii. the case with turning points, Phys. Rev. B 95, 115115 (2017).
  * Ribeiro and Burke [2018] R. F. Ribeiro and K. Burke, Deriving uniform semiclassical approximations for one-dimensional fermionic systems, The Journal of Chemical Physics 148, 10.1063/1.5025628 (2018).
  * Cangi  _et al._ [2013] A. Cangi, E. K. U. Gross, and K. Burke, Potential functionals versus density functionals, Phys. Rev. A 88, 062505 (2013).
  * Cangi  _et al._ [2011] A. Cangi, D. Lee, P. Elliott, K. Burke, and E. K. U. Gross, Electronic structure via potential functional approximations, Phys. Rev. Lett. 106, 236404 (2011).
  * Lindmaa  _et al._ [2014] A. Lindmaa, A. E. Mattsson, and R. Armiento, Quantum oscillations in the kinetic energy density: Gradient corrections from the airy gas, Phys. Rev. B 90, 075139 (2014).
  * Okun and Burke [2022] P. Okun and K. Burke, Asymptotics of eigenvalue sums when some turning points are complex, Journal of Physics A: Mathematical and Theoretical 55, 394003 (2022).
  * Argaman  _et al._ [2022] N. Argaman, J. Redd, A. C. Cancio, and K. Burke, Leading correction to the local density approximation for exchange in large-$z$ atoms, Phys. Rev. Lett. 129, 153001 (2022).
  * Elliott  _et al._ [2015] P. Elliott, A. Cangi, S. Pittalis, E. K. U. Gross, and K. Burke, Almost exact exchange at almost no computational cost in electronic structure, Phys. Rev. A 92, 022513 (2015).

## Appendix: Tables of energies

This appendix contains tables that supplement those in the main text.

$\begin{array}[]{|c|r|r|r|c|c|c|c|c|c|c|}\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Errors
(mH)}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Potential
Functionals}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Density
Functionals}\hfil\lx@intercol\vrule\lx@intercol\\\ \hline\cr
M&\vrule\lx@intercol\hfil
D\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil
N\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil
E/N\hfil\lx@intercol\vrule\lx@intercol&\text{TF}&\text{AEA2}&\text{AEA4}&\text{TF}&\text{GEA2}&\text{MGE2}&\text{GEA4}\\\
\hline\cr 1&12.685&1.293&4.312&-192&9.2&0.003781&-156&-41&-8&-2\\\
2&36.000&6.525&12.189&-190&3.1&0.000452&-159&-35&1&-6\\\
3&70.971&18.318&24.005&-189&1.6&0.000115&-162&-31&7&-7\\\
4&117.599&39.293&39.759&-189&0.9&0.000042&-164&-28&11&-6\\\
5&175.883&72.075&59.451&-189&0.6&0.000019&-165&-26&14&-6\\\
6&245.824&119.288&83.082&-189&0.4&0.000010&-166&-25&16&-6\\\
7&327.422&183.555&110.651&-189&0.3&0.000005&-167&-24&18&-6\\\
8&420.677&267.500&142.159&-189&0.3&0.000003&-168&-23&19&-6\\\
9&525.589&373.746&177.605&-189&0.2&0.000002&-168&-22&21&-5\\\
10&642.157&504.918&216.990&-189&0.2&0.000001&-169&-21&22&-5\\\
\hline\cr\end{array}$

Table 4: Same as Table 1 of the main text, but for the total energy.

Table 4 shows the total energies (not just the kinetic energies) of the
calculations in Table 1 of the main text. In this case, any functionals
evaluated on the exact density include the exact potential energy by
construction. Just as in the self-consistent TF calculation, we expect errors
on (some version of) self-consistent densities to be larger.

$\begin{array}[]{|c|c|c|c|r|r|r|c|c|c|c|}\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Errors
(mH)}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Potential
Functionals}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Density
Functionals}\hfil\lx@intercol\vrule\lx@intercol\\\ \hline\cr
R/R_{c}&T&\text{TF}&\text{GEA2}&\lx@intercol\hfil\text{AEA2'}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{AEA2}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{AEA4}\hfil\lx@intercol\vrule\lx@intercol&\text{TF}&\text{GEA2}&\text{MGE2}&\text{GEA4}\\\
\hline\cr 0&1.890&34&-19&11.8&3.1&-0.038&-77&-22&-5.5&-12\\\
0.1&1.882&35&-17&11.7&3.1&-0.041&-76&-21&-5.5&-12\\\
0.2&1.858&38&-13&11.4&3.2&-0.047&-74&-21&-5.6&-12\\\
0.25&1.841&40&-10&11.1&3.2&-0.052&-72&-21&-5.6&-12\\\
0.3&1.821&41&-6&10.7&3.2&-0.058&-70&-20&-5.6&-11\\\
0.4&1.771&45&1&9.6&3.1&-0.072&-65&-19&-5.4&-11\\\
0.5&1.711&48&8&8.1&2.6&-0.080&-59&-17&-5.0&-10\\\
0.6&1.645&50&13&6.2&2.0&-0.065&-53&-15&-4.3&-9\\\
0.7&1.575&49&17&4.1&1.1&-0.002&-47&-13&-3.4&-7\\\
0.75&1.540&48&17&2.9&0.6&0.056&-44&-12&-2.8&-6\\\
0.8&1.504&47&17&1.8&0.2&0.138&-41&-11&-2.3&-6\\\
0.9&1.435&42&16&-0.3&-0.6&0.390&-36&-9&-1.3&-4\\\
1&1.369&37&12&-2.0&-1.0&0.783&-32&-8&-0.6&-3\\\
1.1&1.309&31&7&-3.2&-0.9&1.330&-29&-7&-0.1&-2\\\
1.2&1.255&25&1&-3.6&-0.2&2.014&-28&-6&-0.2&-1\\\
1.25&1.231&22&-2&-3.6&0.3&2.392&-28&-6&-0.3&-1\\\ \hline\cr\end{array}$

Table 5: Same as Table 3 of the main text, but showing the total kinetic
energies, not the kinetic binding energies.

$\begin{array}[]{|c|r|c|r|r|r|c|c|c|}\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Errors
(mH)}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Potential
Functionals}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Density
Functionals}\hfil\lx@intercol\vrule\lx@intercol\\\ \hline\cr
R/R_{c}&\lx@intercol\hfil
E-2E_{A}\hfil\lx@intercol\vrule\lx@intercol&\text{TF}&\text{AEA2}&\lx@intercol\hfil\text{AEA4}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{TF}\hfil\lx@intercol\vrule\lx@intercol&\text{GEA2}&\text{MGE2}&\text{GEA4}\\\
\hline\cr 0&1.174&-39&-8.1&0.14&-23.8&-11&-7.1&-12\\\
0.1&1.162&-39&-8.2&0.14&-23.0&-11&-7.1&-12\\\
0.2&1.127&-36&-8.5&0.14&-20.6&-10&-7.2&-12\\\
0.25&1.100&-34&-8.6&0.14&-18.8&-10&-7.2&-12\\\
0.3&1.068&-32&-8.7&0.14&-16.8&-9&-7.2&-12\\\
0.4&0.988&-27&-8.7&0.13&-11.8&-8&-7.0&-11\\\
0.5&0.887&-21&-8.3&0.10&-6.0&-6&-6.6&-10\\\
0.6&0.767&-13&-7.4&0.04&0.2&-5&-5.9&-9\\\
0.7&0.630&-6&-6.0&-0.04&6.5&-2&-5.0&-8\\\
0.75&0.556&-2&-5.1&-0.09&9.4&-1&-4.5&-7\\\
0.8&0.479&2&-4.2&-0.14&12.3&0&-3.9&-6\\\
0.9&0.315&9&-2.3&-0.23&17.3&2&-3.0&-4\\\
1&0.143&15&-0.5&-0.28&21.3&3&-2.2&-3\\\
1.1&-0.036&19&0.7&-0.29&24.0&4&-1.8&-2\\\
1.2&-0.217&22&1.3&-0.21&25.3&4&-1.8&-1\\\
1.25&-0.309&22&1.2&-0.12&25.5&4&-1.9&-1\\\ \hline\cr\end{array}$

Table 6: Same as Table 3 of the main text, but with the total binding energy.

$\begin{array}[]{|c|c|c|r|r|c|c|c|c|}\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Errors
(mH)}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Potential
Functionals}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Density
Functionals}\hfil\lx@intercol\vrule\lx@intercol\\\ \hline\cr
R/R_{c}&E&\text{TF}&\lx@intercol\hfil\text{AEA2}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{AEA4}\hfil\lx@intercol\vrule\lx@intercol&\text{TF}&\text{GEA2}&\text{MGE2}&\text{GEA4}\\\
\hline\cr 0&2.845&-99&-7.2&0.15&-77&-22&-5.5&-12\\\
0.1&2.833&-98&-7.3&0.16&-76&-21&-5.5&-12\\\
0.2&2.798&-95&-7.6&0.16&-74&-21&-5.6&-12\\\
0.25&2.772&-94&-7.7&0.16&-72&-21&-5.6&-12\\\
0.3&2.740&-92&-7.8&0.16&-70&-20&-5.6&-11\\\
0.4&2.659&-86&-7.8&0.15&-65&-19&-5.4&-11\\\
0.5&2.558&-80&-7.4&0.12&-59&-17&-5.0&-10\\\
0.6&2.438&-73&-6.5&0.06&-53&-15&-4.3&-9\\\
0.7&2.301&-65&-5.1&-0.03&-47&-13&-3.4&-7\\\
0.75&2.227&-61&-4.3&-0.07&-44&-12&-2.8&-6\\\
0.8&2.150&-57&-3.3&-0.12&-41&-11&-2.3&-6\\\
0.9&1.986&-51&-1.4&-0.21&-36&-9&-1.3&-4\\\
1&1.814&-45&0.4&-0.26&-32&-8&-0.6&-3\\\
1.1&1.636&-40&1.6&-0.27&-29&-7&-0.1&-2\\\
1.2&1.454&-38&2.1&-0.19&-28&-6&-0.2&-1\\\
1.25&1.362&-37&2.1&-0.10&-28&-6&-0.3&-1\\\ \hline\cr\end{array}$

Table 7: Same as Table 3 of the main text, but with the total energy. We fix
the center of each well at 0: $v(0)=0$.

Tables 5-7 supplement Table 3 of the main text, showing total kinetic energies
of the PT dimer slabs, not just binding energies, so that approximations
cannot benefit from cancellation of errors between the PT dimer slabs and the
separated ‘atomic’ slabs. We also give the corresponding total energy and
binding energy.