#  Next-order correction to the Dirac exchange energy
of the free electron gas in the thermodynamic limit
and generalized gradient approximations

Thiago Carvalho Corso and Gero Friesecke Email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>Zentrum Mathematik, Technische
Universität München, Germany

###### Abstract

We derive the next order correction to the Dirac exchange energy for the free
electron gas in a box with zero boundary conditions in the thermodynamic
limit. The correction is of the order of the surface area of the box, and
comes from three different contributions: (i) a real-space boundary layer,
(ii) a boundary-condition-induced small shift of Fermi momentum and bulk
density, and (iii) a long-range electrostatic finite-size correction. Moreover
we show that the LDA, in addition to capturing the bulk term exactly, also
produces a correction of the correct order but not the correct size. GGA
corrections are found to be capable of capturing the surface term exactly,
provided the gradient enhancement factor satisfies a simple explicit integral
constraint. For current GGAs such as B88 and PBE we find that the new
constraint is not satisfied and the size of the surface correction is
overestimated by about ten percent. The new constraint might thus be of
interest for the design of future exchange functionals.

##  1 Introduction

We derive the next order correction to the Dirac exchange energy for the free
electron gas in a box with zero boundary conditions in the thermodynamic
limit. Like Dirac exchange, the correction is of significant interest for
density functional theory (DFT). In particular, it yields a novel exact
constraint on generalized gradient approximations (GGAs).

Motivation from DFT and main result. Original Kohn-Sham DFT [17] is based on
the paradigm that exchange and correlation effects in real many-electron
systems in solid-state and molecular physics are reasonably similar to those
in the three-dimensional uniform electron gas (UEG). This paradigm is
implemented in the form of the celebrated local density approximation (LDA):
one applies an analytical parametrization of the exchange-correlation energy
of the UEG as a function of the density111obtained from high- and low-density
asymptotics as well as numerical quantum Monte Carlo results locally to the
non-uniform density of the real system, thereby obtaining an approximation of
the system’s exact exchange-correlation energy. Thus the LDA corresponds to
the following approximate model for the exchange-correlation energy density at
a point $r$: $e_{xc}(r)=f(\rho(r))$, where $f(\bar{\rho})$ is the exchange-
correlation energy density of the UEG with density $\bar{\rho}$ and $\rho(r)$
is the density at $r$ of the inhomogeneous system.

A milestone in elevating DFT from a useful first approximation to a widely
applicable quantitative method was the development of generalized gradient
approximations (GGAs) [1, 19]. These model the effect of density
inhomogeneities by allowing the exchange-correlation energy density at a point
$r$ to depend not just on the local density $\rho(r)$ but also the local
density gradient, $e_{xc}(r)=f(\rho(r),|\nabla\rho(r)|)$. Deceptively simple
but highly judicious semi-empirical choices for $f$ such as those by Becke [1]
and Perdew, Burke, and Ernzerhof [19] improved the accuracy of total DFT
energies for small atoms and molecules from 1 eV to about 0.2 eV.

The dominating part of the exchange-correlation energy for real molecular and
solid-state systems consists of the exchange energy, on which we focus in this
paper. The LDA exchange energy density is given by Dirac exchange, i.e. the
exchange energy density of the free electron gas with density $\bar{\rho}$
(which is known explicitly, see below) applied to the local density. But what
about gradient corrections? Successful GGA models for exchange, such as the
widely used Becke-88 [1] and PBE [19] functionals, are not obtained from
uniform electron gas theory alone. Instead one makes a low-dimensional (1-
respectively 2-parameter) semi-empirical nonlinear ansatz for the gradient
dependence and fits it to ab initio computations for noble gas atoms (Becke)
respectively the lowest-order term in the gradient expansion of the weakly
inhomogeneous electron gas222in which the electron gas is subjected to a weak
and slowly-varying external potential and its energy density is expanded in
increasingly higher-order density derivatives and the Lieb-Oxford inequality
(PBE). We emphasize that the gradient expansion which partially informs PBE
requires small density gradients and provides no information on the regime
typical for real systems where the density gradient is of order $1$ (in atomic
units).

It seems to us that a very natural reference system which provides just such a
regime is the free electron gas in a box which underlies the celebrated Dirac
exchange, but with an important difference: one replaces periodic boundary
conditions by zero boundary conditions. While the former yield a uniform
density, the latter yield density gradients of order $1$ near the boundary,
see Figure 1. This is the system studied in the present paper. Note that its
exchange energy is captured correctly by the exact Hohenberg-Kohn functional
since the system is a high-density limit of the interacting electron gas in a
box (see eq. (1.8) below). The free electron gas with zero boundary conditions
can be solved exactly for $N$ electrons, just as in the periodic case. One can
then study its asymptotic behaviour in the thermodynamic limit where $N$ and
the sidelength $L$ of the box tend to infinity with the number of electrons
per unit volume, $N/L^{3}=\bar{\rho}$, remaining constant.

Figure 1: Density (left) and density gradient squared (right) of the free
electron gas with 1110 electrons in a three-dimensional box with zero boundary
conditions. The picture shows a two-dimensional cross-section through the
center of the box, and the number of electrons per unit volume was normalized
to $1$. The density gradients of order $1$ near the boundary persist in the
thermodynamic limit.

By careful asymptotic analysis, we are able to determine not just the bulk
contribution to the exchange energy, which is just the familiar Dirac exchange
regardless of the imposed boundary conditions (as has been shown previously
[7]), but also the next-order (surface) contribution to the exchange energy,
see Theorem 1.1 below. The next-order term is to our knowledge new and
captures the inhomogeneous boundary layer depicted in Figure 1. It also
captures two additional effects: a boundary-condition-induced small shift of
Fermi momentum and bulk density, and a long-range electrostatic finite-size
correction which would also be present for periodic boundary conditions (i.e.
unform density). Our asymptotic methods also yield the next-order (surface)
contribution to the exchange energy for GGA exchange functionals with general
$f$. Requiring these contributions to match yields a novel exact constraint on
GGAs (see eq. (1.14) below).

Main result in more detail. The free electron gas in a box consists of $N$
electrons moving freely in a three-dimensional box $Q_{L}=[0,L]^{3}$ of
sidelength $L$ and volume $V=L^{3}$ in the thermodynamic limit $N\to\infty$,
$V\to\infty$, with the number of electrons per unit volume, $\bar{\rho}=N/V$,
remaining constant. Mathematically, ground states of the finite system are
defined as minimizers of kinetic energy

$T[\Psi]=\sum_{s_{1},...,s_{N}\in\mathbb{Z}_{2}}\int_{Q_{L}^{N}}\frac{1}{2}\sum_{i=1}^{N}|\nabla_{r_{i}}\Psi(r_{1},s_{1},...,r_{N},s_{N})|^{2}dr_{1}...dr_{N}$
(1.1)

over square-integrable $N$-electron wave functions $\Psi$ with finite kinetic
energy (i.e., functions in the Sobolev space
$H^{1}((Q_{L}\times\mathbb{Z}_{2})^{N};\mathbb{C})$) subject to the following
constraints: normalization, $||\Psi||_{L^{2}}=1$; antisymmetry,
$\Psi(...,r_{i},s_{i},...,r_{j},s_{j},...)=-\Psi(..,r_{j},s_{j},...,r_{i},s_{i},...)$
for $i\neq j$ (where $(r_{i},s_{i})\in Q_{L}\times\mathbb{Z}_{2}$ are space-
spin coordinates for the ith electron); and one of the boundary conditions

$\displaystyle\Psi(r_{1},s_{1},...,r_{N},s_{N})=\Psi(r_{1}^{\prime},s_{1},...,r_{N}^{\prime},s_{N})\;\mbox{
if }r^{\prime}\\!-r\in L\mathbb{Z}^{3N}\quad\quad$ (periodic case) (1.2)
$\displaystyle\Psi(r_{1},s_{1},...,r_{N},s_{N})=0\;\mbox{ if any
}r_{i}\in\partial[0,L]^{3}\quad\quad$ (Dirichlet case) (1.3)
$\displaystyle\nabla_{r_{i}}\Psi(r_{1},s_{1},...,r_{N},s_{N})\cdot\nu(r_{i})=0\;\mbox{
if any }r_{i}\in\partial[0,L]^{3}\quad\quad$ (Neumann case). (1.4)

Here $\nu(r_{i})$ denotes the outward unit normal to $\partial[0,L]^{3}$ at
$r_{i}$. (Of course, in the Neumann case no boundary conditions are imposed on
the admissible functions; instead, minimizers then automatically satisfy
Neumann conditions.)

For ground states of non-interacting systems, such as the one above, the
exchange energy is defined as the difference between the quantum-mechanical
electron-electron interaction energy and the mean field energy,

$E_{x}[\Psi]=V_{ee}[\Psi]-\frac{1}{2}\int_{Q_{L}^{2}}\frac{\rho(r)\rho(r^{\prime})}{|r-r^{\prime}|}dr\,dr^{\prime}$
(1.5)

with

$V_{ee}[\Psi]=\sum_{s_{1},...,s_{N}\in\mathbb{Z}^{2}}\int_{Q_{L}^{N}}\sum_{1\leq
i<j\leq
N}\frac{1}{|r_{i}-r_{j}|}|\Psi(r_{1},s_{1},...,r_{N},s_{N})|^{2}dr_{1}...dr_{N}$
(1.6)

(interaction energy) and

$\rho(r)=N\sum_{s_{1},...,s_{N}\in\mathbb{Z}^{2}}\int_{Q_{L}^{N-1}}|\Psi(r,s_{1},r_{2},s_{2},...,r_{N},s_{N})|^{2}dr_{2}...dr_{N}$
(1.7)

(single-particle density of the system). For explicit expressions of the
exchange energy in terms of the single-particle orbitals (Laplace
eigenfunctions in the box) see Section 2.

We note that the exchange energy (1.5) is of fundamental interest as a
correction term appearing in the asymptotic expansion of the interacting
ground state energy in the high-density limit. More precisely, denoting the
interacting ground state energy, i.e. that of $T+V_{ee}$, in $Q_{L}$ by
$E^{\rm T+V_{ee}}_{L}$ and considering the re-scaling

$\Psi^{\gamma}(r_{1},s_{1},...,r_{N},s_{N})=\gamma^{3N/2}\Psi(\gamma
r_{1},s_{1},...,\gamma r_{N},s_{N})\;\;\;\;(\gamma>0),$

an elementary change of variables shows that

$E^{T+V_{ee}}_{L/\gamma}=\gamma^{2}E^{T+\gamma^{-1}V_{ee}}_{L}.$

Hence in the high-density limit $\gamma\to\infty$, the interaction is a small
perturbation and thus by standard perturbation theory, provided the ground
state $\Psi$ of $T$ in $Q_{L}$ is non-degenerate,

$\tfrac{1}{\gamma^{2}}E^{T+V_{ee}}_{L/\gamma}=T[\Psi]+\tfrac{1}{\gamma}\Bigl{(}J[\rho]+E_{x}[\Psi]\Bigr{)}+O\bigl{(}\tfrac{1}{\gamma^{2}}\bigr{)}\;\;\;\mbox{
as }\gamma\to\infty$ (1.8)

where $J[\rho]$ denotes the mean-field energy (second term in (1.5)). See the
review [18] for more information on the interacting electron gas and the
recent review [8] for information about the opposite (low-density) limit in
which electrons become strictly correlated.

We are not just interested in the exact exchange energy functional (1.5), but
also want to compare it to two important types of simpler functionals defined
only in terms of the single-particle density:

  * •

The Local Density Approximation (LDA) [17]:

$\displaystyle E_{x}^{\rm LDA}[\rho]=\int_{\mathbb{R}^{3}}e_{x}^{\rm
LDA}(\rho(r))\mathrm{d}r$ (1.9)

where the exchange energy density per unit volume is given by the Dirac-Bloch
formula [2, 6] $e^{\rm LDA}_{x}(\rho)=-c_{x}\rho^{4/3}$ with
$c_{x}=\frac{3}{4}(\frac{3}{\pi})^{\frac{1}{3}}$.

  * •

The GGA functionals [1, 19, 20]:

$\displaystyle E_{x}^{\rm GGA}[\rho]=E^{\rm
LDA}_{x}[\rho]+\underbrace{\int_{\mathbb{R}^{3}}g^{\rm
GGA}(\rho(r),|\nabla\rho(r)|)\mathrm{d}r}_{\coloneqq\Delta E_{x}^{\rm
GGA}[\rho]}$ (1.10)

with the assumption that $g^{\rm GGA}(\bar{\rho},0)=0$ for all $\bar{\rho}\geq
0$ (i.e. the functional reduces to the LDA for the homogeneous density).

In the physics literature [1, 19, 20], GGAs are commonly expressed in terms of
the density and the dimensionless gradient $s=|\nabla\rho|/\rho^{4/3}$. This
has the advantage that, by a scaling argument, one arrives at the simpler
ansatz

$g^{\rm GGA}(\rho,|\nabla\rho|)=e_{x}^{\rm LDA}(\rho)G(s),$ (1.11)

with different GGAs differing only by the choice of $G$.333I.e., by the
“gradient enhancement factor” $F=1+G$ of the overall integrand
$f(\rho,|\nabla\rho|)=e_{x}(\rho)+g^{\rm
GGA}(\rho,|\nabla\rho|)=e_{x}(\rho)F(s)$ The reason we prefer to work with the
density and the density gradient instead is because $s(r)\rightarrow\infty$ as
$r$ approaches the boundary (for Dirichlet boundary conditions) while
$\nabla\rho(r)$ remains bounded, making the mathematical analysis simpler. Our
assumptions on $g^{\rm GGA}$ required in Theorem 1.1 below are as follows:

$g^{\rm GGA}\in C^{0}([0,\infty)^{2})\cap
C^{1}((0,\infty)\times[0,\infty))\mbox{ with }g^{\rm
GGA}(\bar{\rho},0)=0\mbox{ for all }\bar{\rho}.$ (1.12)

These are satisfied for typical GGAs of form (1.11) such as those in [1, 19,
20] (see Appendix B for a proof).

With the functionals (1.5), (1.9), (1.10) in mind, the main result of this
paper can be stated as follows.

######  Theorem 1.1 (Asymptotic expansion of exchange functionals).

Let $N\in\mathbb{N},L>0$, and let $\Omega\subset\mathbb{R}^{3}$ be a
rectangular box. Let $\Psi_{N,L}$ be any determinantal ground state wave
function of the free $N$-electron gas in
$\Omega_{L}=\\{x\in\mathbb{R}^{3}:x/L\in\Omega\\}$ under either Dirichlet,
Neumann, or periodic boundary conditions, and let $\rho_{N,L}$ denote the
associated single-particle density. Moreover, assume that the GGA functional
(1.10) satisfies (1.12). Then in the thermodynamic limit, i.e., for
$N,L\rightarrow\infty$ and $\bar{\rho}=N/(|\Omega|L^{3})=constant$, one has:

  * •

Under periodic boundary conditions:

$\displaystyle E_{x}[\Psi^{\textnormal{Per}}_{N,L}]$
$\displaystyle=-c_{x}\bar{\rho}^{4/3}|\Omega|L^{3}+c_{x,2}^{\textnormal{Per}}\bar{\rho}|\partial\Omega|L^{2}+\mathcal{O}(L^{\frac{45}{23}+\epsilon})$
$\displaystyle E^{\rm LDA}_{x}[\rho^{\textnormal{Per}}_{N,L}]$
$\displaystyle=-c_{x}\bar{\rho}^{4/3}|\Omega|L^{3}+\mathcal{O}(L^{\frac{34}{23}+\epsilon})$
$\displaystyle\Delta E_{x}^{\rm GGA}[\rho^{\textnormal{Per}}_{N,L}]$
$\displaystyle=\mathcal{O}(L^{\frac{34}{23}+\epsilon})$

  * •

Under Dirichlet boundary condtions:

$\displaystyle E_{x}[\Psi^{\textnormal{Dir}}_{N,L}]$
$\displaystyle=-c_{x}\bar{\rho}^{\frac{4}{3}}|\Omega|L^{3}-c_{x,2}^{\textnormal{Dir}}\bar{\rho}|\partial\Omega|L^{2}+\mathcal{O}(L^{\frac{45}{23}+\epsilon})$
$\displaystyle E_{x}^{\rm LDA}[\rho^{\textnormal{Dir}}_{N,L}]$
$\displaystyle=-c_{x}\bar{\rho}^{\frac{4}{3}}|\Omega|L^{3}-c_{\rm
LDA}^{\textnormal{Dir}}\bar{\rho}|\partial\Omega|L^{2}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2})$
$\displaystyle\Delta E_{x}^{\rm GGA}[\rho^{\textnormal{Dir}}_{N,L}]$
$\displaystyle=c_{\rm
GGA}^{\textnormal{Dir}}(\bar{\rho})|\partial\Omega|L^{2}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2})$

  * •

Under Neumann boundary conditions:

$\displaystyle E_{x}[\Psi^{\textnormal{Neu}}_{N,L}]$
$\displaystyle=-c_{x}\bar{\rho}^{\frac{4}{3}}|\Omega|L^{3}-c_{x,2}^{\textnormal{Neu}}\bar{\rho}|\partial\Omega|L^{2}+\mathcal{O}(L^{\frac{45}{23}+\epsilon})$
$\displaystyle E_{x}^{\rm LDA}[\rho^{\textnormal{Neu}}_{N,L}]$
$\displaystyle=-c_{x}\bar{\rho}^{\frac{4}{3}}|\Omega|L^{3}-c_{\rm
LDA}^{\textnormal{Neu}}\bar{\rho}|\partial\Omega|L^{2}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2})$
$\displaystyle\Delta E_{x}^{\rm GGA}[\rho^{\textnormal{Neu}}_{N,L}]$
$\displaystyle=c_{\rm
GGA}^{\textnormal{Neu}}(\bar{\rho})|\partial\Omega|L^{2}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2})$

where $|\Omega|$ and $|\partial\Omega|$ denotes the volume and surface area of
the domain $\Omega$, $h$ is the explicit function $h(t)=3(\sin t-t\cos
t)/t^{3}$, $p_{F}=(3\pi^{2}\bar{\rho})^{1/3}$ (Fermi momentum), and the
constants are given by

$\displaystyle c_{x,2}^{\textnormal{Per}}=\frac{1}{8},\quad
c_{x,2}^{\textnormal{Dir}}=\frac{1-\log 2}{4}\approx 0.0767,\quad c_{\rm
LDA}^{\textnormal{Dir}}=\frac{3}{8\pi}\int_{0}^{\infty}(1-h(t))^{\frac{4}{3}}-1\mathrm{d}t+\frac{3}{8}\approx
0.0673,$ $\displaystyle c_{x,2}^{\textnormal{Neu}}=\frac{3\log 2-2}{4}\approx
0.0199,\quad c_{\rm
LDA}^{\textnormal{Neu}}=\frac{3}{8\pi}\int_{0}^{\infty}(1+h(t))^{\frac{4}{3}}-1\mathrm{d}t-\frac{3}{8}\approx
0.0430,$ $\displaystyle c_{\rm
GGA}^{\textnormal{Dir}}(\bar{\rho})=\frac{1}{2p_{F}}\int_{0}^{\infty}g^{\rm
GGA}\biggr{(}\bar{\rho}(1-h(t)),2\bar{\rho}p_{F}|\dot{h}(t)|)\biggl{)}\mathrm{d}t,$
$\displaystyle c_{\rm
GGA}^{\textnormal{Neu}}(\bar{\rho})=\frac{1}{2p_{F}}\int_{0}^{\infty}g^{\rm
GGA}\biggr{(}\bar{\rho}(1+h(t)),2\bar{\rho}p_{F}|\dot{h}(t)|\biggl{)}\mathrm{d}t.$

This result extends that of a previous work by one of the authors [7] as we
determine not just the leading but also the next-order terms (of order
$L^{2}$) and include the GGA functionals. Also to further illustrate the role
of the boundary conditions we have included the Neumann case.

An immediate corollary of Theorem 1.1 is the following simple exact constraint
on GGAs. The next-order correction to Dirac exchange for the free electron gas
with zero boundary conditions is captured exactly, i.e.

$\displaystyle\frac{E_{x}[\Psi_{N,L}]-E_{x}^{\rm
GGA}[\rho_{N,L}]}{L^{2}}\rightarrow 0\quad\mbox{ as
}N,L\rightarrow\infty\mbox{ with }\bar{\rho}=\frac{N}{L^{3}}=constant$ (1.13)

for all values of the average density $\bar{\rho}$, if and only if the
gradient enhancement factor $F(s)=1+G(s)$ defined by (1.11) satisfies

$\displaystyle\frac{3}{8\pi}\int_{0}^{\infty}(1-h(t))^{\frac{4}{3}}G\biggr{(}2(3\pi^{2})^{\frac{1}{3}}\frac{|\dot{h}(t)|}{(1-h(t))^{\frac{4}{3}}}\biggr{)}\mathrm{d}t=c_{x,2}^{\textnormal{Dir}}-c^{\textnormal{Dir}}_{\rm
LDA}$ (1.14)

