The Local Density Approximation in Density Functional Theory

TL;DR

This paper provides a rigorous mathematical justification for the Local Density Approximation in Density Functional Theory, quantifying its accuracy and conditions under which it is valid.

Contribution

It offers the first rigorous proof and quantitative estimates for the Local Density Approximation's validity in DFT.

Findings

Quantitative error bounds involving gradient terms

Justification of LDA for flat electron densities

Conditions for the approximation's accuracy

Abstract

We give the first mathematically rigorous justification of the Local Density Approximation in Density Functional Theory. We provide a quantitative estimate on the difference between the grand-canonical Levy-Lieb energy of a given density (the lowest possible energy of all quantum states having this density) and the integral over the Uniform Electron Gas energy of this density. The error involves gradient terms and justifies the use of the Local Density Approximation in the situation where the density is very flat on sufficiently large regions in space.

Equations

$$F_{\rm LL}(\rho)\approx\frac{1}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{\rho(x)\rho(y)}{|x-y|}\,dx\,dy+\int_{\mathbb{R}^{3}}e_{\rm UEG}\big{(}\rho(x)\big{)}\,dx.$$

$$\left|F_{\rm LL}(\rho)-\frac{1}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{\rho(x)\rho(y)}{|x-y|}\,dx\,dy-\int_{\mathbb{R}^{3}}e_{\rm UEG}\big{(}\rho(x)\big{)}\,dx\right|\\ \leqslant\varepsilon\int_{\mathbb{R}^{3}}\big{(}\rho(x)+\rho(x)^{2}\big{)}\,dx+\frac{C(1+\varepsilon)}{\varepsilon}\int_{\mathbb{R}^{3}}|\nabla\sqrt{\rho}(x)|^{2}\,dx\\ +\frac{C}{\varepsilon^{4p-1}}\int_{\mathbb{R}^{3}}|\nabla\rho^{\theta}(x)|^{p}\,dx$$

$$\begin{cases}\displaystyle\int_{\mathbb{R}^{3}}|\nabla\sqrt{\rho}(x)|^{2}\,dx\ll\int_{\mathbb{R}^{3}}\big{(}\rho(x)+\rho(x)^{2}\big{)}\,dx,\\[11.38092pt] \displaystyle\int_{\mathbb{R}^{3}}|\nabla\rho^{\theta}(x)|^{p}\,dx\ll\int_{\mathbb{R}^{3}}\big{(}\rho(x)+\rho(x)^{2}\big{)}\,dx.\end{cases}$$

$$\rho_{N}(x):=\rho\big{(}N^{-1/3}x\big{)}$$

$$\left|F_{\rm LL}(\rho_{N})-\frac{N^{\frac{5}{3}}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{\rho(x)\rho(y)}{|x-y|}\,dx\,dy-N\int_{\mathbb{R}^{3}}e_{\rm UEG}\big{(}\rho(x)\big{)}\,dx\right|\leqslant CN^{\frac{11}{12}}.$$

$$F_{\rm LL}(\rho):=\\ \inf_{\begin{subarray}{c}\Gamma_{n}=\Gamma_{n}^{*}\geqslant 0\\ \sum_{n=0}^{\infty}{\rm Tr}\,(\Gamma_{n})=1\\ \sum_{n=1}^{\infty}\rho_{\Gamma_{n}}=\rho\end{subarray}}\Bigg{\{}\sum_{n=1}^{\infty}{\rm Tr}\,_{\mathfrak{H}^{n}}\Bigg{(}-\sum_{j=1}^{n}\Delta_{x_{j}}+\sum_{1\leqslant j<k\leqslant n}\frac{1}{|x_{j}-x_{k}|}\Bigg{)}\Gamma_{n}\Bigg{\}}.$$

$${\mathfrak{H}}^n:=L_a^2(({\mathbb{R}}^3\times \{1,...,q\}{)}^n,\mathbb{C})$$

$$\begin{gathered} {\rho }_{{\Gamma }_n}(x)= \\ n\sum _{\begin{array}{c}{\sigma }_1,...,{\sigma }_n\\ \in \{1,...,q\}\end{array}}{\int }_{{\mathbb{R}}^{3(n-1)}}{\Gamma }_n(x,{\sigma }_1,x_2,...,x_n,{\sigma }_n;x,{\sigma }_1,x_2,...,x_n,{\sigma }_n)d{x}_2\cdots d{x}_n \end{gathered}$$

$${\Gamma }_n(x_{\tau (1)},{\sigma }_{\tau (1)},...,x_{\tau (N)},{\sigma }_{\tau (N)};x_1^{\prime },{\sigma }_1^{\prime },...,x_N^{\prime },{\sigma }_N^{\prime })$$

$$={\Gamma }_n(x_1,{\sigma }_1,...,x_N,{\sigma }_N;x_{\tau (1)}^{\prime },{\sigma }_{\tau (1)}^{\prime }...,x_{\tau (N)}^{\prime },{\sigma }_{\tau (N)}^{\prime })$$

$$=\epsilon (\tau ){\Gamma }_n(x_1,{\sigma }_1,...,x_N,{\sigma }_N;x_1^{\prime },{\sigma }_1^{\prime }...,x_N^{\prime },{\sigma }_N^{\prime })$$

$$\boxed{E(\rho ):=F_{\mathrm{LL}}(\rho )-\frac{1}{2}{\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}\frac{\rho (x)\rho (y)}{|x-y|}dxdy.}$$

$$T(\rho):=\inf_{\begin{subarray}{c}\Gamma_{n}=\Gamma_{n}^{*}\geqslant 0\\ \sum_{n=0}^{\infty}{\rm Tr}\,(\Gamma_{n})=1\\ \sum_{n=1}^{\infty}\rho_{\Gamma_{n}}=\rho\end{subarray}}\Bigg{\{}\sum_{n=1}^{\infty}{\rm Tr}\,_{\mathfrak{H}^{n}}\Bigg{(}-\sum_{j=1}^{n}\Delta_{x_{j}}\Bigg{)}\Gamma_{n}\Bigg{\}}=\inf_{\begin{subarray}{c}0\leqslant\gamma=\gamma^{*}\leqslant 1\\ \rho_{\gamma}=\rho\end{subarray}}{\rm Tr}\,(-\Delta)\gamma$$

$$|\partial {\Omega }_N+B_r|⩽Cr|{\Omega }_N{|}^{2/3},\text{for all r⩽∣ΩN∣1/3/C,}$$

$$\boxed{\lim_{N\to\infty}\frac{E\big{(}\rho_{0}{\mathds{1}}_{\Omega_{N}}\ast\chi_{\delta_{N}}\big{)}}{|\Omega_{N}|}=e_{\rm UEG}\big{(}\rho_{0}\big{)}}$$

$$2⩽p\theta ⩽1+\frac{p}{2}.$$

$$\left|E(\rho)-\int_{\mathbb{R}^{3}}e_{\rm UEG}\big{(}\rho(x)\big{)}\,dx\right|\leqslant\varepsilon\int_{\mathbb{R}^{3}}\big{(}\rho(x)+\rho(x)^{2}\big{)}\,dx\\ +\frac{C(1+\varepsilon)}{\varepsilon}\int_{\mathbb{R}^{3}}|\nabla\sqrt{\rho}(x)|^{2}\,dx+\frac{C}{\varepsilon^{4p-1}}\int_{\mathbb{R}^{3}}|\nabla\rho^{\theta}(x)|^{p}\,dx$$

$${\rho }_N(x)=\rho (N^{-1/3}x).$$

$$\left|E(\rho_{N})-N\int_{\mathbb{R}^{3}}e_{\rm UEG}\big{(}\rho(x)\big{)}\,dx\right|\\ \leqslant CN^{\frac{5}{12}}\int_{\mathbb{R}^{3}}|\nabla\sqrt{\rho}(x)|^{2}\,dx+CN^{\frac{11}{12}}\int_{\mathbb{R}^{3}}\left(\rho(x)+\rho(x)^{2}+|\nabla\rho^{\theta}(x)|^{p}\right)\,dx.$$

$${\int }_{{\mathbb{R}}^3}{\left|\nabla \sqrt{{\rho }_N}(x)\right|}^{2}dx=N^{\frac{1}{3}}{\int }_{{\mathbb{R}}^3}{\left|\nabla \sqrt{\rho }(x)\right|}^{2}dx$$

$${\int }_{{\mathbb{R}}^3}{\left|\nabla {\rho }_N^{1/3}(x)\right|}^{2}dx=N^{\frac{1}{3}}{\int }_{{\mathbb{R}}^3}{\left|\nabla {\rho }^{1/3}(x)\right|}^{2}dx.$$

$$E_{\mathrm{xc}}(\rho )=E(\rho )-T(\rho ),$$

$$\left|T(\rho )-q^{-\frac{2}{3}}{c}_{\mathrm{TF}}{\int }_{{\mathbb{R}}^3}\rho (x{)}^{\frac{5}{3}}dx\right|⩽\epsilon q^{-\frac{2}{3}}{\int }_{{\mathbb{R}}^3}\rho (x{)}^{\frac{5}{3}}dx+\frac{C}{{\epsilon }^{\frac{13}{3}}}{\int }_{{\mathbb{R}}^3}|\nabla \sqrt{\rho }(x){|}^{2}dx$$

$$\left|E_{\rm xc}(\rho)-\int_{\mathbb{R}^{3}}e_{\rm UEG}\big{(}\rho(x)\big{)}\,dx+q^{-\frac{2}{3}}c_{\rm TF}\int_{\mathbb{R}^{3}}\rho(x)^{\frac{5}{3}}\,dx\right|\\ \leqslant\varepsilon\int_{\mathbb{R}^{3}}\big{(}\rho(x)+\rho(x)^{2}\big{)}\,dx+\frac{C(1+\varepsilon)}{\varepsilon}\int_{\mathbb{R}^{3}}|\nabla\sqrt{\rho}(x)|^{2}\,dx\\ +\frac{C}{\varepsilon^{4p-1}}\int_{\mathbb{R}^{3}}|\nabla\rho^{\theta}(x)|^{p}\,dx.$$

$$\begin{gathered} {\rho }_{\sigma }(x)= \\ \sum _{n⩾1}n\sum _{\begin{array}{c}{\sigma }_2,...,{\sigma }_n\\ \in \{1,...,q\}\end{array}}{\int }_{{\mathbb{R}}^{3(n-1)}}{\Gamma }_n(x,\sigma ,x_2,...,x_n,{\sigma }_n;x,\sigma ,x_2,...,x_n,{\sigma }_n)d{x}_2\cdots d{x}_n \end{gathered}$$

$$\left|E(\rho_{1},...,\rho_{q})-\int_{\mathbb{R}^{3}}e_{\rm UEG}\big{(}\rho_{1}(x),...,\rho_{q}(x)\big{)}\,dx\right|\leqslant\varepsilon\int_{\mathbb{R}^{3}}\big{(}\rho(x)+\rho(x)^{2}\big{)}\,dx\\ +\frac{C(1+\varepsilon)}{\varepsilon}\int_{\mathbb{R}^{3}}|\nabla\sqrt{\rho}(x)|^{2}\,dx+\frac{C}{\varepsilon^{4p-1}}\int_{\mathbb{R}^{3}}|\nabla\rho^{\theta}(x)|^{p}\,dx.$$

$$\mathrm{Tr}(-\Delta )\gamma ⩾q^{-\frac{2}{d}}{c}_{\mathrm{LT}}{\int }_{{\mathbb{R}}^d}{\rho }_{\gamma }(x{)}^{1+\frac{2}{d}}dx$$

$$c_{\mathrm{TF}}=\frac{4{\pi }^{2}d}{(d+2)}{\left(\frac{d}{|{\mathbb{S}}^{d-1}|}\right)}^{\frac{2}{d}}$$

$$\mathrm{Tr}(-\Delta )\gamma ⩾q^{-\frac{2}{d}}{c}_{\mathrm{TF}}\left(1-\epsilon \right){\int }_{{\mathbb{R}}^d}{\rho }_{\gamma }(x{)}^{1+\frac{2}{d}}dx-\frac{\kappa }{{\epsilon }^{3+\frac{4}{d}}}{\int }_{{\mathbb{R}}^d}|\nabla \sqrt{{\rho }_{\gamma }}(x){|}^{2}dx$$

$$\mathrm{Tr}(-\Delta )\gamma ⩾{\int }_{{\mathbb{R}}^d}{\left|\nabla \sqrt{{\rho }_{\gamma }}(x)\right|}^{2}dx$$

$$\sum_{n=1}^{\infty}{\rm Tr}\,_{\mathfrak{H}^{n}}\Bigg{(}\sum_{1\leqslant j<k\leqslant n}\frac{1}{|x_{j}-x_{k}|}\Bigg{)}\Gamma_{n}\Bigg{\}}\\ \geqslant\frac{1}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{\rho(x)\rho(y)}{|x-y|}dx\,dy-1.64\int_{\mathbb{R}^{3}}\rho(x)^{\frac{4}{3}}\,dx$$

