PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of acssdACS PublicationsThis JournalSearchSubmit a manuscript
Journal of Chemical Theory and Computation
 
J Chem Theory Comput. 2017 December 12; 13(12): 6089–6100.
Published online 2017 November 7. doi:  10.1021/acs.jctc.7b00998
PMCID: PMC5729548

Fermionic Statistics in the Strongly Correlated Limit of Density Functional Theory

Abstract

An external file that holds a picture, illustration, etc.
Object name is ct-2017-00998p_0010.jpg

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–√λ.

1. Introduction

Density functional theory (DFT)1 and the Kohn–Sham (KS) formalism2 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 density510 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.1214 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,1517 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.610,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 2, we briefly review the theory of SCE and zero point oscillations (ZPO) in the strong coupling limit; then we outline in section 3 the numerical method used to calculate the exact Levy functional for two electrons in 1D. In section 4, we compare the theoretical predictions with the numerical data obtained via the method described in section 3, and in section 5, 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, Vee, with the numerical data in section 5.3. Last, we give our conclusions and outline future steps in section 6.

2. Theoretical Background

The exchange-correlation energy in Kohn–Sham DFT can be expressed exactly in terms of an integral over the coupling strength λ,

equation image
1

of the adiabatic connection integrand, Wλ[ρ],

equation image
2

where Vee is the operator for the electron–electron repulsion,

equation image
3

and U[ρ] is the Hartree functional. In eqs 13, An external file that holds a picture, illustration, etc.
Object name is ct-2017-00998p_m004.gif. 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 2.1.2. The wave function appearing in eq 2, ψλ[ρ], is the Fermionic wave function, which minimizes the generalized Hohenberg–Kohn functional in the constrained-search Levy formulation,22

equation image
4

with T the kinetic energy operator.

If the density ρ is both N and V representable for every λ, ψλ is the ground state of the λ-dependent Hamiltonian

equation image
5

where Vλext[ρ] = ∑i = 1Nvλext[ρ](ri) is the one body operator for the external potential providing the density ρ(r).

2.1. Strictly Correlated Electrons (SCE)

In the limit λ → ∞, the adiabatic connection integrand approaches a finite value,12,14,23,24

equation image
6

given by

equation image
7

where in the last step we used the fact that the external potential v(r) is the Lagrange multiplier for the constraint ψ → ρ.14,2527 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 [left angle bracket]ψλ|Tλ[right angle bracket]O(√λ)12,14,21 (see also ref (19) for a rigorous proof); only an external potential VλextO(λ) can compensate the infinitely strong electronic repulsion in eq 5. Hence, we expect that the asymptotic large-λ expansion of the external potential provides the finite limit

equation image
8

corresponding to the potential needed to counteract exactly the Coulomb repulsion in this semiclassical limit14 (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 (2831), 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 8, the leading order of eq 5 can be written as

equation image
9

The Hamiltonian in eq 9 describes a N particle classical system; minimization in eq 7 requires the associated probability density (a distribution that plays the role of |ψ|2) to be nonzero only on the set Ω0 of configurations r [equivalent] (r1, ..., rN) for which the classical potential energy function,

equation image
10

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

equation image
11

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 properties14

equation image
12

and satisfy

equation image
13

Defining |ψSCE[ρ]|2 [equivalent]λ→∞[ρ]|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

equation image
14

An external file that holds a picture, illustration, etc.
Object name is ct-2017-00998p_m016.gif being any permutation of N particles. Thus, among the set of all functions fi(s) satisfying eqs 12 and 13, the co-motion functions are the minimizers of the electron–electron repulsion, leading to a corresponding SCE potential14,32,33

equation image
15a
equation image
15b

In the rest of section 2, 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).

2.1.1. 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) [equivalent] s, f2(s) [equivalent] f(s), it reads

equation image
16

where

equation image
17

Accordingly, eq 10 reads

equation image
18

where 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, An external file that holds a picture, illustration, etc.
Object name is ct-2017-00998p_m022.gif, with f(s) given by eq 16. As an example, in Figure Figure11, 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 36a. The 1D set Ω0 is shown as a pair of red curves in the contourplot.

2.1.2. 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 |x1x2|–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,

equation image
19

which is not convex in the region An external file that holds a picture, illustration, etc.
Object name is ct-2017-00998p_m024.gif. However, in 1D convexity of the interaction vee(|x|) is a necessary condition34 to prove that Ω0 is determined by the co-motion function of eq 16.

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 3537).). Referring to Figure Figure22, we shall briefly discuss a soft Coulomb interaction with a = 4. We define