where the constants and the function $h$ are those from Theorem 1.1. In
contrast with previous exact conditions on $G$ which refer to small-$s$
asymptotics [19] (for the weakly inhomogeneous electron gas) respectively
large-$s$ asymptotics [1] (for atomic densities), the above condition is an
integral constraint which sees the whole profile of $G$. Note that as $t$
varies from $0$ to $\infty$, the argument $s$ of $G$ (which corresponds to the
reduced density gradient of the Dirichlet free electron gas along a ray moving
from the boundary in perpendicular direction into the interior, see below)
traces out all possible $s$ values from $\infty$ to $0$. The extent to which
current GGAs fail to satisfy (1.14) is discussed in Section 7.

Finally we remark that, due to (1.8), our results correspond to first taking a
high-density limit and then a thermodynamic limit of the interacting
$N$-electron gas in a box. The opposite order of limits was studied in [9]
(note that when first taking a thermodynamic limit the surface correction
terms disappear). In the periodic case a coupled high-density/thermodynamic
limit in which very different corrections (of RPA type) appear is studied in
[5].

Strategy of the proof. We follow the overall strategy introduced in [7] of
deriving an accurate continuum approximation to the ground state density
matrix (see Theorem 4.2 below, or Theorems 4.1 and 4.2 in [7]) and analyzing
the ensuing interior and boundary contributions to the exchange energy. While
the continuum approximation is the same already introduced in [7], the main
advance, and most involved part of our work, is an improved error estimate
(see Theorem 4.2) which shows that it is accurate enough to infer the next-
order contributions to the exchange functionals which are of the order of the
surface area of the box. This is achieved by leveraging, on top of Fourier
analysis techniques [22] as already used in [7], the theory of exponential
sums [10, 14]. The main step is the proof of the following technical lemma
which relies on recent work by Heath-Brown [14].

######  Lemma 1.1.

Let $\alpha\in\mathbb{N}_{0}^{3}$ and $D\in\mathbb{R}^{3\times 3}$ be a
positive diagonal matrix. Then there exists $c=c(\alpha,D)>0$ such that

$\displaystyle\biggr{|}\sum_{k\in\mathbb{Z}^{3}\cap
B^{D}_{R}}(ik)^{\alpha}e^{ik\bm{\cdot}z}-\int_{B^{D}_{R}}(ik)^{\alpha}e^{ik\bm{\cdot}z}\mathrm{d}k\biggl{|}\leq
c(1+R^{|\alpha|+\frac{34}{23}+\epsilon}),$ (1.15)

for all $R>0$ and $z$ with $|z|_{max}\coloneqq\max_{j\leq
3}\\{|z_{j}|\\}\leq\pi$, where $\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$ and
$B^{D}_{R}\coloneqq\\{k\in\mathbb{R}^{3}:|D^{-1}k|\leq R\\}$.

The exponent $\frac{34}{23}$ may seem peculiar at first, and we do not claim
it to be optimal, but the main point is that it improves over the
$\frac{3}{2}$ exponent obtained in [7]. This improvement is necessary for
rigorous derivation of the asymptotic terms of the order of magnitude of the
surface area of the box. This can be quickly seen by integrating the square of
an error proportional to $L^{-\frac{3}{2}}$ (like in Theorem 4.1 and 4.2 of
[7]) against the Coulomb potential in the double box $[0,L]^{6}$, which yields
an error proportional to $L^{2}$ and is therefore not enough for our purpose.

Estimates of this kind were originally motivated by analytic number theory. In
particular, by setting $z=0$ and $\alpha=0$ one recognizes the famous lattice
point counting problem in $\mathbb{R}^{3}$ (also known as the sphere problem)
[23, 4, 3, 13]. The remarkable difference here is that the estimate holds
uniformly in $z$, even though the integrand oscillates, for typical $k\in
B_{R}^{D}$, on the length scale of the lattice.

Figure 2: Boundary layer of the free electron gas with 10014 electrons in a
three-dimensional box with zero boundary conditions. The exact density is
given by (2.10), the leading order continuum approximation by its
thermodynamic limit (1.18), and the next-order continuum approximation by eqs.
(1.16)–(1.17). The number of electrons per unit volume was normalized to the
conduction electron density of copper.

The continuum approximation of the density matrix which we justify with the
help of the above lemma (see Theorem 4.2) entails, in particular, the
following accurate approximation to the boundary layer for zero boundary
conditions and the box $[0,L]^{3}$:

$\rho_{N,L}(r_{0}+r^{\prime})=\bar{\rho}\Bigl{(}\frac{p_{N,L}}{p_{F}}\Bigr{)}^{3}\bigl{(}1-h(2p_{N,L}|r^{\prime}|)\bigr{)}+O(L^{-\frac{35}{23}+\varepsilon})$
(1.16)

whenever $r_{0}$ belongs to the boundary of the box, its distance from the
edges is of order $L$, and $r^{\prime}$ points in normal direction to
$\partial[0,L]^{3}$ into the interior. Here $p_{N,L}$ is the Fermi momentum of
the finite system, which is found (see Lemma 3.2) to differ from its
thermodynamic limit $p_{F}=(3\pi^{2}\bar{\rho})^{1/3}$ by an order $L^{-1}$
shift,

$p_{N,L}=p_{F}+\frac{3\pi}{4}L^{-1}+O(L^{\frac{35}{23}+\epsilon}).$ (1.17)

An important and rather subtle consequence of (1.16)–(1.17) is the following.
The – at first sight reasonable – leading order approximation (thermodynamic
limit) of the boundary layer profile in normal direction,

$\lim_{\begin{subarray}{c}N,\,L\to\infty\\\
N/L^{3}=\bar{\rho}\end{subarray}}\rho_{N,L}(r_{0}+r^{\prime})=\bar{\rho}\bigr{(}1-h(2p_{F}|r^{\prime}|)\bigr{)},$
(1.18)

while enough to resolve the surface contribution for GGA exchange444note that
the integrand in (1.14) is precisely the GGA correction to the exchange energy
density for this profile after making the change of variables
$t=2p_{F}|r^{\prime}|$ and setting $\bar{\rho}$=1, yields a surface
contribution of the wrong size for exact exchange, missing e.g. the $\sim
L^{-1}$ rearrangement of the bulk density implied by (1.16)–(1.17). The much
better accuracy of the next-order approximation (1.16)–(1.17) can be seen in
Figure 2.

Structure of the paper. We start with a small subsection to introduce the
notation used throughout the paper. In Section 2 we begin by recalling some
basic facts about the ground state of the free electron gas in the box under
different boundary conditions. In Section 3 we discuss the control of open
shell effects and the Fermi momentum asymptotics in the thermodynamic limit.
Section 4 contains the proof of Lemma 1.1 and the derivation of the continuum
approximation of the density matrix. In Section 5 we present the proof of
Theorem 1.1 by splitting it into treating semi-local functionals (Theorem 5.1)
and exact exchange (Theorem 5.2). Section 6 briefly discusses the asymptotics
of the kinetic energy, which can easily be extracted with our methods. Section
7 compares the asymptotic behaviour of different exchange functionals (exact
exchange, LDA, B88, PBE, PBEsol) when applied to the free electron gas in a
box with zero boundary conditions and up to 30 000 electrons. We find good
agreement between asymptotics and numerics. Physics-minded readers may want to
skip Sections 2–6 and move forward directly to Section 7.

###  1.1 Notation

The following notation will be used throughout the text.

  * •

We use the standard big-O and small-o notation: for functions
$f:(0,\infty)\rightarrow\mathbb{R}$ and $g:\mathbb{R}\rightarrow(0,\infty)$,
we say that $f=\mathcal{O}(g)$ respectively
$f=\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(g)$
if

$\displaystyle\limsup_{L\rightarrow\infty}\frac{|f(L)|}{g(L)}<\infty\quad\mbox{
respectively }\quad\limsup_{L\rightarrow\infty}\frac{|f(L)|}{g(L)}=0.$

Moreover, for functions $f,g:\mathbb{R}\mapsto\mathbb{R}$, we say that
$f(L)\lesssim g(L)$ or $f(L)\sim g(L)$ to indicate, respectively, the
existence of a constant $C>0$ which does not depend on $L$ such that

$\displaystyle|f(L)|\leq C|g(L)|\quad\mbox{ or }\quad
C^{-1}|f(L)|\leq|g(L)|\leq C|f(L)|$

for all sufficiently large values of $L$. Sometimes we will also use the
notation $\lesssim_{\epsilon}$ to indicate dependence of the implicit constant
on an additional parameter ($\epsilon$ in this case).

  * •

Throughout the text, $D\in\mathbb{R}_{+}^{3\times 3}$ will always denote a
diagonal matrix with entries $d_{1},d_{2},d_{3}>0$, and $|D|\coloneqq\det D$
stands for its determinant. The balls of radius $R$ and $D$-radius $R$ are
denoted by

$\displaystyle B_{R}\coloneqq\\{r\in\mathbb{R}^{3}:|r|\leq R\\}\quad\mbox{ and
}\quad B^{D}_{R}\coloneqq\\{r\in\mathbb{R}^{3}:|D^{-1}r|\leq R\\}.$

The cubic box, the $D$-rectangular box, and their re-scaled versions are
denoted by

$\displaystyle Q\coloneqq[0,1]^{3},\;\;\;$ $\displaystyle
Q^{D}\coloneqq[0,d_{1}]\times[0,d_{2}]\times[0,d_{3}],$ $\displaystyle
Q_{L}\coloneqq[0,L]^{3},\;\;\;$ $\displaystyle
Q_{L}^{D}\coloneqq[0,L\,d_{1}]\times[0,L\,d_{2}]\times[0,L\,d_{3}].$

  * •

For the Fourier transform of a function
$f:\mathbb{R}^{n}\rightarrow\mathbb{C}$, we use the normalization convention

$\displaystyle\widehat{f}(k)=\int_{\mathbb{R}^{n}}f(r)e^{-ik\bm{\cdot}r}\mathrm{d}r$
(1.19)

where $k\bm{\cdot}r\coloneqq\sum_{j=1}^{n}k_{j}r_{j}$ is the standard
Euclidean scalar product. We denote the inverse Fourier transform of $f$ by
$\widecheck{f}$.

  * •

For a set $\Omega\subset\mathbb{R}^{n}$, we use $\chi_{\Omega}$ for its
characteristic function. In particular, with the above convention for the
Fourier transform, in $\mathbb{R}^{3}$ one has that
$\widehat{\chi}_{B_{1}}=\frac{4\pi}{3}h(|k|)$, where the function
$h:\mathbb{R}\mapsto\mathbb{R}$ will appear many times in the sequel and is
given by

$\displaystyle h(t)=3\frac{\sin(t)-t\cos(t)}{t^{3}}.$ (1.20)

For an elementary derivation of this formula see e.g. [7] Lemma 6.1.

  * •

The group generated by reflections at coordinates hyperplanes of
$\mathbb{R}^{3}$ is denoted by $G$ and its elements by $\sigma$, i.e.,

$\displaystyle G=\\{\sigma\in\mathbb{R}^{3\times 3}:\sigma\mbox{ diagonal and
}\sigma_{jj}=\pm 1\mbox{ for any }j=1,2,3\\}.$ (1.21)

  * •

The projection on the $i^{th}$ coordinate hyperplane is denoted by
$\pi_{i}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{2}$, e.g.
$\pi_{1}(r_{1},r_{2},r_{3})=(r_{2},r_{3})$. Moreover, for any
$z\in\mathbb{R}^{3}$, we define $|z|_{\max}\coloneqq\max_{j\leq 3}|z_{j}|$ .

  * •

The positive reals, non-negative integers and the additive group of order 2
are denoted, respectively, by $\mathbb{R}_{+}=(0,\infty)$,
$\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$ and $\mathbb{Z}_{2}=\\{0,1\\}$.

  * •

For any set in $\mathbb{R}^{3}$, we use $|\cdot|$ for either its volume,
surface area or cardinality depending on whether the set has dimension 3, 2,
or 0, e.g, $|B_{1}|=\frac{4\pi}{3}$, $|\partial B_{1}|=4\pi$, and
$|B_{1}\cap\mathbb{Z}^{3}|$ is the number of elements in $\mathbb{Z}^{3}$ with
Euclidean norm smaller than $1$.

##  2 Ground state of the free electron gas : closed shell formulas

We now recall some basic facts and formulas for the ground state of the free
$N$-electron gas in the box subject to Dirichlet, Periodic or Neumann boundary
conditions. (The reader can find a detailed account of the eigenstates in [21]
and of the density matrices and exchange energy in [7].)

It is well known that the Laplacian in $Q^{D}_{L}$ under any of the discussed
boundary conditions (BCs) is diagonalizable in the sense that there exists an
orthonormal basis of $L^{2}(Q^{D}_{L})$ consisting of eigenvectors.
Furthermore, the eigenvectors and eigenvalues can be labelled by

  * •

vectors $k\in\mathbb{Z}^{3}$ in case of periodic boundary conditions:

$\displaystyle\phi^{L}_{k}(r)=\frac{1}{\sqrt{|D|L^{3}}}\,e^{i\frac{2\pi}{L}D^{-1}k\bm{\cdot}r},\quad\lambda_{k}=\frac{4\pi^{2}|D^{-1}k|^{2}}{L^{2}}.$
(2.1)

  * •

vectors $k\in\mathbb{N}^{3}$ in case of Dirichlet boundary condition:

$\displaystyle\phi^{L}_{k}(r)=\frac{1}{\sqrt{|D|L^{3}}}\prod_{i=1}^{3}\sqrt{2}\sin\biggr{(}k_{i}\frac{\pi}{d_{i}L}r_{i}\biggr{)},\quad\lambda_{k}=\frac{\pi^{2}|D^{-1}k|^{2}}{L^{2}}.$
(2.2)

  * •

vectors $k\in\mathbb{N}_{0}^{3}$ in case of Neumann boundary conditions:

$\displaystyle\phi^{L}_{k}(r)=\frac{1}{\sqrt{|D|L^{3}}}\prod_{\begin{subarray}{c}i=1\\\
k_{i}\neq
0\end{subarray}}^{3}\sqrt{2}\cos\biggr{(}k_{i}\frac{\pi}{d_{i}L}r_{i}\biggr{)},\quad\lambda_{k}=\frac{\pi^{2}|D^{-1}k|^{2}}{L^{2}}.$
(2.3)

As a consequence, one possible ground state for the $N$-electron gas (i.e. a
normalized anti-symmetric minimizer of (1.1) under one of the BCs (1.2)–(1.4))
with $N$ even is given by the determinantal wave function (or Slater
determinant)

$\displaystyle\Psi_{N,L}(x_{1},...,x_{N})=\frac{1}{\sqrt{N!}}\det\begin{pmatrix}\psi_{1}(x_{1})&\ldots&\psi_{1}(x_{N})\\\
\vdots&&\vdots\\\ \psi_{N}(x_{1})&\ldots&\psi_{N}(x_{N})\end{pmatrix},$ (2.4)

where $x_{\ell}=(r_{\ell},s_{\ell})\in Q^{D}_{L}\times\mathbb{Z}_{2}$ are the
space-spin variables and $\psi_{i}$ are the space-spin orbitals given by

$\displaystyle\psi_{2i-1}(x)=\phi^{L}_{k_{i}}(r)\chi_{1}(s),\quad\quad\psi_{2i}(x)=\phi^{L}_{k_{i}}(r)\chi_{0}(s)\quad\mbox{
for }i\in\biggr{\\{}1,...,\frac{N}{2}\biggr{\\}},$ (2.5)

where $\phi^{L}_{k_{i}}$ are the eigenfunctions defined in (2.1)–(2.3), and
$\\{k_{i}\\}_{i\leq N/2}$ is any subset of distinct vectors in
$\mathbb{N}^{3}$ (in the Dirichlet case) such that
$B_{\max_{i}|k_{i}|-\epsilon}\cap\mathbb{N}^{3}\subset\\{k_{i}\\}_{i\leq
N/2}$, for all $\epsilon>0$. In fact, the collection of all such Slater
determinants forms a basis for the ground state eigenspace of the free
$N$-electron gas in $Q^{D}_{L}$.

Let us also introduce the Fermi radius555Note that the Fermi radius also
depends on $D$. However, as $D$ will be fixed and $N,L$ will vary, we will not
exhibit this dependence in our notation.

$\displaystyle
R_{N}\coloneqq\begin{dcases}\min\\{R>0:\frac{N}{2}\leq|B^{D}_{R}\cap\mathbb{Z}^{3}|\\},&\mbox{
for periodic BCs,}\\\
\min\\{R>0:\frac{N}{2}\leq|B^{D}_{R}\cap\mathbb{N}^{3}|\\},&\mbox{ for
Dirichlet BCs,}\\\
\min\\{R>0:\frac{N}{2}\leq|B^{D}_{R}\cap\mathbb{N}_{0}^{3}|\\},&\mbox{ for
Neumann BCs.}\end{dcases}$ (2.6)

In the following, for boundary-condition-dependent quantities like $R_{N}$ we
will use superscripts Per, Dir, and Neu. The ground state of the free-electron
gas is unique for any $N$ satisfying the closed shell condition

$\displaystyle\frac{N}{2}=\begin{dcases}|B^{D}_{R^{\textnormal{Per}}_{N}}\cap\mathbb{Z}^{3}|,&\mbox{for
periodic BCs,}\\\
|B^{D}_{R^{\textnormal{Dir}}_{N}}\cap\mathbb{N}^{3}|,&\mbox{for Dirichlet
BCs,}\\\ |B^{D}_{R^{\textnormal{Neu}}_{N}}\cap\mathbb{N}_{0}^{3}|,&\mbox{for
Neumann BCs.}\end{dcases}$ (2.7)

In particular, by recalling that the spinless single-particle density matrix
of a state $\Psi$ is defined as

$\displaystyle\gamma_{\Psi}(r,\tilde{r})\coloneqq
N\\!\\!\\!\\!\\!\sum_{s_{1},...,s_{N}\in\mathbb{Z}^{2}}\int_{(\mathbb{R}^{3})^{N-1}}\\!\\!\Psi(r,s_{1},r_{2},s_{2},...r_{N},s_{N})\overline{\Psi(\tilde{r},s_{1},r_{2},s_{2},...,r_{N},s_{N})}\mathrm{d}r_{2}...\mathrm{d}r_{N},$
(2.8)

for any $N$ satisfying the closed shell condition (2.7), the spinless density
matrix of the (unique) aforementioned ground state is given by

$\displaystyle\gamma_{N,L}(r,\tilde{r})=\begin{dcases}\frac{2}{|D|L^{3}}\sum_{k\in\mathbb{Z}^{3}\cap
B^{D}_{R_{N}^{\textnormal{Per}}}}e^{i\frac{2\pi}{L}D^{-1}k\bm{\cdot}(r-\tilde{r})},&\mbox{
for periodic BCs,}\\\ \frac{1}{4|D|L^{3}}\sum_{\sigma\in
G}\det\sigma\sum_{k\in\mathbb{Z}^{3}\cap
B^{D}_{R_{N}^{\textnormal{Dir}}}}e^{i\frac{\pi}{L}D^{-1}k\bm{\cdot}(r-\sigma\tilde{r})},&\mbox{
for Dirichlet BCs,}\\\ \frac{1}{4|D|L^{3}}\sum_{\sigma\in
G}\sum_{k\in\mathbb{Z}^{3}\cap
B^{D}_{R_{N}^{\textnormal{Neu}}}}e^{i\frac{\pi}{L}D^{-1}k\bm{\cdot}(r-\sigma\tilde{r})},&\mbox{
for Neumann BCs,}\end{dcases}$ (2.9)

where $G$ is the reflection group defined in (1.21). Therefore, the single-
particle (or one-body) density can be written as

$\displaystyle\rho_{N,L}(r)=\begin{dcases}\gamma^{\textnormal{Per}}_{N,L}(r,r)=\frac{N}{|D|L^{3}},&\mbox{
for periodic BCs,}\\\
\gamma^{\textnormal{Dir}}_{N,L}(r,r)=\frac{1}{4|D|L^{3}}\sum_{\sigma\in
G}\det\sigma\sum_{k\in\mathbb{Z}^{3}\cap
B^{D}_{R_{N}^{\textnormal{Dir}}}}e^{i\frac{\pi}{L}D^{-1}k\bm{\cdot}(r-\sigma
r)},&\mbox{ for Dirichlet BCs,}\\\
\gamma^{\textnormal{Neu}}_{N,L}(r,r)=\frac{1}{4|D|L^{3}}\sum_{\sigma\in
G}\sum_{k\in\mathbb{Z}^{3}\cap
B^{D}_{R_{N}^{\textnormal{Neu}}}}e^{i\frac{\pi}{L}D^{-1}k\bm{\cdot}(r-\sigma
r)},&\mbox{ for Neumann BCs,}\end{dcases}$ (2.10)

and the exact exchange energy can be rewritten666The equivalence of (1.5) and
(2.11) with $\gamma_{\Psi}$ defined by (2.8) is in fact valid for any Slater
determinant $\Psi$ of doubly-occupied spatial orbitals. It follows from
(2.4)-(2.5) by straightforward calculation. as