Full text

© 2019 by the authors. This paper may be reproduced, in its entirety, for non-commercial purposes. Final version to appear in Pure and Applied Analysis.

1. 1 Introduction

2. 2 Main result

3. 2.1 The grand-canonical Levy-Lieb functional

4. 2.2 The Uniform Electron Gas

5. 2.3 The Local Density Approximation

6. 3 A priori estimates on $T(\rho)$ and $E(\rho)$

7. 3.1 Known lower bounds

8. 3.2 Upper bound on the best kinetic energy

9. 3.3 Upper bound on $E(\rho)$

10. 4 Proof of Theorems 1 and 2

11. 4.1 Upper bound in terms of local densities

12. 4.2 Lower bound in terms of local densities

13. 4.3 A convergence rate for tetrahedra

14. 4.4 Proof of Theorem 1

15. 4.5 Replacing the local density by a constant density

16. 4.6 Lipschitz regularity of $e_{\rm UEG}$

17. 4.7 Proof of Theorem 2

18. A Classical case

1. Introduction

Density Functional Theory (DFT) [10, 44, 11, 4, 47] is the most efficient approximation of the many-body Schrödinger equation for electrons. It is used in several areas of physics and chemistry and its success in predicting the electronic properties of atoms, molecules and materials is unprecedented. Among the many functionals that have been developed over the years [42], the Local Density Approximation (LDA) is the standard and simplest scheme [18, 22, 10, 44, 45]. It is not as accurate as its successors involving gradient corrections, but it is considered as “the mother of all approximations” [46] and it is still one of the methods of choice in solid state physics.

In the orbital-free formulation of Density Functional Theory [25, 35], the Local Density Approximation consists in replacing the full ground state energy by a local functional as follows:

[TABLE]

Here $\rho$ is the given one-particle density of the system and $F_{\rm LL}(\rho)$ is the Levy-Lieb functional [25, 35], the main object of interest in DFT. This is the lowest possible Schrödinger energy of all quantum states having the prescribed density $\rho$. The first term on the right side is called the direct or Hartree term. It is the classical electrostatic interaction energy of the density $\rho$ and it is the only nonlocal term in the LDA. The second term is the energy of the Uniform Electron Gas (UEG) [10, 44, 13, 30], containing all of the kinetic energy and the exchange-correlation energy in our convention. That is, $e_{\rm UEG}(\rho_{0})$ is the ground state energy per unit volume of the infinite electron gas with the prescribed constant density $\rho_{0}$ over the whole space (from which the direct term has been dropped). The rationale for the approximation (1) is to assume that the density is almost constant locally (in little boxes of volume $dx$), and to replace the local energy per unit volume by that of the infinite gas at that density $\rho(x)$.111It is often more convenient to fix the densities $\rho^{\uparrow}(x)$ and $\rho^{\downarrow}(x)$ of, respectively, spin-up and spin-down electrons instead of the total density $\rho(x)=\rho^{\uparrow}(x)+\rho^{\downarrow}(x)$. All our results apply similarly to this situation, as explained below in Remark 4.

Our goal in this paper is to justify the approximation (1) in the appropriate regime where $\rho$ is flat in sufficiently large regions of $\mathbb{R}^{3}$. We will prove the following quantitative estimate

[TABLE]

for all $\varepsilon>0$, where $F_{\rm LL}(\rho)$ is the grand canonical version of the Levy-Lieb functional. The parameters $p>3$ and $0<\theta<1$ should satisfy some conditions which will be explained below. For instance, $p=4$ and $\theta=1/2$ is allowed. After optimizing over $\varepsilon$, this justifies the LDA when the two gradient terms are much smaller than the local term

[TABLE]

For instance for a rescaled density in the form

[TABLE]

with $\int_{\mathbb{R}^{3}}\rho=1$, we obtain after taking $\varepsilon=N^{-1/12}$

[TABLE]

The bound (2) is, to our knowledge, the first estimate of this kind on the fundamental functional $F_{\rm LL}$. Although it should be possible to extract a definite value of the constant $C$ from our proof, it is probably very large and we have not tried to do it. The factor $1/\varepsilon^{4p-1}$ is also quite large and it is an open problem to improve it. We hope that our work will stimulate more results on the functional $F_{\rm LL}$ in the regime of slowly varying densities.

In physics and chemistry, the exchange-correlation energy is defined by subtracting a kinetic energy term $T(\rho)$ from $F_{\rm LL}(\rho)$. In this paper we also derive a bound on $T(\rho)$ which, when combined with (2), provides a bound on the exchange-correlation energy similar to (2). This is explained below in Remark 3.

In the next section we provide the precise mathematical definition of $F_{\rm LL}$ and $e_{\rm UEG}$, and we state our main theorem containing the estimate (2). In Section 3 we review some known a priori estimates on $F_{\rm LL}$ and prove a new upper bound on the kinetic energy. Section 4 contains the proof of our main results. Finally, in Appendix A we discuss a similar bound in the classical case where the kinetic energy is dropped, extending thereby our previous result in [30].

Acknowledgments.

The authors thank the Institut Henri Poincaré for its hospitality. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreements AQUAMS No 694227 of R.S. and MDFT No 725528 of M.L.).

2. Main result

2.1. The grand-canonical Levy-Lieb functional

Let us consider a density $\rho\in L^{1}(\mathbb{R}^{3},\mathbb{R}_{+})$ such that $\sqrt{\rho}\in H^{1}(\mathbb{R}^{3})$. Naturally we should assume in addition that $\int_{\mathbb{R}^{3}}\rho=N$ is an integer, but here we will work in the grand canonical ensemble where this is not needed. The grand-canonical Levy-Lieb functional [25, 35, 30] is defined by

[TABLE]

Here

[TABLE]

is the $n$-particle space of antisymmetric square-integrable functions on $(\mathbb{R}^{3}\times\{1,...,q\})^{n}$, with $q$ spin states (for electrons $q=2$). The family of operators $\Gamma=\{\Gamma_{n}\}_{n\geqslant 0}$ forms a grand-canonical mixed quantum state, that is, a state over the fermionic Fock space (commuting with the particle number operator). The density of each $\Gamma_{n}$ is defined by

[TABLE]

where $\Gamma_{n}(x_{1},...,\sigma_{n};x_{1}^{\prime},...,\sigma^{\prime}_{n})$ is the kernel of the trace-class operator $\Gamma_{n}$. This kernel is such that

[TABLE]

for every permutation $\tau\in\mathfrak{S}^{N}$ with signature $\varepsilon(\tau)\in\{\pm 1\}$.

If $\int_{\mathbb{R}^{3}}\rho=N\in\mathbb{N}$ and we restrict ourselves to mixed states $\Gamma$ where only $\Gamma_{N}$ is non-zero, we obtain Lieb’s functional [35]. If we further assume that $\Gamma_{N}=|\Psi\rangle\langle\Psi|$ is a rank-one projection, then we find the original Levy or Hohenberg-Kohn functional [18, 25]. It is well known [35] that working with mixed states has several advantages, in particular we obtain a convex function of $\rho$.

