Floating Wigner crystal with no boundary charge fluctuations

TL;DR

This paper introduces a modified trial state for the classical Homogeneous Electron Gas that suppresses boundary charge fluctuations, demonstrating the equivalence of three different ground state energy definitions.

Contribution

It proposes a novel boundary treatment in the trial state for Jellium, unifying various definitions of the ground state energy.

Findings

Boundary charge fluctuations are effectively suppressed.

The three definitions of the ground state energy are shown to coincide.

The modified trial state improves understanding of boundary effects in Jellium.

Abstract

We modify the "floating crystal" trial state for the classical Homogeneous Electron Gas (also known as Jellium), in order to suppress the boundary charge fluctuations that are known to lead to a macroscopic increase of the energy. The argument is to melt a thin layer of the crystal close to the boundary and consequently replace it by an incompressible fluid. With the aid of this trial state we show that three different definitions of the ground state energy of Jellium coincide. In the first point of view the electrons are placed in a neutralizing uniform background. In the second definition there is no background but the electrons are submitted to the constraint that their density is constant, as is appropriate in Density Functional Theory. Finally, in the third system each electron interacts with a periodic image of itself, that is, periodic boundary conditions are imposed on the interaction potential.

Equations

$$\begin{gathered} {\mathcal{E}}_{\mathrm{Jel}}({\Omega }_N,{\mathbf{x}}_1,...,{\mathbf{x}}_N)=\sum _{1⩽j<k⩽N}\frac{1}{|{\mathbf{x}}_j-{\mathbf{x}}_k|} \\ -\sum _{j=1}^N{\int }_{{\Omega }_N}\frac{\mathrm{d}\mathbf{y}}{|{\mathbf{x}}_j-\mathbf{y}|}+\frac{1}{2}{\iint }_{{\Omega }_N\times {\Omega }_N}\frac{\mathrm{d}\mathbf{y}\mathrm{d}\mathbf{z}}{|\mathbf{y}-\mathbf{z}|}. \end{gathered}$$

$$e_{\mathrm{Jel}}=\underset{{\Omega }_N\nearrow {\mathbb{R}}^3}{\lim }\underset{{\mathbf{x}}_1,...,{\mathbf{x}}_N\in {\mathbb{R}}^3}{\min }\frac{{\mathcal{E}}_{\mathrm{Jel}}({\Omega }_N,{\mathbf{x}}_1,...,{\mathbf{x}}_N)}{|{\Omega }_N|}.$$

$$e_{\mathrm{Jel}}={\zeta }_{\mathrm{BCC}}(1)\simeq -1.4442$$

$${\mathcal{E}}_{\mathrm{per},L}({\mathbf{x}}_1,...,{\mathbf{x}}_N)=\sum _{1⩽j<k⩽N}{G}_L({\mathbf{x}}_j-{\mathbf{x}}_k)+\frac{N}{2L}M$$

$$G_L(\mathbf{x})=\frac{G_1(\mathbf{x}/L)}{L}=\frac{4\pi }{L^3}\sum _{\begin{array}{c}\mathbf{k}\in (2\pi /L){\mathbb{Z}}^3\\ \mathbf{k}\ne 0\end{array}}\frac{e^{i\mathbf{k}\cdot \mathbf{x}}}{k^2}$$

$$M=\underset{\mathbf{r}\to 0}{\lim }\left(G_1(\mathbf{r})-\frac{1}{r}\right).$$

$$e_{\mathrm{per}}=\underset{L\to \infty }{\lim }\underset{{\mathbf{x}}_1,...,{\mathbf{x}}_N\in C_L}{\min }\frac{{\mathcal{E}}_{\mathrm{per},L}({\mathbf{x}}_1,...,{\mathbf{x}}_N)}{L^3}$$

$$\begin{gathered} {\mathcal{E}}_{\mathrm{Ind}}(\rho ):=\underset{{\rho }_{\mathbb{P}}=\rho }{\min }{\int }_{{\mathbb{R}}^{3N}}\sum _{1⩽j<k⩽N}\frac{\mathrm{d}\mathbb{P}({\mathbf{x}}_1,...,{\mathbf{x}}_N)}{|{\mathbf{x}}_j-{\mathbf{x}}_k|} \\ -\frac{1}{2}{\iint }_{{\mathbb{R}}^3\times {\mathbb{R}}^3}\frac{\rho (\mathbf{x})\rho (\mathbf{y})}{|\mathbf{x}-\mathbf{y}|}\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{y} \end{gathered}$$

$$e_{\mathrm{UEG}}=\underset{{\Omega }_N\nearrow {\mathbb{R}}^3}{\lim }\frac{{\mathcal{E}}_{\mathrm{Ind}}({\text{mathds}1}_{{\Omega }_N})}{|{\Omega }_N|}.$$

$$\lim_{N\to\infty}\frac{\mathcal{E}_{\rm Ind}\big{(}\rho(\cdot/N^{1/3})\big{)}}{N}=e_{\rm UEG}\int_{{\mathbb{R}}^{3}}\rho(\mathbf{x})^{\frac{4}{3}}\,{\rm d}\mathbf{x}.$$

$${\int }_{{\mathbb{R}}^{3N}}\sum _{1⩽j<k⩽N}\frac{\mathrm{d}\mathbb{P}({\mathbf{x}}_1,...,{\mathbf{x}}_N)}{|{\mathbf{x}}_j-{\mathbf{x}}_k|}-\frac{1}{2}{\iint }_{{\Omega }_N\times {\Omega }_N}\frac{\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{y}}{|\mathbf{x}-\mathbf{y}|}$$

$$={\int }_{{\mathbb{R}}^{3N}}{\mathcal{E}}_{\mathrm{Jel}}({\Omega }_N,{\mathbf{x}}_1,...,{\mathbf{x}}_N)\mathrm{d}\mathbb{P}({\mathbf{x}}_1,...,{\mathbf{x}}_N)$$

$$⩾\underset{{\mathbf{x}}_1,...,{\mathbf{x}}_N\in {\mathbb{R}}^3}{\min }{\mathcal{E}}_{\mathrm{Jel}}({\Omega }_N,{\mathbf{x}}_1,...,{\mathbf{x}}_N).$$