$\displaystyle E_{x}[\Psi_{N,L}]=-\frac{1}{4}\int_{Q^{D}_{L}\times
Q^{D}_{L}}\frac{|\gamma_{N,L}(r,\tilde{r})|^{2}}{|r-\tilde{r}|}\mathrm{d}r\mathrm{d}\tilde{r},$
(2.11)

with $\gamma_{N,L}$ from equation (2.9).

######  Remark.

The above expression of density and density matrix as an average of an
exponential sum over the reflection group was introduced and derived in [7].
For the Neumann case the derivation is analogous.

The above expressions have a few simple but important symmetries that we state
as a lemma for further reference. (The proof is a straightforward
verification.)

######  Lemma 2.1 (Symmetries of $\rho_{N,L}$).

Let $\rho_{N,L}$ be defined by equation (2.10), then $|\nabla\rho_{N,L}|$ and
$\rho_{N,L}$ are unchanged under the following reflections:

$r_{i}\mapsto d_{i}L-r_{i}.$

##  3 Open shell effects and Fermi momentum asymptotics in the thermodynamic
limit

In this section we consider two important aspects of the FEG in the
thermodynamic limit: (i) we justify the use of formulas (2.9) for
determinantal ground states with a general number of particles
$N\in\mathbb{N}$ (Lemma 3.1), and (ii) we derive a two-term asymptotic formula
for the finite-size Fermi momentum (Lemma 3.2).

###  3.1 Open shell effects

Our goal now is to show that, in the thermodynamic limit, the single-particle
density matrix (and its derivatives) for any determinantal ground state is
pointwise close to the unique closed shell formulas in (2.9).

To this end, let us introduce the previous and the current shells for some
$N\in\mathbb{N}$ as

$\displaystyle N_{-}\coloneqq\max\\{n\in\mathbb{N}:n<N\mbox{ and
}n/2=|B^{D}_{R}\cap\mathbb{N}^{3}|\mbox{ for some
}R>0\\}=\max_{\epsilon>0}2|B^{D}_{R_{N}-\epsilon}\cap\mathbb{N}^{3}|,$
$\displaystyle N_{+}\coloneqq\min\\{n\in\mathbb{N}:n\geq N\mbox{ and
}n/2=|B^{D}_{R}\cap\mathbb{N}^{3}|\mbox{ for some
}R>0\\}=2|B^{D}_{R_{N}}\cap\mathbb{N}^{3}|,$

where again $\mathbb{N}$ is replaced by $\mathbb{N}_{0}$ and $\mathbb{Z}$ for
Neumann and periodic boundary conditions, respectively. Then according to [7,
Section 3], one can work with the unique (closed shell) ground state density
matrix $\gamma_{N_{-},L}$ up to a pointwise error proportional to $N^{-1/2}$
(or $L^{-\frac{3}{2}}$). However, as previously remarked, these estimates are
not enough to justify the use of the exact closed shell formulas on the
analysis of the next-to-leading order term in the asymptotic expansion for the
exact exchange. To improve on this error estimate and include the more general
rectangular box case, we use our own estimate in Lemma 1.1. More precisely,
setting $\alpha=0$ and $z=0$ in Lemma 1.1, we obtain777Note that for the case
of the cubic box, finding the optimal (algebraic) coefficient on the remainder
goes under the name of sphere problem and has been studied by many authors
[23, 4, 3, 13]. In particular, better estimates (with smaller exponent) than
(3.1) are available in this case.

$\displaystyle|\mathbb{Z}^{3}\cap
B^{D}_{R}|-\frac{4\pi}{3}|D|R^{3}=\mathcal{O}(R^{\frac{34}{23}+\epsilon}).$
(3.1)

As a consequence, by adapting the arguments in [7] we can prove the following
lemma.

######  Lemma 3.1 (Open shell control).

Let $\alpha,\beta\in\mathbb{N}_{0}^{3}$ and $\epsilon>0$. Then, there exists a
constant $c=c(\alpha,\beta,\epsilon)>0$ independent of $N$ and $L$ such that
for any determinantal ground state of the free $N$-electron gas in $Q^{D}_{L}$
(under either Dirichlet, Neumann or periodic boundary conditions) we have

$\displaystyle|\partial^{\alpha}_{r}\partial^{\beta}_{\tilde{r}}\gamma_{N,L}(r,\tilde{r})-\partial^{\alpha}_{r}\partial^{\beta}_{\tilde{r}}\gamma_{N_{-},L}(r,\tilde{r})|\leq
c\frac{N^{\frac{|\alpha|+|\beta|}{3}+\frac{34}{69}+\epsilon}}{L^{3+|\alpha|+|\beta|}}.$
(3.2)

In particular, if $\bar{\rho}=\frac{N}{|Q^{D}|L^{3}}$ is constant, one has

$\displaystyle\partial^{\alpha}_{r}\partial^{\beta}_{\tilde{r}}\gamma_{N,L}(r,\tilde{r})=\partial^{\alpha}_{r}\partial^{\beta}_{\tilde{r}}\gamma_{N_{-},L}(r,\tilde{r})+\mathcal{O}(L^{-\frac{35}{23}+\epsilon}).$

###### Proof.

For simplicity we disregard spin here. First, we note that by eq. (3.1), the
degeneracy of the open shell can be controlled by

$\displaystyle
d(N)=N_{+}-N_{-}\lesssim_{\epsilon}R_{N}^{\frac{34}{23}+\epsilon}\lesssim_{\epsilon}N^{\frac{34}{69}+\epsilon},$
(3.3)

where the $\epsilon$ in the second inequality is different from the first by a
factor of $1/3$. Next, if $\Psi$ is a determinantal ground state of the free
$N$-electron gas in $Q^{D}_{L}$, then we know (see section 2) that, up to a
phase factor,

$\displaystyle\Psi=\underbrace{\biggr{(}\bigwedge_{k\in
B_{R_{N_{-}}}}\phi^{L}_{k}\biggr{)}}_{=\Psi_{N_{-},L}}\wedge\widetilde{\phi}_{1}\wedge...\wedge\widetilde{\phi}_{N-N_{-}},$

where $\wedge$ is the anti-symmetric tensor product (see (2.4)) and
$\\{\tilde{\phi}_{i}\\}_{i\leq N-N_{-}}$ is a set of orthonormal functions
given by linear combinations of the orbitals in the open shell, i.e.

$\displaystyle\widetilde{\phi}_{i}=\sum_{\begin{subarray}{c}k\in\mathbb{N}^{3}\\\
|k|=R_{N}\end{subarray}}c_{ik}\phi_{k}^{L}\quad\mbox{ for some
$\\{c_{ik}\\}_{i,j}\subset\mathbb{C}$ satisfying
$\sum_{k}c_{ik}\overline{c_{jk}}=\delta_{ij}$.}$

In particular, $C^{\ast}=\\{\overline{c_{ki}}\\}$ is unitary (from
$\mathbb{C}^{N-N_{-}}$ to $\mathbb{C}^{d(N)}$) and we have

$\displaystyle\biggr{|}\sum_{i,j,k}c_{ij}\overline{c_{ik}}a_{j}b_{k}\biggr{|}\leq\biggr{(}\sum_{j}|a_{j}|^{2}\biggr{)}^{\frac{1}{2}}\biggr{(}\sum_{k}|b_{k}|\biggr{)}^{\frac{1}{2}}$
(3.4)

for any $(a_{1},..a_{d(N)}),(b_{1},...,b_{d(N)})\in\mathbb{C}^{d(N)}$. We can
now use the formula

$\displaystyle\gamma_{\Psi}(r,\tilde{r})=\underbrace{\sum_{k\in
B_{R_{N_{-}}}\cap\mathbb{N}^{3}}\phi_{k}^{L}(r)\overline{\phi_{k}^{L}}(\tilde{r})}_{=\gamma_{N_{-},L}(r,\tilde{r})}+\sum_{i=1}^{N-N_{-}}\widetilde{\phi}_{i}(r)\overline{\widetilde{\phi}_{i}}(\tilde{r})$

and the estimates (3.3),(3.4), and $|\partial^{\alpha}\phi_{k}^{L}(r)|\leq
c_{\alpha}(R/L)^{|\alpha|}L^{-\frac{3}{2}}\leq
c_{\alpha}N^{\frac{|\alpha|}{3}}L^{-\frac{3}{2}-|\alpha|}$ (see (2.1)–(2.3))
to conclude that

$\displaystyle|\partial^{\alpha}_{r}\partial^{\beta}_{\tilde{r}}\gamma_{\Psi}(r,\tilde{r})-\partial^{\alpha}_{r}\partial^{\beta}_{\tilde{r}}\gamma_{N_{-},L}(r,\tilde{r})|$
$\displaystyle=\biggr{|}\sum_{i=1}^{N-N_{-}}\partial^{\alpha}\widetilde{\phi_{i}}(r)\partial^{\beta}\overline{\widetilde{\phi}_{i}}(\tilde{r})\biggr{|}=\biggr{|}\sum_{i,j,k}c_{ij}\overline{c}_{ik}\partial^{\alpha}\phi_{j}^{L}(r)\overline{\partial^{\beta}\phi_{k}^{L}}(\tilde{r})\biggr{|}$
$\displaystyle\leq\biggr{(}\sum_{|j|=R_{N}}|\partial^{\alpha}\phi_{j}^{L}(r)|^{2}\bigg{)}^{\frac{1}{2}}\biggr{(}\sum_{|k|=R_{N}}|\overline{\partial^{\beta}\phi_{k}^{L}}(\tilde{r})|^{2}\biggr{)}^{\frac{1}{2}}$
$\displaystyle\lesssim\frac{N^{\frac{|\alpha|+|\beta|}{3}+\frac{34}{69}+\epsilon}}{L^{3+|\alpha|+|\beta|}}.$

∎

######  Remark.

One could equally well consider the density matrix of the current closed shell
$\gamma_{N_{+},L}$ in Lemma 3.1.

###  3.2 Fermi momentum asymptotics

Here we derive the asymptotics of the finite-size Fermi momentum appearing in
(1.17).

For the non-interacting electron gas model, the finite-size Fermi momentum is
defined as the momentum of the highest occupied orbital of the ground state
wave function. For the free $N$-electron gas in $Q^{D}_{L}$, it is simply
given (in atomic units) by

$\displaystyle p_{N,L}\coloneqq\begin{dcases}\frac{2\pi
R^{\textnormal{Per}}_{N}}{L},\quad&\mbox{ for periodic BCs,}\\\ \frac{\pi
R^{\textnormal{Dir}}_{N}}{L},\quad&\mbox{ for Dirichlet BCs,}\\\ \frac{\pi
R^{\textnormal{Neu}}_{N}}{L},\quad&\mbox{ for Neumann BCs,}\end{dcases}$ (3.5)

where $R_{N}$ is the Fermi radius. It is well known that in the thermodynamic
limit, the finite-size Fermi momentum converges to the (continuum) Fermi
momentum, defined as

$\displaystyle p_{F}=(3\pi^{2}\bar{\rho})^{\frac{1}{3}}.$

The next lemma presents the next-order correction of the finite-size Fermi
momentum, which is crucial for deriving the next-order corrections from
Theorem 1.1.

######  Lemma 3.2 (Fermi momentum asymptotics).

Let $\bar{\rho}=\frac{N}{|Q^{D}|L^{3}}$ be constant, $p_{F}$ be the Fermi
momentum and $p_{N,L}$ be the finite-size Fermi momentum. Then

$\displaystyle
p_{N,L}=\begin{dcases}p_{F}+\mathcal{O}(L^{-\frac{35}{23}+\epsilon}),&\mbox{
for periodic BCs,}\\\ p_{F}+\frac{\pi|\partial
Q^{D}|}{8|Q^{D}|}L^{-1}+\mathcal{O}(L^{-\frac{35}{23}+\epsilon}),&\mbox{ for
Dirichlet BCs,}\\\ p_{F}-\frac{\pi|\partial
Q^{D}|}{8|Q^{D}|}L^{-1}+\mathcal{O}(L^{-\frac{35}{23}+\epsilon}),&\mbox{for
Neumann BCs.}\end{dcases}$ (3.6)

###### Proof.

We present the proof only for the Dirichlet case and denote
$R^{\textnormal{Dir}}_{N}$ simply by $R_{N}$. First, from the previous section
we already know that $R_{N}\lesssim N^{\frac{1}{3}}$ and $0\leq N_{+}-N\leq
d(N)\lesssim N^{\frac{34}{69}+\epsilon}$. Since $|Q^{D}|=|D|$, one has

$\displaystyle\frac{\bar{\rho}L^{3}}{2}$
$\displaystyle=\frac{N}{2|D|}=\frac{|\mathbb{N}^{3}\cap
B^{D}_{R_{N}}|}{|D|}+\mathcal{O}(N^{\frac{34}{69}+\epsilon})$
$\displaystyle=\frac{1}{8|D|}(|\mathbb{Z}^{3}\cap
B^{D}_{R_{N}}|-\sum_{j=1}^{3}|\mathbb{Z}^{2}\cap\pi_{j}(B_{R_{N}}^{D})|)+\mathcal{O}(N^{\frac{34}{69}+\epsilon})$
$\displaystyle=\frac{\pi}{6}R_{N}^{3}-\frac{\pi|\partial
Q^{D}|}{16|D|}R_{N}^{2}+\mathcal{O}(N^{\frac{34}{69}+\epsilon}),$

which implies that

$\displaystyle p_{F}^{3}=p_{N,L}^{3}-\frac{3\pi|\partial
Q^{D}|}{8|D|}\frac{1}{L}p_{N,L}^{2}+\frac{\mathcal{O}(N^{\frac{34}{69}+\epsilon})}{L^{3}}.$
(3.7)

Then, with $\bar{\rho}=\frac{N}{|Q^{D}|L^{3}}$ fixed, it is clear that
$p_{N,L}^{3}\rightarrow p_{F}^{3}$. The correction proportional to $1/L$ now
follows from equation (3.7) by Taylor expansion of the (near the solution
$p_{F}$ for $1/L=0$ smooth) solution of the cubic equation for $p_{N,L}$. ∎

######  Remark.

The above lemma can be viewed as a variant of the famous two term Weyl’s law
[16] on the asymptotic behaviour of eigenvalues of the Laplacian (with an
improved remainder).

##  4 Discrete to continuum approximation

In this section, we present the proof of Lemma 1.1 and use it to derive the
continuum approximation for the density matrix with explicit estimates. To
simplify the notation and make the proofs more efficient, we will focus on the
case of the unit cube $Q$ and point out the modifications necessary for the
more general domain $Q^{D}$ at the end of each proof.

###  4.1 Estimate on exponential sums

Here we want to derive non-trivial estimates888the trivial estimate being
$|\sum_{k\in\mathbb{Z}^{n}\cap S}e^{if(k)}|\leq|\mathbb{Z}^{n}\cap S|$, which
holds for any subset $S\subset\mathbb{R}^{n}$ and any real-valued function
$f(k)$ for some weighted exponential sums that appear naturally in the proof
of Lemma 1.1. More precisely, let $\alpha\in\mathbb{N}_{0}^{3}$, $M,R>0$ and
$Q^{h}_{M}$ be the cubic holed box defined by

$\displaystyle
Q^{h}_{M}\coloneqq\\{r\in\mathbb{R}^{3}:M<|r|_{max}\leq\frac{4}{3}M\\},$ (4.1)

then our goal is to find a better than trivial estimate for the sum

$\displaystyle S^{\alpha}_{M}(R,z)\coloneqq\sum_{k\in
Q^{h}_{M}\cap(2\pi\mathbb{Z})^{3}}\frac{(k+z)^{\alpha}}{|k+z|^{|\alpha|+2}}e^{iR|k+z|},$
(4.2)

when $M$ is large and $R\sim M^{2}$. Using the above notation the main
estimate can be stated as follows.

######  Lemma 4.1.

Let $\epsilon>0$, $\alpha\in\mathbb{N}_{0}^{3}$ and $z\in\mathbb{R}^{3}$ with
$|z|_{\max}\leq\pi$. Then

$S_{M}^{\alpha}(R,z)\lesssim_{\epsilon}M^{\epsilon}(R^{\frac{1}{12}}M^{\frac{3}{4}}+M^{\frac{11}{12}}+R^{-\frac{1}{24}}M^{\frac{23}{24}}),$

where the implicit constant depends on $\epsilon$ and $\alpha$, but is
independent of $M,R$ and $z$.

Before proving this result, we recall a recent improvement due to Heath-Brown
[14] of the classical Van der Corput $k^{th}$-derivative estimate [10] on
exponential sums.

######  Theorem 4.1 (Heath-Brown [14]).

Let $k\geq 3$ be an integer and $f\in C^{k}([0,M],\mathbb{R})$. Suppose that

$0<\lambda_{k}\leq f^{(k)}(s)\leq A\lambda_{k},\quad s\in(0,M),$

for some $A>0$. Then

$\sum_{n\leq
M}e^{if(n)}\lesssim_{k,A,\epsilon}M^{1+\epsilon}\Bigl{(}\lambda_{k}^{\frac{1}{k(k-1)}}+M^{-\frac{1}{k(k-1)}}+M^{-\frac{2}{k(k-1)}}\lambda_{k}^{-\frac{2}{k^{2}(k-1)}}\Bigr{)}.$

######  Remark.

Note that for any $c<1$, the estimate above still holds for the sum in the
interval $cM\leq n\leq M$ (with two times the implicit constant above).
Indeed, one can simply consider the sum in this interval as the sum until $M$
minus the sum until $cM$, where both can be controlled by the same factor.

We also recall an elementary partial summation lemma (see e.g. [10, 15] for
the proof) that is used to deal with the weight function in
$S^{\alpha}_{M}(R,z)$.

######  Lemma 4.2 (Partial summation [15]).

Let $g\in C^{1}([a,b])$ and denote the total variation
$\int_{a}^{b}|g^{\prime}|$ of $g$ by $V_{g}[a,b]$. Then for any sequence
$a_{n}$, one has

$\biggr{|}\sum_{n\geq
a}^{b}g(n)a_{n}\biggl{|}\leq(V_{g}[a,b]+|g(a)|)\max_{\gamma\leq
b}\biggr{|}\sum_{\gamma}^{b}a_{n}\biggl{|}.$

###### Proof of Lemma 4.1.

First note that if $M$ is small (compared to $R$), then the estimate is
trivial. Therefore, we need to consider only the case where $M$ is big. Next,
defining

$\displaystyle f(k)\coloneqq R|k+z|\quad\mbox{ and }\quad
g(k)\coloneqq\frac{(k+z)^{\alpha}}{|k+z|^{2+|\alpha|}},$ (4.3)

the idea is to see $f$ and $g$ as functions of one coordinate, estimate the
inner sum by Theorem 4.1 and Lemma 4.2, and then use the trivial estimate for
the outer sums. To verify the assumptions of Theorem 4.1, let us consider the
sets

$\displaystyle U^{j}_{M}\coloneqq\\{k\in\mathbb{R}^{3}:|k_{j}|\geq M\mbox{ and
}\max_{\ell\neq j}\\{|k_{\ell}|\\}\geq M/10\\}\cap Q_{M}^{h},$ (4.4)
$\displaystyle L^{j}_{M}\coloneqq\\{k\in\mathbb{R}^{3}:|k_{j}|\geq M\mbox{ and
}\max_{\ell\neq j}\\{|k_{\ell}|\\}<M/10\\}\cap Q_{M}^{h},$ (4.5)

and observe that a straightforward calculation yields

$\displaystyle\partial_{k_{j}}^{4}f(k)=-3R\frac{|\pi_{j}(k+z)|^{2}(|\pi_{j}(k+z)|^{2}-4(k_{j}+z_{j})^{2})}{|k+z|^{7}}.$
(4.6)

Now the reason for our choice of $\tfrac{4}{3}$ in the definition of
$Q^{h}_{M}$ (any number between $1$ and $\sqrt{2}$ would be enough) is that,
since $|z|_{max}\leq\pi$ is small compared to $M$, inside of $U_{M}^{j}$ we
have

$\displaystyle 4(k_{j}+z_{j})^{2}-|\pi_{j}(k+z)|^{2}\sim
M^{2},\quad\quad|k+z|\sim M,\quad\mbox{ and }\quad|\pi_{j}(k+z)|^{2}\sim
M^{2}$

(with implicit constants independent of $M$ and $R$). In particular, from
(4.6), we find that

$\displaystyle|\partial^{4}_{k_{j}}f(k)|\sim\frac{R}{M^{3}}$ (4.7)

inside $U^{j}_{M}$. Moreover, note that $|g(k)|\lesssim\frac{1}{|k|^{2}}$ and
$|\partial_{k_{j}}g(k)|\lesssim\frac{1}{|k|^{3}}$. Hence, as $|k|\geq M$ for
any $k\in Q_{M}^{h}$, if we consider $g$ as a function of only one coordinate,
say $k_{3}$, then for any $\\{k_{1}\\}\times\\{k_{2}\\}\times I\subset
Q_{M}^{h}$ we have

