Energy Density Functionals From the Strong-Coupling Limit Applied to the Anions of the He Isoelectronic Series

# Energy Density Functionals From the Strong-Coupling Limit Applied to the Anions of the He Isoelectronic Series

## Abstract

Anions and radicals are important for many applications including environmental chemistry, semiconductors, and charge transfer, but are poorly described by the available approximate energy density functionals. Here we test an approximate exchange-correlation functional based on the exact strong-coupling limit of the Hohenberg-Kohn functional on the prototypical case of the He isoelectronic series with varying nuclear charge , which includes weakly bound negative ions and a quantum phase transition at a critical value of , representing a big challenge for density functional theory. We use accurate wavefunction calculations to validate our results, comparing energies and Kohn-Sham potentials, thus also providing useful reference data close to and at the quantum phase transition. We show that our functional is able to bind H and to capture in general the physics of loosely bound anions, with a tendency to strongly overbind that can be proven mathematically. We also include corrections based on the uniform electron gas which improve the results.

Density Functional Theory, Energy Density Functional, Exchange-Correlation, Strong Interaction Limit, Strong Coupling Limit, Strictly Correlated Electrons, Helium Isoelectronic Series

## I Introduction

Density functional theory (DFT),Hohenberg and Kohn (1964) in its Kohn-Sham (KS) formulation,Kohn and Sham (1965) has been a real breakthrough for electronic structure calculations. The key idea of KS DFT is an exact mappingKohn and Sham (1965) between the physical, interacting, many-electron system and a model system of non-interacting fermions with the same density, allowing for a realistic treatment of the electronic kinetic energy. All the complicated many-body effects are embedded in the so-called exchange-correlation (xc) energy functional. Although, in principle, the exact xc functional is unique (or “universal”), in practice a large number of approximations has been developed in the last thirty years, often targeting different systems, different properties, and different phenomena. Common practice for DFT users is nowadays to consult the (rather extensive) benchmark literature to choose the approximate xc functional most suitable for the problem at hand. This reflects the intrinsic difficulty of building a general approximation able to recognize and capture, for each class of systems or process, the many-body effects relevant for its description.

Even in this “specialized-functional” world, there are still important cases in which state-of-the-art KS DFT encounters severe problems, which is why the quest for better xc functionals continues to be a very active research field (for a recent review, see, e.g., Ref. Cohen, Mori-Sánchez, and Yang, 2012). The most notable example is the treatment of near-degeneracy and strong-correlation effects, which involve rearrangement of electrons within partially filled levels. These effects appear in bond dissociation but also at equilibrium geometries, particularly for systems with and unsaturated shells, such as transition metals and actinides. Mott insulators and low-density nanodevices are other examples of strongly-correlated systems whose physics is not captured by the standard approximations. A key problem when dealing with strong (or “static”) correlation is that, similarly to unrestricted Hartree-Fock, approximate KS DFT tries to mimic the physics of strong correlation and near degeneracy with spin and spatial symmetry breaking, which in complex systems may occur erratically and can be very sensitive to the choice of functional.Cramer and Truhlar (2009) This easily leads to a wrong characterization of several properties and to discontinuous potential energy surfaces.Cramer and Truhlar (2009) Being able to capture strong electronic correlation within KS DFT without resorting to symmetry breaking is arguably one of the most important problems of electronic structure theory.Becke (2013a, b, c); Cohen, Mori-Sánchez, and Yang (2012); Cramer and Truhlar (2009)

The mainstream strategies to construct approximate functionals consist of making an ansatz for the dependence of the xc functional on the relevant “ingredients” such as the local density, the local density gradients, the KS kinetic orbital energy, the KS orbitals, etc.Perdew et al. (2005) The ansatz can be constructed in order to fulfill as many exact constraints as possible given the ingredients used.Perdew et al. (2005) Some authors also introduce a (sometimes very large) number of parameters to be fitted to a specific data set (for recent reviews, see, e.g., Refs. Cramer and Truhlar, 2009; Cohen, Mori-Sánchez, and Yang, 2012; Peverati and Truhlar, 2014).

In recent years, an exact piece of information on the exact exchange-correlation functional, namely the limit of infinite correlation,Seidl, Gori-Giorgi, and Savin (2007); Gori-Giorgi, Vignale, and Seidl (2009); Gori-Giorgi, Seidl, and Vignale (2009) has become available. The “strictly-correlated-electrons” (SCE) functional, that utilizes this information, has a highly non-local dependence on the density, but its functional derivative (yielding the KS potential) can be easily constructed via a rigorous and physically transparent shortcut.Malet and Gori-Giorgi (2012); Malet et al. (2013) The SCE functional becomes exact in the limit in which the electron-electron interaction dominates over the electronic kinetic energy, and it has been successfully applied to model low-density quantum wiresMalet and Gori-Giorgi (2012); Malet et al. (2013) and quantum dots.Mendl, Malet, and Gori-Giorgi (2014) In those systems, the SCE functional has been shown capable of capturing the physics of charge localization without introducing magnetic order or any other symmetry breaking. In other words, the SCE functional achieved what was often regarded as practically impossible: making non-interacting electrons behave as strongly-correlated ones, showing that restricted KS DFT with the appropriate functionals can yield results beyond mean-field theory.

Using the SCE functional to address chemical problems also seems very attractive. It provides a new, well defined, starting point to build approximate functionals, deeply different from mainstream approaches. The new ingredient here is the non-locality encoded in the SCE functional and potential, which can capture the physics that is missed by standard approximations. Chemistry, however, is more challenging for the SCE functional than low-density nanostructures, because the kinetic energy and the electron-electron repulsion often have similar importance. For example, in a stretched bond only the bonding electrons are strongly-correlated, while the others are not. Indeed, in a recent paper,Malet et al. (2014) it has been shown that KS SCE dissociates properly a single chemical bond without introducing symmetry breaking, but it overcorrelates in all other aspects. This evidently requires corrections to the SCE functional, which can be built either by including higher-order terms in the expansion at infinite coupling strengthGori-Giorgi, Vignale, and Seidl (2009) or by considering rigorous local and semilocal approximations.Malet et al. (2013)

