Challenging the Lieb–Oxford bound in a systematic way

ABSTRACT

The Lieb–Oxford bound, a nontrivial inequality for the indirect part of the many-body Coulomb repulsion in an electronic system, plays an important role in the construction of approximations in density functional theory. Using the wave function for strictly correlated electrons of a given density, we turn the search over wave functions appearing in the original bound into a more manageable search over electron densities. This allows us to challenge the bound in a systematic way. We find that a maximising density for the bound, if it exists, must have compact support. We also find that, at least for particle numbers N ≤ 60, a uniform density profile is not the most challenging for the bound. With our construction, we improve the bound for N = 2 electrons that was originally found by Lieb and Oxford, we give a new lower bound to the constant appearing in the Lieb–Oxford inequality valid for any N, and we provide an improved upper bound for the low-density uniform electron gas indirect energy.

GRAPHICAL ABSTRACT

KEYWORDS:

• Density functional theory
• Lieb–Oxford bound
• exchange-correlation functionals
• electronic correlation

1. Introduction

Lieb and Oxford (LO) [1,2] proved a nontrivial inequality for the indirect part of the electron–electron interaction energy (total expectation of the interaction minus the Hartree term) with respect to the local-density approximation (LDA) exchange functional. This inequality has been recently extended to include the gradient of the density [3]. The LO bound has played and continues to play a very important role in the construction of approximate exchange-correlation (xc) density functionals [4–13]. While traditionally only the more general LO bound, valid for any number of particles N (and corresponding to N → ∞), has been taken into account in the construction of xc approximations, it has been shown very recently that the bound for N = 1 and N = 2 is important in the context of metaGGA functionals [11,12], and can be imposed as an additional exact condition.

The bound for N = 1 was first given in [14], and proved rigorously in [2]. For N = 2, Lieb and Oxford [2] could only provide a non-optimal estimate of the constant appearing in the bound. In this work, we develop a strategy to systematically challenge the original LO bound for a given number of electrons N. We use optimal trial wave functions for a given density, and we then vary the density in order to challenge the bound as much as possible. After showing that a density that maximally challenges the bound, if it exists, must have compact support, we follow the functional derivative of the bound to challenge it as much as possible without violating N-representability also for densities whose support is the whole space. As the first application of this procedure, we improve the lower bound for N = 2 given by Lieb and Oxford (see Equation (60(60) $${\bar{\lambda }}_3\left(2\right)>1.701052.$$(60) )). Our construction also provides an improved lower bound for the constant appearing in the Lieb–Oxford inequality valid for any N (see Equations (63(63) $$1.91175\le {\bar{\lambda }}_3\le 2.215,$$(63) ) and (64(64) $$1.4119\le C_3\le 1.6358.$$(64) )), and an improved upper bound for the indirect energy of the low-density uniform electron gas (see Equation (61(61) $$-1.45\le w\le -1.343,$$(61) )).

1.1. Notation

In electronic density functional theory (DFT), one is interested in finding the ground-state energy and density of N-electron systems with Hamiltonian (1) $$\widehat{H}=\widehat{T}+{\widehat{V}}_{\mathrm{ee}}+{\widehat{V}}_{\mathrm{ext}},{\widehat{V}}_{\mathrm{ext}}=\sum _{i=1}^{N}v\left({\mathbf{r}}_i\right).$$(1) $\widehat{T}$ and ${\widehat{V}}_{\mathrm{ee}}$ are, respectively, the universal operators of the kinetic energy (in Hartree atomic units used throughout the paper), (2) $$\widehat{T}=-\frac{1}{2}\sum _{i=1}^N\frac{{\partial }^2}{\partial {\mathbf{r}}_i^2},$$(2) and of the interaction (Coulomb repulsion) energy between the N electrons, (3) $${\widehat{V}}_{\mathrm{ee}}=\frac{1}{2}\sum _{i,j=1}^N\frac{1-{\delta }_{ij}}{|{\mathbf{r}}_i-{\mathbf{r}}_j|}.$$(3) The function v(r), in contrast, is a non-universal but arbitrary attractive external potential required to bind the repulsive electrons. Most of the formalism will be carried out for general spatial dimension D = 2 and 3, $\mathbf{r}\in {\mathbb{R}}^D$, focussing later on D = 3 only.

In the following, Ψ denotes a correctly normalised and antisymmetrised, but otherwise arbitrary N-electron wave function (thus, not necessarily eigenstate of (1(1) $$\widehat{H}=\widehat{T}+{\widehat{V}}_{\mathrm{ee}}+{\widehat{V}}_{\mathrm{ext}},{\widehat{V}}_{\mathrm{ext}}=\sum _{i=1}^{N}v\left({\mathbf{r}}_i\right).$$(1) )), (4) $$\Psi =\Psi ({\mathbf{r}}_1{\sigma }_1,...,{\mathbf{r}}_N{\sigma }_N),$$(4) where σ_n are spin variables. By ${\rho }_{\Psi }$, we denote the particle density associated with Ψ, (5) $$\begin{array}{c}\hfill {\rho }_{\Psi }\left(\mathbf{r}\right)=N\sum _{{\sigma }_i}\int d^D{r}_2\cdots d^D{r}_N|\Psi (\mathbf{r}{\sigma }_1,{\mathbf{r}}_2{\sigma }_2,...,{\mathbf{r}}_N{\sigma }_N){|}^2.\end{array}$$(5)

1.2. Indirect Coulomb energy