$\displaystyle V_{g(k_{1},k_{2},\cdot)}[I]+|g(k)|\lesssim M^{-2}$ (4.8)

with implicit constant independent of $k_{1},k_{2},I,M$ and $R$. Therefore, we
can apply Lemma 4.2 and Theorem 4.1 (see (4.7),(4.8)) to conclude that, for
any $I_{1}\times I_{2}\times I_{3}\subset U_{M}^{j}$,

$\displaystyle\sum_{\begin{subarray}{c}k\in I_{1}\times I_{2}\times I_{3}\\\
k\in(2\pi\mathbb{Z})^{3}\end{subarray}}g(k)e^{if(k)}$
$\displaystyle\lesssim_{\epsilon}\sum_{\begin{subarray}{c}k_{\ell}\in
2\pi\mathbb{Z}\cap I_{\ell}\\\ \ell\neq
j\end{subarray}}M^{-2}M^{1+\epsilon}(R^{\frac{1}{12}}M^{-\frac{1}{4}}+M^{-\frac{1}{12}}+R^{-\frac{1}{24}}M^{-\frac{1}{24}})$
$\displaystyle\lesssim
M^{\epsilon}(R^{\frac{1}{12}}M^{\frac{3}{4}}+M^{\frac{11}{12}}+R^{-\frac{1}{24}}M^{\frac{23}{24}}).$
(4.9)

On the other hand, if $I_{1}\times I_{2}\times I_{3}\subset L_{M}^{j}$, then
from equation (4.6) and the definition of $L_{M}^{j}$ (see (4.5)), one can
show that $\partial_{k_{\ell}}^{4}f(k)\sim\frac{R}{M^{3}}$ for $\ell\neq j$.
In particular, by summing over $I_{\ell}$ first and using the same arguments
as before, we conclude that, as long as $I_{1}\times I_{2}\times I_{3}\subset
L_{M}^{j}$, the sum over $I_{1}\times I_{2}\times I_{3}$ is again bounded by
the right hand side of (4.9).

Finally, note that we can split the summation over $Q_{M}^{h}$ as

$\displaystyle\sum_{k\in
Q_{M}^{h}\cap(2\pi\mathbb{Z})^{3}}g(k)e^{if(k)}=\sum_{p=1}^{P}\biggr{(}\sum_{\begin{subarray}{c}k\in
I_{p}^{1}\times I_{p}^{2}\times I_{p}^{3}\\\
k\in(2\pi\mathbb{Z})^{3}\end{subarray}}g(k)e^{if(k)}\biggl{)},$

where the $I_{p}^{j}$ are intervals such that the product $I_{p}^{1}\times
I_{p}^{2}\times I_{p}^{3}$ is contained in one of the $U^{j}_{M}$ or
$L^{j}_{M}$, and the number $P$ is independent of $M$. Therefore, the result
follows from estimating the sums for each $p$ independently by (4.9) and
summing the estimates up (since $P$ is independent of $M$).∎

######  Remark (on the proof above).

  * •

The important point in the estimate above is that it gives better than trivial
estimates for $M\approx\sqrt{R}$, which is crucial for improving the exponent
3/2 in (1.15).

  * •

The classical Van der Corput $4^{th}$-derivative estimate (see [10]) would
actually be enough for our purposes. However, Theorem 4.1 gives a slightly
better estimate on the error term.

  * •

If we assume that $z\in\mathbb{Q}^{3}$, then one can show that the same
estimate from [3, Lemma 3.1] holds, which would then lead to the Vinogradov,
Chen [23, 4] exponent of $4/3$ in (1.15). Unfortunately, in this case the
implicit constant depends on the least common multiple of the denominators of
$z_{1},z_{2},z_{3}$, and therefore, we lose the uniform control that is
essential for our purposes in the next section.

######  Remark (Generalization to $Q^{D}$).

For general $Q^{D}$ we need an analog of Lemma 4.1 with $f$ and $g$ from
equation (4.3) replaced by $f(k)=R|D(k+z)|$ and
$g(k)=\frac{(D(k+z))^{\alpha}}{|D(k+z)|^{|\alpha|+2}}$. In particular, the
only significant modifications are that equation (4.6) has to change
accordingly, and the sets $U^{j}_{M}$ and $L^{j}_{M}$ defined in (4.4) and
(4.5) have to be replaced by their pre-images under $D$ (as a map in
$\mathbb{R}^{3}$).

###  4.2 Proof of Lemma 1.1

To complete the proof of Lemma 1.1, we shall use two classical results. The
first one is the celebrated Poisson summation formula (see e.g. [22]), already
used in the present context in [7], which we write in a slightly different
form to fit our goal.

######  Lemma 4.3 (Poisson summation formula).

Let $u\in C^{\infty}_{c}(\mathbb{R}^{n})$, then (see the Fourier transform
convention adopted (1.19)) the Poisson summation formula states that

$\displaystyle\sum_{k\in\mathbb{Z}^{n}}\frac{1}{(2\pi)^{n}}\widehat{u}(k)e^{ik\bm{\cdot}z}-u(z)=\sum_{\begin{subarray}{c}k\in(2\pi\mathbb{Z})^{n}\\\
k\neq 0\end{subarray}}u(z+k).$ (4.10)

The second one is an (optimal) estimate on the decay of the Fourier transform
of the characteristic function of the ball and its derivatives. The proof is a
straightforward calculation (at least in $\mathbb{R}^{3}$).

######  Lemma 4.4 (Fourier transform of the ball).

Let $\chi_{B_{1}}$ be the characteristic function of the unit ball in
$\mathbb{R}^{3}$. Then

$\widehat{\chi_{B_{1}}}(k)=\frac{4\pi}{3}h(|k|),$

where $h(s)=3(\sin s-s\cos s)/s^{3}$. Moreover, for any
$\alpha\in\mathbb{N}_{0}^{3}$, there exists $c_{\alpha}>0$ such that

$\displaystyle|\partial^{\alpha}\widehat{\chi_{B_{1}}}(k)|\leq\frac{c_{\alpha}}{(1+|k|)^{2}},\quad\forall
k\in\mathbb{R}^{3}.$ (4.11)

###### Proof of Lemma 1.1.

Let $R>0$, then our goal now is to bound the error term

$\displaystyle\mathcal{E}^{\alpha}(R,z)\coloneqq\sum_{k\in
B_{R}\cap\mathbb{Z}^{3}}(ik)^{\alpha}e^{ik\bm{\cdot}z}-\int_{B_{R}}(ik)^{\alpha}e^{ik\bm{\cdot}z}dk.$

For this, the main idea is to apply the Poisson summation formula to a smooth
version of the characteristic function of the ball of radius $R$, use Lemma
4.1 to control the error on the Fourier side and then estimate the difference
between our smoothed error and the error $\mathcal{E}^{\alpha}$ defined above.
So first, let $H>0$ be a small parameter to be chosen later, pick any non-
positive function $\phi\in C_{c}^{\infty}(0,1)$ with $\int\phi(s)ds=-1$, and
define

$\displaystyle f_{R,H}(r)=\begin{cases}1,&\mbox{if}\quad|r|\leq R,\\\
\displaystyle{1+\int_{0}^{|r|-R}\phi^{H}(s)\mathrm{d}s},&\mbox{otherwise,}\end{cases}$
(4.12)

where $\phi^{H}(s)\coloneqq\frac{1}{H}\phi\bigr{(}\frac{s}{H}\bigr{)}$. One
can check that $f_{R,H}$ is smooth, assumes only values between 0 and 1, and
vanishes outside $B_{R+H}$. Moreover, observing that
$\widecheck{f}_{R,H}(k)=\frac{1}{(2\pi)^{3}}\widehat{f}_{R,H}(k)$ (as
$f_{R,H}$ is symmetric) and applying the Poisson summation formula to
$\widehat{u}=f_{R,H}$ we obtain

$\displaystyle\underbrace{\sum_{k\in\mathbb{Z}^{3}}f_{R,H}(k)e^{ik\bm{\cdot}z}-\widehat{f}_{R,H}(z)}_{\coloneqq\mathcal{E}_{H}(R,z)}=\underbrace{\sum_{\begin{subarray}{c}k\in(2\pi\mathbb{Z})^{3}\\\
k\neq
0\end{subarray}}\widehat{f}_{R,H}(z+k)}_{\coloneqq\widehat{\mathcal{E}}_{H}(R,z)}.$
(4.13)

The next step is to use Lemma 4.1 to estimate the right-hand side of (4.13).
For this, we first use spherical coordinates to obtain

$\displaystyle\widehat{f}_{R,H}(z+k)$
$\displaystyle=\int_{\mathbb{R}^{3}}f_{R,H}(r)e^{-ir\bm{\cdot}(z+k)}\mathrm{d}r$
$\displaystyle=4\pi\int_{0}^{R}\frac{\tau\sin(|k+z|\tau)}{|k+z|}\mathrm{d}\tau+4\pi\int_{R}^{R+H}\frac{\tau\sin(|k+z|\tau)}{|k+z|}\biggr{(}1+\int_{0}^{\tau-R}\phi^{H}(s)ds\biggl{)}\mathrm{d}\tau$
$\displaystyle=-4\pi\int_{R}^{R+H}\frac{\sin(|k+z|\tau)-(|k+z|\tau)\cos(|k+z|\tau)}{|k+z|^{3}}\phi^{H}(\tau-R)\mathrm{d}\tau.$

In particular, by partial integration and re-scaling we have

$\displaystyle\widehat{\mathcal{E}}_{H}(R,z)$
$\displaystyle=-4\pi\int_{R}^{R+H}\biggr{(}3\dot{\phi}^{H}(\tau-R)+\tau\ddot{\phi}^{H}(\tau-R)\biggr{)}\sum_{\begin{subarray}{c}k\in(2\pi\mathbb{Z})^{3}\\\
k\neq 0\end{subarray}}\frac{\cos(|k+z|\tau)}{|k+z|^{4}}\mathrm{d}\tau$
$\displaystyle=-4\pi\int_{0}^{1}\biggr{(}\overbrace{\frac{3}{H}\dot{\phi}(\tau)+\frac{\tau}{H}\ddot{\phi}(\tau)}^{=\mathcal{O}(\frac{1}{H})}+\frac{R}{H^{2}}\ddot{\phi}(\tau)\biggr{)}\sum_{\begin{subarray}{c}k\in(2\pi\mathbb{Z})^{3}\\\
k\neq
0\end{subarray}}\frac{\cos\bigr{(}(R+H\tau)|k+z|\bigr{)}}{|k+z|^{4}}\mathrm{d}\tau$
(4.14) $\displaystyle=\mathcal{O}\biggr{(}\frac{1}{H}\biggr{)}-\frac{4\pi
R}{H^{2}}\int_{0}^{1}\ddot{\phi}(\tau)\sum_{\begin{subarray}{c}k\in(2\pi\mathbb{Z})^{3}\\\
k\neq
0\end{subarray}}\frac{\cos\bigr{(}(R+H\tau)|k+z|\bigr{)}}{|k+z|^{4}}\mathrm{d}\tau,$

where for the last equality we used that
$\sum_{k\in(2\pi\mathbb{Z})^{3}\setminus\\{0\\}}|k+z|^{-4}\leq c$ (since
$|z|_{\max}\leq\pi$). Here we slightly abused the big-O notation to denote
terms bounded in absolute value by a constant times $H^{-1}$ by
$\mathcal{O}(H^{-1})$. Now let $M>0$, then by integrating by parts $n$ times,
throwing out the terms outside the box $Q_{M}\setminus
Q_{1}=[-M,M]^{3}\setminus[-1,1]^{3}$, and integrating back by parts (recall
that $\phi\in C^{\infty}_{c}(0,1)$), we find that

$\displaystyle\widehat{\mathcal{E}}_{H}(R,z)$ $\displaystyle=4\pi
R\int_{0}^{1}\phi(\tau)\sum_{\begin{subarray}{c}k\in(2\pi\mathbb{Z})^{3}\\\
\cap Q_{M}\setminus
Q_{1}\end{subarray}}\frac{\cos\bigr{(}(R+H\tau)|(k+z)|\bigr{)}}{|k+z|^{2}}\mathrm{d}\tau+\mathcal{O}\biggr{(}\frac{1}{H}+\frac{R}{H^{n+1}M^{n}}\biggr{)}.$

Next, after possibly replacing the implicit constant above by a factor
independent of $M$, $H$, and $R$, we can assume that $M=(\frac{4}{3})^{q}$ for
some $q\in\mathbb{N}$. Thus by using the decomposition $Q_{M}\setminus
Q_{1}=\bigcup_{k=1}^{q-1}Q_{(\frac{4}{3})^{k}}^{h}$, the facts that $H<1$ and
$|z|_{max}\leq\pi$ (by periodicity), and Lemma 4.1, we find that

$\displaystyle\widehat{\mathcal{E}}_{H}(R,z)\lesssim_{\epsilon,n}$
$\displaystyle\frac{1}{H}+\frac{R}{H^{n+1}M^{n}}+\sum_{k=1}^{q-1}R\biggr{(}\frac{4}{3}\biggr{)}^{k\epsilon}\biggr{(}R^{\frac{1}{12}}\biggr{(}\frac{4}{3}\biggr{)}^{\frac{3k}{4}}+\biggr{(}\frac{4}{3}\biggr{)}^{\frac{11k}{12}}+R^{-\frac{1}{24}}\biggr{(}\frac{4}{3}\biggr{)}^{\frac{23k}{24}}\biggr{)}$
$\displaystyle\lesssim_{\epsilon,n}$
$\displaystyle\frac{1}{H}+\frac{R}{H^{n+1}M^{n}}+\log
M\biggr{(}R^{\frac{13}{12}}M^{\frac{3}{4}+\epsilon}+RM^{\frac{11}{12}+\epsilon}+R^{\frac{23}{24}}M^{\frac{23}{24}+\epsilon}\biggr{)},$
(4.15)

where the $\log M$ can be absorbed in $M^{\epsilon}$. Moreover, we can
estimate the error coming from the smoothing procedure as

$\displaystyle|\mathcal{E}^{0}(R,z)-\mathcal{E}_{H}(R,z)|\leq\\!\\!\sum_{\begin{subarray}{c}k\in\mathbb{Z}^{3}\cap
B_{R+H}\\\ k\not\in
B_{R}\end{subarray}}\\!\\!\\!\\!\\!\\!|f_{R,H}(k)|+|\widehat{f}_{R,H}(z)-\widehat{\chi_{B_{R}}}(z)|\lesssim
R^{2}H.$ (4.16)

Hence by summing estimates (4.15) and (4.16), taking into account (4.13),
setting $M=R^{m}$ and $H=R^{-h}$, and minimizing the exponents, one concludes
that

$\displaystyle\mathcal{E}^{0}(R,z)\lesssim_{\epsilon,n}R^{2-h(n)+\epsilon
m(n)}\log R,\quad\forall n\in\mathbb{N},$

where $h(n)\coloneqq\frac{12n+11}{23n+22}$ and
$m(n)\coloneqq\frac{12n+12}{23n+22}$. Thus, since $h\uparrow\frac{12}{23}$ for
$n\rightarrow\infty$, by choosing $n$ big enough, the result for $\alpha=0$
follows.

For the estimates with $\alpha\neq 0$, we repeat the same arguments by taking
$f^{\alpha}_{R,H}(r)=(ir)^{\alpha}f_{R,H}(r)$. In this case, note that by the
Leibniz rule and using that
$\partial^{\alpha}(|r|^{-4})\lesssim|r|^{-4-|\alpha|}$ away from the origin,
we have

$\displaystyle\partial_{k}^{\alpha}\biggr{(}\frac{\cos((R+H\tau)|k+z|)}{|k+z|^{4}}\biggl{)}-\frac{\bigr{(}iR(k+z)\bigr{)}^{\alpha}}{2|k+z|^{|\alpha|+4}}\biggr{(}e^{i(R+H\tau)|k+z|}+(-1)^{|\alpha|}e^{-i(R+H\tau)|k+z|}\biggr{)}$
$\displaystyle\lesssim\frac{R^{|\alpha|-1}}{|k+z|^{4}}.$ (4.17)

We also have, by Poisson summation and the symmetry of $f^{\alpha}_{R,H}$,
that

$\displaystyle\underbrace{\sum_{k\in\mathbb{Z}^{3}}f^{\alpha}_{R,H}(k)e^{ik\bm{\cdot}z}-\widehat{f}^{\alpha}_{R,H}(z)}_{\coloneqq\mathcal{E}^{\alpha}_{H}(R,z)}=\underbrace{\sum_{\begin{subarray}{c}k\in(2\pi\mathbb{Z})^{3}\\\
k\neq
0\end{subarray}}\widehat{f}^{\alpha}_{R,H}(z+k)}_{\coloneqq\widehat{\mathcal{E}}^{\alpha}_{H}(R,z)}.$
(4.18)

Hence using the identity
$\widecheck{(i\cdot)^{\alpha}f}=\partial^{\alpha}\widecheck{f}$ together with
estimate (4.17), and repeating the same steps from before (see (4.14)), we
conclude that

$\displaystyle\mathcal{E}^{\alpha}_{H}(R,z)$
$\displaystyle=-4\pi\int_{0}^{1}\biggr{(}\frac{3}{H}\dot{\phi}(\tau)+\frac{R+H\tau}{H^{2}}\ddot{\phi}(\tau)\biggr{)}\sum_{\begin{subarray}{c}k\in(2\pi\mathbb{Z})^{3}\\\
k\neq
0\end{subarray}}\partial^{\alpha}\biggr{(}\frac{\cos((R+H\tau)|k+z|)}{|k+z|^{4}}\biggl{)}$
$\displaystyle=-\frac{4\pi
R}{H^{2}}\int_{0}^{1}\ddot{\phi}(\tau)\\!\\!\sum_{\begin{subarray}{c}k\in(2\pi\mathbb{Z})^{3}\\\
k\neq
0\end{subarray}}\\!\\!\frac{\bigr{(}iR(k+z)\bigr{)}^{\alpha}}{2|k+z|^{|\alpha|+4}}\biggr{(}e^{i(R+H\tau)|k+z|}+(-1)^{|\alpha|}e^{-i(R+H\tau)|k+z|}\biggr{)}\mathrm{d}\tau$
$\displaystyle+\mathcal{O}\biggr{(}\frac{R^{|\alpha|}}{H^{2}}\biggr{)}.$

Thus we can again integrate by parts $n$ times, remove all terms for
$|k|\not\in Q_{M}\setminus Q_{1}$ and use Lemma 4.1 to show that

$\displaystyle\mathcal{E}^{\alpha}_{H}(R,z)\lesssim_{\epsilon}R^{|\alpha|}\biggr{(}\frac{1}{H^{2}}+\frac{R}{H^{n+1}M^{n}}+R^{\frac{13}{12}}M^{\frac{3}{4}+\epsilon}+RM^{\frac{11}{12}+\epsilon}+R^{\frac{23}{24}}M^{\frac{23}{24}+\epsilon}\biggr{)}.$
(4.19)

Finally, by combining (4.19) with (4.18) and the following simple estimate for
the smoothing error

$|\mathcal{E}^{\alpha}(R,z)-\mathcal{E}^{\alpha}_{H}(R,z)|\lesssim
R^{2+|\alpha|}H$

one obtains an overall bound which is just $R^{|\alpha|}$ times the right hand
sides of (4.15) plus (4.16), and the result follows by optimizing the exponent
and using the same $n\rightarrow\infty$ argument as before. ∎

######  Remark (Generalization to $Q^{D}$).

To deal with the general case, one has to replace $f_{R,H}$ by
$f^{D}_{R,H}(r)\coloneqq f_{R,H}(D^{-1}r)$ and observe that
$\widehat{f}^{D}_{R,H}(k)=|D|\widehat{f}_{R,H}(Dk)$. The rest of the proof
follows from the same steps by using the generalization of Lemma 4.1 discussed
in the previous remark.

###  4.3 Continuum approximation of density matrices

We can now introduce continuum approximations for the density matrices and
prove that they accurately approximate the exact ground state density
matrices.

Recalling formula (2.9) and replacing the sum over $k$ by an integral, we can
define the continuum version of the density matrix as

$\displaystyle\gamma_{N,L}^{\textnormal{ctm}}(r,\tilde{r})=\begin{dcases}\bar{\rho}\frac{p_{N,L}^{3}}{p_{F}^{3}}h(p_{N,L}|r-\tilde{r}|_{D,L}),&\mbox{
for periodic BCs,}\\\ \bar{\rho}\frac{p_{N,L}^{3}}{p_{F}^{3}}\sum_{\sigma\in
G}\det\sigma\,h(p_{N,L}|r-\sigma\tilde{r}|_{D,2L}),&\mbox{ for Dirichlet
BCs,}\\\ \bar{\rho}\frac{p_{N,L}^{3}}{p_{F}^{3}}\sum_{\sigma\in
G}h(p_{N,L}|r-\sigma\tilde{r}|_{D,2L}),&\mbox{ for Neumann BCs ,}\end{dcases}$
(4.20)