The grand-canonical version (3) is less popular but still important physically.222The grand canonical functional $F_{\rm LL}$ is the weak-$\ast$ lower semi-continuous closure of the canonical Lieb functional [27]. Hence it appears naturally in situations where some particles can be lost, e.g. in scattering processes. It is also a convex function of $\rho$. The fact that we can appeal to states with an arbitrary number $n$ of particles (but still a fixed average number $\int_{\mathbb{R}^{3}}\rho$) will considerably simplify several technical parts of our study. We expect that our main result (Theorem 2 below) holds the same for the canonical functionals, for which the energy is minimized over mixed or pure states with $N$ particles.

It is useful to subtract the direct term from $F_{\rm LL}$, hence to consider the energy

[TABLE]

The (grand-canonical) exchange-correlation energy is defined by $E_{\rm xc}(\rho):=E(\rho)-T(\rho)$ where

[TABLE]

is the lowest possible kinetic energy. We will study the functional $T(\rho)$ in Section 3.2 below.

2.2. The Uniform Electron Gas

In [30, Section 5], we have defined the uniform electron gas, which is obtained in the limit when $\rho$ approaches a constant function in the whole space. This is believed to be the same as the ground state energy of Jellium, where the density is not necessarily constant but the electrons instead evolve in a constant background [36]. This has recently been proved in the classical case in [31, 7] and the same is expected in the quantum case.

The following result is a slight improvement of [30, Thm. 5.1].

Theorem 1** (Quantum Uniform Electron Gas).**

Let $\rho_{0}>0$. Let $\{\Omega_{N}\}\subset\mathbb{R}^{3}$ be a sequence of bounded connected domains with $|\Omega_{N}|\to\infty$, such that $\Omega_{N}$ has a uniformly regular boundary in the sense that

[TABLE]

for some constant $C>0$. Let $\delta_{N}>0$ be any sequence such that $\delta_{N}/|\Omega_{N}|^{1/3}\to 0$ and $\delta_{N}|\Omega_{N}|^{1/3}\to\infty$. Let $\chi\in L^{1}(\mathbb{R}^{3})$ be a radial non-negative function of compact support such that $\int_{\mathbb{R}^{3}}\chi=1$ and $\int_{\mathbb{R}^{3}}|\nabla\sqrt{\chi}|^{2}<\infty$. Denote $\chi_{\delta}(x)=\delta^{-3}\chi(x/\delta)$. Then the following thermodynamic limit exists

[TABLE]

where the function $e_{\rm UEG}$ is independent of the sequence $\{\Omega_{N}\}$, of $\delta_{N}$ and of $\chi$.

For more properties of the UEG energy $e_{\rm UEG}$ we refer to  [30] and the references therein. In [30, Thm. 5.1] we rather optimized over the values of $\rho$ in the transition region around $\partial\Omega_{N}$. We were able to prove the simple limit (6) only when $\Omega_{N}$ is a tetrahedron. Using an upper bound on $E(\rho)$ that will be derived later in Proposition 1, we are now able to treat more reasonable limits in the form of (6). The proof of Theorem 1 is provided in Section 4.4 below.

The function $\chi_{\delta_{N}}$ is used to regularize the function $\rho_{0}{\mathds{1}}_{\Omega_{N}}$ which cannot be the density of a quantum state, since its square root is not in $H^{1}(\mathbb{R}^{3})$ [35]. The first condition $\delta_{N}/|\Omega_{N}|^{1/3}\to 0$ implies that the smearing happens in a neighborhood of the boundary $\partial\Omega_{N}$ which has a negligible volume compared to $|\Omega_{N}|$. The second condition $\delta_{N}|\Omega_{N}|^{1/3}\to\infty$ ensures that the kinetic energy in the transition region stays negligible in the thermodynamic limit.

Remark 1**.**

The same result holds under the weaker condition that $\Omega_{N}$ has an $\eta$–regular boundary, which means that $|\partial\Omega_{N}+B_{r}|\leqslant C|\Omega_{N}|\eta\left(r|\Omega_{N}|^{-1/3}\right)$ for all $r\leqslant|\Omega_{N}|^{1/3}/C$, with $\eta(t)\to 0$ when $t\to 0^{+}$. The condition $\delta_{N}|\Omega_{N}|^{1/3}\to\infty$ is then replaced by $\delta_{N}^{-2}\eta(\delta_{N}|\Omega_{N}|^{-1/3})\to 0$.

2.3. The Local Density Approximation

We are now able to state our main result.

Theorem 2** (Local Density Approximation).**

Let $p>3$ and $0<\theta<1$ such that

[TABLE]

There exists a constant $C=C(p,\theta,q)$ such that

[TABLE]

for every $\varepsilon>0$ and every non-negative density $\rho\in L^{1}(\mathbb{R}^{3})\cap L^{2}(\mathbb{R}^{3})$ such that $\nabla\sqrt{\rho}\in L^{2}(\mathbb{R}^{3})$ and $\nabla\rho^{\theta}\in L^{p}(\mathbb{R}^{3})$.

The constant $C=C(p,\theta,q)$ in our estimate (8) depends on the number of spin states $q$ ($q=2$ for electrons), in addition to the parameters $p$ and $\theta$. It diverges when $p\to 3^{+}$. If $p\to 3^{+}$ then we can take $\theta\to 5/6^{-}$. Our estimate therefore applies to densities $\rho$ with compact support, which vanish at the boundary of their support like $\delta(x)^{a}$ with $a>4/5$, where $\delta(x)={\rm d}(x,\partial\rho^{-1}(\{0\}))$. In particular, densities which vanish linearly are allowed. Our proof allows one to consider more singular densities, that is, to relax the constraint that $\theta p\leqslant 1+p/2$, but then the power of $\varepsilon$ deteriorates.

Our estimate (8) is certainly not optimal and it is an interesting challenge to improve it. We conjecture that a similar inequality holds with $\rho+\rho^{2}$ replaced by $\rho^{4/3}+\rho^{5/3}$ which have the scaling of the Coulomb and kinetic energies, respectively. The higher power $\rho^{2}$ arises from the trial state used in our upper bound (Proposition 1) and it is used to control some errors appearing when merging quantum systems with overlapping supports. This is explained in Section 4.1 below. Finally, the last gradient term in (8) is used to control local variations of $\rho$ in $L^{\infty}$. One could expect gradient errors involving only $\int_{\mathbb{R}^{3}}|\nabla\sqrt{\rho}|^{2}$ and $\int_{\mathbb{R}^{3}}|\nabla\rho^{1/3}|^{2}$ which are believed to arise in the gradient expansion of the uniform electron gas for, respectively, the Coulomb and kinetic energies.

One interesting case is when the density is given by a fixed function $\rho$ with $\int_{\mathbb{R}^{3}}\rho=1$, which is rescaled in the manner

[TABLE]

After taking $\varepsilon=N^{-1/12}$ in (8) we obtain the following simple bound

[TABLE]

It is conjectured [21, 24, 23, 26, 11] that the next order in the expansion of $E(\rho_{N})$ should involve the gradient correction to the kinetic energy

[TABLE]

and the gradient correction to the Coulomb energy

[TABLE]

In particular, the next order should be proportional to $N^{1/3}$. It remains an open problem to establish this rigorously.

In the classical case where the kinetic energy is neglected, the limit of $E_{\rm cl}(\rho_{N})/N$ was found in our previous work [30], but without a quantitative estimate on the remainder. We can give an estimate similar to (8) in the classical case, with a lower power of $\varepsilon$ in front of the gradient term. This is just a slight adaptation of the proof in [30], which is much easier than the quantum case. The argument is explained for completeness in Appendix A.

In the classical case the limit for $E_{\rm cl}(\rho_{N})/N$ was later extended to Riesz interactions $|x|^{-s}$ and other dimensions $d\geqslant 1$ in [8]. Although our result (8) in the quantum case can probably be extended to other Riesz interactions by using ideas from [12, 19, 15, 8], we only consider here the physically relevant 3D Coulomb case for shortness.

Remark 2** (Canonical case).**

We expect an inequality similar to (8) for the (mixed) canonical version of $E(\rho)$ where $\int_{\mathbb{R}^{3}}\rho=N\in\mathbb{N}$ and $\Gamma_{n}=0$ for $n\neq N$. However our proof does not adapt in an obvious way to this case.

Remark 3** (Exchange-correlation energy).**

In physics and chemistry, the LDA is usually expressed in terms of the exchange-correlation energy. In the grand-canonical setting it is defined by

[TABLE]

with $T(\rho)$ the lowest possible kinetic energy (5). The functional $T(\rho)$ is studied in Section 3 below, where it is proved that

[TABLE]

with $c_{\rm TF}=3^{5/3}4^{1/3}\pi^{4/3}/5$ the Thomas-Fermi constant. The lower bound was derived by Nam [43] and we prove the missing upper bound (with a better power of $\varepsilon$) in Theorem 3 below. Actually, by following our proof of Theorem 2 (simply discarding the Coulomb interaction) we can also prove a lower bound on $T(\rho)$, with an error similar to the right side of (8) but with a smaller power of $\varepsilon$ in front of $|\nabla\rho^{\theta}|^{p}$. This provides the following estimate on the exchange-correlation energy

[TABLE]

For a rescaled density $\rho_{N}(x)=\rho(x/N^{1/3})$ we obtain the same rate of convergence $N^{11/12}$ as in (9).

Remark 4** (Local Spin Density Approximation).**

In practice, it is often convenient to not fix the total density but, rather, the density of each spin component

[TABLE]

for $\sigma\in\{1,...,q\}$. Similarly as in Theorem 1 one can define the corresponding spin-polarized UEG energy $e_{\rm UEG}(\rho_{1},...,\rho_{q})$ of the uniform electron gas where the electrons of spin $\sigma$ are assumed to have the constant density $\rho_{\sigma}$. By following the arguments in this paper, one can then prove the estimate similar to (8)

[TABLE]

It is only for simplicity of notation that we work with the total density $\rho=\sum_{\sigma=1}^{q}\rho_{\sigma}$.