The electronic interaction energy in the quantum state Ψ, defined as the expectation (6) $$⟨\Psi \left|{\widehat{V}}_{\mathrm{ee}}\right|\Psi ⟩>0,$$(6) excludes the infinite self-energies of the point electrons (see the factor 1 − δ_ij in Equation (3(3) $${\widehat{V}}_{\mathrm{ee}}=\frac{1}{2}\sum _{i,j=1}^N\frac{1-{\delta }_{ij}}{|{\mathbf{r}}_i-{\mathbf{r}}_j|}.$$(3) )). If the electrons were a classical continuous distribution of negative charge with density ${\rho }_{\Psi }\left(\mathbf{r}\right)$, their interaction energy would be $U\left[{\rho }_{\Psi }\right]$, with the Hartree functional (7) $$U\left[\rho \right]=\frac{1}{2}\int d^{D}r\int d^D{r}^{\text{'}}\frac{\rho \left(\mathbf{r}\right)\rho \left({\mathbf{r}}^{\text{'}}\right)}{|\mathbf{r}-{\mathbf{r}}^{\text{'}}|}>0.$$(7) The most severe error introduced by this classical continuum approximation is a spurious finite self-interaction energy included for each electron. This is particularly evident in the case N = 1, since for any normalised one-electron wave function Ψ, we have $⟨\Psi \left|{\widehat{V}}_{\mathrm{ee}}\right|\Psi ⟩=0$, while $U\left[{\rho }_{\Psi }\right]>0$. The indirect interaction energy W[Ψ] is defined as (8) $$W\left[\Psi \right]\equiv ⟨\Psi |{\widehat{V}}_{\mathrm{ee}}|\Psi ⟩-U\left[{\rho }_{\Psi }\right].$$(8) For wave functions that are ground states of an N-electron Hamiltonian (1(1) $$\widehat{H}=\widehat{T}+{\widehat{V}}_{\mathrm{ee}}+{\widehat{V}}_{\mathrm{ext}},{\widehat{V}}_{\mathrm{ext}}=\sum _{i=1}^{N}v\left({\mathbf{r}}_i\right).$$(1) ) (or good trial wave function for it), W[Ψ] is normally negative. However, for a given density ρ, it is possible to construct wave functions Ψ for which W[Ψ] is positive or even infinity [2,15]. We emphasise that U[ρ] is a density functional, while W[Ψ] is a functional in terms of the wave function Ψ.

1.3. Lieb–Oxford bound

The quantity W[Ψ] is limited by the Lieb–Oxford (LO) bound, (9) $$-C_D\int d^{D}r{\rho }_{\Psi }{\left(\mathbf{r}\right)}^{1+1/D}\le W\left[\Psi \right].$$(9) C_D > 0 is the unknown minimum possible number that makes this inequality true for all wave functions Ψ in D = 2 or 3 dimensions. So far, it is rigorously known that C₃ ≤ 1.6358 [16] and C₂ ≤ 481.28 [17], and it has been argued [18], on physical arguments, that the two bounds can be tightened to C₃ ≤ 1.44 and C₂ ≤ 1.96. The assumption behind these latter conjectured values is that the tightest possible bound is provided by the indirect energy of the uniform electron gas in the low-density limit, which, in turn, is commonly identified with the Wigner crystal total energy. This latter assumption has recently been proven wrong for the 3D case by Lewin and Lieb [3]. The study presented in this paper will also raise doubts on the first assumption that a uniform density is really the most challenging case for the LO bound, after a suitable optimal wave function for each given density has been defined (see Section 4).

In terms of the LDA (10) $$-A_D\int d^{D}r\rho {\left(\mathbf{r}\right)}^{1+1/D}=E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]$$(10) to the D-dimensional exchange energy, with the exact constants $A_3=\frac{3}{4}{\left(\frac{3}{\pi }\right)}^{1/3}\approx 0.739$, $A_2=\frac{4}{3}{\left(\frac{2}{\pi }\right)}^{1/2}\approx 0.798$, Equation (9(9) $$-C_D\int d^{D}r{\rho }_{\Psi }{\left(\mathbf{r}\right)}^{1+1/D}\le W\left[\Psi \right].$$(9) ) reads (11) $$\lambda \left[\Psi \right]\le {\bar{\lambda }}_D,$$(11) where we have defined (12) $$\lambda \left[\Psi \right]\equiv \frac{W\left[\Psi \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[{\rho }_{\Psi }\right]},{\bar{\lambda }}_D\equiv \frac{C_D}{A_D}.$$(12) Considering all antisymmetric wave functions Ψ in D dimensions, we may write (13) $${\bar{\lambda }}_D=\underset{\Psi :D}{sup}\lambda \left[\Psi \right].$$(13) The above rigorous upper bounds for C_D correspond to (14) $${\bar{\lambda }}_2\le 603,{\bar{\lambda }}_3\le 2.215.$$(14) Considering wave functions Ψ → N with a given particle number N, we define (15) $${\bar{\lambda }}_D\left(N\right)=\underset{(\Psi :D)\to N}{sup}\lambda \left[\Psi \right].$$(15) Lieb and Oxford [2] have proven that ${\bar{\lambda }}_3\left(N\right)$ is monotonically increasing with its integer variable N, (16) $${\bar{\lambda }}_3\left(N\right)<{\bar{\lambda }}_3(N+1),\underset{N\to \infty }{lim}{\bar{\lambda }}_3\left(N\right)={\bar{\lambda }}_3.$$(16) They have also proven that ${\bar{\lambda }}_3\left(1\right)=1.4786$ (which was given originally by Gadre et al. [14]) and they have found a lower bound for ${\bar{\lambda }}_3\left(2\right)$, (17) $${\bar{\lambda }}_3\left(2\right)>1.67.$$(17) These bounds in D = 3 for N = 1 and N = 2 have been recently used to improve a certain class of exchange-correlation functionals [11,12].

In this paper, we develop a general strategy to find improved lower bounds for ${\bar{\lambda }}_D\left(N\right)$ by challenging the Lieb–Oxford bound, i.e. by evaluating λ[Ψ] with particularly efficient trial wave functions Ψ. Note that this is different from what is usually called tightening the bound, which means finding improved upper bounds to ${\bar{\lambda }}_D\left(N\right)$.

A new lower bound for ${\bar{\lambda }}_D\left(N\right)$ (or, generally, for ${\bar{\lambda }}_D$) is rigorously obtained each time we find a wave function that gives the highest value ever observed for λ[Ψ] (for a given N, or in general). Until very recently, it was believed that a lower bound for ${\bar{\lambda }}_3$ is given by ${\bar{\lambda }}_3\ge 1.444/A_3=1.955$, corresponding to the total energy of the bcc Wigner crystal in the classical jellium model. However, in the jellium model, one can only identify the total energy with the indirect energy if the electronic density is uniform, exactly equal to the one of the positive background. Only in this case, the electronic Hartree term will be exactly cancelled by the electron–background and the background–background contributions to the total energy. Lewin and Lieb [3] have shown that, in the 3D case trying to make this cancellation happen by taking a superposition of all the possible Wigner lattices to have a uniform electronic density, introduces a shift that does not disappear in the thermodynamic limit. Thus, the value 1.955 does not correspond to the indirect energy of any wave function and is not a valid lower bound for ${\bar{\lambda }}_3$. In Section 4, we report a new lower bound for general N, by considering an optimal trial wave function for N = 60, and we also report an improved upper bound to the indirect energy of the low-density uniform gas.

2. The density functional Λ[ρ]

Considering only those wave functions Ψ → ρ (in D dimensions) that are associated with a given particle density ρ = ρ(r), we define the density functional (18) $$\Lambda \left[\rho \right]\equiv \underset{\Psi \to \rho }{max}\lambda \left[\Psi \right].$$(18) Writing $N_{\Psi }=\int d^{D}r{\rho }_{\Psi }\left(\mathbf{r}\right)$ for the electron number in the state Ψ, we then have (19) $$\lambda \left[\Psi \right]\le \Lambda \left[{\rho }_{\Psi }\right]<{\bar{\lambda }}_D\left(N_{\Psi }\right)<{\bar{\lambda }}_D.$$(19)

2.1. SCE interaction energy

More explicitly, (20) $$\begin{array}{c}\hfill \Lambda \left[\rho \right]\equiv \underset{\Psi \to \rho }{max}\frac{⟨\Psi \left|{\widehat{V}}_{\mathrm{ee}}\right|\Psi ⟩-U\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]}=\frac{V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]-U\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]},\end{array}$$(20) with the SCE interaction energy of Appendix A, (21) $$V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]=\underset{\Psi \to \rho }{min}⟨\Psi |{\widehat{V}}_{\mathrm{ee}}|\Psi ⟩.$$(21) The acronym SCE [19–21] stands for ‘strictly correlated electrons’ and defines a state |Ψ_SCE[ρ]|², which is a very accurate trial wave function (actually a distribution) for the maximising one in Equation (20(20) $$\begin{array}{c}\hfill \Lambda \left[\rho \right]\equiv \underset{\Psi \to \rho }{max}\frac{⟨\Psi \left|{\widehat{V}}_{\mathrm{ee}}\right|\Psi ⟩-U\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]}=\frac{V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]-U\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]},\end{array}$$(20) ), being exact in 1D [22] for any N, and in any dimension for N = 2 [23]. The SCE state is detailed in Appendix A. In other words, Out of all antisymmetric wave functions Ψ that are associated with a given density ρ, the one that provides (or is very close to) the strongest challenge to the Lieb–Oxford bound is the SCE state |Ψ_SCE[ρ]|². Consequently, since V^SCE_ee[ρ] can be evaluated rigorously for a wide class of densities (Equation (A5(A5) $$V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]=\int {\mathrm{d}}^{D}s\frac{\rho \left(\mathbf{s}\right)}{N}\sum _{i=1}^{N-1}\sum _{j=i+1}^N\frac{1}{|{\mathbf{f}}_i\left(\mathbf{s}\right)-{\mathbf{f}}_j\left(\mathbf{s}\right)|}.$$(A5) ) in Appendix A), we no longer need to consider different trial wave functions Ψ, but only different trial densities ρ instead, (22) $${\bar{\lambda }}_D=\underset{\rho :D}{sup}\Lambda \left[\rho \right],{\bar{\lambda }}_D\left(N\right)=\underset{(\rho :D)\to N}{sup}\Lambda \left[\rho \right].$$(22)