where $p_{N,L}$ and $p_{F}$ are the finite-size and continuum Fermi momentum
(see (3.5)), the function $h$ is defined in (1.20), and $|r|_{D,2L}=|r\mod
D(2L\mathbb{Z})^{3}|$ (the torus distance). In fact, in the Dirichlet case
this approximation was already introduced in [7] but the error estimates
obtained there were insufficient to conclude that the approximation is
accurate enough to reveal the surface correction to Dirac exchange.

The following new estimate will be sufficient to obtain the surface
correction. It is a direct consequence of Lemma 1.1.

######  Theorem 4.2 (Continuum approximation estimates).

Let $\bar{\rho}=\frac{N}{|Q^{D}|L^{3}}$ be fixed. Let $\Psi_{N,L}$ be any
determinantal ground state of the free $N$-electron gas on $Q^{D}_{L}$, and
let $\gamma_{N,L}$ be its single-particle density matrix introduced in (2.8).
Then for any $\epsilon>0$ we have

$\displaystyle|\gamma_{N,L}(r,\tilde{r})-\gamma_{N,L}^{\textnormal{ctm}}(r,\tilde{r})|\lesssim_{\epsilon}\bar{\rho}N^{-\frac{35}{69}+\epsilon}\lesssim_{\epsilon}\bar{\rho}^{\frac{34}{69}+\epsilon}L^{-\frac{35}{23}+3\epsilon}.$
(4.21)

More generally, for any $\alpha,\beta\in\mathbb{N}_{0}^{3}$,

$\displaystyle|\partial_{r}^{\alpha}\partial_{\tilde{r}}^{\beta}\gamma_{N,L}(r,\tilde{r})-\partial^{\alpha}_{r}\partial_{\tilde{r}}^{\beta}\gamma_{N,L}^{\textnormal{ctm}}(r,\tilde{r})|\lesssim_{\epsilon,\alpha,\beta}\bar{\rho}^{1+\frac{|\alpha|+|\beta|}{3}}N^{-\frac{35}{69}+\epsilon}\lesssim_{\epsilon,\alpha,\beta}\bar{\rho}^{\frac{34}{69}+\frac{|\alpha|+|\beta|}{3}+\epsilon}L^{-\frac{35}{23}+3\epsilon}.$
(4.22)

Moreover, one has

$\displaystyle|\gamma^{\textnormal{Dir/Neu}}_{N,L}(r,\tilde{r})|\lesssim_{\epsilon}\bar{\rho}\biggr{(}N^{-\frac{35}{69}+\epsilon}+(1+|r-\tilde{r}|_{D,2L})^{-2}\biggl{)},$
(4.23)
$\displaystyle|\gamma^{\textnormal{Per}}_{N,L}(r,\tilde{r})|\lesssim_{\epsilon}\bar{\rho}\biggr{(}N^{-\frac{35}{69}+\epsilon}+(1+|r-\tilde{r}|_{D,L})^{-2}\biggl{)}.$
(4.24)

###### Proof.

According to Lemma 3.1, up to an error
$\lesssim\bar{\rho}^{\frac{34}{69}+\epsilon}L^{-\frac{35}{69}+3\epsilon}$, we
can use the closed shell formulas for any determinantal ground state.

Next, we know that (i) $R\sim L$, (ii)
$\widehat{\chi}_{B^{D}_{R}}(z)=|D|R^{3}\widehat{\chi}_{B_{1}}(RDz)$, (iii)
$\widehat{\chi}_{B_{1}}(z)=\frac{4\pi}{3}h(|z|)$, and (iv) the sums
$\sum_{k\in\mathbb{Z}^{3}\cap B_{R}}e^{i\frac{\pi}{L}k\bm{\cdot}D^{-1}w}$ and
$\sum_{k\in\mathbb{Z}^{3}\cap B_{R}}e^{i\frac{2\pi}{L}k\bm{\cdot}D^{-1}w}$ are
periodic (in $w$) with respect to $D((2L\mathbb{Z})^{3})$ and
$D((L\mathbb{Z})^{3})$ respectively. Hence, estimate (4.21) follows directly
from the closed shell formulas (2.9) by applying Lemma 1.1 with
$z=\frac{\pi}{L}D^{-1}(r-\sigma\tilde{r})\mod(2\pi\mathbb{Z})^{3}$ for the
Dirichlet and Neumann case and
$z=\frac{2\pi}{L}D^{-1}(r-\tilde{r})\mod(2\pi\mathbb{Z})^{3}$ for the periodic
case.

For the derivative estimates in (4.22), one can simply use Lemma 1.1 together
with the identity
$\widehat{(-i\cdot)^{\alpha}f}=\partial^{\alpha}\widehat{f}$. Note that each
derivative gives an additional factor of $1/L$ which compensates for the
$R^{|\alpha|}$ gained in (1.15) and accounts for the factor
$\bar{\rho}^{\frac{|\alpha|}{3}}$ in (4.22). The decay estimates (4.23) and
(4.24) follow from estimate (4.21) and the decay of $h$ (see Lemma 4.4). ∎

##  5 Proof of Theorem 1.1

In this section we make the following simplification to the current notation:
as we are dealing with the thermodynamic limit ($N,L\rightarrow\infty$ with
$\bar{\rho}=\frac{N}{|Q^{D}|L^{3}}$ fixed), all functions and constants
depending on both $N$ and $L$ will simply be indexed by $L$ (e.g.,
$\rho_{L}=\rho_{N,L}$).

###  5.1 Semi-local functionals

Our goal now is to prove the following two-term asymptotic expansion for
general semi-local functionals of the density and its gradient.

######  Theorem 5.1 (Asymptotics for general semi-local functionals).

Let $\frac{N}{|Q^{D}|L^{3}}=\bar{\rho}=constant$, $\Psi_{L}$ be any
determinantal ground state of the free $N$-electron gas in $Q^{D}_{L}$, and
$\rho_{L}$ be the associated single-particle density. Suppose that $f(a,b)\in
C^{0}([0,\infty)^{2})\cap C^{1}((0,\infty)\times[0,\infty))$ and let
$F[\rho_{L}]\coloneqq\int_{Q^{D}_{L}}f(\rho_{L}(r),|\nabla\rho_{L}(r)|)\mathrm{d}r$.
Then for any $\epsilon>0$ we have

$\displaystyle
F[\rho_{L}]=\begin{dcases}f(\bar{\rho},0)|Q^{D}|L^{3}+\mathcal{O}(L^{-\frac{34}{23}+\epsilon}),&\mbox{
for periodic BCs,}\\\
f(\bar{\rho},0)|Q^{D}|L^{3}+\biggr{(}c^{\textnormal{Dir}}_{BL}(\bar{\rho})+c_{FM}(\bar{\rho})\biggr{)}|\partial
Q^{D}|L^{2}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2}),&\mbox{
for Dirichlet BCs,}\\\
f(\bar{\rho},0)|Q^{D}|L^{3}+\biggr{(}c^{\textnormal{Neu}}_{BL}(\bar{\rho})-c_{FM}(\bar{\rho})\biggr{)}|\partial
Q^{D}|L^{2}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2}),&\mbox{
for Neumann BCs,}\end{dcases}$

where the boundary layer and Fermi momentum corrections are given by

$\displaystyle
c_{BL}(\bar{\rho})=\begin{dcases}\frac{1}{2p_{F}}\int_{0}^{\infty}f\bigr{(}\bar{\rho}(1-h(s)),2\bar{\rho}p_{F}|\dot{h}(s)|\bigr{)}-f(\bar{\rho},0)\mathrm{d}s,&\mbox{
for Dirichlet BCs,}\\\
\frac{1}{2p_{F}}\int_{0}^{\infty}f\bigr{(}\bar{\rho}(1+h(s)),2\bar{\rho}p_{F}|\dot{h}(s)|\bigr{)}-f(\bar{\rho},0)\mathrm{d}s,&\mbox{
for Neumann BCs,}\end{dcases}$ $\displaystyle
c_{FM}(\bar{\rho})\coloneqq\frac{3\pi\bar{\rho}}{8p_{F}}\partial_{1}f(\bar{\rho},0).$

###### Proof.

In the periodic case, since $f$ is differentiable near $(\bar{\rho},0)$ for
any $\bar{\rho}>0$, the result follows from Theorem 4.2 and a simple Lipschitz
estimate.

As the Neumann and Dirichlet case are similar to each other, we give the
details only for the Dirichlet case. The proof is somewhat technical, due to
the fact that we have not assumed continuity of the derivative of $f$ with
respect to $|\nabla\rho|$ at $\rho=0$ (otherwise the small parameter $\delta$
in the proof below could be set to zero). But allowing a derivative
discontinuity at $\rho=0$ is essential to include the B88 functional, and
shows the robustness of our results.

First note that, since $f(\rho,|\nabla\rho|)$ depends only on $\rho$ and on
the norm of $|\nabla\rho|$, by the symmetries pointed out in Lemma 2.1 we can
reduce the integration domain to

$\displaystyle
F[\rho_{L}]=8\int_{Q^{D}_{\frac{L}{2}}}f(\rho_{L},|\nabla\rho_{L}|)\mathrm{d}r.$

Thus the continuum approximation from Theorem 4.2 reads

$\displaystyle\rho^{\textnormal{ctm}}_{L}(r)=\gamma^{\textnormal{ctm}}_{L}(r,r)=\bar{\rho}\frac{p_{L}^{3}}{p_{F}^{3}}\sum_{\sigma\in
G}\det\sigma h(p_{L}\lVert r-\sigma r\rVert)$ (5.1)

(since $\lVert r-\sigma r\rVert_{D,2L}=\lVert r-\sigma r\rVert$ for any $r\in
Q^{D}_{\frac{L}{2}}$ and $\sigma\in G$). The idea now is to use the continuum
approximation to show that $\rho_{L}$ is only small close to the faces and
edges of the box. For this, first note that $h(r)=1$ if and only if $r=0$.
Therefore, for any $\delta>0$ there exist $c(\delta)>0$ such that
$1-h(2p_{L}|r_{i}|)\geq c(\delta)$ as long as $|r_{i}|>\delta$ and $L$ is big
enough. Moreover, from (5.1) and the fact that $h$ decays at infinity, we see
that we can choose $R(\delta)>0$ such that

$\displaystyle
C\geq\rho^{\textnormal{ctm}}_{L}(r)\geq\frac{c(\delta)}{2}\quad\mbox{for any
$r\in Q_{L,\delta}$ and some $C>0$},$ (5.2)

where $Q_{L,\delta}$ is the region defined (using the projection $\pi_{j}$
from the notation section which removes the coordinate $r_{j}$) by

$\displaystyle Q_{L,\delta}\coloneqq\\{r\in Q^{D}_{\frac{L}{2}}:\min_{j\leq
3}|r_{j}|>\delta\quad\mbox{ and }\quad\min_{j\leq 3}|\pi_{j}r|>R(\delta)\\}.$

Next, note that by the continuum approximation estimate (see equation (4.22))
and estimate (5.2), we find that

$\displaystyle 2C\geq\rho_{L}(r)\geq\frac{c(\delta)}{4}\quad\mbox{for any
$r\in Q_{L,\delta}$ and $L$ big enough}.$

In particular, by the assumptions on $f$ ($C^{1}$ away from $a=0$ and $C^{0}$
up to $a=0$), there exists a Lipschitz constant $C(\delta)>0$ and a uniform
(with respect to $\delta$ and $L$) constant $C_{0}>0$ such that

$\displaystyle|F[\rho_{L}]-F[\rho^{\textnormal{ctm}}_{L}]|$
$\displaystyle\lesssim\int_{Q^{D}_{\frac{L}{2}}\setminus
Q_{L,\delta}}C_{0}+C(\delta)\int_{Q_{L,\delta}}|\rho_{L}-\rho^{\textnormal{ctm}}_{L}|+|\nabla\rho_{L}-\nabla\rho^{\textnormal{ctm}}_{L}|$
$\displaystyle\lesssim\delta
L^{2}+R(\delta)^{2}L+C(\delta)L^{\frac{34}{23}+\epsilon}.$

Therefore, by dividing the above by $L^{2}$, taking the limit as
$L\rightarrow\infty$, and then the limit $\delta\rightarrow 0$, one has

$\displaystyle
F[\rho_{L}]=F[\rho^{\textnormal{ctm}}_{L}]+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2}).$
(5.3)

The next step is to work with the continuum versions on $Q_{L,\delta}$ and use
a Taylor expansion of $f$ together with the decay of $h$ to determine the
asymptotic coefficients. To shorten the notation, let us define

$\displaystyle
f^{L}_{k}(r,t)\coloneqq\partial_{k}f\bigr{(}\bar{\rho}(1-t)+t\rho^{\textnormal{ctm}}_{L}(r),t|\nabla\rho^{\textnormal{ctm}}_{L}(r)|\bigr{)},\quad\quad
k\in\\{1,2\\},$

with the usual notation $\partial_{k}f$ for the partial derivative of $f$ with
respect to the $k$-th argument. Then, note that by estimate (5.2) and the
assumptions on $f$, we have

$\displaystyle|f_{k}^{L}(r,t)|\leq C(\delta),$ (5.4)

for any $(r,t)\in Q_{L,\delta}\times[0,1]$, for any $k\in\\{1,2\\}$, and for
some $C(\delta)>0$. In addition, by the fundamental theorem of calculus we
find

$\displaystyle F[\rho^{\textnormal{ctm}}_{L}]$
$\displaystyle=\underbrace{8\int_{Q_{\frac{L}{2}}^{D}\setminus
Q_{L,\delta}}f(\rho^{\textnormal{ctm}}_{L},|\nabla\rho^{\textnormal{ctm}}_{L}|)\mathrm{d}r}_{\lesssim\delta
L^{2}+R(\delta)^{2}L}+8\int_{Q_{L,\delta}}\biggr{(}f(\bar{\rho},0)+\int_{0}^{1}f^{L}_{1}(r,t)(\rho^{\textnormal{ctm}}_{L}(r)-\bar{\rho})\mathrm{d}t$
$\displaystyle+\int_{0}^{1}\underbrace{f_{2}^{L}(r,t)|\nabla\rho_{L}^{\rm
ctm}(r)|}_{=\frac{f_{2}^{L}(r,t)}{|\nabla\rho^{\textnormal{ctm}}_{L}(r)|}\sum_{j=1}^{3}\bigl{(}\partial_{r_{j}}\rho^{\textnormal{ctm}}_{L}(r)\bigr{)}^{2}}\mathrm{d}t\biggr{)}\mathrm{d}r$
$\displaystyle=\mathcal{O}(\delta
L^{2}+R(\delta)^{2}L)+f(\bar{\rho},0)|Q^{D}|L^{3}+\underbrace{8\bar{\rho}\int_{0}^{1}\int_{Q_{L,\delta}}\\!f^{L}_{1}(r,t)\biggr{(}\frac{p_{L}^{3}}{p_{F}^{3}}-1\biggr{)}\mathrm{d}r\mathrm{d}t}_{=:I(\delta)}$
$\displaystyle+\sum_{\sigma\neq
id}\underbrace{\frac{8\bar{\rho}p_{L}^{3}}{p_{F}^{3}}\int_{0}^{1}\int_{Q_{L,\delta}}\\!f^{L}_{1}(r,t)\det\sigma
h(p_{L}\lVert r-\sigma r\rVert)\mathrm{d}r\mathrm{d}t}_{=:J_{\sigma}(\delta)}$
$\displaystyle+\sum_{j=1}^{3}\underbrace{\frac{8\bar{\rho}^{2}p_{L}^{6}}{p_{F}^{6}}\int_{0}^{1}\int_{Q_{L,\delta}}\frac{f^{L}_{2}(r,t)}{|\nabla\rho^{\textnormal{ctm}}_{L}(r)|}\biggr{(}\sum_{\sigma_{jj}=-1}\det\sigma
p_{L}\dot{h}(p_{L}\lVert r-\sigma r\rVert)\frac{4r_{j}}{\lVert r-\sigma
r\rVert}\biggr{)}^{2}\mathrm{d}r\mathrm{d}t}_{=:K_{j}(\delta)}.$ (5.5)

Therefore, to complete the proof we need to study the ($L$ dependent) terms
$I(\delta)$, $J_{\sigma}(\delta)$ and $K(\delta)$. Physically, $I$ is a
correction coming from the expansion of the Fermi momentum, $J$ comes from the
density-dependence of $f$ (and is already present for the LDA), and $K$ comes
from the dependence of $f$ on the norm of the density gradient (and is absent
for the LDA). Let us start with $I(\delta)$. In this case, we first note that

$\displaystyle\frac{p_{L}^{3}}{p_{F}^{3}}-1=\frac{3\pi|\partial
Q^{D}|}{8p_{F}|Q^{D}|}\frac{1}{L}+\mathcal{O}(L^{-\frac{35}{23}+\epsilon})$

by Lemma 3.2. Hence by scaling out the $L$ in $I(\delta)$ and making the
following observations: (i) $f^{L}_{1}(r,t)$ is bounded in $Q_{\delta,L}$ by a
constant depending on $\delta$, but independent of $L$ (see (5.4)) , (ii)
$\lim_{L\rightarrow\infty}f^{L}_{1}(Lr,t)=\partial_{1}f(\bar{\rho},0)$ for
a.e. $(r,t)\in Q^{D}_{\frac{1}{2}}\times[0,1]$, and (iii)
$\lim_{L\rightarrow\infty}\chi_{Q_{L,\delta}}(Lr)=\chi_{Q^{D}_{\frac{1}{2}}}(r)$
in $L^{1}(\mathbb{R}^{3})$, we conclude that

$\displaystyle I(\delta)=\frac{3\pi\bar{\rho}|\partial
Q^{D}|}{8p_{F}}\partial_{1}f(\bar{\rho},0)L^{2}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2})=c_{FM}(\bar{\rho})|\partial
Q^{D}|L^{2}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2}).$
(5.6)

We consider next the terms $J_{\sigma}$ with $trace(\sigma)\leq-1$. In this
case, note that by the decay of $h$ (see Lemma 4.4), there exists some
$j\in\\{1,2,3\\}$, such that $|h(p_{L}\lVert r-\sigma
r\rVert)|\lesssim(1+\pi_{j}r)^{-2}$. Hence by estimate (5.4) we find that
$J_{\sigma}\lesssim_{\delta}L\log L$, and therefore, we just need to worry
about the terms $K_{j}$ and $J_{\sigma}$ with $trace(\sigma)=1$. For
simplicity let us label the reflection $\sigma$ with $trace(\sigma)=1$ and
$\sigma_{jj}=-1$ by $\sigma_{j}$. Now note that, by the decay of $h$ and
estimate (5.4), we have

$\displaystyle|f_{1}^{L}(r,t)h(2p_{L}r_{j})+f_{2}^{L}(r,t)\frac{|\partial_{r_{j}}\rho^{\textnormal{ctm}}(r)|^{2}}{|\nabla\rho^{\textnormal{ctm}}(r)|}|\lesssim\frac{C(\delta)}{(1+r_{j})^{2}}.$

As a consequence, up to an error bounded by $C(\delta)L$, one can change the
domain of integration of $K_{j}(\delta)+J_{\sigma_{j}}(\delta)$ from
$Q_{L,\delta}$ to

$\displaystyle
Q_{L,\delta}^{j}\coloneqq\\{r:r_{j}\in(0,\infty),\pi_{j}r\in\pi_{j}Q^{D}_{\frac{L}{2}},\quad\min_{j\leq
3}|r_{j}|\geq\delta,\mbox{ and }\min_{j\leq 3}|\pi_{j}r|>R(\delta)\\}.$

In summary, we have

$\displaystyle K_{j}(\delta)+J_{\sigma_{j}}(\delta)$
$\displaystyle=\frac{8\bar{\rho}p_{L}^{3}}{p_{F}^{3}}\int_{0}^{1}\int_{\mathbb{R}^{3}}\chi_{Q_{L,\delta}^{j}}(r)\biggr{(}\frac{\bar{\rho}p_{L}^{3}}{p_{F}^{3}}\frac{f_{2}^{L}(r,t)}{|\nabla\rho^{\textnormal{ctm}}_{L}(r)|}\biggr{(}\sum_{\sigma_{jj}=-1}\frac{p_{L}\det\sigma\dot{h}(p_{L}\lVert
r-\sigma
r\rVert)4r_{j}}{\lVert(1-\sigma)(r_{j},L\pi_{j}r)\rVert}\biggr{)}^{2}$
$\displaystyle+f_{1}^{L}(r,t)(-h(2p_{L}r_{j}))\biggr{)}\mathrm{d}r\mathrm{d}t+\mathcal{O}(C(\delta)L).$
(5.7)

Now note that by the decay of $h$, for a.e.
$(r_{j},\pi_{j}r)\in[0,\infty)\times\pi_{j}(Q^{D}_{\frac{1}{2}})$, the
following holds:

