Fermionic Statistics in the Strongly Correlated Limit of Density Functional Theory.

Exact pieces of information on the adiabatic connection integrand, Wλ[ρ], which allows evaluation of the exchange-correlation energy of Kohn-Sham density functional theory, can be extracted from the leading terms in the strong coupling limit (λ → ∞, where λ is the strength of the electron-electron interaction). In this work, we first compare the theoretical prediction for the two leading terms in the strong coupling limit with data obtained via numerical implementation of the exact Levy functional in the simple case of two electrons confined in one dimension, confirming the asymptotic exactness of these two terms. We then carry out a first study on the incorporation of the Fermionic statistics at large coupling λ, both numerical and theoretical, confirming that spin effects enter at orders ∼e-√λ.

Introduction

Density functional theory (DFT)[1] and the Kohn–Sham (KS) formalism[2] have been a remarkable advancement for electronic structure calculations, allowing the theoretical study of a vast class of processes in natural sciences, from physics to chemistry to biology. In KS DFT, a self-consistent machinery allows one to map the interacting electronic system into a noninteracting model endowed with the same density. Although formally an exact theory, approximations are needed for the exchange-correlation energy functional, Exc[ρ], which encloses all the complicated effects arising from the electron–electron interaction. Despite the improvement of approximate functionals in the last 30 years, several phenomena are still problematic for DFT: among the most striking cases, KS DFT shows problems in dealing with the description of van der Waals interactions, strong correlation causing charge-localization effects (i.e., low density electronic systems, Mott insulators, etc.) and dissociation processes even in simple molecules.[3,4]

In recent years, a new class of functionals that rely on integrals of the density[5−10] rather than on the usual scheme of the “Jacob’s Ladder”[11] have been proposed, inspired by the mathematical structure of what has become known as the strictly correlated electrons (SCE) limit of DFT.[12−14] In this semiclassical limit, the physical system is mapped onto an infinitely strongly interacting one with the same density ρ, where the electron–electron interaction dominates over the kinetic energy, which is suppressed: in this sense, SCE is the counterpart of the noninteracting KS system. Via the adiabatic connection formalism,[15−17] which is based on an integration over the coupling strength λ, these two limits can provide exact information on Exc[ρ], for example, via interpolated forms of the adiabatic connection integrand.[6−10,13,18]

Although it has been very recently rigorously proved that the SCE provides the exact strong-coupling (or low-density or semiclassical) limit of the Levy–Lieb functional,[19,20] the validity of the expression for the next leading term in the expansion at large λ, first conjectured and studied in refs (12) and (21), has not been proved yet and remains for now only a very plausible hypothesis. Moreover, the inclusion of the statistics in the theory is a problem that has not been investigated at all: the intrinsic semiclassical nature of the SCE limit prevents one from taking into account the difference between bosons and Fermions (which is suppressed, as electrons in the SCE limit are always far apart from each other). Nevertheless, the effects due to the statistics of the particles or due to different spin states become important when the electron–electron interaction is large but not infinite: the kinetic energy, which is nonzero as a consequence of zero point oscillations around the SCE minimum, allows electrons to be subject to Pauli’s principle.

The aim of this work is to address these two issues, namely, (i) to probe the validity of the second term in the asymptotic expansion of the adiabatic connection integrand at large λ, and (ii) to study the inclusion of the Fermionic statistics in the large-λ limit. We focus on the easiest case of N = 2 electrons confined in one dimension (1D) because in this case we can also compute accurate numerical results for the exact Levy functional at large λ, which allows us to carefully validate our asymptotic analytic expansions.

This paper is organized as follows. In section , we briefly review the theory of SCE and zero point oscillations (ZPO) in the strong coupling limit; then we outline in section the numerical method used to calculate the exact Levy functional for two electrons in 1D. In section , we compare the theoretical predictions with the numerical data obtained via the method described in section , and in section , we describe how to induce a Fermionic statistics in the ZPO wave function, comparing the singlet–triplet splitting in the expectation of the electron–electron repulsion, V̂ee, with the numerical data in section . Last, we give our conclusions and outline future steps in section .