$${\mathcal{E}}_{\mathrm{Ind}}({\text{mathds}1}_{{\Omega }_N})⩾\underset{{\mathbf{x}}_1,...,{\mathbf{x}}_N\in {\mathbb{R}}^3}{\min }{\mathcal{E}}_{\mathrm{Jel}}({\Omega }_N,{\mathbf{x}}_1,...,{\mathbf{x}}_N).$$

$$e_{\mathrm{UEG}}⩾e_{\mathrm{Jel}}.$$

$$e_{\mathrm{UEG}}⩽{\zeta }_{\mathrm{BCC}}(1)\simeq -\mathrm{1.4442.}$$

$${\Omega }_N=\bigcup _{j=1}^N(Q+{\mathbf{x}}_j),$$

$$\tilde{\mathbb{P}}={\int }_Q{\delta }_{{\mathbf{x}}_1+\mathbf{a}}\otimes \cdots \otimes {\delta }_{{\mathbf{x}}_N+\mathbf{a}}\mathrm{d}\mathbf{a}$$

$$\frac{1}{2N}\sum _{1⩽j<k⩽N}\frac{1}{|{\mathbf{x}}_j-{\mathbf{x}}_k|}-\frac{1}{2N}{\iint }_{{\Omega }_N\times {\Omega }_N}\frac{\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{y}}{|\mathbf{x}-\mathbf{y}|}.$$

$${\zeta }_{\mathrm{BCC}}(1)+\frac{2\pi }{3}{\int }_Q{x}^2\mathrm{d}\mathbf{x}\simeq -\mathrm{0.9507.}$$

$$\frac{1}{N}{\int }_Q{\mathcal{E}}_{\mathrm{Jel}}({\Omega }_N-\mathbf{a},{\mathbf{x}}_1,...,{\mathbf{x}}_N)\mathrm{d}\mathbf{a}.$$

$$\mathbb{P}={\int }_Q{\delta }_{{\mathbf{x}}_1+\mathbf{a}}\otimes \cdots \otimes {\delta }_{{\mathbf{x}}_N+\mathbf{a}}\otimes {\left(\frac{{\text{mathds}1}_{C\setminus ({\Omega }_N+\mathbf{a})}}{M}\right)}^{\otimes M}\mathrm{d}\mathbf{a}.$$

$${\rho }_{\mathbb{P}}$$

$$={\text{mathds}1}_C+{\text{mathds}1}_{{\Omega }_N}-{\text{mathds}1}_{{\Omega }_N}\ast {\text{mathds}1}_Q$$

$$D(f,g):=\frac{1}{2}{\iint }_{{\mathbb{R}}^3\times {\mathbb{R}}^3}\frac{f(\mathbf{x})g(\mathbf{y})}{|\mathbf{x}-\mathbf{y}|}\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{y}$$

$${\int }_{({\mathbb{R}}^3{)}^{N+M}}\sum _{1⩽j<k⩽N+M}\frac{1}{|{\mathbf{z}}_j-{\mathbf{z}}_k|}\mathrm{d}\mathbb{P}({\mathbf{z}}_1,...,{\mathbf{z}}_{N+M})$$

$$=\sum _{1⩽i<j⩽N}\frac{1}{|{\mathbf{x}}_j-{\mathbf{x}}_k|}+\sum _{j=1}^N{\int }_Q{\int }_{C\setminus ({\Omega }_N+\mathbf{a})}\frac{\mathrm{d}\mathbf{a}\mathrm{d}\mathbf{y}}{|{\mathbf{x}}_j+\mathbf{a}-\mathbf{y}|}$$

$$+\left(1-\frac{1}{M}\right){\int }_{Q}D\left({\text{mathds}1}_{C\setminus ({\Omega }_N+\mathbf{a})}\right)\mathrm{d}\mathbf{a}$$

$$=\sum _{1⩽i<j⩽N}\frac{1}{|{\mathbf{x}}_j-{\mathbf{x}}_k|}-\sum _{j=1}^N{\int }_{{\Omega }_N}\frac{\mathrm{d}\mathbf{y}}{|{\mathbf{x}}_j-\mathbf{y}|}+D({\text{mathds}1}_{{\Omega }_N})$$

$$+D({\text{mathds}1}_C)+2D({\text{mathds}1}_C,{\text{mathds}1}_{{\Omega }_N}-{\text{mathds}1}_{{\Omega }_N}\ast {\text{mathds}1}_Q)$$

$$-\frac{1}{M}{\int }_{Q}D\left({\text{mathds}1}_{C\setminus ({\Omega }_N+\mathbf{a})}\right)\mathrm{d}\mathbf{a}.$$

$$\displaystyle\int_{({\mathbb{R}}^{3})^{N+M}}\sum_{1\leqslant j<k\leqslant N+M}\frac{{\rm d}\mathbb{P}(\mathbf{z}_{1},...,\mathbf{z}_{N+M})}{|\mathbf{z}_{j}-\mathbf{z}_{k}|}-D\big{(}\rho_{\mathbb{P}}\big{)}$$

$$\displaystyle\qquad=\mathcal{E}_{\rm Jel}(\Omega_{N},\mathbf{x}_{1},...,\mathbf{x}_{N})-D\big{(}{\mathds{1}}_{\Omega_{N}}-{\mathds{1}}_{\Omega_{N}}\ast{\mathds{1}}_{Q}\big{)}$$

$$-\frac{1}{M}{\int }_{Q}D\left({\text{mathds}1}_{C\setminus ({\Omega }_N+\mathbf{a})}\right)\mathrm{d}\mathbf{a}.$$

$${\mathcal{E}}_{\mathrm{Ind}}({\rho }_{\mathbb{P}})⩽{\mathcal{E}}_{\mathrm{Jel}}({\Omega }_N,{\mathbf{x}}_1,...,{\mathbf{x}}_N).$$

Full text

I Introduction

The Homogeneous Electron Gas, also called Jellium, is a fundamental system in quantum physics and chemistry Parr and Yang (1994); Giuliani and Vignale (2005). In this paper we introduce a modified “floating crystal” trial state and use it to prove that three possible definitions of the Jellium ground state energy coincide in the thermodynamic limit. In particular, we resolve a conundrum originating in Navet et al. (1980); Choquard et al. (1980); Borwein et al. (1988, 1989, 2014) and raised again in Lewin and Lieb (2015), where it was observed that the usual floating crystal trial state fails for Coulomb interactions.