As a preliminary step, we have used simple analytical trial spherical densities to evaluate Λ[ρ] for D = 3 and N = 2, reporting the results in Table 1. We see that the lower bound (17(17) $${\bar{\lambda }}_3\left(2\right)>1.67.$$(17) ) is readily improved to (23) $${\bar{\lambda }}_3\left(2\right)>1.70097.$$(23) There is no need for considering scaled densities ρ_ξ(r) ≡ ξ^Dρ(ξr), with various values of ξ > 0, since Λ[ρ_ξ] = Λ[ρ] (see Equation (A14(A14) $$V_{\mathrm{ee}}^{\mathrm{SCE}}\left[{\rho }_{\xi }\right]=\xi V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right].$$(A14) ) in Appendix A).

It is interesting to note that, once the most challenging wave function for each given ρ(r) is used, the densities that give the highest values of Λ[ρ] are quite surprising. For example, a density proportional to $e^{-50{(r-1)}^2}$, consisting of a thin spherical shell, is similar to the one of the strongly correlated limit of the Hooke’s atom series. Yet, it gives a value of Λ[ρ] which is much lower than the one obtained from the exponential density. Indeed, the strong correlation limit of the Hooke’s series is known to give λ[ρ] = 1.489 [8], again much less than what we obtain for exponential-like densities. The point is that previous works which analysed numerically the LO bound [5–8] focused on physical Hamiltonians of the kind (1(1) $$\widehat{H}=\widehat{T}+{\widehat{V}}_{\mathrm{ee}}+{\widehat{V}}_{\mathrm{ext}},{\widehat{V}}_{\mathrm{ext}}=\sum _{i=1}^{N}v\left({\mathbf{r}}_i\right).$$(1) ), choosing v(r) that could be particularly challenging for the bound. In that context, exponential-like densities would correspond to the large nuclear-charge limit of the He isoelectronic series, which is a weakly correlated system. With our construction, instead, we use the most challenging wave function for any given density, finding the unexpected trends of Table 1. We also see that the density of a uniform sphere (droplet) is not particularly challenging for the bound, a feature that will be further analysed in Section IV for larger N.

Table 1. Values Λ[ρ] for some simple spherical two-electron trial densities ρ(r) in three dimensions (N = 2, D = 3), obtained numerically from Equations (A9(A9) $$N_e\left(s\right)\equiv {\int }_{\left|\mathbf{r}\right|\le s}{\mathrm{d}}^{D}r\rho \left(r\right),$$(A9) )–(A11(A11) $$V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]=\frac{1}{2}\int {\mathrm{d}}^{D}r\frac{\rho \left(r\right)}{r+f\left(r\right)}.$$(A11) ) of Appendix A. In the last two rows, we consider densities with compact support: ‘droplet’ corresponds to the case of a sphere of uniform density [24], and the density proportional to r⁻³ [11] has been evaluated for R₁ = 10³ and R₂ = 10⁵. [Atomic units are used, where r is a dimensionless radial coordinate.]

ρ(r)∝	Λ[ρ]	ρ(r)∝	Λ[ρ]
$e^{-10{(r-1)}^2}$	1.499	$e^{-50{(r-1)}^2}$	1.262
(1 + r)^−4	1.562	$e^{-r^2}$	1.689
(1 + r)^−5	1.637	e^−r	1.69905
(1 + r)^−6	1.662	r e^−r	1.69866
(1 + r)^−7	1.674	r^1/2e^−r	1.70097
(1 + r)^−10	1.687	r^1/3e^−r	1.70095
1 − r,  r ≤ 1	1.638	r^−3,  r ∈ [R1, R2]	1.145
Droplet	1.498	$cos\left(r\right),r\le \frac{\pi }{2}$	1.627

2.2. Absence of a maximising density without compact support