Both low-density nanostructures and stretched bonds involve charge localization due to strong spatial correlations, which is, by definition, the case in which SCE tends asymptotically to the exact xc functional. To gain insight into the performance of the SCE functional for other classes of chemical systems, we consider here a conceptually simple problem in which electronic correlation plays a crucial role. Despite its simplicity the problem nonetheless is very challenging for both DFT and other approaches. This is the anions of the He isoelectronic series, described by the Hamiltonian (in Hartree atomic units used throughout the paper)

 ^H=−12∇21−12∇22−Zr1−Zr2+1r12, (1)

with . Accurate wavefunction calculationsBaker et al. (1990) have shown that when the nuclear charge is lowered and crosses a critical value, , a quantum phase transition occurs from a bound to an unbound two-electron system. Thus, with this simple hamiltonian we can explore a whole class of very loosely bound anions, including the quantum phase transition at .

As is well known, anions are problematic for state-of-the-art KS DFT. Standard approximations often yield a positive eigenvalue for the highest occupied molecular orbital (HOMO), corresponding to a quasi-bound state (or resonance) instead of a properly bound system. Often, in practice, anions are tackled by finite basis sets within approximate DFT. Estimates of electron affinities are then obtained by the energy difference ( being the number of electrons), ignoring the fact that the HOMO has a positive eigenvalue.Galbraith and Schäfer III (1996) In the complete basis set limit (here by inclusion of plane waves), a positive orbital eigenvalue would lead to an unbound electron extending over the entire space and, consequently, within the finite basis set, orbitals with positive orbital energies should not be occupied. In practice, however, convergence of these calculations can only be achieved if the orbital with positive eigenvalue is occupied, corresponding to an electron artificially bound by the finite basis, a procedure which has been criticized.Rösch and Trickey (1997)

It is worth mentioning that the failure of standard DFT approximations to bind anions properly is often attributed to the self-interaction error (SIE). However, despite being self-interaction free, the Hartree-Fock (HF) method fails for H, yielding a negative binding energy for the second electron in contradiction to experiment.Hotop and Lineberger (1975) Thus in this case, it is correlation that stabilizes the system, so SIE is not the only problem.

In this work we test the KS SCE functional for the hamiltonian of Eq. (1), focusing on the anions close to the quantum phase transition. We compare our results (including energies, densities and KS potentials) with those from a very accurate wavefunction treatment and from standard approximate xc functionals. We also consider local corrections to KS SCE, and we analyze some exact properties of the density and of the KS potential at .

## Ii Theoretical Methods

### ii.1 Variational Calculations from Accurate Wavefunctions

Very accurate energies of the He isoelectronic series with nuclear charge between one and ten,Freund, Huxtable, and Morgan (1984) and for weakly bound anions close to and at the quantum phase transition,Baker et al. (1990) have been obtained using basis functions that depend explicitly on the interelectronic coordinates. For the latter the wavefunction was a linear combination of 476 basis functions consisting of 244 modified Freund, Huxtable, and Morgan (1984) Frankowski-PekerisFrankowski and Pekeris (1966) basis functions , where

 ϕFPn,l,m,j(s,t,u)=sntlum(lns)je−s/2, (2)

and 232 FrankowskiFrankowski (1967) basis functions , where

 ϕFn,l,m,j(s,t,u)=sntlum(lns)j(ect±e−ct)e−s/2. (3)

Here and are flexible scaling parameters and , , and are the Hylleraas coordinates

 s=r1+r2,t=r2−r1,u=r12. (4)

The sign depends on whether is even or odd to assure the proper symmetry of the basis functions under exchange of the two electrons (). The powers , , , and are chosen to duplicate the first several leading terms in the behavior of the exact wavefunction of a helium-like ion near the 3-particle coalescence, which is given by the Fock expansion.Fock (1954, 1958); Morgan III (1986) This composite basis was used in Ref. Baker et al., 1990 to obtain compact and highly accurate representations of the wavefunction of the sole bound state of the helium isoelectronic sequence for values of between and 1. Table 1 shows the approximately optimal values of and used for several values of in the present paper.

### ii.2 Restricted KS DFT with local, semilocal, and hybrid functionals

We quickly review some basic aspects of Kohn-Sham density functional theory, as this helps in clarifying the concepts behind the less familiar SCE functional, introduced in the next subsection.

For any -electron system in the external potential , Hohenberg and Kohn (HK) have provenHohenberg and Kohn (1964) the existence of a “universal” density functional , which in Levy’s constrained minimization formalismLevy (1979) is

 F[ρ]=minΨ→ρ⟨Ψ|^T+^Vee|Ψ⟩, (5)

where “” means that the minimization is carried over all fermionic wavefunctions yielding the same one-electron density , so that the ground-state energy can be obtained by minimizing the energy functional

 E0=minρ{F[ρ]+∫vext(r)ρ(r)dr}. (6)

Since it is extremely difficult to construct approximations for that encode the fermionic nature of the electrons, Kohn and Sham have introduced another functional, the non-interacting KS kinetic energy functional,

 Ts[ρ]=minΨ→ρ⟨Ψ|^T|Ψ⟩, (7)

which defines a non-interacting system of fermions with the same density of the physical, interacting, one. The HK functional is then partitioned as

 F[ρ]=Ts[ρ]+U[ρ]+Exc[ρ], (8)