In its original formulation, due to Wigner Wigner (1934), Jellium is defined as an infinite gas of electrons placed in a positively-charged uniform background. The thermodynamic limit of this system has been rigorously established in Lieb and Narnhofer (1975). This model provides a good description of the deep interior of white dwarfs Salpeter (1961); Baus and Hansen (1980) (where the point charges are the fully ionized atoms evolving in a uniform background of negatively charged electrons). It has also been shown to be of high relevance for valence electrons in alkaline metals, for instance in solid sodium Huotari et al. (2010).

A similar system appears in the Local Density Approximation of Density Functional Theory (DFT) Parr and Yang (1994), where it plays a central role for deriving functionals Perdew (1991); Perdew and Wang (1992); Becke (1993); Perdew et al. (1996). In DFT the density is fixed and there is no background. The natural system arising in this situation is an infinite gas of electrons submitted to the constraint that its density is constant over the whole space. This model was called the Uniform Electron Gas (UEG) in Lewin and Lieb (2015); Lewin et al. (2018, 2019) to avoid any possible confusion with Jellium.

At low density, the electrons in Jellium are believed to form a BCC Wigner crystal Wigner (1934, 1938), hence their density is not at all constant. It has nevertheless been assumed by many authors that the two definitions should coincide. The reason for this belief is that the Wigner crystal has no preferred position and orientation, hence one may consider the mixed state obtained by uniformly averaging over the position of the lattice. This state is sometimes called the floating crystal Bishop and Lührmann (1982); Mikhailov and Ziegler (2002); Drummond et al. (2004) and it has a constant density.

With long range potentials such as Coulomb, one should however be very careful, since boundary effects can easily play a decisive role. It was proved in Lewin and Lieb (2015) that, for the classical gas, computing the energy of the floating crystal in a thermodynamic limit leads to a much higher energy than the Jellium energy of the BCC crystal, with a shift of the order of the volume of the sample due to charge fluctuations close to the boundary. This is very specific to the Coulomb case, which is critical as far as the computation of the energy is concerned. No shift arises for potentials decaying slightly faster than $1/r$ at infinity.

In this paper we provide a simple and physically intuitive proof of the equality of the Jellium and UEG ground state energies in the thermodynamic limit, by explaining how to modify the floating crystal trial state. Our argument is to immerse the floating crystal into a thin layer of fluid. The small layer of fluid around the floating crystal is used to compensate the large charge fluctuations at the boundary of the system, which are responsible for the undesired shift of the energy. The trial state suggests that the UEG ground state in infinite volume is indeed the uniform average of the Jellium crystal, as was believed. In a finite system, the particles of the UEG are probably not crystallized in a neighborhood of the boundary, however.

Our argument will use a third definition of Jellium, which has always been of high relevance in practical computations Brush et al. (1966); Hansen (1973); Pollock and Hansen (1973); Ceperley and Chester (1977); Jones and Ceperley (1996). In this third point of view, the electrons are placed on a large torus without any background, whereas the Coulomb potential is replaced by a periodized version without zero mode. The problem of showing that the periodic system has the same thermodynamic limit has a long history for short range potentials Fisher and Lebowitz (1970); Angelescu and Nenciu (1973). For the Coulomb potential, a rigorous proof seems to have been provided only recently, in a series of works Sandier and Serfaty (2015); Petrache and Serfaty (2017); Rougerie and Serfaty (2016); Leblé and Serfaty (2017); Cotar and Petrache (2019). For completeness we will also give a simple argument for this important fact.

During the preparation of this work, the equality of the Jellium and UEG energies was claimed in a preprint (Cotar and Petrache, 2019, version 5). But the argument is long and indirect. Contrary to our trial state approach, it does not seem to provide any insight on the possible form of the UEG ground state.

The paper is organized as follows. In the next section we introduce the three definitions. In Section III we explain how to modify the floating crystal argument to prove an upper bound on the UEG energy. If Jellium was rigorously proved to be crystallized, this would complete the proof of the equality of Jellium and the UEG. In Section IV we apply the modified floating crystal argument in the case of a unit cell of large but fixed side length $L$, containing $N=L^{3}$ electrons. After passing to the thermodynamic and taking $L\to\infty$ in a second step, this gives an upper bound on the UEG energy in terms of the periodic energy. Finally, in Section V we give a simple proof that the periodic problem coincides with the Jellium problem in the thermodynamic limit, which concludes the proof of the equality of the three definitions. Section VI contains a discussion on how our argument can be generalized to other interaction potentials and other space dimensions. We particularly consider the case of Riesz interaction potentials $r^{-s}$.

II Three definitions of the ground state energy

We only discuss here the classical case where the kinetic energy is dropped. We expect that a similar construction should apply to the quantum model but are unable to make this work at the moment. This is due to the Pauli principle which makes it difficult to merge two quantum systems with overlapping supports, as is explained in Lewin et al. (2019) and is needed to add the thin layer of fluid around the Wigner crystal.

By scaling we may assume in the classical case that the density is $\rho=1$. The Jellium energy of $N$ point charges in a background $\Omega_{N}\subset{\mathbb{R}}^{3}$ (a domain with volume $|\Omega_{N}|=N$) is given by

[TABLE]

For any given $\Omega_{N}$ we may minimize the energy over the positions $\mathbf{x}_{j}$. It does not matter whether we constrain the point charges to stay in $\Omega_{N}$ or allow them to visit the whole space ${\mathbb{R}}^{3}$. After minimization they will always all end up in $\Omega_{N}$, since the energy is a harmonic function outside of $\Omega_{N}$ with respect to each $\mathbf{x}_{j}$, when the other particles are fixed. It was also proved by one of us Lieb that the point charges in $\Omega_{N}$ must have a universal positive distance to each other, a theorem that was recently used in Rougerie and Serfaty (2016); Lieb et al. (2019). We define the Jellium ground state energy per unit volume by

[TABLE]

Under some natural technical conditions on $\partial\Omega_{N}$, the limit was proved to exist and to be independent of the sequence $\Omega_{N}$ in Lieb and Narnhofer (1975). The reader may think of $\Omega_{N}=N^{1/3}\Omega$ where $\Omega$ is a fixed open convex set of volume $|\Omega|=1$, for instance a cube or a ball. It is a famous conjecture Wigner (1934, 1938) that the electrons crystallize on a BCC lattice, that is,