Theoretical Background

The exchange-correlation energy in Kohn–Sham DFT can be expressed exactly in terms of an integral over the coupling strength λ,of the adiabatic connection integrand, Wλ[ρ],where V̂ee is the operator for the electron–electron repulsion,and U[ρ] is the Hartree functional. In eqs –3, . While D = 3 is obviously the most interesting case in Chemistry, in Physics it is common practice to consider also low-dimensional effective problems with D = 1 and 2. Accordingly, while in D = 3 or 2 usually vee(x) = 1/x, in 1D people often resort to an effective interaction, which will be discussed in section . The wave function appearing in eq , ψλ[ρ], is the Fermionic wave function, which minimizes the generalized Hohenberg–Kohn functional in the constrained-search Levy formulation,[22]with T̂ the kinetic energy operator.

If the density ρ is both N and V representable for every λ, ψλ is the ground state of the λ-dependent Hamiltonianwhere V̂λext[ρ] = ∑vλext[ρ](r) is the one body operator for the external potential providing the density ρ(r).

Strictly Correlated Electrons (SCE)

In the limit λ → ∞, the adiabatic connection integrand approaches a finite value,[12,14,23,24]given bywhere in the last step we used the fact that the external potential v(r) is the Lagrange multiplier for the constraint ψ → ρ.[14,25−27] The finiteness of W∞[ρ] stems from the fact that the electrons must be confined in a given finite density and thus cannot escape infinitely far from each other.[12,14,23,24]

Since for λ → ∞, we expect ⟨ψλ|T̂|ψλ⟩ ≈ O(√λ)[12,14,21] (see also ref (19) for a rigorous proof); only an external potential V̂λext ≈ O(λ) can compensate the infinitely strong electronic repulsion in eq . Hence, we expect that the asymptotic large-λ expansion of the external potential provides the finite limitcorresponding to the potential needed to counteract exactly the Coulomb repulsion in this semiclassical limit[14] (notice that here we use the same notation as in refs (14) and (21), in which vSCE is minus the functional derivative of VeeSCE[ρ]; in more recent works, for example, in refs (28−31), the notation vSCE has been used with the opposite sign, to denote a potential that represents, rather than compensates, the net electron–electron repulsion force acting on an electron in r).

As a consequence of eq , the leading order of eq can be written asThe Hamiltonian in eq describes a N particle classical system; minimization in eq requires the associated probability density (a distribution that plays the role of |ψ|2) to be nonzero only on the set Ω0 of configurations ≡ (r1, ..., r) for which the classical potential energy function,assumes its global minimum.

The SCE ansatz consists in searching for potentials that make Ω0 a D-dimensional subset of the configuration space, defined by a set of co-motion functions (or optimal maps)[12,14]Co-motion functions provide, after the measurement of the position of any one chosen reference electron, the positions of the remaining N – 1 electrons. They are endowed with group properties[14]and satisfyDefining |ψSCE[ρ]|2 ≡ |ψλ→∞[ρ]|2, in the SCE limit |ψSCE|2 yields a distribution that represents a gas of electrons frozen in strictly correlated positions, nevertheless yielding a smooth density by behaving as a “floating” Wigner crystal in a metric,[14] being any permutation of N particles. Thus, among the set of all functions f̃(s) satisfying eqs and 13, the co-motion functions are the minimizers of the electron–electron repulsion, leading to a corresponding SCE potential[14,32,33]In the rest of section , we shall restrict to the case of two electrons in 1D: this is the simplest case to study both numerically and analytically, as most quantities of interest can be expressed in closed form. Moreover, mathematical simplification of the concepts outlined so far shall suggest a clearer and physically straightforward interpretation. For the general approach, we refer the reader to refs (14) and (21).

SCE for 2 Electrons in 1D