We now demonstrate that a function ρ(r) that maximises the functional Λ[ρ] for a finite N cannot be a physical density, unless it has compact support. The argument is essentially the same used by Lieb and Oxford [2] for N = 1 and N = 2. In terms of the SCE external potential of Appendix A, (24) $$v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)\equiv \frac{\delta V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)},$$(24) and the Hartree potential (25) $$v_{\mathrm{H}}\left[\rho \right]\left(\mathbf{r}\right)\equiv \frac{\delta U\left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)}=\int {\mathrm{d}}^D{r}^{\text{'}}\frac{\rho \left({\mathbf{r}}^{\text{'}}\right)}{|\mathbf{r}-{\mathbf{r}}^{\text{'}}|},$$(25) we consider the Euler equation for maximising Λ[ρ]. By writing ρ(r) = p(r)² to ensure ρ(r) ≥ 0, and by varying p(r), we obtain (26) $$\begin{array}{ccc}& & \left\{\frac{v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)-v_{\mathrm{H}}\left[\rho \right]\left(\mathbf{r}\right)}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]}-\frac{V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]-U\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}{\left[\rho \right]}^2}{v}_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]\left(\mathbf{r}\right)\right\}p\left(\mathbf{r}\right)\hfill \\ & & =\mu p\left(\mathbf{r}\right).\hfill \end{array}$$(26) If p(r) ≠ 0 everywhere, we obtain the Euler equation (27) $$\begin{array}{ccc}\hfill \frac{\delta \Lambda \left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)}& \equiv & \frac{v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)-v_{\mathrm{H}}\left[\rho \right]\left(\mathbf{r}\right)}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]}\hfill \\ & & -\frac{V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]-U\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}{\left[\rho \right]}^2}{v}_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]\left(\mathbf{r}\right)=\mu ,\hfill \end{array}$$(27) where μ is the Lagrange multiplier ensuring fixed particle number N = ∫d^Dr ρ(r), and (28) $$\begin{array}{c}\hfill v_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]\left(\mathbf{r}\right)\equiv \frac{\delta E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)}=-A_D\left(1+\frac{1}{D}\right)\rho {\left(\mathbf{r}\right)}^{1/D}.\end{array}$$(28) In this case, since $v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)\to \frac{N-1}{r}$ and $v_{\mathrm{H}}\left[\rho \right]\left(\mathbf{r}\right)\to \frac{N}{r}$ for r ≡ |r| → ∞, we have asymptotically (29) $$v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)-v_{\mathrm{H}}\left[\rho \right]\left(\mathbf{r}\right)\to -\frac{1}{r}(r\to \infty ).$$(29) Comparing this with Equation (28(28) $$\begin{array}{c}\hfill v_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]\left(\mathbf{r}\right)\equiv \frac{\delta E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)}=-A_D\left(1+\frac{1}{D}\right)\rho {\left(\mathbf{r}\right)}^{1/D}.\end{array}$$(28) ), we see that a solution ρ(r) of Equation (27(27) $$\begin{array}{ccc}\hfill \frac{\delta \Lambda \left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)}& \equiv & \frac{v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)-v_{\mathrm{H}}\left[\rho \right]\left(\mathbf{r}\right)}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]}\hfill \\ & & -\frac{V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]-U\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}{\left[\rho \right]}^2}{v}_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]\left(\mathbf{r}\right)=\mu ,\hfill \end{array}$$(27) ) must display the asymptotic behaviour (30) $$\rho \left(\mathbf{r}\right)\to \frac{k_1}{r^D}(r\to \infty ),$$(30) with some constant k₁. Such a function is evidently not normalisable, since with the D-dimensional volume element d^Dr = k₂r^{D − 1}dr and a radius R > 0 finite but large enough, we have (31) $${\int }_{\left|\mathbf{r}\right|\ge R}{\mathrm{d}}^{D}r\rho \left(\mathbf{r}\right)={\int }_R^{\infty }\mathrm{d}r\frac{k_1{k}_2}{r}=\infty .$$(31) We emphasise that this reasoning also applies to the modified functional (32) $$\tilde{\Lambda }\left[\rho \right]=\frac{E_{\mathrm{xc}}\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]},$$(32) where the indirect SCE interaction energy V^SCE_ee[ρ] − U[ρ] is replaced with the functional E_xc[ρ] of the exchange-correlation energy, since the xc potential for N-electron systems has the same asymptotic behaviour as Equation (29(29) $$v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)-v_{\mathrm{H}}\left[\rho \right]\left(\mathbf{r}\right)\to -\frac{1}{r}(r\to \infty ).$$(29) ), (33) $$v_{\mathrm{xc}}\left[\rho \right]\left(\mathbf{r}\right)\equiv \frac{\delta E_{\mathrm{xc}}\left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)}\to -\frac{1}{r}(r\to \infty ).$$(33)

Quite interestingly, a density of exactly the same form of Equation (30(30) $$\rho \left(\mathbf{r}\right)\to \frac{k_1}{r^D}(r\to \infty ),$$(30) ) for the 3D case, but restricted to a finite region of space (thus set to zero outside some region r ∈ [R₁, R₂]), has been considered by Perdew et al. [11] to study a general feature of GGA approximations related to the LO bound. Note, however, that if we consider this kind of densities, ρ(r)∝r⁻³ in r ∈ [R₁, R₂], even by choosing R₁ and R₂ very large, we get quite low values for Λ[ρ], indicating that the asymptotic condition is anyway not enough to give a large Λ value (see Table 1).

Even more generally, in a fictitious universe where the electron–electron repulsion is multiplied by a factor α ≥ 0, the density functional of their xc energy is given by (34) $$E_{\mathrm{xc},\alpha }\left[\rho \right]={\int }_0^{\alpha }\mathrm{d}\beta \left\{⟨{\Psi }_{\beta }\left[\rho \right]\left|{\widehat{V}}_{\mathrm{ee}}\right|{\Psi }_{\beta }\left[\rho \right]⟩-U\left[\rho \right]\right\}.$$(34) Here, out of all antisymmetric wave functions Ψ that are associated with the same density ρ, Ψ_β[ρ] is the one that minimises the expectation $⟨\Psi |\widehat{T}+\beta {\widehat{V}}_{\mathrm{ee}}|\Psi ⟩$, for any number β ≥ 0. Since the corresponding α-dependent xc potential has the asymptotic behaviour (35) $$v_{\mathrm{xc},\alpha }\left[\rho \right]\left(\mathbf{r}\right)\equiv \frac{\delta E_{\mathrm{xc},\alpha }\left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)}\to -\frac{\alpha }{r}(r\to \infty ),$$(35) we conclude that even for the functional (36) $${\Lambda }_{\alpha }\left[\rho \right]=\frac{\frac{1}{\alpha }{E}_{\mathrm{xc},\alpha }\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]},$$(36) the maximising function ρ(r) must have compact support. Note that ${\Lambda }_{\alpha =1}\left[\rho \right]=\tilde{\Lambda }\left[\rho \right]$ and lim_{α → ∞}Λ_α[ρ] = Λ[ρ].

If p(r) = 0 for |r| ≥ r₀, we see that, in principle, a maximising density in Equation (26(26) $$\begin{array}{ccc}& & \left\{\frac{v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)-v_{\mathrm{H}}\left[\rho \right]\left(\mathbf{r}\right)}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]}-\frac{V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]-U\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}{\left[\rho \right]}^2}{v}_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]\left(\mathbf{r}\right)\right\}p\left(\mathbf{r}\right)\hfill \\ & & =\mu p\left(\mathbf{r}\right).\hfill \end{array}$$(26) ) could exist. However, with our numerical investigation, we have always found larger values of Λ for densities with unbounded support.

3. Following the functional gradient of Λ[ρ]

Although Λ[ρ] has no maximising density ρ without compact support, the functional gradient δΛ/δρ tells us how to increase the value Λ[ρ] (or challenge the Lieb–Oxford bound) systematically. Starting from an N-electron density ρ = ρ(r) with a high value Λ[ρ], we consider a small density variation, (37) $$\rho \left(\mathbf{r}\right)\to \rho \left(\mathbf{r}\right)+ϵ\sigma \left(\mathbf{r}\right),\int d^{3}r\sigma \left(\mathbf{r}\right)=0.$$(37) Provided that εσ(r) is truly ‘small’, which precisely means that (38) $$\int d^{3}r\sigma {\left(\mathbf{r}\right)}^2=1$$(38) and |ε| ≪ 1, we have (39) $$\Lambda [\rho +ϵ\sigma ]-\Lambda \left[\rho \right]\approx ϵ\int d^{3}rG\left[\rho \right]\left(\mathbf{r}\right)\sigma \left(\mathbf{r}\right),$$(39) with the gradient G[ρ](r) ≡ δΛ[ρ]/δρ(r) given by Equation (27(27) $$\begin{array}{ccc}\hfill \frac{\delta \Lambda \left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)}& \equiv & \frac{v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)-v_{\mathrm{H}}\left[\rho \right]\left(\mathbf{r}\right)}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]}\hfill \\ & & -\frac{V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]-U\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}{\left[\rho \right]}^2}{v}_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]\left(\mathbf{r}\right)=\mu ,\hfill \end{array}$$(27) ). Although ∫d³r σ(r) = 0, the right-hand side of Equation (39(39) $$\Lambda [\rho +ϵ\sigma ]-\Lambda \left[\rho \right]\approx ϵ\int d^{3}rG\left[\rho \right]\left(\mathbf{r}\right)\sigma \left(\mathbf{r}\right),$$(39) ) can nevertheless be >0, provided that G[ρ](r), as a function of r, is different from a constant, G[ρ](r) ≠ const.

3.1. Formal optimisation of the increment