$\displaystyle\lim_{L\rightarrow\infty}\rho^{\textnormal{ctm}}_{L}(r_{j},L\pi_{j}r)-\bar{\rho}=-\bar{\rho}h(2p_{F}r_{j}),$
$\displaystyle\lim_{L\rightarrow\infty}|\nabla\rho^{\textnormal{ctm}}_{L}(r_{j},L\pi_{j}r)|=2\bar{\rho}p_{F}|\dot{h}(2p_{F}r_{j})|,$
$\displaystyle\lim_{L\rightarrow\infty}\sum_{\sigma_{jj}=-1}\det\sigma\bar{\rho}p_{L}\dot{h}\bigr{(}p_{L}\lVert(1-\sigma)(r_{j},L\pi_{j}r)\rVert\bigr{)}\frac{4r_{j}}{\lVert(1-\sigma_{j})(r_{j},L\pi_{j}r_{j})\rVert}=-2\bar{\rho}p_{F}\dot{h}(2p_{F}r_{j}),$
$\displaystyle\lim_{L\rightarrow\infty}\chi_{Q_{L,\delta}^{j}}(r_{j},L\pi_{j}r)=\chi_{(\delta,\infty)}(r_{j})\chi_{\pi_{j}(Q^{D}_{\frac{1}{2}})}(\pi_{j}r).$

Therefore, by scaling out the $L$ in the variables $\pi_{j}r$ in
$K_{j}(\delta)$ and noting that the integrand in equation (5.7) is bounded by
$C(\delta)\chi_{Q_{L,\delta}}(r)(1+|r_{j}|)^{-2}$ (by the decay of $h$), we
conclude from dominated convergence in $Q^{j}_{L,\delta}$ and the continuity
of $\nabla f$ in $Q^{j}_{L,\delta}$ that

$\displaystyle K_{j}(\delta)+J_{\sigma_{j}}(\delta)$
$\displaystyle=2\biggr{(}\prod_{\ell\neq
j}d_{\ell}\biggr{)}L^{2}\int_{\delta}^{\infty}\\!\int_{0}^{1}\\!\biggr{(}\partial_{2}f\bigr{(}\bar{\rho}(1-th(2p_{F}r_{j})),t2\bar{\rho}p_{F}|\dot{h}(2p_{F}r_{j})|\bigr{)}2\bar{\rho}p_{F}|\dot{h}(2p_{F}r_{j})|$
$\displaystyle-\partial_{1}f\bigr{(}\bar{\rho}(1-th(2p_{F}r_{j})),t2\bar{\rho}p_{F}|\dot{h}(2p_{F}r_{j})|\bigr{)}\bar{\rho}h(2p_{F}r_{j})\biggr{)}\mathrm{d}t\mathrm{d}r_{j}+\mathcal{O}(\delta
L^{2})+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}_{\delta}(L^{2})$
$\displaystyle=\frac{\prod_{\ell\neq
j}d_{\ell}}{p_{F}}L^{2}\int_{2p_{F}\delta}^{\infty}f\bigr{(}\bar{\rho}(1-h(r_{j}),2\bar{\rho}p_{F}|\dot{h}(r_{j})|\bigr{)}-f(\bar{\rho},0)\mathrm{d}r_{j}+\mathcal{O}(\delta
L^{2})+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}_{\delta}(L^{2}),$
(5.8)

where
$\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}_{\delta}$
emphasizes that the bounds may depend on $\delta$. As a consequence, the proof
follows from equations (5.3),(5.5),(5.6), and (5.8) by first taking the limit
$L\rightarrow\infty$ and then $\delta\rightarrow 0$. ∎

######  Remark (on Theorem 5.1).

  * •

Note that $c_{FM}(\bar{\rho})$ above depends on $\bar{\rho}$ which is in
contrast with the constant $c_{FM}$ defined in Theorem 5.2 below. The reason
is that for the “physically” relevant cases (see the discussion before Theorem
1.1), we have $c_{FM}(\bar{\rho})=-c_{FM}\bar{\rho}$, where
$c_{FM}=\frac{3}{8}$ is precisely the value defined there (see Corollary 5.1
below).

  * •

The same arguments can be used for semi-local functionals of higher order
derivatives by considering generalized variables like
$(\rho,\nabla\rho,...,\partial^{\alpha}\rho)$ and applying estimate (1.15). As
we are only interested in LDA and GGAs for the moment, we leave the
asymptotics of functionals for higher order derivatives for future works

  * •

If $f$ is more regular and one knows the ratio at which the derivatives
diverge when $(a,b)\rightarrow 0$, one can further use the decay of $h$ to
improve the remainder term from
$\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2})$
(for Dirichlet and Neumann cases) to, possibly,
$\mathcal{O}(L^{\frac{34}{23}+\epsilon})$.

As straightforward corollaries of Theorem 5.1, we obtain the two-term
asymptotic expansion for LDA and GGAs from Theorem 1.1.

######  Corollary 5.1 (Asymptotics of LDA).

Let $\rho_{L}$ be the single-particle density of any determinantal ground
state of the free $N$-electron gas in $Q_{L}^{D}$ (under our usual boundary
conditions). Then in the thermodynamic limit we have

$\displaystyle
E_{x}^{LDA}[\rho_{L}]=\begin{dcases}-c_{x}\bar{\rho}^{4/3}|Q^{D}|L^{3}+\mathcal{O}(L^{\frac{34}{23}+\epsilon}),&\mbox{for
periodic BCs,}\\\
-c_{x}\bar{\rho}^{4/3}|Q^{D}|L^{3}-c_{LDA}^{\textnormal{Dir}}\bar{\rho}|\partial
Q^{D}|L^{2}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2}),&\mbox{for
Dirichlet BCs,}\\\
-c_{x}\bar{\rho}^{4/3}|Q^{D}|L^{3}-c_{LDA}^{\textnormal{Neu}}\bar{\rho}|\partial
Q^{D}|L^{2}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2}),&\mbox{for
Neumann BCs,}\end{dcases}$

where the constants are

$\displaystyle
c_{LDA}=\begin{dcases}c_{FM}+\frac{3}{8\pi}\int_{0}^{\infty}(1-h(s))^{\frac{4}{3}}-1\mathrm{d}s,&\mbox{for
Dirichlet BCs,}\\\
-c_{FM}+\frac{3}{8\pi}\int_{0}^{\infty}(1+h(s))^{\frac{4}{3}}-1\mathrm{d}s,&\mbox{for
Neumann BCs,}\end{dcases}$

with $c_{FM}=\frac{3}{8}$. (Compare with the constant in Theorem 5.2.)

######  Corollary 5.2 (Asymptotics for GGA).

Let $g^{GGA}(a,b)\in C^{0}([0,\infty)^{2})\cap
C^{1}((0,\infty)\times[0,\infty))$ such that $g^{GGA}(a,0)=0$, for all $a\geq
0$. Moreover, let $\rho_{L}$ be the single-particle density of any
determinantal ground state of the free $N$-electron gas in $Q_{L}^{D}$ (under
our usual boundary conditions). Then, for $\Delta
E_{x}^{GGA}[\rho_{L}]=\int_{Q^{D}_{L}}g^{GGA}(\rho_{L},|\nabla\rho_{L}|)$, in
the thermodynamic limit we have

$\displaystyle\Delta
E_{x}^{GGA}[\rho_{L}]=\begin{dcases}\mathcal{O}(L^{-\frac{34}{23}+\epsilon})&\mbox{for
periodic BCs,}\\\ c^{\textnormal{Dir}}_{GGA}(\bar{\rho})|\partial
Q^{D}|L^{2}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2}),&\mbox{for
Dirichlet BCs,}\\\ c_{GGA}^{\textnormal{Neu}}(\bar{\rho})|\partial
Q^{D}|L^{2}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(L^{2}),&\mbox{for
Neumann BCs,}\end{dcases}$

where the constants are given by

$\displaystyle
c_{GGA}(\bar{\rho})=\begin{dcases}\frac{1}{2p_{F}}\int_{0}^{\infty}g^{GGA}\biggr{(}\bar{\rho}(1-h(s)),2\bar{\rho}p_{F}|\dot{h}(s)|\biggl{)}\mathrm{d}s,&\mbox{
for Dirichlet BCs,}\\\
\frac{1}{2p_{F}}\int_{0}^{\infty}g^{GGA}\biggr{(}\bar{\rho}(1+h(s)),2\bar{\rho}p_{F}|\dot{h}(s)|\biggl{)}\mathrm{d}s,&\mbox{
for Neumann BCs.}\end{dcases}$

###### Proof of Corollary 5.1 and 5.2.

In the first corollary, just apply Theorem 5.1 to $f(\rho)=c_{x}\rho^{4/3}$.
For the second one, set $f(\rho,|\nabla\rho|)=g^{GGA}(\rho,|\nabla\rho|)$ and
note that $f(a,0)=0$ for all $a\geq 0$ implies that $\partial_{1}f(a,0)=0$ for
all $a>0$. ∎

###  5.2 Exact exchange

Now we turn to the asymptotic expansion of the exact exchange. The analysis
here is different as we are dealing with a functional of the density matrix,
not the density.

######  Theorem 5.2 (Asymptotics of exact exchange).

Let $\bar{\rho}=\frac{N}{|Q^{D}|L^{3}}=constant$, $\Psi_{L}$ be a
determinantal ground state of the free $N$-electron gas in $Q_{L}^{D}$. Then,
we have

$\displaystyle
E_{x}[\Psi_{L}]=\begin{dcases}-c_{x}\bar{\rho}^{\frac{4}{3}}|Q^{D}|L^{3}+c_{FS}\bar{\rho}|\partial
Q^{D}|L^{2}+\mathcal{O}(L^{\frac{45}{23}+\epsilon}),&\mbox{for periodic
BCs,}\\\
-c_{x}\bar{\rho}^{\frac{4}{3}}|Q^{D}|L^{3}-(c_{BL}^{\textnormal{Dir}}+c_{FM}-c_{FS})\bar{\rho}|\partial
Q^{D}|L^{2}+\mathcal{O}(L^{\frac{45}{23}+\epsilon}),&\mbox{for Dirichlet
BCs,}\\\
-c_{x}\bar{\rho}^{\frac{4}{3}}|Q^{D}|L^{3}-(c_{BL}^{\textnormal{Neu}}-c_{FM}-c_{FS})\bar{\rho}|\partial
Q^{D}|L^{2}+\mathcal{O}(L^{\frac{45}{23}+\epsilon}),&\mbox{for Neumann
BCs,}\end{dcases}$

where the finite-size, Fermi momentum and boundary layer corrections are

$\displaystyle c_{FS}=\frac{1}{8},\quad c_{FM}=\frac{3}{8},\quad
c_{BL}^{\textnormal{Dir}}=-\frac{\log 2}{4}\mbox{, and}\quad
c_{BL}^{\textnormal{Neu}}=\frac{3\log 2}{4}.$

###### Proof.

As before, we prove the Dirichlet case in detail and outline the proof for the
other two boundary conditions at the end.

The first step is again to justify the use of the continuum density matrices
(4.22). For this, we use the identity

$|\gamma_{L}(r,\tilde{r})|^{2}-|\gamma^{\textnormal{ctm}}_{L}(r,\tilde{r})|^{2}=Re\\{\overline{(\gamma_{L}(r,\tilde{r})-\gamma^{\textnormal{ctm}}_{L}(r,\tilde{r}))}(\gamma_{N}(r,\tilde{r})+\gamma^{\textnormal{ctm}}_{N}(r,\tilde{r}))\\}$

together with estimates (4.21) and (4.23) from Theorem 4.2 to obtain

$\displaystyle\biggr{|}\int_{Q_{L}^{D}\times
Q_{L}^{D}}\frac{|\gamma_{L}(r,\tilde{r})|^{2}-|\gamma^{\textnormal{ctm}}_{L}(r,\tilde{r})|^{2}}{|r-\tilde{r}|}\biggl{|}\mathrm{d}r\mathrm{d}\tilde{r}$
$\displaystyle\lesssim L^{-\frac{70}{23}+\epsilon}\int_{Q_{L}^{D}\times
Q_{L}^{D}}|r-\tilde{r}|^{-1}\mathrm{d}r\mathrm{d}\tilde{r}$
$\displaystyle+L^{-\frac{35}{23}+\epsilon}\int_{Q_{L}^{D}\times
Q_{L}^{D}}\frac{(1+|r-\tilde{r}|_{D,2L})^{-2}}{|r-\tilde{r}|}\mathrm{d}r\mathrm{d}\tilde{r}$
$\displaystyle\lesssim
L^{\frac{45}{23}+\epsilon}+L^{\frac{34}{23}+\epsilon}\log L.$

Therefore, by the continuum formulas we have

$\displaystyle E_{x}[\rho_{L}]\approx-\sum_{\tau,\sigma\in
G}\det(\sigma\tau)\frac{\bar{\rho}^{2}p_{L}^{6}}{4p_{F}^{6}}\underbrace{\int_{Q_{L}^{D}\times
Q_{L}^{D}}\frac{h(p_{L}|r-\sigma\tilde{r}|_{D,2L})h(p_{L}|r-\tau\tilde{r}|_{D,2L})}{|r-\tilde{r}|}\mathrm{d}r\mathrm{d}\tilde{r}}_{\coloneqq
J_{\sigma,\tau}(L)},$ (5.9)

where $\approx$ will be used throughout this proof to denote equality up to
errors included in the remainder of Theorem 5.2. Now, to estimate the terms
$J_{\sigma,\tau}$ we start with the following lemma.

######  Lemma 5.1.

If $trace(\sigma)\leq-1$ or $trace(\tau)\leq-1$, then
$J_{\sigma,\tau}(L)\lesssim L(\log L)^{2}$. Furthermore, if
$trace(\sigma)=1=trace(\tau)$ and $\sigma\neq\tau$, we have
$J_{\sigma,\tau}(L)\lesssim L(\log L)^{4}$.

###### Proof.

First observe that since $r_{j},\tilde{r}_{j}\in[0,d_{j}L]$, one has

$\displaystyle|r_{j}+\tilde{r}_{j}\mod
2d_{j}L|=\min\\{r_{j}+\tilde{r}_{j},2d_{j}L-r_{j}-\tilde{r}_{j}\\}\geq|r_{j}-\tilde{r}_{j}|.$

In particular, by the decay of $h$, we see that

$\displaystyle|h(p_{L}|r-\sigma\tilde{r}|_{D,2L})|\lesssim\frac{1}{(1+|r-\sigma\tilde{r}|_{D,2L})^{2}}\lesssim\frac{1}{(1+|r-\tilde{r}|)^{2}},$
(5.10)

for any $\sigma\in G$. On the other hand, if $trace(\sigma)\leq-1$, there
exists $j\in\\{1,2,3\\}$ such that

$\displaystyle|r-\sigma\tilde{r}|\geq\min_{\begin{subarray}{c}p\in\pi_{j}(D(2L\mathbb{Z})^{3})\\\
p\in\pi_{j}(Q^{D}_{2L})\end{subarray}}|\pi_{j}(r+\tilde{r})-p|,$

and therefore,

$\displaystyle|h(p_{L}|r-\sigma\tilde{r}|_{D,2L})|\lesssim\sum_{\begin{subarray}{c}p\in\pi_{j}(D(2L\mathbb{Z})^{3})\\\
p\in\pi_{j}(Q^{D}_{2L})\end{subarray}}\frac{1}{(1+|\pi_{j}(r+\tilde{r})-p|)^{2}}.$
(5.11)

As a result, assuming that $trace(\sigma)\leq-1$, we can see from (5.10) and
(5.11) that

$\displaystyle J_{\sigma,\tau}(L)\lesssim\int_{Q^{D}_{L}\times
Q^{D}_{L}}\sum_{\begin{subarray}{c}p\in\pi_{j}(D(2L\mathbb{Z})^{3})\\\
p\in\pi_{j}(Q^{D}_{2L})\end{subarray}}\frac{1}{(1+|\pi_{j}(r+\tilde{r})-p|)^{2}}\frac{1}{(1+|r-\tilde{r}|)^{2}}\frac{1}{|r-\tilde{r}|}\mathrm{d}r\mathrm{d}\tilde{r}\lesssim
L(\log L)^{2}.$

For the terms $J_{\sigma,\tau}$ with $trace(\sigma)=1=trace(\tau)$ and
$\sigma\neq\tau$, note that there exists $j\neq\ell\in\\{1,2,3\\}$ such that