where is the classical Hartree energy, and the exchange-correlation energy is defined as the correction needed to make Eq. (8) exact. With the non-interacting KS system, the full minimization for the many-electron energy becomes equivalent to the solution of the KS one-particle equations

 [−12∇2+vKS(r)]ψi(r)=ϵiψi(r), (9)

with the KS potential given by the functional derivatives

 vKS(r) =vext(r)+δU[ρ]δρ(r)+δExc[ρ]δρ(r) ≡vext(r)+uH(r)+vxc(r), (10)

where is the Hartree potential and is the exchange-correlation potential. The density is obtained as , with the sum running over the occupied orbitals, and the KS equations are solved self-consistently.

The simplest approximation for is the local density approximation (LDA), defining the exchange-correlation energy as a functional of the local density alone. The next levels of refinement are the generalized-gradient approximations (GGA), obtained by including the gradient of the local density , and the meta-GGA functionals which use also the local Laplacian of the density and/or the local kinetic energy density . For the special case of the two-electron systems considered here, the Hartree-Fock method becomes equivalent to KS DFT with the exact exchange functional, as the non-local HF exchange potential reduces to a local one-body potential. For systems with higher electron number the non-local Hartree-Fock exchange can be transformed into a local potential via the optimized effective potential method, yielding a well defined orbital-dependent functional (called exact exchange). Hybrid functionals are obtained by mixing a fraction of single determinant exchange with GGA or metaGGA functionals, and are normally computed using a non-local potential, a treatment outside the KS framework.

### ii.3 Restricted KS DFT with the SCE functional

The HK functional of Eq. (5) and the KS kinetic energy functional of Eq. (7) can be seen as the values at and of a more general functional , in which the electronic interaction is rescaled by a coupling strength parameter ,

 Fλ[ρ]=minΨ→ρ⟨Ψ|^T+λ^Vee|Ψ⟩. (11)

The SCE functional, first introduced in the seminal work of Seidl and coworkers,Seidl (1999); Seidl, Perdew, and Levy (1999) is the strong-interaction limit, , of , which corresponds to minimizing the electron-electron repulsion alone for a given density ,

 VSCEee[ρ]=minΨ→ρ⟨Ψ|^Vee|Ψ⟩. (12)

It is the natural counterpart of the KS kinetic energy functional of Eq. (7). The functional describes the physical situation in which the electrons are perfectly correlated, so that the position of one of them fixes all the other positions (or, to be more precise, all the interparticle distances) via the so-called co-motion functions , .Seidl, Gori-Giorgi, and Savin (2007) The co-motion functions are highly non-local functionals of the density , satisfying the differential equationSeidl, Gori-Giorgi, and Savin (2007)

 ρ(fi(r))dfi(r)=ρ(r)dr, (13)

which can be derived from the constraint “” of Eq. (12).Seidl, Gori-Giorgi, and Savin (2007); Gori-Giorgi, Vignale, and Seidl (2009) They also obey group properties that ensure the indistinguishability of the electrons,

 f1(r)≡r,f2(r)≡f(r),f3(r)=f(f(r)),f4(r)=f(f(f(r))), ⋮f(f(…f(f(r))))N times=r. (14)

The minimizing -electron density in Eq. (12), which becomes a distribution in this limit,Seidl, Gori-Giorgi, and Savin (2007); Gori-Giorgi, Vignale, and Seidl (2009); Buttazzo, De Pascale, and Gori-Giorgi (2012); Cotar, Friesecke, and Klüppelberg (2013) is the strictly-correlated state:

 |ΨSCE(r1,r2,…,rN)|2=1N!∑℘∫drρ(r)Nδ(r1−f℘(1)(r))×δ(r2−f℘(2)(r))⋯δ(rN−f℘(N)(r)), (15)

where denotes a permutation of . Equations (13)-(14) together with the properties of the Dirac -function guarantee that . In terms of the co-motion functions, the SCE functional isSeidl, Gori-Giorgi, and Savin (2007); Mirtschink, Seidl, and Gori-Giorgi (2012)

 VSCEee[ρ]=12∫d3rρ(r)N∑i=21|r−fi(r)|, (16)

and its functional derivative

 vSCE(r)=δVSCEee[ρ]δρ(r) (17)

can be obtained from the equationMalet and Gori-Giorgi (2012); Malet et al. (2013)

 ∇vSCE(r)=−N∑i=2r−fi(r)|r−fi(r)|3, (18)

which has a simple physical meaning: as the position of one electron fixes all the relative distances, the net electron-electron repulsion acting on an electron at becomes a function of alone, and can be represented as the gradient of a one-body potential. Equation (18) is a very powerful shortcut to compute the functional derivative of the highly non-local functional of Eqs. (13)-(16).

The KS SCE approach consists of approximating the constrained minimization in the universal functional of Eq. (5) as the sum of two constrained minima

 F[ρ] ≈minΨ→ρ⟨Ψ|^T|Ψ⟩+minΨ→ρ⟨Ψ|^Vee|Ψ⟩ =Ts[ρ]+VSCEee[ρ]. (19)

Obviously, the minimizing wavefunction is different for and : for the former, it is usually a single Slater determinant, while for the latter it is the strictly-correlated state of Eq. (15). The key point here is that, for a given , is a well defined density functional, whose functional derivative can be easily computed via Eq. (16). Our total energy functional is then

 Extra open brace or missing close brace (20)

By varying with respect to the single-particle orbitals appearing in , one obtains the usual KS equations, so that Eq. (19) is completely equivalent to making the approximation

 Exc[ρ]≈VSCEee[ρ]−U[ρ], (21)

and thus .

The minimization of our energy density functional of Eq. (20) reduces then to solving the standard restricted KS equations self-consistently:

 [−12∇2+vSCE[ρ](r)+vext(r)]ψi(r)=ϵiψi(r). (22)