equation image
20

where 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 7 (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 Figure11) where the interaction is An external file that holds a picture, illustration, etc.
Object name is ct-2017-00998p_m078.gif.

In panel c of Figure Figure22, we report EpotSCE(r) and Epotdual(r) along the negative diagonal x2 = −x1. We see that in this case the manifold described by eq 16 is only a local minimum for EpotSCE(r), which has its global minimum at (0,0). In the energy landscape Epotdual(r), 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 16.

On the other hand, an effective Coulomb interaction in 1D of the form

equation image
21

being always convex, does not suffer from these problems: with this interaction, as it can be seen from Figure Figure33, the manifold Ω0 is parametrized by the co-motion functions of eq 16. 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 21 with a = 1.

Figure 3

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

2.2. 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 9, which is due to zero-point oscillations of the strongly interacting electrons, can be obtained following the treatment in ref (21) and reads as

equation image
22

In 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 Figure11).

Considering only small oscillations around the minimum of Epot allows for a harmonic expansion around the manifold Ω0,

equation image
23

where f1(s) = s, f2(s) = f(s), ESCE = Epot(s, f(s)), and Mμν(s) is the Hessian of Epot evaluated in Ω0:

equation image
24

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 that

equation image
25a
equation image
25b

Since ω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):

equation image
26

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 Figure44 for an illustration.

Figure 4

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

Explicitly, the coordinate transformation reads

equation image
27

Equation 23 becomes diagonal in terms of these local normal modes:

equation image
28

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 by37

equation image
29

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[ρ] reads

equation image
30

which is a particular case of eq 81 in ref (21). The corresponding W[ρ] reads in this case

equation image
31

3. Constrained Search Method for Two Electrons in 1D

The Levy constrained-search functional for a N-representable density is defined as22

equation image
32

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 32.

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 An external file that holds a picture, illustration, etc.
Object name is ct-2017-00998p_m039.gif. 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(x12(x2) – ϕ1(x22(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 17

equation image
33a
equation image
33b

For practical calculations on a finite grid, the orbitals have to overlap at the two grid-points xi = L and xi+1 = 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 equations

equation image
34a
equation image
34b
equation image
34c
equation image
34d

for normalization, density constraint, and zero overlap. The solution is given by

equation image
35

and the other points determined from eqs 34a34c with one negative square root chosen to satisfy eq 34d.

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 32. 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 λ.

4. 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 3 with eqs 22, 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 21.

Our first two densities,

equation image
36a
equation image
36b

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 21) on the interval [−5,5] and has no analytical form.

Using eqs 30, 29,, and 31, we find for the different densities the values of Table 1, 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 32 at λ = ∞, W[ρ] = minΨ→ρ[left angle bracket]Ψ|Vee[right angle bracket]U[ρ], and W[ρ] is calculated by finite difference, W = (W500W)√500. The asymptotic expansion of eq 22, with the values of W[ρ] and W[ρ] obtained from eqs 2931, is also compared to the numerical results for the Levy functional at large λ in Figure Figure55 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 study41 addressed the SCE leading term only, while an exact result was recently obtained for a uniform density defined on a ring31). 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[ρ].

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 22 with the values of W ...

Table 1

W[ρ] and W[ρ] from the Analytical Treatment, Eqs 2931, and the Numerical Constrained Search Method for the Densities Considered

5. 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 5 is, in the curvilinear coordinates system, the equation of a harmonic oscillator whose spring constant depends on s,21

equation image
37

where the term An external file that holds a picture, illustration, etc.
Object name is ct-2017-00998p_m050.gif, 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 37 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 suggested21 that, since the Hamiltonian (eq 37) 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,

equation image
38

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 be21,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 38. This hypothesis is the analogue for a nonuniform density of the one used by Carr43 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 Vee evaluated on the singlet and on the triplet state:

equation image
39

We will check if the hypothesis

equation image
40

is consistent with the results provided both via an explicit construction of an antisymmetric and a symmetric state starting from eq 38 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.

5.1. Explicit Antisymmetrization of the ZPO Wave Function

Being expressed in the (s, q) curvilinear coordinate system, the wave function in the form of eq 38 is not suitable for a straightforward antisymmetrization. In order to do so, we first have to retrieve the Cartesian coordinates, that is, write

equation image
41

by inverting eq 27 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 Figure66, 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). ...

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 38) in Cartesian coordinates with respect to the two different branches

equation image
42

It should be noted that, since