[TABLE]

where $\zeta_{\rm BCC}(s)$ is the Epstein Zeta function of the (density one) BCC lattice, see Coldwell-Horsfall and Maradudin (1960); Blanc and Lewin (2015) and (Giuliani and Vignale, 2005, p. 43).

Next we turn to periodic Jellium, which is formally obtained when we repeat periodically a Jellium configuration in the whole space and compute its energy per unit volume. For simplicity we work with a cube (that is, we place the particles on the torus), but the argument is the same for other tilings. For $N=L^{3}$ we introduce

[TABLE]

where

[TABLE]

with $G_{1}$ the ${\mathbb{Z}}^{3}$–periodic Coulomb potential, that is, the unique solution of the equation $-\Delta G_{1}=4\pi(\sum_{\mathbf{z}\in{\mathbb{Z}}^{3}}\delta_{\mathbf{z}}-1)$ such that $\int_{C_{1}}G_{1}=0$, with $C_{1}=(-1/2,1/2)^{3}$ the unit cube. The constant $M$ appearing in (3) is the Madelung constant of the cubic lattice which may be defined by

[TABLE]

In another point of view, $M/2=\zeta_{{\mathbb{Z}}^{3}}(1)$ is the Jellium energy per unit volume of the cubic lattice, that is, the interaction of each particle with all its periodic images. Except for the unimportant constant $M/(2L)$ which disappears in the thermodynamic limit, one can obtain (3) from (1) by replacing $1/r$ by the periodic function $G_{L}(\mathbf{r})$ whenever $\Omega_{N}$ is a box. This is because $\int_{C_{L}}G_{L}=0$ hence the two background terms disappear. We define the ground state energy by

[TABLE]

with $C_{L}=(-L/2,L/2)^{3}$. The limit on the right clearly exists when $L=2^{n}L_{0}$ because we can use as trial state a $2^{n}L_{0}$–periodic configuration in a cube of size $2^{n+1}L_{0}$, hence the right side is decreasing. The existence of the limit for $L\to\infty$ was proved in Sandier and Serfaty (2015); Petrache and Serfaty (2017); Rougerie and Serfaty (2016); Leblé and Serfaty (2017) but it will also be a consequence of our analysis.

We finally turn to the UEG ground state energy. In this case there is no background but the electrons are assumed to form a constant charge density, say over a given set $\Omega_{N}\subset{\mathbb{R}}^{3}$. The indirect energy of a given density $\rho$ with $\int_{{\mathbb{R}}^{3}}\rho(\mathbf{r}){\rm d}\mathbf{r}=N$ reads

[TABLE]

where the first minimum is taken over all $N$-particle probability measures $\mathbb{P}$ with one-particle density $\rho$. Since the electrons are indistinguishable we should restrict ourselves to symmetric $\mathbb{P}$’s, but the minimum is the same. Note that $\mathcal{E}_{\rm Ind}(\rho)$ can be obtained from the Levy-Lieb functional of DFT Levy (1979); Lieb (1983) by taking $\hbar\to 0$ or, equivalently, scaling the density in the manner $\lambda^{3}\rho(\lambda\mathbf{x})$ with $\lambda\to 0$ Seidl (1999); Cotar et al. (2013); Bindini and De Pascale (2017); Lewin (2018); Cotar et al. (2018); Grossi et al. (2017). The ground state energy per unit volume of the UEG is given by

[TABLE]

It was proved in Lewin et al. (2018) that the limit exists under the same conditions on $\Omega_{N}$ as for (2).

One can replace the characteristic function ${\mathds{1}}_{\Omega_{N}}$ by any sequence of densities $\rho_{N}$ which are equal to 1 well inside $\Omega_{N}$ (at a distance $\ell\ll|\Omega_{N}|^{1/3}$ from the boundary), equal to 0 well outside, and which stay bounded in the transition region. While such a $\rho_{N}$ is not exactly constant, we proved in Lewin et al. (2018) that $\mathcal{E}_{\rm Ind}(\rho_{N})$ has the same thermodynamic limit as in (6). We shall take advantage of this relaxed formulation in the following.

As shown in Lewin et al. (2018, 2019), the constant $e_{\rm UEG}$ naturally arises in the Local Density Approximation of DFT. For instance, we have for a very spread out density in the form $\rho(\mathbf{x}/N^{1/3})$

[TABLE]

The classical UEG has been the object of many recent numerical works, based on methods from optimal transportation Seidl (1999); Seidl et al. (1999, 2007); Gori-Giorgi and Seidl (2010); Seidl et al. (2017). In addition to providing interesting properties of DFT at low density, the classical UEG has been used to get numerical bounds on the best constant in the Lieb-Oxford inequality Lieb (1979); Lieb and Oxford (1980); Lieb and Seiringer (2010); Odashima and Capelle (2007); Räsänen et al. (2009); Lewin and Lieb (2015).

For any $N$-particle probability measure $\mathbb{P}$ such that $\rho_{\mathbb{P}}={\mathds{1}}_{\Omega_{N}}$, we have

[TABLE]

Hence, after optimizing over $\mathbb{P}$ we obtain

[TABLE]

After passing to the thermodynamic limit this yields the lower bound

[TABLE]

The question of equality has been left open. Our main result is the following

Theorem 1**.**

We have $e_{\rm Jel}=e_{\rm per}=e_{\rm UEG}$.

The proof will be given in Sections IV and V.

III The floating crystal

Before showing Theorem 1 and as an illustration of the main idea, we first prove that

[TABLE]

If we had a proof that Jellium is crystallized in a BCC lattice, then this would immediately imply that $e_{\rm UEG}=e_{\rm Jel}$, due to (9). Note that (10) also implies that the best constant in the Lieb-Oxford inequality Lieb (1979); Lieb and Oxford (1980); Lewin and Lieb (2015) is at least as large as $-\zeta_{\rm BCC}(1)\simeq 1.4442$.