In the 1D case, a conjectured solution for the co-motion functions for any number of electrons N was presented in ref (12) and proved to be exact later in ref (34). For N = 2, defining f1(s) ≡ s, f2(s) ≡ f(s), it readswhereAccordingly, eq readswhere vSCE(x) can be obtained by integrating the last line of eq 15. In 1D, the support Ω0 of the minimum of Epot(x1, x2) is just a parametric curve (s, f(s)) on the (x1, x2) plane, , with f(s) given by eq . As an example, in Figure , we report Epot(x1, x2) and the corresponding Ω0 for a simple analytic density (a Lorentzian, see the following for details).

Figure 1

Function Epot(x1, x2) as a 3D plot (top) and as a contour plot (bottom) for the Lorentzian density ρ2(x) of eq . The 1D set Ω0 is shown as a pair of red curves in the contourplot.

On the Convexity of Interaction in 1D

In 1D, it is not suitable to use the interaction 1/|x|, since some key features of the physical model are lost: due to the divergence of |x1 – x2|–1 at x1 = x2, both bosonic and Fermionic wave functions are forced to have the same nodal surface and thus the same energy; moreover, the Hartree energy, U[ρ], is not finite. It is thus customary to resort to an effective 1D interaction, which is finite at the origin. One of the most commonly used ones is the soft Coulomb, that is,which is not convex in the region . However, in 1D convexity of the interaction vee(|x|) is a necessary condition[34] to prove that Ω0 is determined by the co-motion function of eq .

We believe it is important to clarify this with an example, as nonconvex interactions are often used when probing DFT approximations using 1D physics and chemistry models (e.g., see refs (31 and 35−37).). Referring to Figure , we shall briefly discuss a soft Coulomb interaction with a = 4. We definewhere vSCE(r) is obtained via eqs 15 and 16, and vdual(r) is obtained numerically from the dual problem, which basically corresponds to the last line of eq (see refs (27, 33, and 38) for details on the implementation of the numerical dual formulation of the SCE functional).

Figure 2

Case of a 1D Lorentzian density (the density is the same as in Figure ) where the interaction is .

In panel c of Figure , we report EpotSCE() and Epotdual() along the negative diagonal x2 = −x1. We see that in this case the manifold described by eq is only a local minimum for EpotSCE(), which has its global minimum at (0,0). In the energy landscape Epotdual(), instead, the two minima become degenerate. As it can be seen from inset of panel d, the support of the minimum of Epotdual, getting contribution also from x1 = −x2 close to the origin, is not provided by a solution of the kind in eq .

On the other hand, an effective Coulomb interaction in 1D of the formbeing always convex, does not suffer from these problems: with this interaction, as it can be seen from Figure , the manifold Ω0 is parametrized by the co-motion functions of eq . In this case, vSCE(r) and vdual(r) are exactly equal. In order to work in this framework (which correctly models the 3D physics, in which the electrons stay always away from each other in the strong-coupling limit with Coulomb interaction), throughout the rest of this paper, we use eq with a = 1.

Figure 3

Case of a 1D Lorentzian density with vee = (4 + |x|)−1.

Zero Point Oscillations

Equation 15 provides an expression for the leading term of the adiabatic connection integrand in the λ → ∞ limit. An ansatz for the subleading term in eq , which is due to zero-point oscillations of the strongly interacting electrons, can be obtained following the treatment in ref (21) and reads asIn analogy with the expansion of the adiabatic connection at λ = 0, (Wλ→0[ρ] = W0[ρ] + W0′[ρ]λ + ...), the subleading contribution in the large λ limit has been denoted as W∞′[ρ][12] (note that the prime on W∞′ does not denote a derivative).

In the λ → ∞ limit, we expect the electrons to be forced to stay in the vicinity of Ω0, with the (relatively small) kinetic energy due to zero-point oscillations allowing them to explore the part of potential energy landscape Epot(x1, x2) close to this degenerate minimum (i.e., the darker regions around the red curve in Figure ).