3. A priori estimates on $T(\rho)$ and $E(\rho)$

Lower bounds on $E(\rho)$ in (4) are well known and will be recalled below. Upper bounds are somewhat difficult to derive due to the constraint that the quantum states considered need to have the exact given density $\rho$. In this section we prove an upper bound on the best kinetic energy and use it to derive an upper bound on $E(\rho)$. Because our bounds are of independent interest we work in this section in any dimension $d\geqslant 1$. First we quickly recall the known lower bounds.

3.1. Known lower bounds

We recall that the Lieb-Thirring inequality [39, 40, 38] states that there exists a positive constant $c_{\rm LT}=c_{\rm LT}(d)>0$ such that

[TABLE]

for every self-adjoint operator $\gamma$ on $L^{2}(\mathbb{R}^{d},\mathbb{C}^{q})$ such that $0\leqslant\gamma\leqslant 1$. The best constant $c_{\rm LT}$ is unknown but has been conjectured to be the semi-classical constant

[TABLE]

in dimension $d\geqslant 3$. In [43], Nam has proved that

[TABLE]

for every $\varepsilon>0$ and some constant $\kappa=\kappa(d)$, in all space dimensions $d\geqslant 1$.

We also recall the Hoffmann-Ostenhof inequality [17] which states that

[TABLE]

and always imposes that $\sqrt{\rho}\in H^{1}(\mathbb{R}^{d})$. The inequality (16) does not require the fermionic constraint $0\leqslant\gamma\leqslant 1$.

The Lieb-Oxford inequality [32, 37, 20, 38] states that the total Coulomb energy is bounded from below by

[TABLE]

where $\Gamma=\{\Gamma_{n}\}$ is a grand-canonical quantum state satisfying the conditions in (4). Inspired by [2], this bound was recently generalized in [29] to

[TABLE]

but the constant $(3/5)({9\pi}/{2})^{1/3}\simeq 1.4508$ is not expected to be the optimal Lieb-Oxford constant.

Using (13) together with (17), we obtain the following.

Corollary 1** (Lower bound on $E(\rho)$).**

We have

[TABLE]

for every $\rho\geqslant 0$ such that $\sqrt{\rho}\in H^{1}(\mathbb{R}^{d})$.

The constants can be improved at the expense of adding gradient corrections, using (15) and (18).

3.2. Upper bound on the best kinetic energy

Let us recall that the lowest possible kinetic energy of a fixed density $\rho\in L^{1}(\mathbb{R}^{d},\mathbb{R}_{+})$ with $\sqrt{\rho}\in H^{1}(\mathbb{R}^{d})$ reads

[TABLE]

In [30], we have shown that for $\rho\leqslant 1$,

[TABLE]

by using the trial state

[TABLE]

Our goal in this section is to prove a similar bound without the assumption that $\rho\leqslant 1$. Coherent states [33] can usually give good bounds on the kinetic energy but they do not preserve the density. The main difficulty here is to construct a state having the exact given density $\rho$. The next result says that the semi-classical approximation to the kinetic energy is an upper bound to the exact $T(\rho)$, up to some gradient corrections.

Theorem 3** (Upper bound on the best kinetic energy).**

There are two constants $\kappa_{1},\kappa_{2}>0$ depending only on the space dimension $d$ such that

[TABLE]

for every $\rho\in L^{1}(\mathbb{R}^{d},\mathbb{R}_{+})$ with $\sqrt{\rho}\in H^{1}(\mathbb{R}^{d})$ and every $\varepsilon>0$, where $c_{\rm TF}$ is the Thomas-Fermi constant (14).

Note that the gradient correction in (20) has a better behavior in $\varepsilon$ than in Nam’s lower bound (15).

In dimension $d=1$, March and Young have given the proof of a better estimate without the parameter $\varepsilon$, in [41, Eq. (9)]:

[TABLE]

In the same paper they also state a result in 3D (for a constant $c>c_{\rm TF}$) but the proof has a mistake. This was mentioned as a conjecture in [35, Sec. 5.B]. Our result (20) can therefore be seen as a solution to the March-Young problem. We conjecture that a similar bound holds without the parameter $\varepsilon$ in dimension $d=2,3$ as well.

Remark 5** (Explicit constants in 3D).**

In dimension $d=3$ one can take

[TABLE]

These constants are not optimal and they are only displayed for concreteness. Our proof allows to slightly improve the constants under the assumption that $\varepsilon$ is small enough. For instance for $\varepsilon\leqslant 1$, we have the better inequality

[TABLE]

Proof of Theorem 3.

For simplicity we only write the proof in the no-spin case $q=1$. Recall that the free Fermi sea

[TABLE]

has the constant density $\rho_{P_{r}}=r$ and the constant kinetic energy density $c_{\rm TF}r^{1+2/d}$. In particular, if we take

[TABLE]

for some $f\geqslant 0$, we obtain

[TABLE]

See for instance [30, Sec. 5] for details. In addition, we have in the sense of operators $0\leqslant\gamma\leqslant f(x)$, hence $\gamma$ is a fermionic one-particle density matrix under the additional condition that $f\leqslant 1$.

Let now $\eta\geqslant 0$ be a smooth non-negative function such that

[TABLE]

Using the smooth layer cake principle

[TABLE]

we introduce the trial state

[TABLE]

In the sense of operators, we have

[TABLE]

In addition, $\gamma$ has the required density

[TABLE]

Hence $\gamma$ is admissible and it can be used to get an upper bound on $T(\rho)$. From (22), its kinetic energy is

[TABLE]

Note that

[TABLE]

Hence we have proved that

[TABLE]

for all $\eta\geqslant 0$ satisfying the two constraints (23). The smallest constant we can get in front of the $\rho^{1+2/d}$ term is $c_{\rm TF}$, by concentrating $\eta$ at the point $t=1$, but this makes the other term blow up. If we fix

[TABLE]

then the best constant we can get in front of the gradient term is given by the variational problem

[TABLE]

We claim that $C(\varepsilon)\leqslant{\rm const.}(1+\varepsilon^{-1})$ for $\varepsilon$ small enough, which we prove by an appropriate choice of $\eta$.

Let us first take, for instance,

[TABLE]

Then $\int\eta_{\varepsilon}=1$ and $\int\eta_{\varepsilon}(t)/t\,dt\leqslant 1$ since $\eta_{\varepsilon}$ is supported on $[1,\infty)$. Using the simple bounds

[TABLE]

and

[TABLE]

we obtain (after changing $\varepsilon$ into $\varepsilon/2$)

[TABLE]

The behavior of the correction in front of the semi-classical term $c_{\rm TF}$ is not optimal for small $\varepsilon$. It can be replaced by $1+\kappa_{1}\varepsilon^{2}$, for $\varepsilon\leqslant 1$. To see this we slightly translate the function (24) to the left by an amount $-\varepsilon b$ and introduce

[TABLE]

Then we have

[TABLE]

and

[TABLE]

The unique $b_{\varepsilon}$ such that $\int_{0}^{\infty}\eta_{\varepsilon,b_{\varepsilon}}(t)\frac{dt}{t}=1$ satisfies

[TABLE]

and for this $b_{\varepsilon}$ we have

[TABLE]

[TABLE]

This is how we can get (21) for $\varepsilon$ small enough (after replacing $\varepsilon^{2}$ by $\varepsilon$). ∎

Remark 6**.**

In the 3D case we can take for instance $b=1-\varepsilon/10-4\varepsilon^{3}/350$. One can then verify that

[TABLE]

and

[TABLE]

for all $\varepsilon\leqslant 1$. Hence

[TABLE]

Combining with (25), we find the estimate (20) for $\kappa_{1}=1$ and $\kappa_{2}=48$.

3.3. Upper bound on $E(\rho)$

It is well known that any fermionic one-particle density matrix $\gamma$ (i.e., an operator satisfying $0\leqslant\gamma=\gamma^{*}\leqslant 1$) is representable by a quasi-free state $\Gamma_{\gamma}$ in Fock space [1]. The two-particle density matrix of such a state is given by Wick’s formula

[TABLE]

In particular, the corresponding interaction energy with pair potential $w$ is

[TABLE]

From this we immediately obtain the following.

Corollary 2** (Upper bound on $E(\rho)$ in dimension $d\geqslant 1$).**

We have

[TABLE]

for every $\rho\geqslant 0$ such that $\sqrt{\rho}\in H^{1}(\mathbb{R}^{d})$.

This is for the grand-canonical version (4) of the Levy-Lieb functional which is the object of concern in this paper. It was proved in [34] that any fermionic $\gamma$ with integer trace $N={\rm Tr}\,(\gamma)$ is also the one-particle density matrix of an $N$-particle mixed state $\Gamma$ on the fermionic space $\mathfrak{H}^{N}$, such that the corresponding two-particle density matrix satisfies

[TABLE]

in the sense of operators. From the positivity of the Coulomb potential we deduce immediately the following result for mixed canonical states.

Corollary 3** (Upper bound in the mixed canonical case).**

Let $\rho\in L^{1}(\mathbb{R}^{d},\mathbb{R}_{+})$ be such that $\int_{\mathbb{R}^{d}}\rho=N\in\mathbb{N}$, and $\sqrt{\rho}\in H^{1}(\mathbb{R}^{d})$. Then there exists a mixed state $\Gamma$ on the fermionic space $\bigwedge_{1}^{N}L^{2}(\mathbb{R}^{d}\times\{1,...,q\})$ such that $\rho_{\Gamma}=\rho$ and

[TABLE]

4. Proof of Theorems 1 and 2

Our proof is divided into several steps. The first is to show that the energy is essentially local, that is, to prove that

[TABLE]

where $\{\chi_{k}\}$ is a smooth partition of unity, $\sum_{k}\chi_{k}=1$. The precise statement of (30) will involve upper and lower bounds, as well as an average over the translations, rotations and dilations of the partition itself. The lower bound was indeed already shown in [30] using the Graf-Schenker inequality [14]. The upper bound is the main new ingredient of our proof. The two bounds are derived in Sections 4.1 and 4.2. This will allow us to provide a rather simple proof of Theorem 1 in Section 4.4, using the convergence for tetrahedra which will be studied in Section 4.3.