$\displaystyle|h(p_{L}|r-\sigma\tilde{r}|)|\lesssim\sum_{\begin{subarray}{c}p_{j}\in\\{0,2d_{j}L\\}\end{subarray}}\frac{1}{(1+|r_{j}+\tilde{r}_{j}-p_{j}|)(1+|\pi_{j}(r-\tilde{r}|)},$
$\displaystyle|h(p_{L}|r-\tau\tilde{r}|)|\lesssim\sum_{\begin{subarray}{c}p_{\ell}\in\\{0,2d_{\ell}L\\}\end{subarray}}\frac{1}{(1+|r_{\ell}+\tilde{r}_{\ell}-p_{\ell}|)(1+|\pi_{\ell}(r-\tilde{r})|)}.$

The lemma thus follows by integrating the product of the estimates above
against the Coulomb potential in the box $Q_{L}^{D}\times Q_{L}^{D}$. ∎

From the lemma above, it is enough to study
$J_{\textnormal{id},\textnormal{id}}(L)$, $J_{\sigma,\textnormal{id}}(L)$ and
$J_{\sigma,\sigma}(L)$ where $trace(\sigma)=1$.

We start with $J_{\textnormal{id},\textnormal{id}}(L)$. In this case, we first
note that $|r-\tilde{r}|_{D,2L}=|r-\tilde{r}|$ for any $r,\tilde{r}\in
Q^{D}_{L}$. In particular, by the change of variables
$(w(r,\tilde{r}),\tilde{w}(r,\tilde{r}))=(r-\tilde{r},\tilde{r})$, we obtain

$\displaystyle J_{\textnormal{id},\textnormal{id}}(L)$
$\displaystyle=|Q^{D}|L^{3}\int_{Q_{L}^{D}-Q_{L}^{D}}\frac{h(p_{L}|w|)^{2}}{|w|}\mathrm{d}w-L^{2}\sum_{j=1}^{3}\biggr{(}\prod_{\ell\neq
j}d_{\ell}\biggr{)}\int_{Q^{D}_{L}-Q_{L}^{D}}\frac{h(p_{L}|w|)^{2}|w_{j}|}{|w|}\mathrm{d}w$
$\displaystyle+\int_{Q^{D}_{L}-Q^{D}_{L}}\frac{h(p_{L}|w|)^{2}}{|w|}\biggr{(}L\sum_{j=1}^{3}d_{j}\prod_{\ell\neq
j}|w_{\ell}|-\prod_{j=1}^{3}|w_{j}|\biggr{)}\mathrm{d}w.$

By the decay of $h$ and a simple estimate, up to an error $\lesssim L\log L$
the integral on the first two terms can be taken over the whole of
$\mathbb{R}^{3}$, and the third integral can be neglected. In addition, by
scaling out $p_{L}$, using spherical coordinates, and recalling that
$|\partial Q^{D}|=2(d_{1}d_{2}+d_{1}d_{3}+d_{2}d_{3})$, we have

$\displaystyle
J_{\textnormal{id},\textnormal{id}}(L)\approx\frac{4\pi|Q^{D}|L^{3}}{p_{L}^{2}}\underbrace{\int_{0}^{\infty}h(r)^{2}r\mathrm{d}r}_{\coloneqq
I_{0}}-\frac{\pi|\partial
Q^{D}|L^{2}}{p_{L}^{3}}\underbrace{\int_{0}^{\infty}h(r)^{2}r^{2}\mathrm{d}r}_{\coloneqq
I_{1}}.$ (5.12)

Next, for the term $J_{\sigma,\sigma}$ we assume without loss of generality
that $\sigma_{11}=-1$. Then, by invariance of the integrand under the
reflections $r_{i}-\tilde{r}_{i}\mapsto\tilde{r}_{i}-r_{i}$, the change of
variables
$(w(r,\tilde{r}),\tilde{w}(r,\tilde{r}))=(r-\tilde{r},r_{1}+\tilde{r}_{1},\tilde{r}_{2},\tilde{r}_{3})$
(notice
$\frac{\mathrm{d}w\,\mathrm{d}\tilde{w}}{2}=\mathrm{d}r\,\mathrm{d}\tilde{r}$),
and the decay of $h$, we have

$\displaystyle J_{\sigma,\sigma}(L)$
$\displaystyle=\int_{Q^{D}_{L}-Q^{D}_{L}}\int_{|w_{1}|}^{2d_{1}L-|w_{1}|}\frac{h(p_{L}|(\tilde{w}_{1}\mod
2d_{1}L,\pi_{1}w)|)^{2}\prod_{\ell\neq
1}(d_{\ell}L-|w_{\ell}|)}{|w|}\frac{\mathrm{d}\tilde{w}_{1}\mathrm{d}w}{2}$
$\displaystyle=8d_{2}d_{3}L^{2}\int_{Q^{D}_{L}}\int_{w_{1}}^{d_{1}L}\frac{h(p_{L}|(\tilde{w}_{1},\pi_{1}w)|)^{2}}{|w|}\mathrm{d}\tilde{w}_{1}\mathrm{d}w+\mathcal{O}(L\log
L)$
$\displaystyle=8d_{2}d_{3}L^{2}\int_{Q^{D}_{L}}h(p_{L}|(\tilde{w}_{1},\pi_{1}w)|)^{2}\biggr{(}\int_{0}^{\tilde{w}_{1}}\frac{1}{|w|}\mathrm{d}w_{1}\biggr{)}\mathrm{d}\tilde{w}_{1}\mathrm{d}\pi_{1}w+\mathcal{O}(L\log
L)$

(where we inverted the order of integration between $\tilde{w}_{1}$ and
$w_{1}$ in the last step). In addition, since
$\int_{0}^{|\tilde{w}_{1}}\frac{1}{|w|}\mathrm{d}w_{1}=\frac{1}{2}\log\bigr{(}\frac{\tilde{w}_{1}+|(\tilde{w}_{1},\pi_{1}w)|}{|(\tilde{w}_{1},\pi_{1}w)|-\tilde{w}_{1}}\bigr{)}\lesssim\frac{|\tilde{w}_{1}|}{|\pi_{1}w|}$,
by the decay of $h$ one can see that, up to an error $\lesssim L$, we can
replace the domain of the outer integration by $\mathbb{R}_{+}^{3}$. Hence, by
scaling out $p_{L}$, changing to spherical coordinates, and recalling the
definition of $I_{1}$ in (5.12), we find that

$\displaystyle J_{\sigma,\sigma}(L)$ $\displaystyle\approx\frac{2\pi
d_{2}d_{3}}{p_{L}^{3}}L^{2}\int_{0}^{\infty}h(r)^{2}r^{2}\mathrm{d}r\underbrace{\int_{0}^{\frac{\pi}{2}}\log\biggr{(}\frac{1+\cos\theta}{1-\cos\theta}\biggr{)}\sin\theta\mathrm{d}\theta}_{\mathclap{=(\cos\theta-1)\log(1-\cos\theta)-(1+\cos\theta)\log(1+\cos\theta)}}=\frac{4\pi\log(2)}{p_{L}^{3}}d_{2}d_{3}I_{1}L^{2}.$
(5.13)

At last, for $J_{\textnormal{id},\sigma}$ (again assuming that
$\sigma_{11}=-1$ without loss of generality), by the change of variables
$(w(r,\tilde{r}),\tilde{w}(r,\tilde{r}))=(r-\tilde{r},r_{1}+\tilde{r}_{1},\tilde{r}_{2},\tilde{r}_{3})$
and the same arguments as before we conclude that

$\displaystyle J_{\textnormal{id},\sigma}$
$\displaystyle=\int_{Q^{D}_{L}-Q^{D}_{L}}\int_{|w_{1}|}^{d_{1}L}\frac{h(p_{L}|(\tilde{w}_{1},\pi_{1}w)|)h(p_{L}|w|)\prod_{\ell\neq
1}(d_{\ell}L-|w_{\ell}|)}{|w|}\mathrm{d}\tilde{w}_{1}\mathrm{d}w$
$\displaystyle=8d_{2}d_{3}L^{2}\int_{\mathbb{R}_{+}^{3}}\int_{|w_{1}|}^{\infty}\frac{h(p_{L}|(\tilde{w}_{1},\pi_{1}w)|)h(p_{L}|w|)}{|w|}\mathrm{d}\tilde{w}_{1}\mathrm{d}w+\mathcal{O}(L\log
L).$

Hence, scaling out $p_{L}$ and using spherical coordinates for $w$, we have

$\displaystyle J_{\textnormal{id},\sigma}(L)\approx\frac{4\pi
d_{2}d_{3}}{p_{L}^{3}}L^{2}\underbrace{\int_{0}^{\infty}\int_{0}^{\frac{\pi}{2}}\int_{r\cos\theta}^{\infty}h(\sqrt{(r\sin\theta)^{2}+\tilde{w}_{1}^{2}})h(r)r\sin\theta\mathrm{d}\tilde{w}_{1}\mathrm{d}\theta\mathrm{d}r}_{\coloneqq
I_{2}}.$ (5.14)

Thus by plugging (5.12),(5.13) and (5.14) into (5.9), we have

$\displaystyle E_{x}[\rho_{L}]$
$\displaystyle\approx-\frac{\bar{\rho}^{2}p_{L}^{6}}{4p_{F}^{6}}\biggr{(}J_{\textnormal{id},\textnormal{id}}(L)+\sum_{j=1}^{3}J_{\sigma_{j},\sigma_{j}}(L)-\sum_{j=1}^{3}\bigr{(}J_{\textnormal{id},\sigma_{j}}(L)+J_{\sigma_{j},\textnormal{id}}(L)\bigr{)}\biggr{)}$
$\displaystyle\approx-\frac{\pi\bar{\rho}^{2}p_{L}^{4}}{p_{F}^{6}}I_{0}|Q^{D}|L^{3}-\biggr{(}-\frac{\pi\bar{\rho}^{2}p_{L}^{3}}{4p_{F}^{6}}I_{1}+\frac{\pi\log
2\bar{\rho}^{2}p_{L}^{3}}{2p_{F}^{6}}I_{1}-\frac{\pi\bar{\rho}^{2}p_{L}^{3}}{p_{F}^{6}}I_{2}\biggr{)}|\partial
Q^{D}|L^{2}.$

As a result, using Lemma 3.2 to replace $p_{L}$ by $p_{F}$ plus correction, we
conclude that

$\displaystyle
E_{x}[\rho_{L}]\approx-\underbrace{\frac{\pi}{(3\pi^{2})^{\frac{2}{3}}}I_{0}}_{c_{x}}\bar{\rho}^{\frac{4}{3}}|Q^{D}|L^{3}-\biggr{(}\underbrace{\frac{1}{6}I_{0}}_{c_{FM}}-\underbrace{\frac{1}{12\pi}I_{1}}_{c_{FS}}+\underbrace{\frac{\log
2}{6\pi}I_{1}-\frac{1}{3\pi}I_{2}}_{c_{BL}}\biggr{)}\bar{\rho}|\partial
Q^{D}|L^{2}.$

The proof is completed by using the values of the integrals computed in the
next lemma.

######  Lemma 5.2.

$I_{0}=\frac{9}{4}$, $I_{1}=\frac{3\pi}{2}$ and $I_{2}=\frac{3\pi\log 2}{2}$.

###### Proof.

For $I_{0}$, we can use the identity $\int\frac{(\sin s-s\cos
s)^{2}}{s^{5}}\mathrm{d}s=\frac{-2s^{2}+2s\sin(2s)+\cos(2s)-1}{s^{4}}$ (see
e.g., [7, Lemma 6.1] for an elegant evaluation) to obtain

$\displaystyle I_{1}=\int_{0}^{\infty}9\frac{(\sin s-s\cos
s)^{2}}{s^{5}}\mathrm{d}s=-9\lim_{s\rightarrow
0}\frac{-2s^{2}+2s\sin(2s)+\cos(2s)-1}{s^{4}}=\frac{9}{4},$

where the limit can be computed by L’Hôpital’s rule. Next, note that for any
$a>0$, by Plancherel’s theorem, we have

$\displaystyle\int_{0}^{\infty}h(r)h(ar)\,r^{2}\mathrm{d}r$
$\displaystyle=\frac{1}{4\pi|B_{1}|^{2}}\int_{\mathbb{R}^{3}}\widehat{\chi}_{B_{1}}(k)\widehat{\chi}_{B_{1}}(ak)\mathrm{d}k$
$\displaystyle=\frac{(2\pi)^{3}}{4\pi|B_{1}|^{2}a^{3}}\int_{\mathbb{R}^{3}}\chi_{B_{1}}(k)\chi_{B_{1}}\biggr{(}\frac{k}{a}\biggr{)}\mathrm{d}k$
$\displaystyle=\frac{3\pi}{2\max\\{a,1\\}^{3}}.$ (5.15)

In particular, the value of $I_{1}$ follows by setting $a=1$. For $I_{2}$,
first note that by using the inverse of polar coordinates
$(r,\theta)=(\sqrt{x^{2}+y^{2}},\arctan(y/x))$, one has

$\displaystyle
I_{2}=\int_{\mathbb{R}_{+}^{2}}\int_{x}^{\infty}h\bigr{(}\sqrt{y^{2}+\tilde{w}_{1}^{2}}\bigr{)}h\bigr{(}\sqrt{x^{2}+y^{2}}\bigr{)}\frac{y}{\sqrt{x^{2}+y^{2}}}\mathrm{d}\tilde{w}_{1}\mathrm{d}x\mathrm{d}y.$

Next, set
$\mathcal{R}\coloneqq\\{(x,y,\tilde{w}_{1})\in\mathbb{R}^{3}:x>0,y>0,\tilde{w}_{1}>x\\}$
and consider the change of coordinates

$\displaystyle
T:\biggr{\\{}(s,\phi,v)\in\mathbb{R}^{3}:s>0,0<\phi<\frac{\pi}{2},0<v<\log\biggr{(}\frac{\sin\phi+1}{\cos\phi}\biggr{)}\biggr{\\}}\rightarrow\mathcal{R}$
$\displaystyle(s,\phi,v)\mapsto(x,y,\tilde{w}_{1})=T(s,\phi,v)=(s\cos\phi\sinh
v,s\sqrt{\sin^{2}\phi-\cos^{2}\phi\sinh^{2}v},s\cos\phi\cosh v),$

for which

$\displaystyle\det
DT(s,\phi,v)=s^{2}\frac{\sin(\phi)\cos(\phi)}{\sqrt{\sin^{2}\phi-\cos^{2}\phi\sinh^{2}v}}=s^{2}\cos\phi\frac{\sqrt{x^{2}+y^{2}}}{y}.$

Then, by (5.15) and the substituion $\tau=\sin\phi$, we conclude that

$\displaystyle I_{2}$
$\displaystyle=\int_{0}^{\frac{\pi}{2}}\cos\phi\,\log\biggr{(}\frac{\sin\phi+1}{\cos\phi}\biggr{)}\underbrace{\int_{0}^{\infty}h(s)h(s\sin\phi)s^{2}\mathrm{d}s}_{=\frac{3\pi}{2}}\mathrm{d}\phi$
$\displaystyle=\frac{3\pi}{2}\int_{0}^{1}\log\biggr{(}\frac{1+\tau}{\sqrt{1-\tau^{2}}}\biggr{)}\mathrm{d}\tau=\frac{3\pi}{2}\int_{0}^{1}\log\biggr{(}\frac{\sqrt{1+\tau}}{\sqrt{1-\tau}}\biggr{)}\mathrm{d}\tau=\frac{3\pi\log
2}{2},$

which completes the proof of Theorem 5.2 for the Dirichlet case. For the
Neumann case the same arguments work with the proper change of signs (due to
the missing $\det\sigma$ in (2.9)). ∎

For the periodic case, we first note that by the decay of $h$,

$\displaystyle|h(p_{L}|r-r^{\prime}|_{L,D}|)^{2}-\sum_{\begin{subarray}{c}p\in
DL\mathbb{Z}^{3}\\\ |D^{-1}p|_{\max}\leq
L\end{subarray}}h(p_{L}|r-r^{\prime}-p|)^{2}\lesssim L^{-4}.$

Moreover, one can now show that

$\displaystyle\int_{Q^{D}_{L}\times
Q^{D}_{L}}\frac{h(p_{L}|r-\tilde{r}-p|)^{2}}{|r-\tilde{r}|}\mathrm{d}r\mathrm{d}\tilde{r}\lesssim
L\log L\quad\mbox{for any $p\in DL\mathbb{Z}^{3}\setminus\\{0\\}$,}$

which implies that
$E_{x}[\Psi_{L}]\approx-\frac{\bar{\rho}^{2}p_{L}^{6}}{4p_{F}^{6}}J_{id,id}(L)$.
The result now follows from Lemma 3.2, estimate (5.12) and Lemma 5.2. ∎

##  6 Kinetic energy

In this section we use Lemma 1.1 to compute the asymptotic expansion of the
kinetic energy.

######  Theorem 6.1 (Asymptotics of kinetic energy).

Let $D$ fixed, $\bar{\rho}=\frac{N}{|Q^{D}|L^{3}}=constant$ and $\Psi_{N,L}$
be any ground state of the free $N$-electron gas in the box $Q^{D}_{L}$ under
our usual boundary conditions. Let $T$ be the kinetic energy functional
defined in (1.1). Then, we have

$\displaystyle
T[\Psi_{N,L}]=\begin{dcases}c_{TF}\bar{\rho}^{5/3}|Q^{D}|L^{3}+\mathcal{O}(L^{\frac{34}{23}+\epsilon}),&\mbox{for
periodic BCs,}\\\
c_{TF}\bar{\rho}^{5/3}|Q^{D}|L^{3}+c_{K}\bar{\rho}^{4/3}|\partial
Q^{D}|L^{2}+\mathcal{O}(L^{\frac{34}{23}+\epsilon}),&\mbox{for Dirichlet
BCs,}\\\ c_{TF}\bar{\rho}^{5/3}|Q^{D}|L^{3}-c_{K}\bar{\rho}^{4/3}|\partial
Q^{D}|L^{2}+\mathcal{O}(L^{\frac{34}{23}+\epsilon}),&\mbox{for Neumann
BCs,}\end{dcases}$

where $c_{TF}=\frac{3}{10}(3\pi^{2})^{2/3}$ is the Thomas-Fermi constant and
$c_{K}=\frac{3\pi}{32}(3\pi^{2})^{1/3}$.

###### Proof.

For the Dirichlet case, just note that

$\displaystyle T[\Psi_{N,L}]$
$\displaystyle=\langle\Psi_{N,L},-\frac{\Delta}{2}\Psi_{N,L}\rangle=\\!\\!\\!\sum_{k\in\mathbb{N}^{3}\cap
B^{D}_{R_{N}}}\\!\\!\\!\frac{\pi^{2}|D^{-1}k|^{2}}{L^{2}}-(N_{+}-N)\frac{\pi^{2}R_{N}^{2}}{L^{2}}$
$\displaystyle=\frac{\pi^{2}}{L^{2}}\frac{1}{8}\biggr{(}\sum_{k\in\mathbb{Z}^{3}\cap
B^{D}_{R_{N}}}\\!\\!\\!|D^{-1}k|^{2}-\sum_{j=1}^{3}\sum_{\begin{subarray}{c}k\in\mathbb{Z}^{3}\cap
B^{D}_{R_{N}}\\\
k_{j}=0\end{subarray}}|D^{-1}k|^{2}\biggr{)}+\mathcal{O}(L^{\frac{34}{23}+\epsilon}).$

Moreover, by a simple estimate999In fact, by adapting the proof of [7], one
can get a remainder of order $\mathcal{O}(R^{2+\frac{2}{3}})$ in estimate
(6.1). we find that

$\displaystyle\sum_{\begin{subarray}{c}k\in B^{D}_{R_{N}}\cap\mathbb{Z}^{3}\\\
k_{j}=0\end{subarray}}|D^{-1}k|^{2}=\int_{\mathbb{R}^{2}}\chi_{\pi_{j}(B^{D}_{R_{N}})}(k)\sum_{\ell\neq
j}(d_{\ell}^{-1}k_{\ell})^{2}\mathrm{d}k+\mathcal{O}(R_{N}^{3})=\frac{\pi}{4}R_{N}^{4}\prod_{\ell\neq
j}d_{\ell}+\mathcal{O}(R_{N}^{3}).$ (6.1)

Thus using estimate (6.1) and Lemma 1.1, we conclude that

$\displaystyle
T[\Psi_{N,L}]=\frac{1}{10\pi^{2}}p_{N,L}^{5}|Q^{D}|L^{3}-\frac{1}{32\pi}p_{N,L}^{4}|\partial
Q^{D}|L^{2}+\mathcal{O}(L^{\frac{34}{23}+\epsilon}).$

The result now follows from the asymptotics of $p_{L}$ in Lemma 3.2. The
Neumann and periodic cases are entirely analogous. ∎

######  Remark.

Note that we do not assume $\Psi_{N,L}$ to be a determinantal ground state as
the kinetic energy is simply the ground state energy of the Laplacian and
therefore unique (even if the ground state is not).

##  7 Numerical results and final discussion

We now compare our asymptotic results to numerical values of different
exchange functionals for the free electron gas with zero boundary conditions,
for up to $30~{}\\!000$ electrons. Our numerical computations were carried out
in Matlab. All energy functionals other than exact exchange were evaluated by
direct numerical integration of the exact formulas given in section 2. For
exact exchange, accurate direct numerical evaluation of the expression (2.11),
(2.9) is impossible, because of the high-dimensionality of the domain of
integration (6D) and the ${1}/{|r-r^{\prime}|}$ singularity of the integrand.
We tackled these obstructions by reducing the problem to the numerical
computation of a small ($O(N^{1/3})$) number of one-dimensional integrals of
smooth functions (see appendix A for a detailed description). Moreover, we
focus here on the case of a cubic box $[0,L]^{3}$ and $\bar{\rho}=1$.

To begin with, in Figure 3 we have plotted the exact exchange energy per unit
volume, together with the theoretical one-term (just $c_{x}$) and two-term
asymptotics ($c_{x}+c_{x,2}\cdot 6L^{-1}$) from Theorem 1.1. For comparison we
have also included the LDA exchange energy per unit volume.

Figure 3: Exact exchange energy per unit volume of the free electron gas in a
box with zero boundary conditions, compared with the LDA, one-term asymptotics
(Dirac exchange constant) and two-term asymptotics (present work, Theorem
1.1).

Note that even for small $N$, the two-term asymptotics is a much better
approximation than the one-term asymptotics, and also a better one than the
more complicated LDA. Note that the latter requires integration of an
inhomogeneous $N$-electron exchange energy density of the system.

Let us now look in more detail at the next-order contribution. Besides exact
exchange and the LDA, we consider the widely used GGAs introduced by Becke in
1988 (B88) [1] and Perdew, Burke and Ernzerhof in 1996 (PBE) [19], and the
modified version of PBE introduced by Perdew et al. in 2008 (PBEsol) [20]. For
convenience of the reader, we recall the expressions for these functionals
here:

$\displaystyle
g^{B88}(\rho,|\nabla\rho|)=\frac{2^{\frac{1}{3}}\beta\bigr{(}|\nabla\rho|/\rho^{\frac{4}{3}}\bigr{)}^{2}}{1+6\beta
2^{\frac{1}{3}}\bigr{(}|\nabla\rho|/\rho^{\frac{4}{3}}\bigr{)}\sinh^{-1}\bigr{(}2^{\frac{1}{3}}|\nabla\rho|/\rho^{\frac{4}{3}}\bigr{)}}\rho^{\frac{4}{3}}$
(7.1) $\displaystyle
g^{PBE}(\rho,|\nabla\rho|)=c_{x}\frac{\mu\bigr{(}|\nabla\rho|/\rho^{\frac{4}{3}}\bigr{)}^{2}}{4(3\pi^{2})^{\frac{2}{3}}+\frac{\mu}{\kappa}\bigr{(}|\nabla\rho|/\rho^{\frac{4}{3}}\bigr{)}^{2}}\rho^{\frac{4}{3}}$
(7.2)

where $\sinh^{-1}$ is the inverse hyperbolic sine and the constants are
$\beta=0.0042$, $\kappa=0.804$, and $\mu=0.2195$. For PBEsol, one has the same
expression as for PBE in (7.2), but with $\mu=0.1235$.

Figure 4: Exchange energy of the free electron gas in a box with zero boundary
conditions minus leading order term, per unit boundary area, for various
functionals. The number of electrons per unit volume was normalized to 1.
Solid lines: Numerical values. Dashed lines: asymptotic values (second order
coefficients from Theorem 1.1, present work).

Numerical evaluation of the exact one-dimensional integral expression for the
GGA constant in Theorem 1.1 gives the following values:

$c^{\textnormal{Dir}}_{PBE}\approx 0.0157,\quad
c^{\textnormal{Dir}}_{B88}\approx 0.0192,\quad
c^{\textnormal{Dir}}_{PBEsol}\approx 0.0105.$ (7.3)

To numerically verify the next-order asymptotics, we plotted in Figure 4 the
graph of the energy functionals minus the leading order term divided by the
boundary area $|\partial[0,L]^{3}|=6L^{2}$, together with the asymptotic
values predicted by Theorem 1.1 and (7.3). Precisely, since $\bar{\rho}=1$,
the values for the dashed lines in Figure 4 are, respectively,
$c^{\textnormal{Dir}}_{x,2}$, $c^{\textnormal{Dir}}_{\rm LDA}$,
$c^{\textnormal{Dir}}_{\rm LDA}+c^{\textnormal{Dir}}_{B88}$,
$c^{\textnormal{Dir}}_{\rm LDA}+c^{\textnormal{Dir}}_{PBE}$, $c_{\rm
LDA}^{\textnormal{Dir}}+c^{\textnormal{Dir}}_{PBEsol}$.

Overall, there is a good match between numerics and asymptotics as $N$ gets
large. More detailed observations are the following.

  * •

Asymptotically, the LDA underestimates the surface term by 12%, whereas B88
and PBE overestimate it by 13% respectively 8%.

  * •

Asymptotically, only PBEsol is much more accurate than the LDA, exhibiting an
error of just 1.4%. This should not come as a surprise to experts, as certain
surface data (although not the ones considered here) entered into the choice
of the parameters. Thus one may say that the present work provides an
alternative theoretical justification of the PBEsol parameters. The price to
pay is that PBEsol is the least accurate of the GGAs for very small $N$.

  * •

B88 is the most accurate GGA for very small $N$. This is not unexpected given
the fact that the parameter $\beta$ was fitted to data for the first few noble
gas atoms. The price to pay is that B88 does not improve on the LDA beyond a
few hundred electrons.

  * •

In the regime of 20 to 100 electrons, which is certainly relevant in
applications, particularly in chemistry, PBE fares best.

  * •

The slowest convergence to the asymptotic value, and the largest fluctuations,
occur for exact exchange. Neither asymptotics up to second order nor any of
the GGA functionals correctly reproduce these significant finite-$N$
fluctuations. Note that they would be captured exactly by the universal
Hohenberg-Kohn functional.

In the context of our model system, the free electron gas in a box with zero
boundary conditions, our rigorous asymptotic results and the above
observations illustrate both the advances that have been made in the physics
and chemistry literature in designing computationally simple exchange-
correlation functionals, and the immense difficulties in improving on the
current state of the art. For the latter, we hope that the new exact
constraint on GGAs presented here (eq. (1.14)) will in the future turn out to
be useful.

## Acknowledgments

We are grateful to Fernando Chamizo for a helpful email exchange regarding his
results. We would also like to thank Mi-Song Dupuy for useful comments on an
early draft of this manuscript and helpful discussions regarding the numerical
evaluation of exact exchange.

##  Appendix A Numerical scheme for exact exchange

Here we detail our scheme for accurate and efficient evaluation of exact
exchange for the free electron gas in a box. As already explained, the closed-
form expression (2.11), (2.9) cannot be evaluated directly by numerical
integration, because of the high-dimensionality of the domain of integration
(6D) and the Coulomb singularity of the integrand.

Recall that the eigenfunctions of the Laplacian with zero boundary conditions
on $Q=[0,1]^{3}$ are given by the following expression:

$\displaystyle\phi_{k}(x)=\prod_{i=1}^{3}\sqrt{2}\sin(\pi
k_{i}x_{i})\;\;\;(k\in\mathbb{N}^{3}).$ (A.1)

Hence for closed shell $N$, the ground state $\Psi_{N,L}$ of the free
$N$-electron gas in $Q_{L}=[0,L]^{3}$, $N/L^{3}=1$, is unique and the exact
exchange energy is, by rescaling to the fixed domain $Q$,

$\displaystyle E_{x}[\Psi_{N,L}]$ $\displaystyle=-\frac{1}{L}\int_{Q\times
Q}\frac{\left|\sum_{k\in\mathbb{N}^{3},\,|k|\leq
R_{N}}\phi_{k}(x)\overline{\phi_{k}(y)}\right|^{2}}{|x-y|}dx\,dy$
$\displaystyle=-\frac{1}{L}\sum_{\begin{subarray}{c}k,\ell\in\mathbb{N}^{3}\\\
|k|,|\ell|\leq R_{N}\end{subarray}}\underbrace{\int_{Q\times
Q}\frac{\prod\limits_{i=1}^{3}4\sin(\pi
k_{i}x_{i})\sin(\pi\ell_{i}x_{i})\sin(\pi
k_{i}y_{i})\sin(\pi\ell_{i}y_{i})}{|x-y|}dx\,dy}_{=:I_{k,\ell}}$ (A.2)

where $R_{N}$ is the Dirichlet Fermi radius defined in (2.6).

The starting point of our numerical scheme is a simple calculation which
reduces the above six-dimensional integral to the three-dimensional integral
of a separable function times the Coulomb potential over a finite region. For
periodic boundary conditions such a reduction is trivial because the system is
translation invariant, which implies that the exchange integrand depends only
on the relative coordinate $z=x-y$; but the zero boundary condition breaks the
translation invariance. Nevertheless the following holds:

######  Lemma A.1.

For $k$, $\ell\in\mathbb{N}^{3}$, and $I_{k,\ell}$ as defined above,

$\displaystyle
I_{k,\ell}=8\int_{[0,1]^{3}}\frac{\prod_{i=1}^{3}f_{k_{i},\ell_{i}}(z_{i})}{|z|}dz$
(A.3)

where for $a,b\in\mathbb{N}$, $f_{a,b}$ is defined as

$\displaystyle f_{a,b}(\tau)$
$\displaystyle=\frac{1-\tau}{2}\Bigl{(}\cos(\pi(a+b)\tau)+\cos(\pi(a-b)\tau)\Bigr{)}+\biggr{(}\frac{1}{a}+\frac{1}{b}-\frac{1}{(a+b)}\biggr{)}\frac{\sin(\pi(a+b)\tau)}{2\pi}$
$\displaystyle+\begin{cases}\Bigl{(}\frac{1}{a}-\frac{1}{b}-\frac{1}{a-b}\Bigr{)}\frac{\sin(\pi(a-b)\tau)}{2\pi}&\mbox{
if }a\neq b,\\\ \;\;\;\;\frac{1-\tau}{2}&\mbox{otherwise.}\end{cases}$ (A.4)

We remark that there is a well known alternative reduction of any 6D Coulomb
integral of the form $\int u(x)\frac{1}{|x-y|}v(y)\,dx\,dy$ to a 3D integral
over reciprocal space, by using the convolution theorem for the Fourier
transform. But this leads to an integral over an unbounded domain, a stronger
($\sim 1/|k|^{2}$) singularity, and – in our case – a slow decay of the
integrand, making the expression (A.3) numerically much more favourable.

###### Proof.

Using the identity $2\sin(A)\sin(B)=\cos(A-B)-\cos(A+B)$ and the change of
variables $w=x+y,\,z=x-y$, and noting that the volume element becomes
$\mathrm{d}x_{i}\mathrm{d}y_{i}=\frac{1}{2}\mathrm{d}w_{i}\mathrm{d}z_{i}$,

$I_{k,\ell}=\int_{[-1,1]^{3}}\frac{1}{|z|}\prod_{i=1}^{3}\underbrace{\left\\{\frac{1}{2}\int_{|z_{i}|}^{2-|z_{i}|}\bigl{(}\cos\pi
k_{i}z_{i}-\cos\pi
k_{i}w_{i}\bigr{)}\bigl{(}\cos\pi\ell_{i}z_{i}-\cos\pi\ell_{i}w_{i}\bigr{)}dw_{i}\right\\}}_{=:f_{k_{i},\ell_{i}}(z_{i})}dz.$

The integral in the definition of $f_{a,b}(z_{i})$ is elementary to evaluate
by using the identity $2\cos A\cos B=cos(A+B)+\cos(A-B)$ for the term $\cos\pi
aw_{i}\cos\pi bw_{i}$, yielding

$\displaystyle f_{a,b}(z_{i})$ $\displaystyle=(1-|z_{i}|)\cos(\pi
az_{i})\cos(\pi bz_{i})+\cos(\pi az_{i})\frac{\sin(\pi b|z_{i}|)}{\pi
b}+\cos(\pi bz_{i})\frac{\sin(\pi a|z_{i}|)}{\pi a}$
$\displaystyle-\frac{\sin(\pi(a+b)|z_{i}|)}{2\pi(a+b)}\;\;\begin{cases}-\frac{\sin(\pi(a-b)|z_{i}|)}{2\pi(a-b)}&\mbox{if
}a\neq b,\\\ +\frac{1-|z_{i}|}{2}&\mbox{otherwise.}\end{cases}$

Next, we note that since $f_{a,b}(-z_{i})=f_{a,b}(z_{i})$ and the Coulomb
potential $1/|z|$ is invariant under the transformations $z_{i}\mapsto-z_{i}$,
we can replace $[-1,1]^{3}$ by $[0,1]^{3}$ by adding a factor of $8$ in front
of the integral. The final expression for $f_{a,b}$ given in the lemma now
follows from the trigonometric identities $2\cos A\cos B=\cos(A+B)+\cos(A-B)$
and $\alpha\sin A\cos B+\beta\cos A\sin
B=\frac{\alpha+\beta}{2}\sin(A+B)+\frac{\alpha-\beta}{2}\sin(A-B)$. ∎

Since the integrand in (A.3) is a separable function except for the Coulomb
potential, the idea now is to also approximate the latter by separable
functions, therefore reducing the problem to the computation of one
dimensional integrals. Such an approximation is provided by recent advances in
low-rank tensor approximation; more specifically, we use results of Hackbusch
([12, Section 9.8.2]). The Coulomb potential can be very accurately
approximated by a sum of weighted Gaussians,

$\displaystyle\frac{1}{r}\approx\sum_{j=1}^{M}\omega_{j}e^{-\alpha_{j}r^{2}}.$
(A.5)

Plugging this approximation into equation (A.3) and factorizing
$e^{-\alpha_{j}|z|^{2}}=\prod_{i=1}^{3}e^{-\alpha_{j}z_{j}^{2}}$, one obtains

$\displaystyle I_{k,\ell}\approx
8\sum_{j=1}^{M}\omega_{j}\prod_{i=1}^{3}I_{k_{i},\ell_{i},j}\;\mbox{ with
}\;I_{a,b,j}=\int_{0}^{1}f_{a,b}(t)e^{-\alpha_{j}t^{2}}dt,$ (A.6)

which reduces the 3D integral in (A.3) to one-dimensional integrals of
analytic functions.101010In fact, one could represent these integrals exactly
in terms of the error function ${\tt erf}$ and the imaginary error function
${\tt erfi}$, but we do not use this fact on our scheme.

To reduce the overall number of 1D integrals that must be computed, let us
introduce, for $p\in\\{0,1,...,2\lfloor{R_{N}^{\textnormal{Dir}}}\rfloor\\}$
and $j\in\\{1,...,M\\}$ (where $\lfloor\;\rfloor$ denotes the integer part),
the auxiliary integrals

$\displaystyle C_{p,j}\coloneqq\int_{0}^{1}\frac{1-t}{2}\cos(\pi
p\,t)e^{-\alpha_{j}t^{2}}\mathrm{d}t\quad\mbox{ and }\quad
S_{p,j}\coloneqq\int_{0}^{1}\frac{\sin(\pi
p\,t)}{2\pi}e^{-\alpha_{j}t^{2}}\mathrm{d}t.$ (A.7)

It follows from the explicit expression for $f_{a,b}$ in (A.4) that

$\displaystyle I_{a,b,j}$
$\displaystyle=C_{a+b,j}+C_{|a-b|,j}+\Bigl{(}\tfrac{1}{a}+\tfrac{1}{b}-\tfrac{1}{a+b}\Bigr{)}S_{a+b,j}$
$\displaystyle\;\;\;+\begin{cases}\Bigl{(}\tfrac{1}{a}-\tfrac{1}{b}-\tfrac{1}{a-b}\Bigr{)}{\rm
sign}(a-b)S_{|a-b|,j}&\mbox{ if }a\neq b\\\ \;C_{0,j}&\mbox{ if
}a=b.\end{cases}$ (A.8)

Thus in total

$E_{x}[\Psi_{N,L}]\approx-\frac{8}{L}\sum_{\begin{subarray}{c}k,\ell\in\mathbb{N}^{3}\\\
|k|,|\ell|\leq
R_{N}\end{subarray}}\sum_{j=1}^{M}\omega_{j}\prod_{i=1}^{3}I_{k_{i},\ell_{i},j}$
(A.9)

with $I_{a,b,j}$ given by (A.7)–(A.8). In particular, as
$R_{N}^{\textnormal{Dir}}\sim N^{1/3}$, calculating the exchange energy of the
free $N$-electron gas reduces to the problem of evaluating
$\mathcal{O}(N^{1/3}M)$ one-dimensional integrals of analytic functions on the
interval $[0,1]$, and multiplying and summing them according to equation
(A.9).

Next, let us discuss the choice of weights and exponents, and the error, in
(A.5). We used the values $\\{\omega^{H}_{j},\alpha^{H}_{j}\\}_{j=1}^{51}$
given on Hackbusch’s webpage [11] for the (approximately) best approximation
of the Coulomb potential as the sum of $M=51$ Gaussians, which satisfy $\lVert
1/r-v^{H}(r)\rVert_{L^{\infty}([1,10^{9}])}\leq 10^{-9}$ where
$v^{H}(r)=\sum_{j=1}^{51}\omega_{j}^{H}e^{-\alpha_{j}^{H}r^{2}}$. Moreover
since we are interested in a good approximation of $1/r$ on the unit cube, we
rescaled Hackbusch’s parameters by setting

$\displaystyle\omega_{j}=\frac{\omega^{H}_{j}}{r_{0}}\quad\mbox{ and
}\quad\alpha_{j}=\frac{\alpha_{j}^{H}}{r_{0}^{2}},$

which yields a pointwise error of

$\displaystyle\lVert v(r)-1/r\rVert_{L^{\infty}([r_{0},10^{9}r_{0}])}\leq
r_{0}^{-1}\times 10^{-9}$ (A.10)

where $v(r)=\sum_{j=1}^{51}\omega_{j}e^{-\alpha_{j}r^{2}}$. In our numerical
results we chose $r_{0}=10^{-4}$, to achieve good accuracy in (A.3) both in
the region $|z|<r_{0}$ (note that the integral of $1/|z|$ over this region is
$\sim r_{0}^{2}$) and outside it.

Finally, let us discuss evaluation of the 1D integrals (A.7), which requires a
moment’s thought as one needs to resolve both the oscillatory trigonometric
factor and the Gaussian factor. The wavevector $\pi p$ is $\leq\pi\cdot
2R_{N}^{\textnormal{Dir}}$ and hence $\lesssim 200$ for up to $N=30~{}000$
electrons (in which case $R_{N}^{\textnormal{Dir}}\approx 31$), so the
trigonometric oscillations can be accurately resolved by any standard
quadrature method. The Gaussian factor, however, turns out to be more
delicate, as the maximum value of $\alpha_{j}$ is $\approx 8\times 10^{8}$.
For $\alpha_{j}>1$ we therefore used the following alternative expressions
obtained by re-scaling:

$\displaystyle C_{p,j}$
$\displaystyle=\frac{1}{\sqrt{\alpha_{j}}}\int_{0}^{\sqrt{\alpha_{j}}}\biggr{(}\frac{1}{2}-\frac{t}{2\sqrt{\alpha_{j}}}\biggr{)}\cos(\pi
p\frac{t}{\sqrt{\alpha_{j}}})e^{-t^{2}}\mathrm{d}t,$ (A.11) $\displaystyle
S_{p,j}$
$\displaystyle=\frac{1}{\sqrt{\alpha_{j}}}\int_{0}^{\sqrt{\alpha_{j}}}\frac{1}{2\pi}\sin(\pi
p\frac{t}{\sqrt{\alpha_{j}}})e^{-t^{2}}\mathrm{d}t.$ (A.12)

Note that even though the integration interval may be big, for practical
purposes one can truncate at $\min\\{\sqrt{\alpha_{j}},10\\}$ (as
$\int_{10}^{\infty}e^{-t^{2}}\mathrm{d}t\approx 10^{-45}$).

##  Appendix B Assumptions on GGAs

We now show that the expressions for the PBE and B88 functionals (see
equations (7.2) and (7.1)) satisfy our assumptions in Theorem 1.1.

The $C^{1}$ regularity in $(0,\infty)\times\mathbb{R}$ is straightforward, so
we just need to worry about the continuity when $a$ goes to zero and $b$
remains bounded. For this, let us rewrite equations (7.1) and (7.2) as

$\displaystyle g^{PBE}(a,b)=c_{x}\frac{\mu
s^{2}}{4(3\pi^{2})^{\frac{2}{3}}+\frac{\mu}{\kappa}s^{2}}a^{\frac{4}{3}}$
$\displaystyle g^{B88}(a,b)=\frac{2^{\frac{1}{3}}s^{2}}{1+6\beta
2^{\frac{1}{3}}s\sinh^{-1}(2^{\frac{1}{3}}s)}a^{\frac{4}{3}}=\frac{2^{\frac{1}{3}}s}{1+6\beta
2^{\frac{1}{3}}s\sinh^{-1}(2^{\frac{1}{3}}s)}b$

where $s=b/a^{\frac{4}{3}}$. Thus for the PBE, since the (enhancement) factor
in front of $a^{\frac{4}{3}}$ is bounded, we see that $g^{PBE}(a,b)\rightarrow
0$ as $a\rightarrow 0$ regardless of $b$. For the B88, we make two
observations: (i) thanks to the superlinear growth of the denominator in the
enhancement factor, we see that $g^{B88}(a,b)$ goes to zero if
$s\rightarrow\infty$ and $b$ stays bounded, and (ii) if $s$ is bounded and
$b\rightarrow 0$, $g^{B88}(a,b)$ also goes to zero. In particular, as taking
the limit $a\rightarrow 0$ with $b$ bounded falls into one of these two cases,
property (1.12) holds.

## References

  * [1] A. D. Becke.  Density-functional exchange-energy approximation with correct asymptotic behavior.  Physiscal Review A, 38:3098–3100, Sep 1988.
  * [2] F. Bloch.  Bemerkung zur Elektronentheorie des Ferromagnetismus und der elektrischen Leitfähigkeit.  Zeitschrift für Physik, 57(7):545–555, Jul 1929.
  * [3] F. Chamizo and H. Iwaniec.  On the sphere problem.  Revista Matemática Iberoamericana, 11(2):417–429, 1995.
  * [4] J.-R. Chen.  Improvement on the asymptotic formulas for the number of lattice points in a region of three dimensions (ii).  Scientia Sinica, 12(5):739–741, 1963.
  * [5] M. R. Christiansen, C. Hainzl, and P. T. Nam.  The Gell-Mann–Brueckner formula for the correlation energy of the electron gas: A rigorous upper bound in the mean-field regime.  Communications in Mathematical Physics, Feb 2023.
  * [6] P. A. M. Dirac.  Note on exchange phenomena in the Thomas atom.  Mathematical Proceedings of the Cambridge Philosophical Society, 26(3):376–385, 1930.
  * [7] G. Friesecke.  Pair correlations and exchange phenomena in the free electron gas.  Communications in Mathematical Physics, 184(1):143–171, Mar 1997\.
  * [8] G. Friesecke, A. Gerolin, and P. Gori-Giorgi.  The strong-interaction limit of density functional theory.  arXiv preprint arXiv:2202.09760, 2022.
  * [9] G. M. Graf and J. P. Solovej.  A correlation estimate with applications to quantum systems with Coulomb interactions.  Reviews in Mathematical Physics, 6(05a):977–997, 1994.
  * [10] S. W. Graham and G. Kolesnik.  Van der Corput’s Method of Exponential Sums.  London Mathematical Society Lecture Note Series. Cambridge University Press, 1991.
  * [11] W. Hackbusch.  https://www.mis.mpg.de/scicomp/EXP_SUM/1_sqrtx/.  Accessed: 22-02-2022.
  * [12] W. Hackbusch.  Tensor spaces and numerical tensor calculus, volume 42.  Springer, 2012.
  * [13] D. R. Heath-Brown.  Lattice points in the sphere.  In In Number theory in progress, pages 883–892, 1999.
  * [14] D. R. Heath-Brown.  A new kth derivative estimate for exponential sums via Vinogradov’s mean value.  Proceedings of the Steklov Institute of Mathematics, 296(1):88–103, 2017.
  * [15] M. N. Huxley.  Area, lattice points, and exponential sums, volume 13.  Clarendon Press, 1996.
  * [16] V. Ivrii.  100 years of Weyl’s law.  Bulletin of Mathematical Sciences, 6(3):379–452, 2016.
  * [17] W. Kohn and L. J. Sham.  Self-consistent equations including exchange and correlation effects.  Physical Review, 140:A1133–A1138, Nov 1965.
  * [18] P.-F. Loos and P. M. W. Gill.  The uniform electron gas.  WIREs Computational Molecular Science, 6(4):410–429, 2016.
  * [19] J. P. Perdew, K. Burke, and M. Ernzerhof.  Generalized gradient approximation made simple.  Physical Review Letters, 77:3865–3868, Oct 1996.
  * [20] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke.  Restoring the density-gradient expansion for exchange in solids and surfaces.  Physical Review Letters, 100:136406, Apr 2008.
  * [21] M. Reed and B. Simon.  Methods of Modern Mathematical Physics. IV Analysis of Operators.  Academic Press, 1978.
  * [22] C. D. Sogge.  Fourier integrals in classical analysis, volume 210.  Cambridge University Press, 2017.
  * [23] I. M. Vinogradov.  On the number of integer points in a sphere.  Izvestiya Akademii Nauk SSSR. Seriya Matematicheskaya, 27:957–968, 1963.