Considering only small oscillations around the minimum of Epot allows for a harmonic expansion around the manifold Ω0,where f1(s) = s, f2(s) = f(s), ESCE = Epot(s, f(s)), and Mμν(s) is the Hessian of Epot evaluated in Ω0:Diagonalization of Mμν(s) suggests a natural set of coordinates associated with its (non-negative) eigenvalues, ωμ(s)2, which can be labeled in such a way thatSince ω12(s) is proportional to the curvature of Epot along Ω0 (which is flat, as the minimum is degenerate), while ω22(s) is connected to the curvature orthogonal to Ω0, it is possible to introduce a set of curvilinear coordinates in which every point in the configuration space sufficiently close to Ω0 can be described in terms of its closest point to the manifold Ω0 and its distance from it.[21] We therefore introduce a local coordinate transformation, from Cartesian to the coordinates associated with the eigenvectors of the Hessian Mμν(s):The coordinate q gives the distance of point (x1, x2) from the closest manifold branch, while s is the parametric value of the closest point on the manifold Ω0, around which the oscillation takes place, see Figure for an illustration.

Figure 4

Coordinate transformation (x1, x2) → (s, q).

Explicitly, the coordinate transformation readsEquation becomes diagonal in terms of these local normal modes:and we see that ω2(s) can be associated with the zero-point vibrational frequency around the SCE minimum. The only nonzero frequency associated with the Hessian of Epot for 2 electrons in 1D is simply given by[37]The correction due to the zero point oscillations to the adiabatic connection can now be written as a weighted sum of harmonic oscillators’ energies, since the degeneracy with respect to s allows one to weight the energy of each configuration with the density ρ(s): W∞′[ρ] readswhich is a particular case of eq 81 in ref (21). The corresponding W∞[ρ] reads in this case

Constrained Search Method for Two Electrons in 1D

The Levy constrained-search functional for a N-representable density is defined as[22]By restricting the search over spatially symmetric (ΨS) or antisymmetric (ΨT) wave functions, it is possible to define, respectively, FLevyλ, S[ρ] and FLevyλ, T[ρ], finding the corresponding minimizing wave function for a singlet and triplet state associated with the same physical density ρ(x).

In previous work,[39] the Levy constrained search was found for the exact density-matrix functional of the two-site Hubbard model using an analytic formula. However, in this work, the constrained search is carried out via a stochastic minimization of the wave function as in ref (40) to give the exact density functional of eq .

We will focus on the details to carry out a general optimization for two electrons. First, construct an initial wave function that integrates exactly to the density, ρ(x). For the singlet, this is trivial, as . However, for the triplet, one route is to find two orbitals that sum up to the given density, ρ(x) = ϕ12(x) + ϕ22(x) and then an initial wave function can be constructed, ΨinitialT(x1, x2) = {ϕ1(x1)ϕ2(x2) – ϕ1(x2)ϕ2(x1)}/√2. The simplest way to find two orbitals is to use a division of space into two, which is actually done by the inverse cumulant of eq For practical calculations on a finite grid, the orbitals have to overlap at the two grid-points x = L and x = R on the left and right of the point in which the density integrates to 1, L < Ne–1(1) and R > Ne–1(1), and they satisfy the following equationsfor normalization, density constraint, and zero overlap. The solution is given byand the other points determined from eqs –34c with one negative square root chosen to satisfy eq .

With these initial wave functions that integrate to ρ(x), the key to the procedure is to define moves of the spatial part of the wave function that maintain the density. When the density is represented on a grid (we generally use 200 grid points), this can be done based on a move of four points of the wave function at once as outlined in ref (40). These moves are attempted and accepted if they lower the energy of eq . This is then repeated many times to carry out a stochastic optimization of the wave function, and convergence is typically found in 20 000 steps for all values of λ.

Adiabatic Connection at Large λ: Numerical and Analytic Results

The main purpose of this section is to compare the data obtained via the constrained search method outlined in section with eqs , 30, and 31. In order to probe the validity of the ZPO approach, we shall discuss a set of three 1D densities, which integrate to N = 2 particles in a box, interacting via the effective Coulomb interaction of eq .

Our first two densities,share the property of having both an analytical expression as well as analytical co-motion functions, reported in Appendix A. Our third one, ρ3(x), is a numerical density for the 1D He atom with the same interaction (eq ) on the interval [−5,5] and has no analytical form.