In Section 4.5 we will estimate the deviation of the energy when we replace $\rho\chi_{k}$ by a constant function, say $\rho(x_{k})\chi_{k}$ for some $x_{k}$ in the support of $\chi_{k}$. If $\rho$ is essentially constant in the corresponding region, the error will be small, but if $\rho$ is not constant we bound the energy using some gradient terms, utilizing our upper bound (29).

After showing the Lipschitz regularity of $e_{\rm UEG}$ in Section 4.6, we will be able to conclude the proof of Theorem 2 in Section 4.7. We replace $E(\rho(x_{k})\chi_{k})$ by $\int_{\mathbb{R}^{3}}e_{\rm UEG}(\rho(x))\chi_{k}(x)\,dx$ when $\rho$ is large enough on ${\rm supp}(\chi_{k})$, using some quantitative estimates derived for tetrahedra in Section 4.3.

In the rest of the paper we call $C$ a generic constant which can sometimes change from line to line, but which only depends on $q$ (the number of spin states) and $p,\theta$, the two parameters appearing in the statement of Theorem 2.

4.1. Upper bound in terms of local densities

Our main goal here is to give an upper bound on the energy $E(\rho)$ by splitting $\rho$ into a sum of local densities. In the classical case, we have the exact subadditivity property (see [30, Lem. 2.5])

[TABLE]

which considerably simplifies the analysis and was one of the main tools of our previous work [30]. In particular we immediately find an upper bound in the form

[TABLE]

for a partition of unity $\chi_{k}$. In the quantum case this is not as easy. The first difficulty is that we cannot cut sharply and have to use a smooth partition of unity. This has the consequence that neighboring local densities overlap. But then, for two densities $\rho_{1}$ and $\rho_{2}$ with overlapping support, it is not obvious how to relate $E(\rho_{1}+\rho_{2})$ with $E(\rho_{1})$ and $E(\rho_{2})$. This is due to the fermionic nature of the electrons which puts a very strong constraint on trial states. If we take two trial quantum states for $\rho_{1}$ and $\rho_{2}$, we cannot simply take their tensor product and use it as a trial state for $\rho_{1}+\rho_{2}$. The tensor product does not have the fermionic symmetry, and if we anti-symmetrize it the density is not equal to $\rho_{1}+\rho_{2}$ anymore.

For this reason, we will use an incomplete partition of unity with holes, in order to make sure that the local quantum states are not overlapping. Since we want to get the exact density, holes are however in principle not allowed. Instead of filling the holes with electrons, our idea is to rather average over all the possible rotations, translations and dilations of the partition of unity, which will make the holes disappear in average. All the arguments of this section apply the same to a tiling made of cubes, but for a better matching with the lower bound we will consider a tiling made of tetrahedra. Our lower bound relies on the Graf-Schenker inequality [14] which requires the use of tetrahedra.

Let us consider the unit cube $C_{1}=(-1/2,1/2)^{3}$, which is the union of 24 disjoint identical tetrahedra $\mathbf{\Delta}_{1},...,\mathbf{\Delta}_{24}$, all of volume $1/24$. Since the cube can be repeated in the whole space, we obtain a tiling of $\mathbb{R}^{3}$ with tetrahedra:

[TABLE]

and the corresponding partition of unity

[TABLE]

for any fixed size $\ell>0$ of the tiles. Any $\mathbf{\Delta}_{j}$ can be written as $\mathbf{\Delta}_{j}=\mu_{j}\mathbf{\Delta}$ where $\mathbf{\Delta}$ is a reference tetrahedron with [math] as its center of mass. Here $\mu_{j}=(z_{j},R_{j})\in C_{1}\times SO(3)$ is an appropriate translation and rotation, which acts as $\mu_{j}x=R_{j}x-z_{j}$, hence $\mathbf{\Delta}_{j}=R_{j}\mathbf{\Delta}-z_{j}$. In each tetrahedron we now place the regularized characteristic function

[TABLE]

Here $\eta_{\delta}(x)=(10/\delta)^{3}\eta_{1}(10x/\delta)$ where $\eta_{1}$ is a fixed $C^{\infty}_{c}$ non-negative radial function with support in the unit ball and such that $\int_{\mathbb{R}^{3}}\eta_{1}=1$. Assuming that $\delta\leqslant\ell/2$, the function $\chi_{\ell,\delta,j}$ has its support well inside $\ell\mathbf{\Delta}_{j}$, at a distance proportional to $\delta$ from its boundary. The prefactor has been chosen to ensure that

[TABLE]

The function

[TABLE]

is equal to $(1-\delta/\ell)^{-3}>1$ inside the tiles but vanishes in a neighborhood of the boundary of the tiles. It is the incomplete partition of unity which we have mentioned above. We obtain a partition of unity after averaging over the translations of the tiling:

[TABLE]

Here $C_{\ell}=(-\ell/2,\ell/2)^{3}=\ell C_{1}$ is the cube of side length $\ell$. This is because for any $f\in L^{1}(\mathbb{R}^{3})$,

[TABLE]

The main result of this section is the following upper bound.

Proposition 1** (Upper bound in terms of local densities).**

There exists a universal constant $C$ such that for any $\sqrt{\rho}\in H^{1}(\mathbb{R}^{3})$, any $0<\delta<\ell/2$ and any $0<\alpha<1/2$,

[TABLE]

In particular, we can find $\ell^{\prime}\in(\ell(1-\alpha),\ell(1+\alpha))$, $\delta^{\prime}\in(\delta(1-\alpha),\delta(1+\alpha))$ and an isometry $(\tau,R)\in\mathbb{R}^{3}\times SO(3)$ such that

[TABLE]

The right side of (34) involves our incomplete partition of unity with holes, which is rotated (with the rotation $R$), translated (with the translation $\tau$) and dilated (with the dilation parameter $t$). The error is small only when $\delta$ (the size of the holes) is small. However we cannot take $\delta=0$ since that would make the gradient of the densities $\chi_{t\ell,t\delta,j}(R\,\cdot\,-t\ell z-\tau)\rho$ blow up. Nevertheless, the statement is that the energy decouples and the holes can be neglected, at the expense of an error of the order $\delta^{2}\int_{\mathbb{R}^{3}}\rho^{2}$. In (34) we use dilations for purely technical reasons, in order to better control error terms.

Proof of Proposition 1.

Using (33), we write our density $\rho$ as follows

[TABLE]

For every fixed $\tau\in C_{\ell}$, we can construct a grand canonical trial state $\Gamma_{\tau}$ having the density

[TABLE]

For this we pick $\Gamma_{\tau}=\bigotimes_{j=1}^{24}\bigotimes_{z\in\mathbb{Z}^{3}}\Gamma_{\tau,z,j}$ where each $\Gamma_{\tau,z,j}$ has the density

[TABLE]

and minimizes the corresponding energy $E(\chi_{\ell,\delta,j}(\cdot-\ell z-\tau)\rho)$. Since the quantum states $\Gamma_{\tau,z,j}$ have disjoint supports, we can anti-symmetrize the state $\Gamma_{\tau}$ in the standard manner. We denote by $\Gamma_{\tau,a}$ the anti-symmetrized state. The energy of $\Gamma_{\tau,a}$ is equal to that of $\Gamma_{\tau}$ and so is its density $\rho_{\Gamma_{\tau,a}}=\rho_{\Gamma_{\tau}}$. Finally we take as trial state

[TABLE]

which satisfies by construction that $\rho_{\Gamma}=\rho$. We find the upper bound

[TABLE]

Here we employ the usual notation

[TABLE]

In the second line of (37) we have used that the energy of a tensor product of states of disjoint supports is the sum of the energies of the pieces [30]. This is because the cross terms in the direct energy exactly cancel with the many-particle interactions of different states in the tensor product.

The error term in (37) is solely due to the nonlinearity of the direct term and it may be rewritten as

[TABLE]

For every real-valued $(\ell\mathbb{Z}^{3})$–periodic function $f$, we have

[TABLE]

Hence we obtain

[TABLE]

with

[TABLE]

for $\epsilon=\delta/\ell$. We have

[TABLE]

Since all the functions appearing in the sum on the right side of (39) are supported in the unit cube, we also obtain for $k\in 2\pi\mathbb{Z}^{3}\setminus\{0\}$,

[TABLE]

This results in the final formula for the error term

[TABLE]

In order to control the denominator $|p-k/\ell|^{2}$, we are going to average our calculation over all the rotations of the tiling. We also replace $\ell$ and $\delta$ by, respectively, $t\ell$ and $t\delta$ and we average over $t\in(1-\alpha,1+\alpha)$ with a weight $t^{-4}$. Rotating the tiling is the same as rotating $\rho$. In addition, $\epsilon=\delta/\ell$ is independent of $t$. Hence we are left with estimating

[TABLE]

with $p^{\prime}=p\ell/|k|$ and where $A_{\alpha}$ is the annulus

[TABLE]

We will use the following estimate

Lemma 1**.**

We have

[TABLE]

for all $\alpha\leqslant 1/2$.

The proof of (42) is a simple computation which is provided at the very end of the proof. Using (42) we obtain

[TABLE]

and it remains to estimate the sum in the parenthesis. For this we have to bound $\sum_{j=1}^{24}\widehat{{\mathds{1}}_{\mu_{j}(1-\varepsilon)\mathbf{\Delta}}}(k)$.

Lemma 2** (Fourier transform of the reduced tetrahedra).**

We have

[TABLE]

for all $0<\varepsilon<1/2$ and all $k\in 2\pi\mathbb{Z}^{3}\setminus\{0\}$.

Proof.

We recall that $\mathbf{\Delta}_{j}=R_{j}\mathbf{\Delta}_{1}-z_{j}$ with $\mu_{j}=(z_{j},R_{j})\in C_{1}\times SO(3)$. We have

[TABLE]

Since $k\in 2\pi\mathbb{Z}^{3}\setminus\{0\}$, the integral vanishes at $\varepsilon=0$ after summing over $j$. Inserting the derivative at $\varepsilon=0$ yields