We first explain the floating crystal and the problem associated with its use as a trial state for estimating the UEG energy. We use the same notation as in (Lewin and Lieb, 2015, App. B). Let $\mathcal{L}$ be the BCC lattice, with Wigner-Seitz unit cell $Q$ centered at 0, such that $|Q|=1$ Ashcroft and Mermin (1976). We place the particles on the intersection of the lattice $\mathcal{L}$ with a large cube $C^{\prime}$ and call $\mathcal{L}\cap C^{\prime}=\{\mathbf{x}_{1},...,\mathbf{x}_{N}\}$ the corresponding positions of the $N$ particles. We then take

[TABLE]

the union of the cells centered at the particles. The floating crystal Bishop and Lührmann (1982); Mikhailov and Ziegler (2002); Drummond et al. (2004) is obtained by taking the delta function distribution of the $N$ particles, then translating by an amount $\mathbf{a}\in{\mathbb{R}}^{3}$ and integrating $\mathbf{a}$ over the unit cell $Q$. This corresponds to the $N$-particle probability

[TABLE]

which has the constant density $\rho_{\widetilde{\mathbb{P}}}={\mathds{1}}_{\Omega_{N}}$. The indirect energy per particle of this state is

[TABLE]

In the limit $N\to\infty$, it has been shown in (Lewin and Lieb, 2015, App. B) to converge to

[TABLE]

By (7) the indirect energy per particle (12) can also be written in terms of a moving background in the form

[TABLE]

As explained in (Lewin and Lieb, 2015, App. B), the difference ${\mathds{1}}_{\Omega_{N}-\mathbf{a}}-{\mathds{1}}_{\Omega_{N}}$ describes a monopole layer in a neighborhood of the surface which produces an electric potential felt by all the particles in the system. This survives in the thermodynamic limit and gives rise to the positive shift in (13).

We now explain our key new idea to avoid the energy shift. We immerse the crystal in a thin layer of fluid of density one. As we average over the positions of the crystal, the fluid gets displaced as depicted in Figure 1. To this end, we choose a slightly larger cubic container $C$ such that $\Omega_{N}+Q\subset C$ with the volume of the fluid $|C\setminus\Omega_{N}|=M$ being an integer. We will need $M\ll N$, so that the fluid layer around the floating crystal has a negligible volume in the thermodynamic limit. In practice, we choose the cube $C$ to be at a finite distance to the boundary of $\Omega_{N}$, this distance being larger than the diameter of $Q$. Then $M\sim N^{2/3}$. Our new trial state has the $N$ particles on the floating crystal, translated by $\mathbf{a}\in Q$ as before, together with $M$ other particles forming a fluid in $C\setminus(\Omega_{N}+\mathbf{a})$, the set remaining after we have subtracted the moving background:

[TABLE]

Note that the fluid is correlated with the position of the crystal. A sketch of the set-up is depicted in Figure 1. The density of this trial state equals

[TABLE]

with $f\ast g(\mathbf{x})=\int_{{\mathbb{R}}^{3}}f(\mathbf{y})g(\mathbf{x}-\mathbf{y})\,{\rm d}\mathbf{y}$ the convolution between two functions. In order to compute the Coulomb energy of $\mathbb{P}$, it is convenient to denote the Hartree energy by

[TABLE]

and to use the shorthand notation $D(f):=D(f,f)$. Then we find

[TABLE]

The first line is the Jellium energy $\mathcal{E}_{\rm Jel}(\Omega_{N},\mathbf{x}_{1},...,\mathbf{x}_{N})$ of the finite crystal, whereas the second line is equal to $D(\rho_{\mathbb{P}})-D({\mathds{1}}_{\Omega_{N}}-{\mathds{1}}_{\Omega_{N}}\ast{\mathds{1}}_{Q})$. Hence we have shown that the indirect energy of our trial state equals

[TABLE]

Using that $D(f)\geqslant 0$, the last two terms can be neglected for an upper bound. We have therefore proved that

[TABLE]

The function $\rho_{\mathbb{P}}$ is equal to 1 when $\mathbf{x}$ is inside $\Omega_{N}$, at a distance at least equal to the diameter of $Q$ from the boundary $\partial\Omega_{N}$, whereas it is equal to 0 outside of $C$. It varies between 0 and 2 in the intermediate region. Since $M\ll N$ we can make use of the relaxed formulation in Lewin et al. (2018) mentioned after (6) to conclude that

[TABLE]

Therefore, after passing to the limit in (18), we find the claimed upper bound

[TABLE]

IV Upper bound on the Uniform Electron Gas energy

The previous upper bound (19) is not enough to conclude that $e_{\rm UEG}=e_{\rm Jel}$ since we have no rigorous proof that Jellium is crystallized. However, the previous section contains the main idea. Here we show that

[TABLE]

for any fixed $\mathbf{x}_{1},...,\mathbf{x}_{n}$ points in the cube $C_{\ell}$, with $\ell^{3}=n$. The bound (19) simply corresponds to $n=1$ for the BCC lattice, but for simplicity we work here with the cubic lattice. After minimizing over the $\mathbf{x}_{j}$ and passing to the thermodynamic limit $\ell\to\infty$, we obtain the inequality

[TABLE]

It has been shown in Cotar and Petrache (2019) based on results from Sandier and Serfaty (2015); Petrache and Serfaty (2017); Rougerie and Serfaty (2016); Leblé and Serfaty (2017) that $e_{\rm per}=e_{\rm Jel}$, hence this concludes the proof of the theorem. Since the proof in these works is quite long and delicate, we provide a simpler argument based on Lieb and Narnhofer (1975) in the next section, for completeness. But here we concentrate on proving (20) and (21).

Let us consider $n$ distinct points $\mathbf{x}_{1},...,\mathbf{x}_{n}$ inside the cube $C_{\ell}$ and denote $\tau=n^{-1}\sum_{j=1}^{n}\mathbf{x}_{j}$ their center of mass. If we shift the background by $\tau$, we obtain a configuration with no dipole moment:

[TABLE]

Next we repeat our configuration $\mathbf{x}_{j}$ periodically in space and add a layer of fluid as before, in a background shifted by $\tau$. This is done as follows. We define the large cube of side length $\ell(2K+1)$

[TABLE]

with $N=\ell^{3}(2K+1)^{3}$ and pick $C$ to be a slightly larger cube so that $\Omega_{N}+2C_{\ell}\subset C$ and $|C\setminus\Omega_{N}|=M\ll N$. Our trial state is, similarly as in the previous section, given by