equation image
43a
equation image
43b

we also have

equation image
44a
equation image
44b

As a consequence, the exchange of the two particles’ position actually means switching branch in eq 42. In this way, antisymmetrization of eq 38 reads as

equation image
45

where 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 to

equation image
46

as the terms neglected would be of higher order in e–√λ.

In Figure Figure77, 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 Figure88, 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 10. 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, ...

Evaluating the spin splitting in the expectation value of the electron–electron interaction in the singlet and triplet state from our construction yields

equation image
47

an expression that is clearly of orders e–√λ and that will be compared with the numerical results from the constrained-search method in section 5.3.

5.2. Alternative Strategies to Include the Statistics in the λ [dbl greater-than sign] 1 Regime

In this section, we outline some strategies to simplify the procedure of section 5.1, 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 5 becomes

equation image
48

An 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)

equation image
49

Defining M11(s) [equivalent] Ω12(s) and M22(s) [equivalent] Ω22(s) = Ω12(f(s)) and

equation image
50

it is clear that, for every fixed s, a properly antisymmetrized eigenfunction for eq 49 reads

equation image
51

where 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 Vee 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 define

equation image
52a
equation image
52b

where ΨλS, T[ρ], as already mentioned in section 3, 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 theorems

equation image
53

and defining Δλkin[ρ] [equivalent] TS[ρ](λ) – TT[ρ](λ) ≤ 0, we can also obtain the singlet–triplet splitting from

equation image
54

This 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.