Formally, maximising the integral ∫d³r G[ρ](r)σ(r) with respect to σ(r) subject to the two constraints ∫d³r σ(r) = 0 and ∫d³r σ(r)² = 1, (40) $$\begin{array}{ccc}& & \frac{\delta }{\delta \sigma \left(\mathbf{r}\right)}\{\int {\mathrm{d}}^{3}rG\left[\rho \right]\left(\mathbf{r}\right)\sigma \left(\mathbf{r}\right)-{\mu }_1\int {\mathrm{d}}^{3}r\sigma \left(\mathbf{r}\right)\hfill \\ & & -{\mu }_2\int {\mathrm{d}}^{3}r\sigma {\left(\mathbf{r}\right)}^2\}=0,\hfill \end{array}$$(40) yields the Euler Equation G[ρ](r) − μ₁ − 2μ₂σ(r) = 0, with the solution (41) $${\sigma }_0\left(\mathbf{r}\right)=\frac{G\left[\rho \right]\left(\mathbf{r}\right)-{\mu }_1}{2{\mu }_2}.$$(41) The first Lagrange multiplier μ₁ is fixed by the normalisation constraint ∫d³r σ₀(r) = 0. The second one μ₂ is absorbed in the small parameter ε, guaranteeing the validity of the approximation (39(39) $$\Lambda [\rho +ϵ\sigma ]-\Lambda \left[\rho \right]\approx ϵ\int d^{3}rG\left[\rho \right]\left(\mathbf{r}\right)\sigma \left(\mathbf{r}\right),$$(39) ). [Independently, smallness of εσ(r) is necessary (but not sufficient) for the resulting density to be non-negative, ρ(r) + εσ(r) ≥ 0 for all r.]

For a N-electron (finite) density ρ, Equations (27(27) $$\begin{array}{ccc}\hfill \frac{\delta \Lambda \left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)}& \equiv & \frac{v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)-v_{\mathrm{H}}\left[\rho \right]\left(\mathbf{r}\right)}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]}\hfill \\ & & -\frac{V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]-U\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}{\left[\rho \right]}^2}{v}_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]\left(\mathbf{r}\right)=\mu ,\hfill \end{array}$$(27) )–(29(29) $$v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)-v_{\mathrm{H}}\left[\rho \right]\left(\mathbf{r}\right)\to -\frac{1}{r}(r\to \infty ).$$(29) ) imply the large-r behaviour (r → ∞) (42) $$\begin{array}{c}\hfill G\left[\rho \right]\left(\mathbf{r}\right)\to \frac{1}{|E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]|}\left[\frac{1}{r}-\Lambda \left[\rho \right]A_D\left(1+\frac{1}{D}\right)\rho {\left(\mathbf{r}\right)}^{1/D}\right].\end{array}$$(42) Necessarily, σ₀(r) → 0 for r → ∞, implying μ₁ = 0 in Equation (41(41) $${\sigma }_0\left(\mathbf{r}\right)=\frac{G\left[\rho \right]\left(\mathbf{r}\right)-{\mu }_1}{2{\mu }_2}.$$(41) ). Consequently, due to the term 1/r in Equation (42(42) $$\begin{array}{c}\hfill G\left[\rho \right]\left(\mathbf{r}\right)\to \frac{1}{|E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]|}\left[\frac{1}{r}-\Lambda \left[\rho \right]A_D\left(1+\frac{1}{D}\right)\rho {\left(\mathbf{r}\right)}^{1/D}\right].\end{array}$$(42) ), ∫d³r σ₀(r) cannot be zero (or even finite). In other words, ρ(r) + εσ₀(r), with ε ≠ 0, must, again, yield a density with compact support. In the following, we give an analytical example for the case of a density with compact support.

3.2. Analytical example for densities with compact support