[TABLE]

and it has the density

[TABLE]

This density is again equal to one well inside $\Omega_{N}$ and 0 outside of $C$. Defining $\mathbf{x}_{n+1},...,\mathbf{x}_{N}$ to be the $\mathbf{x}_{j}+\ell\mathbf{k}$ with $\mathbf{k}\neq 0$ (ordered in any chosen fashion), the exact same calculations as in the previous section give

[TABLE]

Hence we obtain

[TABLE]

As before, when $N\to\infty$ and $M/N\to 0$, the left side converges to $e_{\rm UEG}$. Since the repeated configuration has no dipole moment, it is a well-known fact that

[TABLE]

This concludes the proof of (20), hence of (21).

For completeness, we briefly explain how to derive the limit (23). We start with the upper bound, which turns out to be sufficient for our purpose. Since the points $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ are strictly inside $C_{\ell}$, the periodically repeated points are located at a positive distance from each other, independent of $N$. Then we replace the point charges by small uniform balls of radius $\eta$. By Newton’s theorem, this does not change the interaction between the point charges, whereas the interaction with the background is increased. With $\chi_{\eta}(\mathbf{r})=\eta^{-3}{\mathds{1}}(4\pi r^{3}/3\leqslant\eta^{3})$ we obtain

[TABLE]

The density in the parenthesis equals

[TABLE]

with $f=\sum_{j=1}^{n}\chi_{\eta}(\cdot-\mathbf{x}_{j})-{\mathds{1}}_{C_{\ell}+\tau}$. Passing to Fourier variables we can write

[TABLE]

Note that $\int_{C_{\ell}}f(\mathbf{x})\,{\rm d}\mathbf{x}=0$ and that $f$ has no dipole moment:

[TABLE]

Hence the function $|\widehat{f}(\mathbf{p})|^{2}/p^{2}$ is continuous and vanishes at $p=0$. On the other hand, the square of the Dirichlet kernel converges weakly to the Dirac comb of the dual lattice

[TABLE]

as $K\to\infty$. Going back to configuration space, this gives the convergence

[TABLE]

Passing then to the limit $\eta\to 0$ using that

[TABLE]

we obtain, as was claimed, the upper bound in (23).

The proof of the lower bound in (23) is similar. It requires to estimate the error made in the interaction with the background when we replace the point particles by uniform balls. By Newton’s theorem, there is no error when the ball is outside of the background. In case of an intersection the error can be bounded by

[TABLE]

a term which disappears in the thermodynamic limit since we take $\eta\to 0$.

V Upper bound on the periodic energy

We give here a short proof of the inequality

[TABLE]

Our strategy follows the one of Lieb and Narnhofer (1975), which is based on the earlier work in Lieb and Lebowitz (1972) and is also described in Lieb and Seiringer (2010). In combination with (9) and (21) this completes the proof of Theorem 1.

We are going to use the important fact that Newton’s theorem holds in the periodic cell, for neutral systems. More precisely, if we have a radial charge distribution $\rho$ compactly supported in a ball $B_{R}$ and such that $\int_{B_{R}}\rho(\mathbf{r}){\rm d}\mathbf{r}=0$, then for $L$ large enough so that $B_{R}\subset C_{L}$ we claim that $\rho\ast G_{L}=\rho\ast r^{-1}$ in $C_{L}$. This is because $V=\rho\ast 1/r$ vanishes outside of $B_{R}$, by Newton’s theorem, hence the periodized potential $\sum_{\mathbf{k}\in{\mathbb{Z}}^{d}}V(\mathbf{r}+\ell\mathbf{k})$ solves the same equation as $\rho\ast G_{L}$ and we must have $\rho\ast G_{L}=\rho\ast 1/r+K$ in $C_{L}$. The constant $K$ is found to vanish after integration over $C_{L}$. We also infer that $\rho\ast G_{L}$ vanishes on $C_{L}\setminus B_{R}$.

We use the Swiss cheese theorem (Lieb and Seiringer, 2010, Sec. 14.5) to cover the cube $C_{L}$ with many balls (all of integer volume) of sizes ranging from some fixed $\ell_{0}$ to the largest one of order $\ell$. The volume not covered by the balls is small compared to $L^{3}$ if $\ell$ is large, and in particular goes to zero relative to $L^{3}$ if $\ell\to\infty$ after $L\to\infty$ (for fixed $\ell_{0}$ or, more generally, if $\ell_{0}\ll\ell$). In each ball $B_{n}$, we place $N_{n}=|B_{n}|$ particles in the optimal Jellium configuration of the ball. The remaining $M=N-\sum_{n}|B_{n}|$ particles are placed uniformly in the left-over cheese $S=C_{L}\setminus\bigcup_{n}B_{n}$. We obtain an upper bound on the minimal energy in the box

[TABLE]

of the form

[TABLE]

where $\mathbf{x}_{1},...,\mathbf{x}_{N-M}$ denote the positions of the $N-M$ particles in $\cup_{n}B_{n}$. We then use that $S=C_{L}\setminus\cup_{n}B_{n}$ and the fact that $\int_{C_{L}}G_{L}=0$. Discarding the term of order $1/M$ for an upper bound, we find that the right side is bounded above by

[TABLE]

This is exactly the energy obtained by putting point particles together with a uniform background over $\cup_{n}B_{n}$. Next we can average the particle configurations in each ball $B_{n}$ over rotations. Due to Newton’s theorem recalled above, this cancels the interactions between the systems in different balls. We obtain the upper bound

[TABLE]

where $\mathbf{x}_{j}^{(n)}$ denote the point charges in the ball $B_{n}$ and $\widetilde{G_{L}}$ denotes the average of $G_{L}$ over rotations of the ball $B_{n}$. As an upper bound, we thus obtain the sum of the jellium energy in each ball $B_{n}$, with interaction kernel $\widetilde{G_{L}}(\mathbf{x}-\mathbf{y})$ in place of $|\mathbf{x}-\mathbf{y}|^{-1}$. As $L\to\infty$, the former converges to the latter, however. Hence dividing by $L^{3}$ and taking the successive limits $L\to\infty$, $\ell\to\infty$ and $\ell_{0}\to\infty$, we arrive at the desired result.