5.3. 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 40, we compare in Figure Figure99 the splitting from eq 47 with data from numerical constrained search method, which numerically proves the ansatz of eq 40. The bottom panel of Figure Figure99 shows in fact that log Δλ[ρ] is linear in √λ both for the constrained search method (blue) and the calculation from eq 47 (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 α[ρ ...

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 Figure88, 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.

6. 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 Figure55), 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 order19O(√λ), 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 kind

equation image
55

The 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.

Acknowledgments

Financial support was provided by the European Research Council under H2020/ERC Consolidator Grant corr-DFT [Grant Number 648932]. P.M.S. acknowledges funding from MINECO Grant No. FIS2015-64886-C5-5-P.

Appendix A. Co-motion Functions for the Analytical Densities

The Density ρ1(x)

Let’s consider ρ1(x) = sech(x)/[2 arctan(tanh(5))]. From eq 17, we have

equation image
56a
equation image
56b

and using eq 16, we find

equation image
57

with the Gudermannian function, gd(s) = arcsin(tanh(s)).

The Density ρ2(x)

Let’s consider ρ2(x) = [(1 + x2) arctan(10)]−1. From eq 17, we have

equation image
58a
equation image
58b

and using eq 16, we find

equation image
59

Notes

The authors declare no competing financial interest.

References

  • Kohn W. Electronic Structure of Matter - Wave Functions and Density Functionals. Rev. Mod. Phys. 1999, 71, 1253..10.1103/RevModPhys.71.1253 [Cross Ref]
  • Kohn W.; Sham L. J. Self-consistent equations including exchange and correlation effects. Phys. Rev. 1965, 140, A1133..10.1103/PhysRev.140.A1133 [Cross Ref]
  • Cohen A.; Mori-Sánchez P.; Yang W. Insights into current limitations of density functional theory. Science 2008, 321, 792–794.10.1126/science.1158722 [PubMed] [Cross Ref]
  • Cohen A.; Mori-Sánchez P.; Yang W. Challenges for density functional theory. Chem. Rev. 2012, 112, 289–320.10.1021/cr200107z [PubMed] [Cross Ref]
  • Wagner L. O.; Gori-Giorgi P. Electron avoidance: A nonlocal radius for strong correlation. Phys. Rev. A: At., Mol., Opt. Phys. 2014, 90, 052512..10.1103/PhysRevA.90.052512 [Cross Ref]
  • Zhou Y.; Bahmann H.; Ernzerhof M. Construction of exchange-correlation functionals through interpolation between the non-interacting and the strong-correlation limit. J. Chem. Phys. 2015, 143, 124103..10.1063/1.4931160 [PubMed] [Cross Ref]
  • Bahmann H.; Zhou Y.; Ernzerhof M. The shell model for the exchange-correlation hole in the strong-correlation limit. J. Chem. Phys. 2016, 145, 124104..10.1063/1.4962738 [PubMed] [Cross Ref]
  • Vuckovic S.; Irons T.; Savin A.; Teale A.; Gori-Giorgi P. Exchange-correlation functionals via local interpolation along the adiabatic connection. J. Chem. Theory Comput. 2016, 12, 2598–2610.10.1021/acs.jctc.6b00177 [PubMed] [Cross Ref]
  • Vuckovic S.; Irons T.; Wagner L.; Teale A.; Gori-Giorgi P. Interpolated energy densities, correlation indicators and lower bounds from approximations to the strong coupling limit of DFT. Phys. Chem. Chem. Phys. 2017, 19, 6169–6183.10.1039/C6CP08704C [PubMed] [Cross Ref]
  • Vuckovic S.; Gori-Giorgi P. Simple fully non-local density functionals for the electronic repulsion energy. J. Phys. Chem. Lett. 2017, 8, 2799..10.1021/acs.jpclett.7b01113 [PubMed] [Cross Ref]
  • Perdew J. P.; Schmidt K. In Density Functional Theory and Its Application to Materials; Van Doren V., editor. , et al. , Eds.; AIP Press: Melville, NY, 2001.
  • Seidl M. Strong-interaction limit of density-functional theory. Phys. Rev. A: At., Mol., Opt. Phys. 1999, 60, 4387..10.1103/PhysRevA.60.4387 [Cross Ref]
  • Seidl M.; Perdew J. P.; Levy M. Strictly correlated electrons in density-functional theory. Phys. Rev. A: At., Mol., Opt. Phys. 1999, 59, 51..10.1103/PhysRevA.59.51 [Cross Ref]
  • Seidl M.; Gori-Giorgi P.; Savin A. Strictly correlated electrons in density-functional theory: A general formulation with applications to spherical densities. Phys. Rev. A: At., Mol., Opt. Phys. 2007, 75, 042511..10.1103/PhysRevA.75.042511 [Cross Ref]
  • Harris J. Adiabatic-connection approach to Kohn-Sham theory. Phys. Rev. A: At., Mol., Opt. Phys. 1984, 29, 1648..10.1103/PhysRevA.29.1648 [Cross Ref]
  • Langreth D. C.; Perdew J. P. Exchange-correlation energy of a metallic surface: Wave-vector analysis. Phys. Rev. B 1977, 15, 2884..10.1103/PhysRevB.15.2884 [Cross Ref]
  • Gunnarsson O.; Lundqvist B. I. Exchange and correlation in atoms, molecules, and solids by the spin-density-functional formalism. Phys. Rev. B 1976, 13, 4274..10.1103/PhysRevB.13.4274 [Cross Ref]
  • Seidl M.; Perdew J. P.; Kurth S. Simulation of all-order density-functional perturbation theory, using the second order and the strong-correlation limit. Phys. Rev. Lett. 2000, 84, 5070..10.1103/PhysRevLett.84.5070 [PubMed] [Cross Ref]
  • Lewin M.. Semi-classical limit of the Levy-Lieb functional in Density Functional Theory. arXiv preprint arXiv:1706.02199, 2017, https://arxiv.org/abs/1706.02199.
  • Cotar C.; Friesecke G.; Klüppelberg C.Smoothing of transport plans with fixed marginals and rigorous semiclassical limit of the Hohenberg-Kohn functional. arXiv preprint arXiv:1706.05676, 2017, https://arxiv.org/abs/1706.05676.
  • Gori-Giorgi P.; Vignale G.; Seidl M. Electronic zero-point oscillations in the strong-interaction limit of density functional theory. J. Chem. Theory Comput. 2009, 5, 743–753.10.1021/ct8005248 [PubMed] [Cross Ref]
  • Levy M. Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem. Proc. Natl. Acad. Sci. U. S. A. 1979, 76, 6062–6065.10.1073/pnas.76.12.6062 [PubMed] [Cross Ref]
  • Lieb E. H. A lower bound for Coulomb energies. Phys. Lett. A 1979, 70, 444–446.10.1016/0375-9601(79)90358-X [Cross Ref]
  • Lieb E. H.; Oxford S. Improved lower bound on the indirect Coulomb energy. Int. J. Quantum Chem. 1981, 19, 427–439.10.1002/qua.560190306 [Cross Ref]
  • Lieb E. H. Density functionals for Coulomb systems. Int. J. Quantum Chem. 1983, 24, 243–277.10.1002/qua.560240302 [Cross Ref]
  • Buttazzo G.; De Pascale L.; Gori-Giorgi P. Optimal-transport formulation of electronic density-functional theory. Phys. Rev. A: At., Mol., Opt. Phys. 2012, 85, 062502..10.1103/PhysRevA.85.062502 [Cross Ref]
  • Mendl C. B.; Lin L. Kantorovich dual solution for strictly correlated electrons in atoms and molecules. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 125106..10.1103/PhysRevB.87.125106 [Cross Ref]
  • Malet F.; Gori-Giorgi P. Strong correlation in Kohn-Sham density functional theory. Phys. Rev. Lett. 2012, 109, 246402..10.1103/PhysRevLett.109.246402 [PubMed] [Cross Ref]
  • Mendl C. B.; Malet F.; Gori-Giorgi P. Wigner localization in quantum dots from Kohn-Sham density functional theory without symmetry breaking. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 125106..10.1103/PhysRevB.89.125106 [Cross Ref]
  • Lani G.; Di Marino S.; Gerolin A.; van Leeuwen R.; Gori-Giorgi P. The adiabatic strictly-correlated-electrons functional: kernel and exact properties. Phys. Chem. Chem. Phys. 2016, 18, 21092–21101.10.1039/C6CP00339G [PubMed] [Cross Ref]
  • Cort L.; Karlsson D.; Lani G.; van Leeuwen R. Time-dependent density-functional theory for strongly interacting electrons. Phys. Rev. A: At., Mol., Opt. Phys. 2017, 95, 042505..10.1103/PhysRevA.95.042505 [Cross Ref]
  • Colombo M.; Di Marino S.Annali di Matematica Pura ad Applicata; Springer: Berlin Heidelberg, 2013; pp 1–14.
  • Seidl M.; Di Marino S.; Gerolin A.; Nenna L.; Giesbertz K. J.; Gori-Giorgi P.The strictly-correlated electron functional for spherically symmetric systems revisited. arXiv preprint arXiv:1702.05022, 2017, https://arxiv.org/abs/1702.05022.
  • Colombo M.; De Pascale L.; Di Marino S. Multimarginal optimal transport maps for 1-dimensional repulsive costs. Canad. J. Math 2015, 67, 350–368.10.4153/CJM-2014-011-x [Cross Ref]
  • Wagner L. O.; Stoudenmire E. M.; Burke K.; White S. R. Phys. Chem. Chem. Phys. 2012, 14, 8581..10.1039/c2cp24118h [PubMed] [Cross Ref]
  • Helbig N.; Fuks J. I.; Casula M.; Verstraete M. J.; Marques M.; Tokatly I.; Rubio A. Density functional theory beyond the linear regime: Validating an adiabatic local density approximation. Phys. Rev. A: At., Mol., Opt. Phys. 2011, 83, 032503..10.1103/PhysRevA.83.032503 [Cross Ref]
  • Malet F.; Mirtschink A.; Giesbertz K.; Wagner L.; Gori-Giorgi P. Exchange-correlation functionals from the strong interaction limit of DFT: applications to model chemical systems. Phys. Chem. Chem. Phys. 2014, 16, 14551–14558.10.1039/C4CP00407H [PubMed] [Cross Ref]
  • Vuckovic S.; Wagner L. O.; Mirtschink A.; Gori-Giorgi P. Hydrogen Molecule Dissociation Curve with Functionals Based on the Strictly Correlated Regime. J. Chem. Theory Comput. 2015, 11, 3153–3162.10.1021/acs.jctc.5b00387 [PubMed] [Cross Ref]
  • Cohen A. J.; Mori-Sánchez P. Landscape of an exact energy functional. Phys. Rev. A: At., Mol., Opt. Phys. 2016, 93, 042511..10.1103/PhysRevA.93.042511 [Cross Ref]
  • Mori-Sánchez P.; Cohen A. J.Exact density functional obtained via the Levy constrained search. arXiv preprint arXiv:1709.10284, 2017, https://arxiv.org/abs/1709.10284.
  • Chen H.; Friesecke G. Pair densities in density functional theory. Multiscale Model. Simul. 2015, 13, 1259–1289.10.1137/15M1014024 [Cross Ref]
  • Gori-Giorgi P.; Seidl M.; Vignale G. Density-functional theory for strongly interacting electrons. Phys. Rev. Lett. 2009, 103, 166402..10.1103/PhysRevLett.103.166402 [PubMed] [Cross Ref]
  • Carr W. Jr Energy, specific heat, and magnetic properties of the low-density electron gas. Phys. Rev. 1961, 122, 1437..10.1103/PhysRev.122.1437 [Cross Ref]

Articles from ACS AuthorChoice are provided here courtesy of American Chemical Society