As an example, we evaluate Equation (41(41) $${\sigma }_0\left(\mathbf{r}\right)=\frac{G\left[\rho \right]\left(\mathbf{r}\right)-{\mu }_1}{2{\mu }_2}.$$(41) ) for the spherical two-electron density [24] (43) $$\begin{array}{c}\hfill \rho \left(\mathbf{r}\right)=\left\{\begin{array}{cc}{\rho }_0& (r\le R),\\ 0& (r& >R)\end{array}\right.{\rho }_0=\frac{3}{2\pi R^3}.\end{array}$$(43) This density corresponds to a uniformly charged sphere (droplet) with radius R and total charge 2, (44) $$\begin{array}{c}\hfill U\left[\rho \right]=\frac{12}{5R},v_{\mathrm{H}}\left[\rho \right]\left(r\right)=\frac{3R^2-r^2}{R^3}(r\le R).\end{array}$$(44) The exact exchange energy E_x[ρ] is given by $-E_x\left[\rho \right]=\frac{1}{2}U\left[\rho \right]=\frac{1.2}{R}$, while [24] (45) $$\begin{array}{c}\hfill -E_x^{\mathrm{LDA}}\left[\rho \right]=\frac{1.1545}{R},\end{array}$$(45) (46) $$\begin{array}{c}\hfill v_x^{\mathrm{LDA}}\left[\rho \right]\left(r\right)=-{\left(\frac{9}{2{\pi }^2}\right)}^{1/3}\frac{1}{R}(r\le R).\end{array}$$(46) From [24], we have Λ[ρ] = 1.498 and the SCE co-motion function (see Appendix A) is (47) $$f\left(r\right)=R{\left(1-\frac{r^3}{R^3}\right)}^{1/3}.$$(47) The resulting SCE external potential is given by (48) $$\begin{array}{ccc}& & v^{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)=v^{\mathrm{SCE}}\left[\rho \right]\left(0\right)-{\int }_0^r\frac{du}{{[u+f\left(u\right)]}^2}\hfill \\ & & =v^{\mathrm{SCE}}\left[\rho \right]\left(0\right)-\frac{1}{R}{\int }_0^{r/R}\frac{dx}{{[x+{(1-x^3)}^{1/3}]}^2}.\hfill \end{array}$$(48) Eventually, Equation (41(41) $${\sigma }_0\left(\mathbf{r}\right)=\frac{G\left[\rho \right]\left(\mathbf{r}\right)-{\mu }_1}{2{\mu }_2}.$$(41) ) reads (49) $$\sigma \left(r\right)=\frac{-v^{\mathrm{SCE}}\left[\rho \right]\left(r\right)+v_{\mathrm{H}}\left[\rho \right]\left(r\right)-{\tilde{\mu }}_1}{2{\tilde{\mu }}_2},$$(49) where ${\tilde{\mu }}_1={\mu }_1-\Lambda \left[\rho \right]v_x^{\mathrm{LDA}}$ (note that v^LDA_x does not depend on r in the present example) and $2{\tilde{\mu }}_2=-E_x^{\mathrm{LDA}}\left[\rho \right]2{\mu }_2>0$. The constant v^SCE[ρ](0) can also be absorbed by the multiplier ${\tilde{\mu }}_1$ which is fixed by the condition ∫^R₀dr 4πr²σ(r) = 0. Then, we have (50) $$\begin{array}{c}\hfill \sigma \left(r\right)=\frac{1}{2{\tilde{\mu }}_{2}R}\left\{3-\frac{r^2}{R^2}-{\tilde{\mu }}_1+{\int }_0^{r/R}\frac{dx}{{[x+{(1-x^3)}^{1/3}]}^2}\right\}.\end{array}$$(50) A simple but accurate approximation to this function (for R = 1) is (51) $$\begin{array}{c}\hfill {\sigma }_{\mathrm{appr}}\left(r\right)=\frac{1}{2{\tilde{\mu }}_2}\left[0.4r^3-1.85r^2+r+0.16\right]r\le 1.\end{array}$$(51) We therefore consider the densities (for r ≤ 1) (52) $$\begin{array}{c}\hfill {\rho }_a\left(r\right)={\rho }_0+a\left[0.4r^3-1.85r^2+r+0.16\right](a\ge 0)\end{array}$$(52) to obtain the values Λ[ρ_0.2] = 1.521, Λ[ρ_0.5] = 1.551, Λ[ρ_1.0] = 1.590, Λ[ρ_1.5] = 1.611, Λ[ρ_1.6] = 1.612. [For a > 1.6, the density ρ_a(r) becomes negative.]

3.3. Compromise for N-representability for densities with unbounded support

Since we observe higher values of Λ[ρ] for densities for which p(r) ≠ 0 everywhere, we consider here a compromise to follow the gradient of Λ[ρ] without violating N-representability. We perturb the density with some function σ(r) in Equation (37(37) $$\rho \left(\mathbf{r}\right)\to \rho \left(\mathbf{r}\right)+ϵ\sigma \left(\mathbf{r}\right),\int d^{3}r\sigma \left(\mathbf{r}\right)=0.$$(37) ) that depends on a certain number of parameters and satisfies Equations (37(37) $$\rho \left(\mathbf{r}\right)\to \rho \left(\mathbf{r}\right)+ϵ\sigma \left(\mathbf{r}\right),\int d^{3}r\sigma \left(\mathbf{r}\right)=0.$$(37) ) and (38(38) $$\int d^{3}r\sigma {\left(\mathbf{r}\right)}^2=1$$(38) ), keeping the perturbed density N-representable with suitable constraints. One can then choose the parameter values in order to maximise the overlap with the gradient [the right-hand-side of Equation (39(39) $$\Lambda [\rho +ϵ\sigma ]-\Lambda \left[\rho \right]\approx ϵ\int d^{3}rG\left[\rho \right]\left(\mathbf{r}\right)\sigma \left(\mathbf{r}\right),$$(39) )].

As an example, we start from the 3D exponential two-electron density (D = 3,  N = 2) (53) $$\rho \left(\mathbf{r}\right)=\frac{{\mathrm{e}}^{-r}}{4\pi },$$(53) which already gives the high value Λ[ρ] = 1.69905 (see Table 1). We choose for σ(r) the parametrised form (54) $${\sigma }_a\left(\mathbf{r}\right)=\sqrt{\frac{3a^3}{\pi }}\left(1-\frac{ar}{3}\right){\mathrm{e}}^{-ar}\left(a>0\right),$$(54) which obeys the conditions ∫d³r σ_a(r) = 0 and ∫d³r σ_a(r)² = 1 for all values of the parameter a > 0, so that the function ρ_{a, ε}(r) = ρ(r) + εσ_a(r) is always correctly normalised. In addition, N-representability requires that ρ_{a, ε}(r) ≥ 0 for all r ≥ 0. For any value of a, this is fulfilled for ε_min(a) ≤ ε ≤ ε_max(a), where ε_min(a) ≤ 0 and ε_max(a) ≥ 0 are given by (55) $$\begin{array}{c}\hfill {ϵ}_{\min }\left(a\right)=\left\{\begin{array}{ccc}\hfill {\displaystyle \frac{3(a-1)}{4a\sqrt{3a^3\pi }}{\mathrm{e}}^{-(3-4a)/a}}& \hfill (0& \hfill {\displaystyle <a\le \frac{3}{4}),}\\ \hfill {\displaystyle -\frac{1}{4\sqrt{3a^3\pi }}}& \hfill {\displaystyle (a\ge \frac{3}{4}),}\end{array}\right.\end{array}$$(55) (56) $$\begin{array}{c}\hfill {ϵ}_{\max }\left(a\right)=\left\{\begin{array}{ccc}\hfill 0& \hfill (0& \hfill <a\le 1),\\ \hfill {\displaystyle \frac{3(a-1)}{4a\sqrt{3a^3\pi }}{\mathrm{e}}^{-(3-4a)/a}}& \hfill (a\ge 1).\end{array}\right.\end{array}$$(56)

Evaluating numerically the functional gradient G[ρ](r) ≡ δΛ[ρ]/δρ(r) of Equation (27(27) $$\begin{array}{ccc}\hfill \frac{\delta \Lambda \left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)}& \equiv & \frac{v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)-v_{\mathrm{H}}\left[\rho \right]\left(\mathbf{r}\right)}{E_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]}\hfill \\ & & -\frac{V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]-U\left[\rho \right]}{E_{\mathrm{x}}^{\mathrm{LDA}}{\left[\rho \right]}^2}{v}_{\mathrm{x}}^{\mathrm{LDA}}\left[\rho \right]\left(\mathbf{r}\right)=\mu ,\hfill \end{array}$$(27) ) for the density ρ(r) = ρ(r) of Equation (53(53) $$\rho \left(\mathbf{r}\right)=\frac{{\mathrm{e}}^{-r}}{4\pi },$$(53) ), we consider, as a function of a, the overlap integral (57) $$I\left(a\right)={\int }_0^{\infty }\mathrm{d}r\left(4\pi r^2\right)G\left[\rho \right]\left(r\right){\sigma }_a\left(r\right).$$(57) For any value of a > 0, the maximum possible value in Equation (39(39) $$\Lambda [\rho +ϵ\sigma ]-\Lambda \left[\rho \right]\approx ϵ\int d^{3}rG\left[\rho \right]\left(\mathbf{r}\right)\sigma \left(\mathbf{r}\right),$$(39) ) is approximately (if the first-order expansion holds) (58) $$\Lambda \left[\rho +ϵ\left(a\right){\sigma }_a\right]-\Lambda \left[\rho \right]\approx ϵ\left(a\right)I\left(a\right),$$(58) where ε(a) = ε_max(a) ≥ 0 for I(a) ≥ 0 and ε(a) = ε_min(a) ≤ 0 for I(a) ≤ 0.

Numerically, I(a) > 0 for 0 < a < 1 and I(a) < 0 for a > 1, with a strong maximum I(a₁) ≈ 2.9267 · 4π at a₁ ≈ 0.079 and a weak minimum I(a₂) ≈ −0.01479 · 4π = −0.2302 at a₂ ≈ 2.49. While ε(a₁) = 0, we have ε(a₂) = −0.0207, and Equation (58(58) $$\Lambda \left[\rho +ϵ\left(a\right){\sigma }_a\right]-\Lambda \left[\rho \right]\approx ϵ\left(a\right)I\left(a\right),$$(58) ) for a = a₂ gives (59) $$\Lambda \left[\rho +ϵ\left(a_2\right){\sigma }_{a_2}\right]-\Lambda \left[\rho \right]\approx 0.004765.$$(59) In Table 2, we report the values of $\Lambda [\rho +ϵ{\sigma }_{a_2}]$ as a function of ε and compare them with the ones from the first-order expansion. As predicted, we see that Λ[ρ] increases for small ε. However, the first-order expansion breaks down before ε(a₂), so that the maximum value of Λ[ρ] that we obtain is less than the one predicted by Equation (59(59) $$\Lambda \left[\rho +ϵ\left(a_2\right){\sigma }_{a_2}\right]-\Lambda \left[\rho \right]\approx 0.004765.$$(59) ). The improvement in this case is very small, but we suspect that this is due to the fact that for N = 2, the exact ${\bar{\lambda }}_3\left(2\right)$ is very close to 1.701, so that we are really hitting the boundary. In fact, in the previous example of Section 3.2, we have seen that when we start from a much less optimal density, the improvement in Λ[ρ] with our procedure is much larger.

Table 2. Exact values $\Lambda [\rho +ϵ{\sigma }_{a_2}]$ for various values of ε < 0, compared with the first-order expansion (see Section 3.3).