Using eqs , 29,, and 31, we find for the different densities the values of Table , where we also report the values extracted from the numerical data obtained via the constrained search method. The numerical W∞[ρ] is the value of eq at λ = ∞, W∞[ρ] = minΨ→ρ⟨Ψ|Vee|Ψ⟩ – U[ρ], and W∞′[ρ] is calculated by finite difference, W∞′ = (W500 – W∞)√500. The asymptotic expansion of eq , with the values of W∞[ρ] and W∞′[ρ] obtained from eqs –31, is also compared to the numerical results for the Levy functional at large λ in Figure for the three densities. We see that the agreement is excellent. This provides the first numerical evidence that the zero point term should be exact for arbitrary symmetrical density, at least for one-dimensional systems (a related numerical study[41] addressed the SCE leading term only, while an exact result was recently obtained for a uniform density defined on a ring[31]). We hope that this result will trigger, similarly to what has been done recently for the leading term W∞[ρ],[19,20] works on a rigorous proof for the subleading term W∞′[ρ].

Table 1

W∞[ρ] and W∞′[ρ] from the Analytical Treatment, Eqs –31, and the Numerical Constrained Search Method for the Densities Considered

	W∞[ρ] + U[ρ]		W∞′[ρ]
	analytic	numerical	analytic	numerical
ρ1(x)	0.31229	0.31237	0.12209	0.12076
ρ2(x)	0.27282	0.27291	0.11635	0.11573
ρ3(x)	0.40208	0.40212	0.17223	0.17521

Figure 5

Exchange-correlation energy in the strongly correlated limit of DFT for different densities. Insets: plot of the related density. Blue dots: Numerical results from the constrained search method. Red curve: the expansion of eq with the values of W∞[ρ] and W∞′[ρ] computed from eqs and 31.

The Effects of the Spin State at Large λ

The Schrödinger equation corresponding to the order O(√λ) in the asymptotic expansion of the density fixed λ-dependent Hamiltonian of eq is, in the curvilinear coordinates system, the equation of a harmonic oscillator whose spring constant depends on s,[21]where the term , denoted in ref (21) as V(0), is the correction to the external potential of order √λ computed on the manifold.[21] Its role is to keep the energy E(0) in the right-hand-side of eq independent of s (otherwise the wave function would collapse in one particular value of s, the one with lowest energy, and the density constraint would be lost; see ref (21) for details).

It has been suggested[21] that, since the Hamiltonian (eq ) describes an uncoupled set of harmonic oscillators, the leading order in the wave function ψλ factorizes into a product of Gaussians, with amplitude depending on √λ and on s through the curvature of the manifold,J(s, q) being the Jacobian of the transformation from Cartesian to curvilinear coordinates. As a consequence, the effect on the energy of the introduction of statistics has been conjectured to be[21,42] to the leading order in the λ → ∞ limit, ∼e–√λ, this being the order of magnitude of the overlap between two Gaussians centered in different positions having the form of eq . This hypothesis is the analogue for a nonuniform density of the one used by Carr[43] for the uniform electron gas at low density.

The purpose of this section is hence to investigate the splitting in energy between the expectation value of V̂ee evaluated on the singlet and on the triplet state:We will check if the hypothesisis consistent with the results provided both via an explicit construction of an antisymmetric and a symmetric state starting from eq and via the accurate results from the constrained search method. We will also discuss possible routes to simplify the inclusion of spin starting from the large-λ expansion.

Explicit Antisymmetrization of the ZPO Wave Function

Being expressed in the (s, q) curvilinear coordinate system, the wave function in the form of eq is not suitable for a straightforward antisymmetrization. In order to do so, we first have to retrieve the Cartesian coordinates, that is, writeby inverting eq and only then proceed to construct a symmetric (singlet) and an antisymmetric (triplet) state.