Notice that the self-consistent KS SCE total energy that we obtain in this way is always a lower bound to the exact one. In fact, since the minimum of a sum is always larger or equal than the sum of the minima, for the exact ground-state density we have

 F[ρ]+∫ρvext≥Ts[ρ]+VSCEee[ρ]+∫ρvext. (23)

This inequality becomes even stronger when we minimize the right-hand side by solving self-consistently the KS SCE equations. The SCE functional is also self-interaction free, as for any one-electron density. Thus, as and , the self-consistent KS SCE method will certainly bind all the anions of the He isoelectronic series that are physically bound, and its error will always be towards overbinding, providing a lower bound for , which, however, turns out to be not very tight (see Sec. III).

It is also useful to represent graphically the approximation made in KS SCE in terms of the standard adiabatic connection of KS DFTLangreth and Perdew (1975) (see Fig. 1). By denoting the minimizing wave function in Eq. (11), and by defining the indirect part of the electron-electron repulsion at coupling strength ,

 Wλ[ρ]=⟨Ψλ[ρ]|^Vee|Ψλ[ρ]⟩−U[ρ], (24)

one obtains the well-known exact formulaLangreth and Perdew (1975) for ,

 Exc[ρ]=∫10Wλ[ρ]dλ. (25)

In Fig. 1 we show, schematically, as a function of for a weakly- and a strongly-correlated system. The area between , the -axis and the vertical lines corresponding to and gives the exchange-correlation energy . The KS SCE approximates with its value at for all , and thus the exchange-correlation energy with the area of the rectangle limited by , the -axis, and the vertical lines corresponding to and . This is evidently a good approximation only when the system is very correlated.

Evaluating the co-motion functions in the general case is still an open problem, although progress has been madeMendl and Lin (2013) by using the dual Kantorovich formulation,Buttazzo, De Pascale, and Gori-Giorgi (2012) which allows one to evaluate and its functional derivative in a different way, bypassing the co-motion functions. In the special case of spherically symmetric densities, like the ones considered here, an explicit solution is knownSeidl, Gori-Giorgi, and Savin (2007) in terms of the function ,

 Ne(r)=∫r04πx2ρ(x)dx (26)

and its inverse . For a system, the two electrons in the SCE solution are always opposite to each other with respect to the nucleus (maximum angular correlation), at a relative angle . Their distances from the nucleus, and , are related by the single co-motion function

 f(r)=N−1e[2−Ne(r)]. (27)

Equations (26)-(27) clearly show the non-local dependence of on the density. The SCE potential is then simply obtained by integrating the spherically-symmetric equivalent of Eq. (18),

 v′SCE(r)=−1[r+f(r)]2, (28)

with boundary condition . Notice that has the correct asymptotic behavior of the Hartree plus xc potential, , since . This is true for the general -electron case also, since the correct asymptotic leading term can be similarly derivedSeidl, Gori-Giorgi, and Savin (2007) from Eq. (18).

### ii.4 Local corrections to the SCE functional

As can be expected from Fig. 1, KS SCE can capture correlation effects at all correlation regimes, but good quantitative accuracy is obtained only when correlation becomes very strong.Malet and Gori-Giorgi (2012); Malet et al. (2013); Mendl, Malet, and Gori-Giorgi (2014) In practice, however, the systems of interest in chemistry are in between the weak- and strong correlation regimes and it is thus desirable to improve the approximation of Eq. (19). Therefore, we consider the more general decomposition of the universal functional

 F[ρ]=Ts[ρ]+VSCEee[ρ]+Tc[ρ]+Vdee[ρ], (29)

where (kinetic correlation energy) is the difference between the true kinetic energy and the KS one,

 Tc[ρ]=⟨Ψλ=1[ρ]|^T|Ψλ=1[ρ]⟩−Ts[ρ], (30)

and (decorrelation energyGori-Giorgi, Seidl, and Vignale (2009); Liu and Burke (2009)) is the difference between the true electron-electron repulsion energy and the SCE value,

 Vdee[ρ]=⟨Ψλ=1[ρ]|^Vee|Ψλ=1[ρ]⟩−VSCEee[ρ]. (31)