ε	$\Lambda [\rho +ϵ{\sigma }_{a_2}]$	Λ[ρ] + εI(a2)
0	1.699052	1.699052
−0.01	1.700487	1.701354
−0.015	1.700833	1.702505
−0.02	1.700868	1.703656
−0.0207	1.700843	1.703817

We have also repeated the procedure using as a starting density the one corresponding to ε = −0.02 in Table 2, but we could only slightly improve the result obtaining Λ[ρ] = 1.701052, which is, so far, our best value, (60) $${\bar{\lambda }}_3\left(2\right)>1.701052.$$(60)

4. Is a uniform density the most challenging for the Lieb–Oxford bound?

In [18], it has been argued that the tightest bound should correspond to the case of the uniform electron gas at extremely low density (equivalent to the SCE limit for a uniform density). This suggestion was made by considering electronic Hamiltonians of the form (1(1) $$\widehat{H}=\widehat{T}+{\widehat{V}}_{\mathrm{ee}}+{\widehat{V}}_{\mathrm{ext}},{\widehat{V}}_{\mathrm{ext}}=\sum _{i=1}^{N}v\left({\mathbf{r}}_i\right).$$(1) ) with particularly challenging v(r), keeping in mind that the bound increases [2] with the number of electrons N (see Equation (16(16) $${\bar{\lambda }}_3\left(N\right)<{\bar{\lambda }}_3(N+1),\underset{N\to \infty }{lim}{\bar{\lambda }}_3\left(N\right)={\bar{\lambda }}_3.$$(16) )).

With our formalism, we directly consider the most challenging wave function (or one which is very close to it, thus providing anyway a lower bound for Λ[ρ]) for each given density, and we can thus question whether a uniform density profile is really the most challenging for the bound. Already by putting together existing data, we can compare, in Table 3, the values of Λ[ρ] obtained from the (sphericalised) atomic densities of Li, Be, C, B, and Ne [21], with the ones obtained from spheres of uniform density (droplets) [24]: we clearly see that the atomic densities yield significantly higher values of Λ[ρ], as already observed for N = 2 in Table 1.

Table 3. For different values of N, we compare the Λ[ρ] obtained from atomic densities (values from [21]) with the ones obtained from spherical droplets of uniform density (values from [24]).

N	Λ[ρ] atomic	Λ[ρ] droplet
3	1.713	1.550
4	1.731	1.603
5	1.747	1.627
6	1.767	1.657
10	1.816	1.708

We have also performed calculations with the fixed spherical density profile ρ(r)∝r^1/2e^−r, which was particularly promising for N = 2 (see Table 1), for particle numbers N ≤ 60, and compared the values with the ones for spheres of uniform density, extending the calculations of [24] up to N = 60. The results are reported in Figure 1, where we clearly see that the uniform droplets give values significantly lower for Λ[ρ]. This suggests that a similar behaviour may arise in the limit N → ∞: a density with particular modulations might challenge the bound more than the uniform one.

Figure 1. Values of Λ[ρ] for the fixed density profile ρ(r)∝r^1/2e^−r compared with those for spheres of uniform density (droplets) as a function of the particle number N. The values for ρ(r)∝r^1/2e^−r are significantly higher than those for uniform densities. The size extrapolation for uniform densities is also shown, where the fitting parameters are a = 1.918, b = −0.3253, c = −0.2791.

Our new value for the uniform sphere at N = 60, Λ = 1.818 sets an improved upper bound [3], equal to −1.343, for the low-density uniform electron gas indirect energy per particle w, which, then, must be between (61) $$-1.45\le w\le -1.343,$$(61) where the lower bound −1.45 has been proven in [25]. We have also performed a size extrapolation of our Λ[ρ] for the droplets of uniform density of Figure 1, by fitting our data to a liquid-drop model expansion (62) $${\Lambda }_{\mathrm{unif}}\left[N\right]=a+b{N}^{-1/3}+c{N}^{-2/3},$$(62) finding a = 1.918, b = −0.3253, and c = −0.2791. The fitting function is also shown in Figure 1. The value of the fitted parameter a gives our N → ∞ extrapolation for Λ[ρ] in the uniform electron gas, Λ_unif[N → ∞] = 1.918. This value can be compared with the one obtained by taking the r_s → ∞ limit of popular LDA parametrisations: for example, the PW92 [26] parametrisation yields 1.947 at zero polarisation and 1.977 for the fully polarised case, while the VWN [27] at zero polarisations gives 1.9043.

After Lewin and Lieb [3] showed that the value 1.955 = 1.4442/A₃ does not correspond to an indirect energy, our value Λ[ρ] = 1.91175 for N = 60 and spherically symmetric density profile ρ(r)∝r^1/2e^−r is the highest value of λ₃[Ψ] ever observed, setting a new lower bound for ${\bar{\lambda }}_3\left(N\right)$ for any N, so that, rigorously (63) $$1.91175\le {\bar{\lambda }}_3\le 2.215,$$(63) or, in terms of the constant C₃ in Equation (9(9) $$-C_D\int d^{D}r{\rho }_{\Psi }{\left(\mathbf{r}\right)}^{1+1/D}\le W\left[\Psi \right].$$(9) ) (64) $$1.4119\le C_3\le 1.6358.$$(64)

5. Conclusions and perspectives

We have developed a method to maximally challenge the Lieb–Oxford bound, using optimal (or nearly optimal) trial wave functions that can be constructed from a given density. This allows us to rewrite the most challenging bound for a given number of particles directly as a density functional. As a first application of the method,

• we improved (see Equation (60(60) $${\bar{\lambda }}_3\left(2\right)>1.701052.$$(60) )) the constant in the LO bound for N = 2, which provides a constraint to develop new metaGGA functionals [12];

• we have given an improved lower bound for the constant appearing in the LO inequality valid for all particle numbers N (see Equations (63(63) $$1.91175\le {\bar{\lambda }}_3\le 2.215,$$(63) ) and (64(64) $$1.4119\le C_3\le 1.6358.$$(64) ));

• we have obtained an improved upper bound for the indirect energy per particle of the low-density uniform electron gas (see Equation (61(61) $$-1.45\le w\le -1.343,$$(61) )).

In future works, we will analyse systematically the bound for larger particle numbers N, trying to give improved lower bounds for ${\bar{\lambda }}_D\left(N\right)$ and for ${\bar{\lambda }}_D$.

More generally, from this study, we have learned that it is quite difficult to predict which densities will maximally challenge the bound (see for example, Table 1: the trends reported there seem totally unpredictable). For sure, we observe that, for finite N, a uniform density is not the one that challenges the bound the most, suggesting that the indirect energy of the uniform gas at low density may not provide the tightest bound, contrary to what was previously suggested.

Acknowledgments

We are really happy to dedicate this paper to Andreas Savin. His curiosity, deep knowledge, integrity, creativity and generosity are an infinite source of inspiration. We are very grateful to Lukas Schimmer, Mathieu Lewin, Kieron Burke and John Perdew for a critical reading of the manuscript, providing several suggestions to improve the paper.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

Financial support was provided by the European Research Council under H2020/ERC Consolidator Grant “corr-DFT” [grant number 648932]; and the Netherlands Organization for Scientific Research (NWO) through an ECHO grant [717.013.004].

Appendix 1. Strictly correlated electrons