[TABLE]

We claim that the second term is uniformly bounded with respect to $k$. Indeed, one integration by parts gives

[TABLE]

Here $n_{j}(x)$ is the normalized vector perpendicular to $\partial\mathbf{\Delta}_{j}$ pointing outwards. Integrating once more in the same manner (involving the edges of the faces of $\partial\mathbf{\Delta}_{j}$ for the first term), we see that (45) is bounded uniformly in $k$, hence we obtain (44). ∎

Inserting (44) in (43), we obtain the two error terms

[TABLE]

(using here that $\left\|\widehat{\eta_{1}}\right\|_{L^{\infty}}\leqslant(2\pi)^{3/2}$) and

[TABLE]

Recalling that $\varepsilon=\delta/\ell$, our final estimate on the averaged error is proportional to

[TABLE]

In order to conclude the proof of Proposition 1, it remains to provide the

Proof of Lemma 1.

We have

[TABLE]

For $|x|\leqslant 1/5$ and $0<\alpha<1/2$, the integrand is bounded on the corresponding interval and we obtain

[TABLE]

Similarly, for $|x|\geqslant 4$ and $0<\alpha<1/2$ the integrand can be estimated by $r^{2}$, which gives again

[TABLE]

Finally, for $1/5\leqslant|x|\leqslant 4$, we have

[TABLE]

The last integral is over an interval of length

[TABLE]

The integral is maximum when the interval is placed at the divergence point $r=1$. So we have

[TABLE]

∎

This concludes the proof of Proposition 1. ∎

4.2. Lower bound in terms of local densities

Next we turn to the lower bound. We are going to use the same tiling made of tetrahedra, with the difference that we do not insert any hole. Similarly to (32), we introduce

[TABLE]

which forms a smooth partition of unity, without holes,

[TABLE]

Proposition 2** (Lower bound in terms of local densities).**

There exists a universal constant $C$ such that for any $\sqrt{\rho}\in H^{1}(\mathbb{R}^{3})$ and any $\delta>0$ with $0<\delta/\ell<1/C$, we have

[TABLE]

In particular, we can find an isometry $(\tau,R)\in\mathbb{R}^{3}\times SO(3)$ such that

[TABLE]

Proof.

For a state $\Gamma=\bigoplus_{n\geqslant 0}\Gamma_{n}$ on Fock space (commuting with the particle number operator) and an interaction potential $w$, we introduce the simplified notation

[TABLE]

and

[TABLE]

For the Coulomb potential $w(x)=|x|^{-1}$ we simply use the notation $\mathcal{C}(\Gamma)$ and $D(\rho)$. For the kinetic energy, we write

[TABLE]

and finally denote by

[TABLE]

the total energy, with the direct term subtracted.

The proof uses the well known fact that, for any interaction potential $w$ and any state $\Gamma$ on Fock space, we have

[TABLE]

For the first bound see for instance [30, Eq. (4.8)]. The second bound uses only that $\mathcal{C}_{w}(\Gamma)\geqslant 0$ and $D_{w}(\rho)\leqslant\left\|w\right\|_{L^{1}}\left\|\rho\right\|_{L^{2}}^{2}/2$, by Young’s inequality.

We are now ready to prove (48). The smeared Graf-Schenker inequality from [14, Lemma 6] states that the potential

[TABLE]

has a positive Fourier transform for all $\ell>\kappa\delta$, with $\widetilde{w}_{\ell}(0)=-\kappa\ell^{-1}\widetilde{h}_{\ell,\delta}(0)$, where

[TABLE]

Here $\mathbf{\Delta}$ is a tetrahedron and $\kappa>0$ is a large enough constant. In addition, we have from [30, Proof of Lem. 5.5] that the potential

[TABLE]

is positive and has positive Fourier transform, with $\widetilde{W}_{\delta}(0)=\sqrt{2}\delta^{-2}$.

Arguing exactly as in [30, Lem. 5.5] using (52), we find that for any fermionic grand-canonical mixed state $\Gamma=\oplus_{n\geqslant 0}\Gamma_{n}$ with density $\rho_{\Gamma}$, we have

[TABLE]

Here $\Gamma_{|f}$ is the geometrically $f$–localized state on Fock space [9, 16, 27], that is, the unique state which has the $k$-particle reduced density matrices $f^{\otimes k}\Gamma^{(k)}f^{\otimes k}$. The last term proportional to $\delta^{3}/\ell$ comes from the $L^{1}$ norm of $\ell^{-1}\delta e^{-\frac{\sqrt{2}}{\delta}|x|}|x|^{-1}$.

For the kinetic energy we use the IMS formula as in [14, 16] and [30, Lem. 5.6], which yields an error in the form

[TABLE]

where $N=\int_{R^{3}}\rho_{\Gamma}$. For the total energy we obtain

[TABLE]

which yields the result. ∎

4.3. A convergence rate for tetrahedra

In this section we study the convergence of the energy per unit volume for tetrahedra and find a convergence rate. We introduce the energy per unit volume of a tetrahedron at constant density $\rho_{0}>0$

[TABLE]

where $\eta_{\delta}(x)=(10/\delta)^{3}\eta_{1}(10x/\delta)$ with $\eta_{1}$ a fixed $C^{\infty}_{c}$ non-negative radial function with support in the unit ball and such that $\int_{\mathbb{R}^{3}}\eta_{1}=1$. We prove the following

Proposition 3** (Thermodynamic limit for tetrahedra).**

For every fixed $\rho_{0}>0$, we have

[TABLE]

For $\delta\leqslant\ell/C$ and $0<\alpha<1/2$, we have the upper bound

[TABLE]

and the averaged lower bound

[TABLE]

If in addition $\rho_{0}^{1/3}\ell\geqslant C$, we have the pointwise lower bound

[TABLE]

The constant $C$ only depends on the chosen regularizing function $\eta_{1}$. It is independent of $\rho_{0},\ell,\delta,\alpha$.

We will later see that the condition $\delta^{3}/\ell\to 0$ is actually not needed in the limit (56). It is an interesting problem to replace the error term in the lower bound (59) by an error similar to the upper bound (57). Note that the error term in (58) goes to zero only when $\delta\to 0$ whereas (59) does not require $\delta\to 0$.

Proof.

For fixed $\rho_{0}>0$ and $\delta>0$, the existence of the limit (56) for $\ell\to\infty$ was proved in [30], using a lower bound similar to (48).

We consider a large tetrahedron $\ell^{\prime}\mathbf{\Delta}$, smeared at a scale $\delta^{\prime}$ and a tiling of smaller tetrahedra of size $\ell\ll\ell^{\prime}$, smeared at scale $\delta$. Applying our lower bound (47), we find

[TABLE]

where for the error term we have used that

[TABLE]

For all the tetrahedra such that $R(\ell\mu_{j}\mathbf{\Delta}-\ell z-\tau)+B_{\delta/10}\subset(\ell^{\prime}-\delta^{\prime})\mathbf{\Delta}$, we obtain exactly $|\ell\mathbf{\Delta}|\,e_{\mathbf{\Delta}}(\rho_{0},\ell,\delta)$ in the integral. The other tetrahedra are at a distance proportional to $\ell+\delta+\delta^{\prime}$ from the boundary of $\ell^{\prime}\mathbf{\Delta}$. Hence, using our lower bound (19) on the energy, they give rise to an error term of the order $\rho_{0}^{4/3}(\ell+\delta+\delta^{\prime})/\ell^{\prime}$. We obtain

[TABLE]

Here $\sigma=1$ if $e_{\mathbf{\Delta}}(\rho_{0},\ell,\delta)\geqslant 0$ and $\sigma=0$ otherwise.

After taking the limit $\ell^{\prime}\to\infty$ at fixed $\ell,\delta,\delta^{\prime},\rho_{0}$, we obtain

[TABLE]

It follows from our upper bound (29) (see also [30, Rmk. 5.4]) that

[TABLE]

for all $\rho_{0}>0$. Hence after dividing by $1-C\sigma\delta/\ell$ we have shown the claimed upper bound

[TABLE]

We may use exactly the same argument using our upper bound (34) in place of (48) and we obtain the lower bound (58).

Next we replace $\ell$ by $t\ell$ and $\delta$ by $t\delta$ in our lower bound (60) on the energy of the large simplex of size $\ell^{\prime}$, and average over $t\in(1/2,3/2)$ with the measure $t^{-4}$. We then insert our lower bound (58) and, after collecting the different error terms, we obtain

[TABLE]

It is natural to choose $\delta=\ell^{-1/3}(\rho_{0})^{-4/9}$ which provides the estimate

[TABLE]

and then $\ell=(\ell^{\prime})^{3/5}\rho_{0}^{-2/15}$ which gives

[TABLE]

under the assumption that $\ell^{\prime}(\rho_{0})^{1/3}\geqslant C$. This is exactly (59). The two bounds (57) and (59) give the limit (56). ∎

4.4. Proof of Theorem 1

Let $\Omega_{N}$ be a sequence of domains as in the statement, that is, such that $|\Omega_{N}|\to\infty$ and $|\partial\Omega_{N}+B_{r}|\leqslant Cr|\Omega_{N}|^{2/3}$ for all $r\leqslant|\Omega_{N}|^{1/3}/C$. Assume also that $\delta_{N}|\Omega_{N}|^{-1/3}\to 0$.

By following the proof of (60) we see that a similar inequality holds with the large tetrahedron $\ell^{\prime}\mathbf{\Delta}$ replaced by $\Omega_{N}$. This gives

[TABLE]

Under the sole condition that $\delta_{N}|\Omega_{N}|^{-1/3}\to 0$, the right side tends to $e_{\rm UEG}(\rho_{0})$ if we take for instance $\delta$ fixed and $\ell=|\Omega_{N}|^{1/6}$.

We then use the upper bound (34) with $\alpha=1/2$ as well as the fact that

[TABLE]

by (29), for the tetrahedra close to the boundary. We find

[TABLE]

We have used here that

[TABLE]