Both corrections are evidently always positive. A simple way to construct the correcting term is to make a local density approximation, which can be defined as the correction that makes Eq. (29) exact when the density becomes uniform,

 TLDAc[ρ]+Vd,LDAee[ρ]=∫d3rρ(r){(tc[ρ(r)]+vdee[ρ(r)]}, (32)

where and are the kinetic correlation energy per particle and the electron-electron decorrelation energy per particle of the homogeneous electron gas (HEG) of density . They can be easily obtained as

 tc(ρ)+vdee(ρ)=ϵxc(ρ)−ϵSCE(ρ), (33)

where and are, respectively, the exchange-correlation energy per particle and the indirect part of the SCE interaction energy per particle for the HEG. The latter can be obtained by considering that in the external potential due to an infinite uniform background with positive charge density the minimum possible electron-electron repulsion is attained with the electrons localized at the sites of the bcc crystal with lattice parameter . A uniform electronic density is constructed by taking a linear superposition of all the possible origins and orientations of the crystal. In other words, in the simple uniform-density case, the co-motion functions are just the lattice vectors of the bcc crystal with origin in the reference electron, whose position is distributed uniformly. This means that for all values of the density parameter the SCE energy of the uniform electron gas is equal to the low-density leading term of the HEG energy,

 ϵSCE(ρ)=−d0rs(ρ). (34)

At high densities, the SCE energy is very far from the exact one, and at low densities it becomes asymptotically exact. This is also true, more generally, for the functional , which approaches the exact Hartree plus exchange correlation functional when the density is scaled as and . Here we have set , which is the value from the Perdew-Wang-92 LDA parametrization.Perdew and Wang (1992) We denote this method KS SCE+LDA.

It is also possible to consider the local correction only for the electron-electron repulsion part, assuming that the error made by the KS kinetic energy is, for these systems, less serious than the one made by the SCE functional, so that the correction needs to rebalance the two terms. This corresponds to taking as correction only

 Vd,LDAee[ρ]=∫d3rρ(r)vdee[ρ(r)], (35)

where is obtained by subtracting from Eq. (33) the kinetic correlation contribution . We call this approximation KS SCE+L.

## Iii Results

### iii.1 Accurate solution and exact properties at Zcrit

Before presenting and discussing the KS SCE results, we extend the work of Umrigar and GonzeUmrigar and Gonze (1994) by studying the accurate densities and KS potentials obtained from the wavefunctions of Sec. II.1 close to the quantum phase transition. The densities and KS potentials (obtained by inversion of the KS equations)Umrigar and Gonze (1994) for selected values of are shown in Fig. 2.

The exact density of an atomic or molecular system is known to decay (with exceptions when the ground-state of the ion is not asymptotically accessible by symmetry) asKatriel and Davidson (1980); Levy, Perdew, and Sahni (1984); Almbladh and von Barth (1985) , an expansion which is valid for , where is the ionization energy. When () the density remains compact, in agreement with the rigorous result of Ref. Hoffmann-Ostenhof, Hoffmann-Ostenhof, and Simon, 1983, where it has been proven that the density at satisfies

 C−(δ)r−3/2−δe−2[8(1−Zcrit)r]1/2≤ρ(r)≤C+(δ)r−3/2+δe−2[8(1−Zcrit)r]1/2, (36)

where is an arbitrary small positive number and are constants depending on .

We can further understand the asymptotic decay of the density at by studying the corresponding differential equationLevy, Perdew, and Sahni (1984) for (which for a singlet coincides with the KS equation). At the quantum phase transition with the asymptotic potential to fourth orderAlmbladh and von Barth (1985); Umrigar and Gonze (1994) this equation is

 [−12∇2r−Z−N+1r+O(1r4)]√ρ(r)=0. (37)

By solving Eq. (37) asymptotically (), we obtain, order by order, a solution for the leading terms to order ,

 ρ(r→∞)∼ e−4a√rr3/2(1+38ar1/2−3128a2r+151024a3r3/2−40532768a4r2), (38)

with . This decay agrees to leading order with Eq. (36). The accurate density at the quantum phase transition together with the decays from Eqs. (36) and (38) are displayed in Fig. 3, where in both cases the proportionality constant has been adjusted to match the accurate density at the end of the radial grid (). Notice that Eq. (37) implies that for the exact KS system (which yields the exact ground-state density) the equality also holds at , when .

From Fig. 2 we see that the KS potentials have a bump at intermediate length scale. This bump increases for smaller as can be expected from the asymptotic first order contribution at large , that will be positive for . The bump is present also for the Hydrogen anion, where this first order contribution vanishes.

In Fig. 4 we show the correlation potentials for selected values of . We see that, as was found in Ref. Umrigar and Gonze, 1994, the accurate correlation potential close to the nucleus has a nearly quadratic behavior. In Refs. Qian, 2007; Qian and Sahni, 2007 it has been shown that the linear term in the correlation potential is due to the kinetic contribution, which, thus, turns out to be very small.

### iii.2 Zcrit from KS-DFT with standard functionals and from KS SCE

For the He isoelectronic series with we solved self-consistently the restricted KS equations with various approximate functionals. Calculations for the HF (or exact exchange) method, KS LDA, KS SCE, and KS SCE with the two local corrections of Sec. II.4 were performed with a numerical code developed in our group. We chose the Perdew-Wang-92 functional (PW92)Perdew and Wang (1992) LDA parametrization. To compare our calculations with the available standard approximations we have further performed restricted KS-DFT calculations with the Amsterdam Density Functional package (ADF).te Velde et al. (2001); ?; ? From the GGA class of functionals we chose PBE,Perdew, Burke, and Ernzerhof (1996) from the metaGGA class the revTPSSTao et al. (2003) functional and for the hybrid functional we chose B3LYP.Becke (1993); Stephens et al. (1994); Lee, Yang, and Parr (1988); Vosko, Wilk, and Nusair (1980) If not mentioned otherwise, all ADF calculations were carried out in the even-tempered (ET) QZ3P basis supported by 3 diffuse s-functions with the parametrization of Hydrogen.Chong et al. (2004) To assess the quality of the basis set we also performed KS-LDA (PW92 functional) calculations with the ADF package and compare them to our numerical solution of the KS equations. To assess the quality of the basis set we also performed KS-LDA calculations with the ADF package (PW92 functional) and compare them to our numerical solution of the KS equations.

We define the critical nuclear charge for the various DFT approximations to be the value of at which either the ionization energy becomes smaller than or the HOMO eigenvalue becomes positive, whichever is larger. Although the equality does not hold in general for approximate functionals, we invoke the HOMO eigenvalue criterion to avoid the conceptual and numerical issue of occupying orbitals with a positive eigenvalue already discussed.

Table 2 shows the predicted for the quantum phase transition together with the corresponding ionization energy and the HOMO energies for the various approximations. Of the DFT approximations considered only the SCE functionals (SCE and SCE with local corrections) and the hybrid functional are able to bind the Hydrogen anion. The hybrid functional however, yields an unphysical description of the bound anion as we will further discuss below. Remarkably, all the standard functionals at different levels of approximation yield a similar value of . This shows that the nonlocality encoded in the SCE functional is able to capture different many-body effects than the standard approximations.

As already discussed in Sec. II.3, the KS SCE self-consistent results yield a lower bound to the total energy, which for these systems is not very tight as can be seen from the underestimation of versus the actual value . This is due to the inherent strong correlation nature of the electrons in the SCE formulation that results in underestimating the electron-electron repulsion energy. Self-consistently thus, the KS-SCE densities become quite compact until the kinetic energy starts to dominate in Eq. (19). This is manifested in Fig. 5 where the density of H is displayed for several methods, the KS-SCE yielding the most compact density. Physically, this is due to the fact that the two electrons, being perfectly correlated, can avoid each other as much as possible and can get much closer to the nucleus to lower the total energy.

The local corrections to the SCE functional improve considerably the predicted and give a more realistic description of the electronic interactions. We observe that the KS-SCE+LDA density is too spread out compared to the accurate data (e.g. Fig. 5). This can be attributed to the self-interaction error that is introduced by the LDA correction which is obvious from Eq. (33) – the energy densities do not vanish for a density integrating to 1.

In Fig. 5 we display also the Hartree-Fock density. It is possible to do this because the HOMO eigenvalue is negative even though . (If the electron number were treated as a variational parameter, the minimum energy would be attained for .) We see that the HF density resembles the accurate density more closely than the density from other functionals considered. This supports the point of view of Ref. Kim, Sim, and Burke, 2013, and the general idea of using HF densities as input for DFT energies in the case of negative ions,Lee, Furche, and Burke (2010); Kim, Sim, and Burke (2011) even when HF does not bind the last electron.

For the hybrid functional in the ET-QZ3P+3diffuse basis we obtain a negative and for H. Formally the hybrid thus binds the Hydrogen anion. When inspecting the density however, one observes that it escapes partially from the nucleus, as shown in Fig. 7. When removing the 3 diffuse basis functions from the basis set to prevent the density accumulation in the outside regions, we obtain a value of in between that from HF and conventional DFT, as expected.

We now discuss the Kohn-Sham and exchange-correlation potentials for the self-consistent densities, displayed in Figs. 6 and 8.

We see that the SCE total Kohn-Sham potential does not develop the bump for , but only for smaller nuclear charges when the interelectronic repulsion dominates over the weaker nuclear attraction (Fig. 6). As already observed in Fig. 5, this corresponds to a very compact density. For larger distances, the SCE potential is in good agreement with the accurate one, as expected from the absence of the self-interaction error in the SCE. From Fig. 8 we also see that the SCE potential is quadratic close to the nucleus, as can be easily proven analytically from Eq. (28), since when we have , so that . This is in agreement with the findings of Refs. Qian, 2007; Qian and Sahni, 2007, as there is no kinetic contribution in the SCE potential.

Although the SCE functional approximates exchange and correlation together, in Fig. 9 we show the SCE correlation potential alone, obtained by subtracting from the xc SCE potential the exchange potential constructed from the self-consistent KS SCE densities. We see that the SCE correlation potential is always negative, in contrast to the exact one. The positive part of the exact correlation potential is mainly due to kinetic correlation effectsBuijse, Baerends, and Snijders (1989); Gritsenko, van Leeuwen, and Baerends (1997) that are missed in the bare SCE.

At least qualitatively, the bump for H in the total KS potential is captured by the KS-SCE with the two local corrections, though the bump is too pronounced, particularly in KS SCE+LDA. This is also responsible for the overestimation of and can be partially attributed to the self-interaction error. However, the self-interaction error present in the KS-SCE+LDA approach is substantially different from the self-interaction error in standard KS-LDA or KS-GGA. In KS-LDA and GGA the self-interaction error manifests in the wrong asymptotic decay of the KS potential ( instead of ). KS-SCE has the correct decay and this is not altered by the exponentially vanishing LDA contribution upon going from KS-SCE to KS-SCE-LDA. The KS SCE+L is more attractive at short distance than the exact KS potential, achieving error compensation with the overestimation of the bump (less severe than in the KS-SCE+LDA method), which results in a good estimate for . Of the methods studied, the KS-SCE approach with the local corrections is the one in which the HOMO energy deviates the least from the corresponding (see Tab. 2).

The HF (or exact exchange) potential is also shown in Fig. 8, although, once more, we have to keep in mind that in this case , so that the system is not really physically bound.

### iii.3 Fractional Electron Numbers at Z=1

We complete our analysis by also allowing for fractional electron numbers , with , in the Hydrogen nuclear potential, which is often considered a paradigmatic model for a Mott insulator.Mori-Sanchez, Cohen, and Yang (2009) In exact KS DFT, it is known that the HOMO eigenvalue should be constant between any two adjacent integer electron numbers (say, and ), equal to the negative of the exact, interacting, ionization energy , and should jump whenever an integer electron number is crossed.Perdew et al. (1982); Sagvolden and Perdew (2008)

The KS DFT results with the standard functionals at fractional electron numbers can be easily obtained by giving fractional occupation to the HOMO.Vydrov, Scuseria, and Perdew (2007); Gaiduk, Firaha, and Staroverov (2012) As discussed in the Introduction, here we consider the challenging case of the restricted KS method, where, for singlet systems, as we increase the occupancy of the HOMO orbital we should observe a jump in its energy at . Notice that in the restricted KS method the conditions regarding the spin degree of freedomMori-Sanchez, Cohen, and Yang (2009) are automatically fulfilled, so that the gap at is the same as the “Mott gap” for spin-up and spin-down electrons.

The KS SCE method needs, additionally, the construction of the SCE functional for fractional electron numbers. This has been rigorously done in Ref. Mirtschink, Seidl, and Gori-Giorgi, 2013, and, in this case, corresponds to setting the co-motion function of Sec. II.3 equal to

 f(r)={N−1e[2−Ne(r)]r>N−1e(2−Q)∞otherwise. (39)

The physical meaning of Eq. (39) is very simple: the two electronic positions are always separated by a radial distance such that the density integrates to 1 (total suppression of fluctuations),

 ∫f(r)r4πx2ρ(x)dx=1, (40)

so that for densities integrating to less than 2 there are values of for which the second electron “cannot enter” in the density.Mirtschink, Seidl, and Gori-Giorgi (2013)

Figure 10 displays for various approaches in the restricted KS scheme. As observed before,Mirtschink, Seidl, and Gori-Giorgi (2013) the SCE functional shows a vertical change in the HOMO energy even in the restricted KS approach. A sharp step, however, is only obtained with KS SCE in the extremely strong correlation (or low-density) limit,Mirtschink, Seidl, and Gori-Giorgi (2013) from which H is still far. KS-SCE with the two local corrections exhibits the same smoothed step, but the self-interaction error leads to a non-constant HOMO energy between . The HF curve we report here has been obtained by keeping the occupancies of the two electrons equal at all . This is what it should be compared in the restricted case, and it is the situation encountered in restricted HF when stretching a bond or expanding a lattice.Mori-Sanchez, Cohen, and Yang (2009)

Finally, Fig. 10 allows for a determination of the maximum number of electrons bound by the conventional DFT approaches. The results are compiled in Tab. 3. We observe, similarly to , that the predicted value of is insensitive to the level of approximation of the standard functionals, further supporting the idea behind the model potential of Ref. Gaiduk, Firaha, and Staroverov, 2012.

## Iv Conclusions and Perspectives

We have applied functionals based on the exact strong-coupling limit of DFT to the loosely bound negative ions of the He isoelectronic series, which are a prototypical case for the delicate physics of anions and radicals. Whereas standard DFT functionals either do not bind anions or bind them with unphysical long-range features in the charge density, the functionals based on the strictly-correlated-electrons have a rigorous tendency to overbind that can be mitigated by local corrections. This shows that the SCE functional and its corrections are able to capture many-body effects radically different than the ones described by the standard functionals, although improvements are still needed. In particular, one should aim at building corrections based on correlation kinetic energy effects Malet et al. (2014) and/or on exact exchange.Malet et al. (2014)

Besides improving the accuracy of the functionals based on SCE, the challenge for the future is also to implement SCE physics into routinely applicable approximations. This can be done by either developing algorithms to evaluate the exact SCE functional exploiting its formal similarity to an optimal transport problem,Buttazzo, De Pascale, and Gori-Giorgi (2012); Cotar, Friesecke, and Klüppelberg (2013) as in the pilot implementation of Ref. Mendl and Lin, 2013, or by constructing new approximations based on the idea of co-motion functions, i.e., by trying to build approximate and simplified co-motion functions. These, in turn, could be used in a local interpolation along the adiabatic connection that preserves size consistency.Mirtschink, Seidl, and Gori-Giorgi (2012)

Finally, our study also provides reference data for the anions of the He isoelectronic series close to and at the quantum phase transition that can be valuable to test the accuracy of new DFT approximations (see, e.g., Ref. Bleiziffer et al., 2013 which presents correlation potentials from RPA approaches that are good approximations to the true correlation potential).

###### Acknowledgements.
We are very grateful to Erik van Lenthe for his assistance with the ADF package, and to A. J. Cohen and P. Mori-Sanchez for pointing out an error in Fig. 10 concerning the Hartree-Fock curve. This work was supported by the Netherlands Organization for Scientific Research (NWO) through a Vidi grant and by the NSF through grant DMR-0908653.

### Footnotes

1. Not supported by ADFte Velde et al. (2001); ?; ?
2. ET-QZ3P basis
3. by integrating over the inside region () in Fig. 7.

### References

1. P. Hohenberg and W. Kohn, Phys. Rev. 136, B 864 (1964).
2. W. Kohn and L. J. Sham, Phys. Rev. 140, A 1133 (1965).
3. A. J. Cohen, P. Mori-Sánchez,  and W. Yang, Chem. Rev. 112, 289 (2012).
4. C. J. Cramer and D. G. Truhlar, Phys. Chem. Chem. Phys. 11, 10757 (2009).
5. A. D. Becke, J. Chem. Phys. 138, 074109 (2013a).
6. A. D. Becke, J. Chem. Phys. 138, 161101 (2013b).
7. A. D. Becke, J. Chem. Phys. 139, 021104 (2013c).
8. J. P. Perdew, A. Ruzsinszky, J. Tao, V. N. Staroverov, G. E. Scuseria,  and G. I. Csonka, J. Chem. Phys. 123, 062201 (2005).
9. R. Peverati and D. G. Truhlar, Phil. Trans. R. Soc. A 372, 20120476 (2014).
10. M. Seidl, P. Gori-Giorgi,  and A. Savin, Phys. Rev. A 75, 042511 (2007).
11. P. Gori-Giorgi, G. Vignale,  and M. Seidl, J. Chem. Theory Comput. 5, 743 (2009).
12. P. Gori-Giorgi, M. Seidl,  and G. Vignale, Phys. Rev. Lett. 103, 166402 (2009).
13. F. Malet and P. Gori-Giorgi, Phys. Rev. Lett. 109, 246402 (2012).
14. F. Malet, A. Mirtschink, J. C. Cremon, S. M. Reimann,  and P. Gori-Giorgi, Phys. Rev. B 87, 115146 (2013).
15. C. B. Mendl, F. Malet,  and P. Gori-Giorgi, Phys. Rev. B 89, 125106 (2014).
16. F. Malet, A. Mirtschink, K. J. H. Giesbertz, L. O. Wagner,  and P. Gori-Giorgi, Phys. Chem. Chem. Phys. , accepted (2014).
17. J. D. Baker, D. E. Freund, R. N. Hill,  and J. D. Morgan III, Phys. Rev. A 41, 1247 (1990).
18. J. M. Galbraith and H. F. Schäfer III, J. Chem. Phys. 105, 862 (1996).
19. N. Rösch and S. B. Trickey, J. Chem. Phys. 106, 8940 (1997).
20. H. Hotop and W. C. Lineberger, J. Chem. Phys. Ref. Data 4, 539 (1975).
21. D. E. Freund, B. D. Huxtable,  and J. D. Morgan, Phys. Rev. A 29, 980 (1984).
22. K. Frankowski and C. L. Pekeris, Phys. Rev. 146, 46 (1966).
23. K. Frankowski, Phys. Rev. 160, 1 (1967).
24. V. A. Fock, Izv. Akad. Nauk SSSR, Ser. Fiz. 18, 161 (1954).
25. V. A. Fock, Det Koneglige Norske Videnskabers Selskabs Forhandlinger 31, 138 (1958).
26. J. D. Morgan III, Theor. Chim. Acta 69, 181 (1986).
27. M. Levy, Proc. Natl. Acad. Sci. U.S.A. 76, 6062 (1979).
28. M. Seidl, Phys. Rev. A 60, 4387 (1999).
29. M. Seidl, J. P. Perdew,  and M. Levy, Phys. Rev. A 59, 51 (1999).
30. G. Buttazzo, L. De Pascale,  and P. Gori-Giorgi, Phys. Rev. A 85, 062502 (2012).
31. C. Cotar, G. Friesecke,  and C. Klüppelberg, Comm. Pure Appl. Math. 66, 548 (2013).
32. A. Mirtschink, M. Seidl,  and P. Gori-Giorgi, J. Chem. Theory Comput. 8, 3097 (2012).
33. D. C. Langreth and J. P. Perdew, Solid State Commun. 17, 1425 (1975).
34. C. B. Mendl and L. Lin, Phys. Rev. B 87, 125106 (2013).
35. Z. F. Liu and K. Burke, J. Chem. Phys. 131, 124124 (2009).
36. J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
37. C. J. Umrigar and X. Gonze, Phys. Rev. A 50, 3827 (1994).
38. J. Katriel and E. R. Davidson, Proc. Natl. Acad. Sci. USA 77, 4403 (1980).
39. M. Levy, J. P. Perdew,  and V. Sahni, Phys. Rev. A 30, 2745 (1984).
40. C.-O. Almbladh and U. von Barth, Phys. Rev. B 31, 3231 (1985).
41. M. Hoffmann-Ostenhof, T. Hoffmann-Ostenhof,  and B. Simon, J. Phys. A 16, 1125 (1983).
42. Z. Qian, Phys. Rev. B 75, 193104 (2007).
43. Z. Qian and V. Sahni, Phys. Rev. A 75, 032517 (2007).
44. G. te Velde, F. M. Bickelhaupt, S. J. A. van Gisbergen, C. F. Guerra, E. J. Baerends, J. G. Snijders,  and T. Ziegler, J. Comp. Chem. 22, 931 (2001).
45. C. F. Guerra, J. G. Snijders, G. te Velde,  and E. J. Baerends, Theor. Chem. Acc. 99, 391 (1998).
46. ADF2013, SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam,  and T. Netherlands, http://www.scm.com.
47. J. P. Perdew, K. Burke,  and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
48. J. Tao, J. P. Perdew, V. N. Staroverov,  and G. E. Scuseria, Phys. Rev. Lett. 91, 146401 (2003).
49. A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
50. P. J. Stephens, F. J. Devlin, C. F. Chabalowski,  and M. J. Frisch, J. Phys. Chem. 98, 11623 (1994).
51. C. Lee, W. Yang,  and R. G. Parr, Phys. Rev. B 37, 785 (1988).
52. S. J. Vosko, L. Wilk,  and M. Nusair, Can. J. Phys. 58, 1200 (1980).
53. D. P. Chong, E. van Lenthe, S. J. A. van Gisbergen,  and E. J. Baerends, J. Comp. Chem. 25, 1030 (2004).
54. M.-C. Kim, E. Sim,  and K. Burke, Phys. Rev. Lett. 111, 073003 (2013).
55. D. Lee, F. Furche,  and K. Burke, J. Phys. Chem. Lett. 1, 2124 (2010).
56. M.-C. Kim, E. Sim,  and K. Burke, J. Chem. Phys. 134, 171103 (2011).
57. M. A. Buijse, E. J. Baerends,  and J. G. Snijders, Phys. Rev. A 40, 4190 (1989).
58. O. V. Gritsenko, R. van Leeuwen,  and E. J. Baerends, J. Chem. Phys. 104, 8535 (1997).
59. P. Mori-Sanchez, A. J. Cohen,  and W. Yang, Phys. Rev. Lett. 102, 066403 (2009).
60. J. P. Perdew, R. G. Parr, M. Levy,  and J. L. Balduz, Phys. Rev. Lett. 49, 1691 (1982).
61. E. Sagvolden and J. P. Perdew, Phys. Rev. A 77, 012517 (2008).
62. O. A. Vydrov, G. E. Scuseria,  and J. P. Perdew, J. Chem. Phys. 126, 154109 (2007).
63. A. P. Gaiduk, D. S. Firaha,  and V. N. Staroverov, Phys. Rev. Lett. 108, 253005 (2012).
64. A. Mirtschink, M. Seidl,  and P. Gori-Giorgi, Phys. Rev. Lett. 111, 126402 (2013).
65. P. Bleiziffer, A. Heßelmann, C. J. Umrigar,  and A. Görling, Phys. Rev. A 88, 042513 (2013).
100752