The minimising antisymmetric wave function Ψ = Φ_KS[ρ] in the definition of the density functional of the non-interacting kinetic energy [28], (A1) $$T_{\mathrm{s}}\left[\rho \right]=\underset{\Psi \to \rho }{min}⟨\Psi |\widehat{T}|\Psi ⟩,$$(A1) is usually a Slater determinant of Kohn–Sham orbitals. In contrast, the minimising one in (A2) $$V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]=\underset{\Psi \to \rho }{min}⟨\Psi |{\widehat{V}}_{\mathrm{ee}}|\Psi ⟩$$(A2) is (or it is very close to) a state Ψ = Ψ_SCE[ρ] with strictly correlated electrons (SCE). Ψ_SCE[ρ] is not a regular wave function but a Dirac-type distribution. Its position representation is singular, (A3) $$\begin{array}{c}\hfill {\Psi }^{\mathrm{SCE}}\left[\rho \right]({\mathbf{r}}_1{\sigma }_1,...,{\mathbf{r}}_N{\sigma }_N)=0\text{for}({\mathbf{r}}_1,...,{\mathbf{r}}_N)\notin {\Omega }_0\left[\rho \right].\end{array}$$(A3) Here, Ω₀[ρ] is a D-dimensional subspace of the (N × D)-dimensional configuration space of the N-electron system, given by (A4) $${\Omega }_0\left[\rho \right]=\left\{(\mathbf{s},{\mathbf{f}}_2\left(\mathbf{s}\right),...,{\mathbf{f}}_N\left(\mathbf{s}\right))\right|\mathbf{s}\in S_{\rho }\subset {\mathbb{R}}^D\},$$(A4) where $S_{\rho }=\{\mathbf{r}\in {\mathbb{R}}^D|\rho \left(\mathbf{r}\right)\ne 0\}$ is the spatial region of non-zero density. f₁(s) ≡ s, f₂(s),… , f_N(s) are co-motion functions: In an SCE state, a configuration (r₁,… , r_N) is observable only when its N positions obey the relations, r_n = f_n(r₁), n = 2,… , N. Then, the distance between electrons i and j is |f_i(s) − f_j(s)|, fixed by the position r₁ = s of electron 1. Consequently, we have (A5) $$V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]=\int {\mathrm{d}}^{D}s\frac{\rho \left(\mathbf{s}\right)}{N}\sum _{i=1}^{N-1}\sum _{j=i+1}^N\frac{1}{|{\mathbf{f}}_i\left(\mathbf{s}\right)-{\mathbf{f}}_j\left(\mathbf{s}\right)|}.$$(A5) This is truly a density functional, since the co-motion functions are fixed by the density, f_n(s) = f_n[ρ](s). For a large class of densities ρ, the functions f_n[ρ](s) and thus the functional V^SCE_ee[ρ] can be evaluated rigorously [19,21–23]. Its functional derivative turns out to be [29,30] (A6) $$\frac{\delta V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]}{\delta \rho \left(\mathbf{r}\right)}=v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right),$$(A6) with the SCE external potential v_SCE[ρ](r), fixed by (A7) $$\nabla v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)=-\sum _{n=2}^N\frac{\mathbf{r}-{\mathbf{f}}_n\left(\mathbf{r}\right)}{|\mathbf{r}-{\mathbf{f}}_n{\left(\mathbf{r}\right)|}^3}.$$(A7) As usual, the functional derivative is determined up to a constant, which for finite systems we fix by requiring that the potential vanishes at infinity.

For example, a spherical two-electron density ρ(r) = ρ(r) in D dimensional space has the co-motion function [19,23] (A8) $${\mathbf{f}}_2\left(\mathbf{s}\right)=-f\left(s\right)\frac{\mathbf{s}}{\left|\mathbf{s}\right|}.$$(A8) In terms of the invertible function (A9) $$N_e\left(s\right)\equiv {\int }_{\left|\mathbf{r}\right|\le s}{\mathrm{d}}^{D}r\rho \left(r\right),$$(A9) the radial co-motion function is given by (A10) $$f\left(s\right)=N_e^{-1}\left(2-N_e\left(s\right)\right).$$(A10) Equation (A8(A8) $${\mathbf{f}}_2\left(\mathbf{s}\right)=-f\left(s\right)\frac{\mathbf{s}}{\left|\mathbf{s}\right|}.$$(A8) ) implies that |r − f₂(r)| = r + f(r) and, due to Equation (A5(A5) $$V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]=\int {\mathrm{d}}^{D}s\frac{\rho \left(\mathbf{s}\right)}{N}\sum _{i=1}^{N-1}\sum _{j=i+1}^N\frac{1}{|{\mathbf{f}}_i\left(\mathbf{s}\right)-{\mathbf{f}}_j\left(\mathbf{s}\right)|}.$$(A5) ), (A11) $$V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]=\frac{1}{2}\int {\mathrm{d}}^{D}r\frac{\rho \left(r\right)}{r+f\left(r\right)}.$$(A11) Due to Equation (A7(A7) $$\nabla v_{\mathrm{SCE}}\left[\rho \right]\left(\mathbf{r}\right)=-\sum _{n=2}^N\frac{\mathbf{r}-{\mathbf{f}}_n\left(\mathbf{r}\right)}{|\mathbf{r}-{\mathbf{f}}_n{\left(\mathbf{r}\right)|}^3}.$$(A7) ), the SCE external potential, with v_SCE → 0 for r → ∞, is (A12) $$v_{\mathrm{SCE}}\left[\rho \right]\left(r\right)={\int }_r^{\infty }\frac{\mathrm{d}s}{{[s+f\left(s\right)]}^2}.$$(A12)

For any N-electron density ρ(r) with co-motion functions f_n(r) (n = 2,… , N), we may consider the continuous series of scaled N-electron densities ρ_ξ(r) = ξ^Dρ(ξr), with ξ > 0 and ∫d^Dr ρ_ξ(r) = ∫d^Dr ρ(r) = N. The co-motion functions f^(ξ)_n(r) of ρ_ξ(r) are given by (A13) $${\mathbf{f}}_n^{\left(\xi \right)}\left(\mathbf{r}\right)=\frac{1}{\xi }{\mathbf{f}}_n\left(\xi \mathbf{r}\right).$$(A13) Therefore, the functional (A5(A5) $$V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]=\int {\mathrm{d}}^{D}s\frac{\rho \left(\mathbf{s}\right)}{N}\sum _{i=1}^{N-1}\sum _{j=i+1}^N\frac{1}{|{\mathbf{f}}_i\left(\mathbf{s}\right)-{\mathbf{f}}_j\left(\mathbf{s}\right)|}.$$(A5) ) has the simple scaling property (A14) $$V_{\mathrm{ee}}^{\mathrm{SCE}}\left[{\rho }_{\xi }\right]=\xi V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right].$$(A14)

We should remark that the SCE wave function as a minimiser for the electron–electron interaction energy has been first conjectured on physical grounds [19,21]. In recent years, it was recognised that the problem posed by the minimisation (A2(A2) $$V_{\mathrm{ee}}^{\mathrm{SCE}}\left[\rho \right]=\underset{\Psi \to \rho }{min}⟨\Psi |{\widehat{V}}_{\mathrm{ee}}|\Psi ⟩$$(A2) ) is equivalent to an optimal transport problem with Coulomb cost [23,31]. Since then, the optimal transport community has produced several rigorous results. In particular, the SCE state has been proven to be the true minimiser for any N in 1D [22] and in any dimension for N = 2 [23]. For more general cases, it has been shown that the minimiser might not be of the SCE form [32]. Even in that case, however, SCE-like solutions seem to be able to go very close to the true minimum [33].