First, a remark is in order: as it can be seen from Figure , there are regions where the (s, q) coordinates are ill-defined (a cone in the second and fourth quadrants, respectively, symmetric with respect to the diagonal x2 = −x1). Nevertheless, as the Fermionic statistics affects particles mostly on the diagonal x2 = x1, the contributions from these regions should be negligible for our purposes.

Figure 6

Top: (sA, B, qA, B) describe the position of a particle as a function of its distance from the branch of the manifold (A = red, B = orange). Bottom: a generic point (x1, x2) can be written as a function of (sA, qA) (red) or (sB, qB) (orange). When we exchange the position of the particles, the roles of the curvilinear coordinate exchange accordingly.

Given the set of positions (x1, x2), the curvilinear frame we used in the ZPO regime prescribes to choose the closest branch of the manifold Ω0: labeling these branches “A” and “B”, this means choosing among two possible coordinates, namely (sA, qA) and (sB, qB), taking the one with the smallest q. However, if we want to describe spin effects, we must take into consideration the overlap of the ZPO wave functions centered on the two different branches, since swapping positions between two electrons amounts to swap the point (s, f(s)) around which the oscillation in curvilinear coordinates takes place with respect to the diagonal x1 = x2.

This means actually writing the ZPO wave function (eq ) in Cartesian coordinates with respect to the two different branchesIt should be noted that, sincewe also haveAs a consequence, the exchange of the two particles’ position actually means switching branch in eq . In this way, antisymmetrization of eq reads aswhere we have labeled with A and B the two branches of the co-motion function and approximated the λ-dependent normalization constant to 1/√2, according toas the terms neglected would be of higher order in e–√λ.

In Figure , we show the singlet and triplet wave functions obtained in this way from the density ρ2(x) for λ = 100. We see that the two wave functions are both concentrated around the manifold Ω0, with the triplet having the expected node at x1 = x2. In Figure , we compare our singlet and triplet wave functions with the ones obtained via the constrained search method for the density ρ2(x) and λ = 500. We see that the singlet and triplet ZPO wave functions agree very well with the accurate ones for the constrained search method. In particular, in panels c and f, we report the difference between the ZPO and constrained-search singlet and triplet, respectively, which appears to be rather small.

Figure 7

3D plot of singlet and triplet wave function associated with density ρ2(x), with coupling constant λ = 100, over the contour plot of Epot(x1, x2) as from eq . Top: singlet wave function. Bottom: triplet wave function.

Figure 8

Comparison of the ZPO wave function for singlet (a) and triplet (d) state with the wave function provided by the constrained search method for the density ρ2(x) with λ = 500 (respectively, panels b and e). Panels c and f show, respectively, the difference between panels a and b and the difference between panels d and e.

Evaluating the spin splitting in the expectation value of the electron–electron interaction in the singlet and triplet state from our construction yieldsan expression that is clearly of orders e–√λ and that will be compared with the numerical results from the constrained-search method in section .

Alternative Strategies to Include the Statistics in the λ ≫ 1 Regime

In this section, we outline some strategies to simplify the procedure of section , namely, disentangling the oscillations of the two electrons around their equilibrium positions and using the Hellman–Feynman theorem to provide an exact relation for the singlet–triplet splitting uniquely in terms of the kinetic energy operator.

With the use of 23, eq becomesAn uncoupled approximation is justified when the off-diagonal elements of the Hessian are small compared to the diagonal ones. In our picture, this is equivalent to removing the dependence of the s coordinate from (x1, x2), leaving us with a Hamiltonian that depends parametrically on s and that describes uncoupled oscillations around their equilibrium positions s and f(s)Defining M11(s) ≡ Ω12(s) and M22(s) ≡ Ω22(s) = Ω12(f(s)) andit is clear that, for every fixed s, a properly antisymmetrized eigenfunction for eq readswhere Nλ± is just the normalization factor. However, in our case this approximation is hardly going to hold: the off-diagonal elements of Mμν in the basis of Cartesian coordinates are of the same order of magnitude of the diagonal ones, and such approximations typically largely overshoot the V̂ee expectation value. However, this approximation might be used to construct a basis to expand the full ZPO wave function, which will be explored in future works.