VI Extension to Riesz potentials in all space dimensions

Our argument in Sections III and IV applies to more general potentials in any dimension, since we have essentially only used that the interaction has a positive Fourier transform, so that $D(f)\geqslant 0$. Here we quickly describe how to generalize our findings to Riesz potentials, which are defined by

[TABLE]

For instance, $s=1$ is the 3D Coulomb case which can also be considered in dimensions $d=1,2$. The case $s=0$ plays a central role in many situations. This is the natural interaction arising in random matrix theory for $d=1,2$ Forrester (2010). It also arises in the study of star polymer solutions, at least at short distances, see Witten and Pincus (1986) and (Likos, 2001, Sec. 5). It is very convenient to enclose all these important physical situations in the one-parameter family of Riesz interactions. This has been useful to better understand how the decay of correlations Alastuey and Martin (1985); Martin (1988) and sum rules Gruber et al. (1980); Martin and Yalcin (1980); Gruber et al. (1981); Fontaine and Martin (1984); Martin (1988) depend on the decay of the potential, that is, the parameter $s$. It does not seem to be a well known fact that adding the background in the spirit of Wigner is a very robust method which, as we will demonstrate, works for all $-2\leqslant s<d$ in any dimension, and not only in the Coulomb case.

We define the Jellium energy of $N$ point particles by

[TABLE]

and always assume $s<d$ to ensure the finiteness of the last two terms. Of course, no background is necessary in the short range case $s>d$. The following says that the system is thermodynamically stable for all $-2\leqslant s<d$.

Lemma 1** (Stability for Riesz potentials).**

Let $d\geqslant 1$ and $-2\leqslant s<d$. We have, for a universal constant $C(d,s)$,

[TABLE]

for every $\mathbf{x}_{1},...,\mathbf{x}_{N}\in{\mathbb{R}}^{d}$ and every bounded open set $\Omega\subset{\mathbb{R}}^{d}$.

Many authors work under the constraint that $s\geqslant d-2$, but the Coulomb case $s=d-2$ is not a natural threshold in the family of Riesz potentials. Since the previous result does not seem to be well known, we provide a proof in Appendix A. A similar lower bound was previously derived in (Cotar and Petrache, 2019, App. B.2) for $0<s<d$.

We can now define the lowest Jellium energy in a given background $\Omega_{N}$ with $|\Omega_{N}|=N$ by

[TABLE]

This function is uniformly bounded from below, due to Lemma 1. We claim that it is also bounded from above for “reasonable” sets. To prove this we have to construct one trial state with an energy of order $N$. Taking the uniform average for all the points in $\Omega_{N}$, we find

[TABLE]

For $s>0$ the right side is negative, proving that $e_{{\rm Jel},d,s}(\Omega_{N})\leqslant 0$. For $s<0$ the right side diverges to $+\infty$ as $N\to\infty$, and the uniform average is not a good trial state. Let us instead consider $\Omega_{N}$ to be the union of $N$ smaller cubes of size one ($\Omega_{N}$ can be made a cube if $N=K^{3}$ with $K\in{\mathbb{N}}$). In each of the small cube we put exactly one particle, which we average uniformly over its small cube only. This cancels exactly the background and we are just left with the self-energies of the small cubes:

[TABLE]

This is of order one as claimed. This argument applies to all $s<d$ and any $\Omega_{N}$ which can be partitioned into $N$ sets of volume one and uniformly bounded diameter. This leads us to conjecture that the Jellium model with Riesz interaction has a thermodynamic limit for all $-2\leqslant s<d$ in any dimension.

We would like to consider the corresponding energy

[TABLE]

where $\Omega_{N}$ is any reasonable sequence of domains like cubes or balls, with $|\Omega_{N}|\in{\mathbb{N}}$. The existence of this limit has been obtained for $s=d-2$ in any dimension in Kunz (1974); Lieb and Narnhofer (1975); Sari and Merlini (1976) and for $d-2<s<d$ (resp. $0\leqslant s<d$ for $d=1,2$) in Petrache and Serfaty (2017); Leblé and Serfaty (2017); Cotar and Petrache (2019). To our knowledge no proof has yet been given for smaller values of $s$. In those cases we define $e_{\rm Jel}(d,s)$ by a $\liminf$ instead of a limit.

We then consider the indirect energy

[TABLE]

which satisfies as in (8)

[TABLE]

for every domain $\Omega_{N}$. Our upper bound (28) applies to the UEG as well, showing that $|\Omega_{N}|^{-1}\mathcal{E}_{{\rm Ind},d,s}({\mathds{1}}_{\Omega_{N}})$ is uniformly bounded for “reasonable” sets. By following the proof of (Lewin et al., 2018, Thm. 2.6), based on the subadditivity of the indirect energy, one can show that the limit

[TABLE]

exists and does not depend on $\Omega_{N}$, for regular-enough sequences $\Omega_{N}\nearrow{\mathbb{R}}^{d}$ and for all $-2\leqslant s<d$.

Finally, we consider the periodic problem

[TABLE]

where $N=L^{d}$ and the periodic function $G_{d,s,L}$ has its Fourier coefficients equal to $c_{d,s}k^{s-d}$ for $0\neq\mathbf{k}\in(2\pi/L){\mathbb{Z}}^{d}$ and the appropriate constant $c_{d,s}$. Here $M_{d,s}=2{\rm sgn}(s)\zeta_{{\mathbb{Z}}^{d}}(s)$ is the corresponding Madelung constant. At $s=0$ the term $M_{d,s}/(2L^{s})$ is replaced by $\zeta_{{\mathbb{Z}}^{d}}^{\prime}(0)-\zeta_{{\mathbb{Z}}^{d}}(0)\ln L$. The limit

[TABLE]

exists for subsequences in the form $L=2^{n}L_{0}$, due to the monotonicity of the energy. In the absence of a proof for general sequences we define $e_{\rm per}(d,s)$ by a liminf instead of a limit.

The Coulomb case $s=-1$ in dimension $d=1$ is completely understood. It is well known that Jellium is crystallized Baxter (1963); Choquard (1975); Kunz (1974) and a calculation furnishes

[TABLE]