Under the additional assumption that $\delta_{N}|\Omega_{N}|^{1/3}\to\infty$, we may choose for instance $\ell=|\Omega_{N}|^{1/6}$ and $\delta=|\Omega_{N}|^{-1/12}$, which yields the result.∎

4.5. Replacing the local density by a constant density

The goal of this section is to provide estimates on the variation of the energy in a (smeared) tetrahedron, when we replace the local density by a constant, chosen to be either the minimum or the maximum of the density in the tetrahedron.

Proposition 4** (Replacing $\rho$ by a constant locally).**

Let $p>3$ and $0<\theta<1$ such that

[TABLE]

There exists a constant $C=C(p,\theta,q)$ such that, for $\ell\geqslant C$ and $\delta\leqslant\ell/C$, we have

[TABLE]

and

[TABLE]

for all $0<\varepsilon\leqslant 1/2$, where

[TABLE]

are respectively the minimum and maximum value of $\rho$ on the support of ${\mathds{1}}_{\ell\mathbf{\Delta}}\ast\eta_{\delta}$.

Under the assumption that $\int_{\ell\mathbf{\Delta}+B_{\delta}}|\nabla\rho^{\theta}|^{p}$ is finite, the density $\rho$ is continuous on $\overline{\ell\mathbf{\Delta}+B_{\delta}}$, so that $\underline{\rho}$ and $\overline{\rho}$ are well defined.

We have already discussed in the beginning of Section 4.1 the difficulty of deriving a subadditivity-type estimate relating $E(\rho_{1}+\rho_{2})$ to $E(\rho_{1})$ and $E(\rho_{2})$. The following lemma provides a rather rough inequality, which however will be sufficient for our purposes.

Lemma 3** (Rough subadditivity estimate).**

Let $\rho_{1},\rho_{2}\in L^{1}(\mathbb{R}^{3},\mathbb{R}_{+})$ be two densities such that $\sqrt{\rho_{1}},\sqrt{\rho_{2}}\in H^{1}(\mathbb{R}^{3})$. Then

[TABLE]

for all $0<\varepsilon\leqslant 1$.

Here we have in mind that $\rho_{2}$ is small compared to $\rho_{1}$ and we estimate $E(\rho_{1}+\rho_{2})$ in terms of $E(\rho_{1})$ plus some error terms. The worse error in the estimate (65) is $D(\rho_{2})/\varepsilon$, because it grows much faster than the volume. Later we will only use (65) locally and this bad term will not be too large. But it will be responsible for the large power of $\varepsilon$ in front of the gradient correction in our main estimate (8). We conjecture that there is an inequality similar to (65) without the term $D(\rho_{2})/\varepsilon$.

Note that we can estimate

[TABLE]

by the convexity of $\rho\mapsto|\nabla\sqrt{\rho}|^{2}$.

Proof of Lemma 3.

Fix an $\varepsilon\in(0,1]$ and consider two optimal states $\Gamma_{1}$ and $\Gamma_{2}$ in Fock space, for $\rho_{1}$ and $\rho_{2}/\varepsilon+\rho_{1}$, respectively. Then

[TABLE]

is a proper quantum state which has the density

[TABLE]

Inserting this trial state and using (29) for $E(\rho_{2}/\varepsilon+\rho_{1})$, we deduce that

[TABLE]

We have

[TABLE]

By the Lieb-Oxford inequality $E(\rho_{1})\geqslant-C\int_{\mathbb{R}^{3}}\rho_{1}^{4/3}$, and the result follows. ∎

We are now able to provide the

Proof of Proposition 4.

We write $\rho=\underline{\rho}+(\rho-\underline{\rho})$ and apply (65). We obtain

[TABLE]

In the first line we have used that $\underline{\rho}\leqslant\rho$ on the support of ${\mathds{1}}_{\ell\mathbf{\Delta}}\ast\eta_{\delta}$ and that ${\mathds{1}}_{\ell\mathbf{\Delta}}\ast\eta_{\delta}\leqslant 1$. First we can bound $\rho^{4/3}+\rho^{5/3}$ by $\rho+\rho^{2}$. Next, using

[TABLE]

and $\nabla(\rho-(1-\varepsilon)\underline{\rho})=\nabla\rho=2\sqrt{\rho}\nabla\sqrt{\rho}$, we can bound the gradient term pointwise by

[TABLE]

Since $\rho\geqslant\underline{\rho}$, we have $\varepsilon\rho\leqslant\rho-(1-\varepsilon)\underline{\rho}$ and hence

[TABLE]

This gives the estimate on the gradient term

[TABLE]

Next we estimate the terms involving $\rho-\underline{\rho}$ in terms of the gradient of $\rho^{\theta}$. We use the Sobolev inequality in the bounded set $\ell\mathbf{\Delta}+B_{\delta}$

[TABLE]

for $p>3$ and every continuous $u$ which vanishes at least at one point in $\ell\mathbf{\Delta}+B_{\delta}$ (we always assume $\delta\leqslant\ell/C$ so that $\ell\mathbf{\Delta}+B_{\delta}$ is included in a ball of radius proportional to $\ell$). By the Hardy-Littlewood-Sobolev inequality, this gives

[TABLE]

In the second estimate we have used that

[TABLE]

since $\theta\leqslant 1$. In the third estimate we have used Hölder’s inequality to obtain an integral to the power $1-2/p$. This yields some power of $\ell$ which has been taken into account in the first factor. In order to bound $\rho^{\frac{2p}{p-2}(1-\theta)}$ by $\rho+\rho^{2}$ we need that

[TABLE]

which is equivalent to our assumption (62).

Similarly, we can bound the other error term as follows

[TABLE]

where $0<a\leqslant 1$ is a parameter to be chosen. As before we need the condition

[TABLE]

in order to bound the last term by $\rho+\rho^{2}$. This is equivalent to

[TABLE]

where the left inequality is always satisfied under our assumption (62). If we choose $a=1$ then the upper bound on $p\theta$ is stronger than (62). Hence we rather choose $a=4/5$ and obtain (63).

The argument for (64) is similar. This time we write $\overline{\rho}=\rho+(\overline{\rho}-\rho)$ and obtain from (65)

[TABLE]

The gradient term can be bounded above by

[TABLE]

The other terms are estimated as before, using that $\overline{\rho}-\rho\leqslant\overline{\rho}$. ∎

4.6. Lipschitz regularity of $e_{\rm UEG}$

In this section we prove that the UEG energy $e_{\rm UEG}$ is locally Lipschitz. The main result is the following.

Proposition 5** (Lipschitz regularity of $e_{\rm UEG}$).**

There exists a universal constant $C$ so that

[TABLE]

for every $0\leqslant\rho^{\prime}\leqslant\rho$. In particular, we have

[TABLE]

Proof.

By scaling we have

[TABLE]

for $\alpha\leqslant 1$. This proves that

[TABLE]

Passing to the limit using Proposition 3, we find

[TABLE]

for every $0\leqslant\alpha=1-\varepsilon\leqslant 1$, hence

[TABLE]

Here we have used the notation $x_{-}=\max(-x,0)$ for the negative part. Using that $e_{\rm UEG}(\rho_{0})\geqslant-c_{\rm LO}\rho_{0}^{4/3}$ by the Lieb-Oxford inequality (with $c_{\rm LO}\leqslant 1.64$), we obtain

[TABLE]

for all $0\leqslant\varepsilon\leqslant 1$. This proves the upper bound in (69).

Similarly, we can write (still for $0\leqslant\alpha\leqslant 1$)

[TABLE]

This gives as before

[TABLE]

Using this time $e_{\rm UEG}(\rho_{0})\leqslant C\rho_{0}^{5/3}$, we obtain

[TABLE]

for all $0\leqslant\varepsilon\leqslant 1$. ∎

4.7. Proof of Theorem 2

We have derived all the estimates we need to prove the main inequality (8) in Theorem 2.

Let $\rho\in L^{1}(\mathbb{R}^{3},\mathbb{R}_{+})\cap L^{2}(\mathbb{R}^{3},\mathbb{R}_{+})$ be any density so that $\nabla\sqrt{\rho}\in L^{2}(\mathbb{R}^{3})$ and $\nabla\rho^{\theta}\in L^{p}(\mathbb{R}^{3})$. First we recall from our upper bound (29) and the lower bound (19) that

[TABLE]

Similarly, we have

[TABLE]

In particular, the inequality (8) is obvious for large $\varepsilon$ and we only have to consider small $\varepsilon$.

In our upper bound (34) and our lower bound (47), the worse coefficient involving $\ell$ and $\delta$ in front of $\rho+\rho^{2}$ is $\delta^{2}+1/(\ell\delta)$. This suggests to take

[TABLE]

which we do for the rest of the proof. In fact, in our proof we will replace $\ell$ and $\delta$ by $t\ell$ and $t\delta$ and average over $t\in[1/2,3/2]$.

Step 1. Upper bound

Let us first take $1/4\leqslant\varepsilon^{3/2}\ell\leqslant 2$ and $1/4\leqslant\delta\varepsilon^{-1/2}\leqslant 2$ and derive an upper bound on $E(\rho(1_{\ell\mathbf{\Delta}}\ast\eta_{\delta}))$. We recall that $\mathbf{\Delta}$ is a tetrahedron of volume $1/24$ as described in Section 4.1 and that $\eta_{\delta}(x)=(10/\delta)^{3}\eta_{1}(10x/\delta)$ with $\eta_{1}$ a fixed $C^{\infty}_{c}$ non-negative radial function with support in the unit ball and such that $\int_{\mathbb{R}^{3}}\eta_{1}=1$. We denote by

[TABLE]

the maximal and minimal values of $\rho$ on the support of ${\mathds{1}}_{\ell\mathbf{\Delta}}\ast\eta_{\delta}$, as in Proposition 4. We use the upper bound (63) from Proposition 4 which quantifies the error made when replacing $E(\rho(1_{\ell\mathbf{\Delta}}\ast\eta_{\delta}))$ by $E(\underline{\rho}(1_{\ell\mathbf{\Delta}}\ast\eta_{\delta}))$. For the latter we then use our estimate (57) in Proposition 3 on the energy of a smeared tetrahedron. With our choice (72) of $\ell$ and $\delta$ in terms of $\varepsilon$, this leads to