Finally, another way to compute Δλ[ρ] is by making use of the Hellman–Feynman theorem. We definewhere ΨλS, T[ρ], as already mentioned in section , is the wave function minimizing FλS, T[ρ] when the search is constrained to the corresponding symmetry sector. Since both singlet and triplet wave functions are required to be stationary, we will have two separate Hellmann–Feynman theoremsand defining Δλkin[ρ] ≡ TS[ρ](λ) – TT[ρ](λ) ≤ 0, we can also obtain the singlet–triplet splitting fromThis approach should bypass the numerical difficulties arising from evaluation of integrals involving 2-body operators, and it might be, at a later stage, more suitable for implementing in realistic models the ideas explained in this paper and will be subject of future works.

Results for the Singlet–triplet Splitting

In this section, we compare the results of our analysis on the ZPO wave function with the data obtained via constrained search method. In particular, to check the validity of eq , we compare in Figure the splitting from eq with data from numerical constrained search method, which numerically proves the ansatz of eq . The bottom panel of Figure shows in fact that log Δλ[ρ] is linear in √λ both for the constrained search method (blue) and the calculation from eq (red).

Figure 9

Splitting in the Vee expectation energy between singlet and triplet state. Inset: plot of the related density. Numerical fit provides α[ρ2] = 0.293, β[ρ2] = 0.978 for constrained search method (blue) and α[ρ2] = 0.361, β[ρ2] = 0.725 for eq (red).

Although our results show qualitative agreement with the data, quantitative discrepancy is evident. Since the agreement between the two different wave functions used, as shown in Figure , is quite good, this discrepancy could be due to either the numerical noise arising from the smallness of the numbers involved, or the fact that, the effect being small, the differences between the two wave functions are still relevant.

Conclusions and Perspectives

We have investigated the validity of the expansion of the adiabatic connection integrand in the strong coupling limit as proposed in ref (21) for three 1D densities with N = 2 electrons, by comparing the theoretical prediction with numerical data for the Levy functional (see Figure ), finding excellent agreement and thus providing the first numerical evidence of the exactness of this term for nonuniform densities.

We have implemented the Fermionic statistics in the strong-interaction limit of DFT by retrieving the zero-point wave function in Cartesian coordinates, and we have used it to evaluate the singlet–triplet splitting, comparing the results with numerical data. In this case, we had qualitative but not quantitative agreement. The main result is the confirmation that spin effects enter at orders ∼e–√λ when λ → ∞.

We expect the order at which these effects appear when λ → ∞ to be the same also for the three-dimensional case with full Coulomb interaction. In fact, it has been recently proved in refs (19) and (20) that the SCE limit is the same regardless of the statistics of the particles and that the next leading term should be order[19]O(√λ), which suggests that the relevant physics in the λ → ∞ limit is well captured by the simple 1D model considered here. In other words, even if in 3D the nodal surface of the wave function, which dictates how the spin state affects the energy, has topological features that cannot be captured by 1 or 2D models, the order at which spin effects enter is likely to remain the same (∼e–√λ), since when λ → ∞ the physics should be that of zero-point motion around the SCE minimum, with the exchange integrals exponentially vanishing with √λ. For the special case of the 3D uniform electron gas, this was already conjectured by Carr.[43] It would be very interesting to have a rigorous proof of this conjecture for the general nonuniform 3D case.

In future work, we aim to find a more explicit (approximate) expression for spin effects in terms of spin densities, namely, to provide an expression of the kindThe main challenge will then be to build approximations for the general three-dimensional case, based on quantities that can be computed routinely. A possible way to do that, could be the generalization to spin-densities of the functionals as described in refs (5, 7, and 10), which are inspired to the SCE limit and use as key ingredient the integral of the spherically averaged density around a given position r.

Another promising reasearch line is the study of the next leading term of the large-λ expansion, which could provide an improvement in the correction of the density to the required order in the ZPO wave function and could give better estimates of the electron–electron interaction. Work along all these directions is in progress.