On the other hand, it has been proved in Colombo et al. (2015) that the floating crystal, defined similarly as in (11), is the exact ground state of the indirect energy $\mathcal{E}_{\rm Ind}({\mathds{1}}_{[-N/2,N/2]})$. Hence the computations in Lewin and Lieb (2015) imply that there is an energy shift:

[TABLE]

Jellium and the UEG differ at $s=-1$ in 1D. Other values of $s$ are considered in Leblé (2015); Di Marino (2019).

Next we discuss the adaptation of the argument in Section IV to Riesz potentials. Our result is the following.

Theorem 2**.**

In space dimension $d\geqslant 1$ we have

[TABLE]

for all $\max\left(0,d-4\right)<s<d$. There is equality for $0<s<d$ in dimensions $d=1,2$ and for $d-2\leqslant s<d$ in dimensions $d\geqslant 3$.

The proof goes as follows. The first inequality is an immediate consequence of (30). For the second inequality we follow the argument in Section IV. The computation (22) continues to hold for the potential $V_{s}$. The last two Hartree terms in this equation continue to be negative when $s>0$ because the Fourier transform of $V_{s}$ is positive. The first Hartree term on the second line of (22) is negative even for $-2\leqslant s\leqslant 0$ because the function in the argument has a vanishing integral (see Appendix A). But the Hartree term on the third line is positive for $s<0$ and it has no particular sign for $s=0$. This term is a $O(MN^{-s/d})$, or $O(M\ln N)$ for $s=0$. In our case where $M\sim N^{1-1/d}$, the last term in (22) is a $o(N)$ under the condition that $s>-1$.

The convergence to the periodic Jellium energy in (23) continues to hold under the condition that $d-s<4$ for configurations which have no dipole moment. The convergence of $\mathcal{E}_{\rm Ind}(\rho_{\mathbb{P}})/N$ to $e_{\rm UEG}(d,s)$ was established in Cotar and Petrache (2019) for all $s>0$. Hence we obtain the inequality for $\max\left(0,d-4\right)<s<d$, as claimed. Should the last convergence hold for all $s\geqslant-2$, as we believe, then the same theorem would hold under the weaker assumption that $\max\left(-1,d-4\right)<s<d$.

Finally, the equality $e_{\rm per}(d,s)=e_{\rm Jel}(d,s)$ is shown in Petrache and Serfaty (2017); Leblé and Serfaty (2017); Cotar and Petrache (2019) for $\max(0,d-2)<s<d$ (resp. $0\leqslant s<d$ in $d=1,2$). Note that our proof that $e_{\rm per}(d,s)\leqslant e_{\rm Jel}(d,s)$ in Section V easily extends to $s=d-2$ in all dimensions, but not to other values of $s$, because it relies on Newton’s theorem.

VII Conclusion

In this paper we have compared definitions of the minimum energy of the uniform electron gas and of jellium. For many years it has been an open problem to prove rigorously that they are the same to leading order in the volume since it was known that the obvious method for constructing the uniform gas definitely did not lead to the desired equivalence. We have succeeded in proving the equivalence, and thus provide a firm foundation for some aspects of density functional theory.

Appendix A Proof of Lemma 1 on the stability of Jellium with Riesz potentials

We start with $s>0$. Our idea is to replace $V_{s}$ by a smaller potential $0\leqslant V_{s,M}\leqslant V_{s}$ which is continuous at the origin and still has a positive Fourier transform. Then we use that, for this potential,

[TABLE]

with $\mu=\sum_{j=1}^{N}\delta_{\mathbf{x}_{j}}-{\mathds{1}}_{\Omega}$. To define the potential $V_{s,M}$ we follow Fefferman and de la Llave (1986); Hainzl and Seiringer (2002) and first remark that for any radial function $\chi\geqslant 0$ with $\int_{{\mathbb{R}}^{d}}\chi=1$,

[TABLE]

where

[TABLE]

This suggests to introduce the truncated potential

[TABLE]

which satisfies $V_{s,M}(0)=\frac{c(s)M^{s}}{s}\int_{{\mathbb{R}}^{d}}\chi^{2}$ and

[TABLE]

Using $V_{s,M}\leqslant V_{s}$ for the self energies of the particles and of the background, as well as

[TABLE]

for their mutual interaction, we find

[TABLE]

Optimizing over $M$ we obtain

[TABLE]

Our conclusion is that

[TABLE]

where

[TABLE]

The best bound is obtained after optimizing over $\chi$ but here we keep it fixed for simplicity. In the limit $s\to 0$, we have

[TABLE]

whereas

[TABLE]

When $N=|\Omega|$, the term $-N/2$ cancels on both sides of (34). Hence we can divide by $s$ and pass to the limit $s\to 0$. We get the claimed estimate for $s=0$, with the appropriate definition of $C_{\chi}(d,0)$.

Finally, for $s<0$ we can write directly

[TABLE]

The Fourier transform of $\mu$ has the behavior at the origin

[TABLE]

where $\mathbf{P}=\sum_{j=1}^{N}\mathbf{x}_{j}-\int_{\Omega}\mathbf{x}\,{\rm d}\mathbf{x}$ is the corresponding dipole moment. Under the assumption that $N=|\Omega|$ and $-2<s<0$, we obtain that $k^{-d-|s|}|\widehat{\mu}(\mathbf{k})|^{2}$ is integrable at the origin, and therefore we obtain

[TABLE]

where $c_{d,s}=-(2\pi)^{\frac{d}{2}}2^{\frac{d}{2}-1-s}\Gamma\left(\frac{d-s}{2}\right)\Gamma\left(\frac{s}{2}\right)^{-1}$ is positive for $-2<s<0$ (but negative for $-4<s<-2$). At $s=-2$ we can compute directly $\mathcal{E}_{{\rm Jel},d,-2}(\Omega,\mathbf{x}_{1},...,\mathbf{x}_{N})=|P|^{2}$ and we conclude, as we have claimed, that $\mathcal{E}_{{\rm Jel},d,s}(\Omega,\mathbf{x}_{1},...,\mathbf{x}_{N})\geqslant 0$ for all $-2\leqslant s<0$ when $|\Omega|=N$.

Acknowledgments. 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.). Part of this work was done when the authors were visiting the Institut Henri Poincaré in Paris and the Erwin Schrödinger Institute in Vienna.