[TABLE]

In the first line we have bounded the error terms in (57) by

[TABLE]

in order to to simplify our final bound. In the support of ${\mathds{1}}_{\ell\mathbf{\Delta}}\ast\eta_{\delta}$ we have by (69)

[TABLE]

hence

[TABLE]

Similarly as we did for (68), we can bound for $0<a\leqslant 1$

[TABLE]

Again we need

[TABLE]

which is equivalent to

[TABLE]

where the left side is automatically satisfied under our main assumption (7). In order to get an error controlled by the other gradient terms, we need

[TABLE]

which requires $a\geqslant 2/5$. Taking $a=2/3$ provides the smallest power of $\varepsilon$. Collecting our estimates, we have proved the following upper bound on the energy in a tetrahedron

[TABLE]

Here we have considered a tetrahedron placed at the origin for simplicity, but we of course get a similar inequality for any tetrahedron, by translating and rotating $\rho$.

Next we recall our upper bound (34) on the total energy $E(\rho)$

[TABLE]

with $\chi_{\ell,\delta,j}:=(1-\varepsilon^{2})^{-3}{\mathds{1}}_{\ell\mu_{j}(1-\varepsilon^{2})\mathbf{\Delta}}\ast\eta_{\delta}$. We also recall from Section 4.1 that $\delta/\ell=(t\delta)/(t\ell)=\varepsilon^{2}$. Inserting (75) into (76) and using the fact (33) that $\chi_{t\ell,t\delta,j}$ forms a partition of unity after averaging over translations and rotations, we obtain

[TABLE]

Note that when we sum over the tiling, the sets $t\ell\mu_{j}\mathbf{\Delta}+B_{t\delta}$ have finitely many intersections, which just results in a bigger constant in front of $|\nabla\rho^{\theta}|^{p}$. We have also used that

[TABLE]

From Proposition 5 and (71), we have

[TABLE]

hence we obtain the desired upper bound

[TABLE]

for $\varepsilon$ small enough.

Step 2. Lower bound

The lower bound is slightly more tedious since all our lower estimates involve $\overline{\rho}$ which can in general not be bounded by $\rho$. We shall argue as follows. First we average our lower bound (47) over $t$. This gives

[TABLE]

We recall that here $\xi_{\ell,\delta,j}={\mathds{1}}_{\ell\mu_{j}\mathbf{\Delta}}\ast\eta_{\delta}$, see Section 4.2. In order to use the same argument as for the upper bound, we are going to prove the estimate

[TABLE]

That the last integral is over the larger set $2\ell\mathbf{\Delta}+B_{2\delta}$ will only affect the multiplicative constant $C$. Inserting (80) into (79) gives a bound as in (78) but in the opposite direction. This concludes the proof of the theorem and it therefore only remains to prove (80).

With an abuse of notation we consider the minimal and maximal values over the larger set $2(\ell\mathbf{\Delta}+B_{\delta})$,

[TABLE]

instead of the corresponding definitions on the smaller set $\ell\mathbf{\Delta}+B_{\delta}$. First we again recall that, by (29) and (19), we have

[TABLE]

Hence there is nothing to prove when

[TABLE]

This is the case if $\overline{\rho}^{1/3}\leqslant C\varepsilon$, for instance. Hence we may assume in the following that $\overline{\rho}\geqslant C\varepsilon^{3}$ and that

[TABLE]

By (66), this implies

[TABLE]

that is,

[TABLE]

Under our assumption (7) on $p$ and $\theta$ the exponent is positive, hence we deduce that for $\varepsilon$ small enough

[TABLE]

on $2\ell\mathbf{\Delta}+B_{2\delta}$. With this additional information we can use our previous estimates.

By arguing exactly as in the proof of Proposition 4 with $\overline{\rho}$ the maximum over $2\ell\mathbf{\Delta}+B_{2\delta}$ instead of the support of ${\mathds{1}}_{t\ell\mathbf{\Delta}}\ast\eta_{t\delta}$, we get the estimate similar to (64)

[TABLE]

From the fact that $\overline{\rho}\leqslant C\rho$, the second term on the right side can be bounded by

[TABLE]

Then we average over $t$ and use our lower estimate (58) on the averaged energy of a tetrahedron. This gives

[TABLE]

The last term can again be bounded by

[TABLE]

and included into the average over $t$. Finally, using (69) and $\overline{\rho}\leqslant C\rho$, we infer that

[TABLE]

on the support of ${\mathds{1}}_{t\ell\mathbf{\Delta}}\ast\eta_{t\delta}$. To conclude the proof of (80) we can proceed in the same way as for the upper bound (75). This concludes the proof of Theorem 2. ∎

Appendix A Classical case

In the classical case where the kinetic energy is neglected, the grand canonical energy functional is defined [30] by

[TABLE]

where each $\mathbb{P}_{n}$ is a symmetric probability measure on $(\mathbb{R}^{3})^{n}$ with density

[TABLE]

This classical energy (83) is obtained from the quantum energy in the limit

[TABLE]

see [5, 3, 28, 6]. When $\int_{\mathbb{R}^{3}}\rho=N\in\mathbb{N}$, the canonical version of $E_{\rm cl}$ reads

[TABLE]

In [30] we have shown that for $\rho_{N}(x)=\rho_{1}(N^{-1/3}x)$ with $\int_{\mathbb{R}^{3}}\rho_{1}=1$,

[TABLE]

where

[TABLE]

is the energy per unit volume of the classical uniform electron gas at density $1$. In this appendix we quickly explain how to derive the following quantitative estimate on the convergence rate in (85).

Theorem 4** (Estimate in the (grand-canonical) classical case).**

Let $p>3$ and $0<\theta<1$ such that $\theta p\geqslant 4/3$. There exists a universal constant $C=C(p,\theta)$ such that

[TABLE]

with

[TABLE]

for every $\varepsilon>0$ and every non-negative density $\rho\in L^{1}(\mathbb{R}^{3})\cap L^{4/3}(\mathbb{R}^{3})$ such that $\nabla\rho^{\theta}\in L^{p}(\mathbb{R}^{3})$.

Under the condition that $\theta p\leqslant 1+p/3$, which is slightly more restrictive than in the quantum case (7), we get the much smaller power $2p-1$ of $\varepsilon$ in front of the gradient term. Then, after optimizing (86) in $\varepsilon$, we obtain the quantitative estimate

[TABLE]

for every $\rho_{N}(x)=\rho_{1}(N^{-1/3}x)$ and for $4/3\leqslant\theta p\leqslant 1+p/3$. The rate $N^{5/6}$ is better than the $N^{11/12}$ obtained in the quantum case, but still far from the expected rate $N^{1/3}$.

Proof.

The estimate (86) follows from the Lieb-Oxford inequality (17) (without the gradient term) when $\varepsilon$ is large, so we only have to consider the case where $\varepsilon$ is small.

We use again the tiling (31) and, as in the proof in [30], the upper bound

[TABLE]

where $\underline{\rho}=\min_{\ell\mu_{j}\mathbf{\Delta}+\ell z}\rho$. The inequality (88) is a consequence of the subadditivity and the negativity of $E_{\rm cl}$. Now it follows from [30, Cor. 3.4] and from the Graf-Schenker inequality as in Section 4.3, that in a tetrahedron

[TABLE]

This provides the upper bound

[TABLE]

For the rest of the argument we use the notation $\varepsilon=1/\ell$. In each tetrahedron we can follow the argument in (74) and estimate

[TABLE]

with $0<a\leqslant 1$. In order to estimate the second term by $\rho+\rho^{4/3}$, we need that

[TABLE]

which is equivalent to

[TABLE]

and which we assume for the rest of the proof.

We finally turn to the lower bound. Up to an appropriate rotation and translation of the tiling (or, equivalently, of the density $\rho$), the Graf-Schenker inequality gives the following lower bound [30, p. 100]

[TABLE]

We recall that

[TABLE]

by the Lieb-Oxford inequality (17). We have nothing to prove in any tetrahedron $\ell\mu_{j}\mathbf{\Delta}+\ell z$ such that

[TABLE]

Here $A$ is a large constant to be chosen later. This is in particular the case when $\overline{\rho}^{1/3}=\max_{\ell\mu_{j}\mathbf{\Delta}+\ell z}\rho\leqslant\varepsilon/c_{\rm LO}$. So we may assume that

[TABLE]

and that $\overline{\rho}>(\varepsilon/c_{\rm LO})^{3}$. This implies

[TABLE]

The power of $\varepsilon$ is non-negative when, again,

[TABLE]

For $A$ large enough (or $\varepsilon$ small enough) this gives $\overline{\rho}\leqslant C\underline{\rho}\leqslant C\rho$ in the simplex $\ell\mu_{j}\mathbf{\Delta}+\ell z$. The rest of the argument is then exactly the same as for the upper bound.

As a conclusion we obtain the bound (86) with the error term

[TABLE]

and the restrictions that $p>3$, $0<\theta\leqslant 1$, $0<a\leqslant 1$ and

[TABLE]

In order to minimize the power of $\varepsilon$ we want to take $a$ as large as possible, that is,

[TABLE]

For $\theta p\leqslant 1+p/3$ we take $a=1$ whereas for $\theta p>1+p/3$ we choose the other value and get the stated inequality (86). ∎

Remark 7** (Canonical case).**

As was mentioned in (85), in [30] we could also handle the canonical case. We would easily obtain a quantitative estimate on the canonical energy $E^{\rm can}_{\rm cl}(\rho)$ if we knew the speed of convergence of $\ell^{-3}E_{\rm cl}^{\rm can}(\rho_{0}{\mathds{1}}_{\ell\mathbf{\Delta}})$ to its limit $c_{\rm UEG}|\mathbf{\Delta}|$. Unfortunately, the argument used in [30, Lem. 3.2] to prove that the limit coincides with the grand canonical one does not seem to produce a quantitative bound.