Communication: Strong-interaction limit of an adiabatic connection in Hartree-Fock theory

The adiabatic connection formalism has always been a powerful theoretical tool to build approximations for the exchange-correlation (XC) energy in density functional theory (DFT), in particular, for hybrid^1–3 and double-hybrid functionals,^4–7 but also for other kinds of XC functionals.^8–15 In the Kohn-Sham (KS) framework, the density fixed adiabatic connection can be defined via the λ-dependent hamiltonian^16,17

where $T^$ is the kinetic energy operator for the N electrons, $V^ee$ is their mutual Coulomb repulsion, and $V^λ[ρ]=∑i=1Nvλ(ri,[ρ])$ is the one body potential that makes the ground-state wavefunction of Eq. (1), $ΨλDFT$⁠, yield the density ρ(r) ≡ ρ_λ=1(r) for all values of λ. From Eq. (1), one can derive an exact formula for the KS DFT XC energy^16,17

where

with U[ρ] being the Hartree energy. The coupling-constant integrand of Eq. (3) has the known small-¹⁸ and large-λ expansions¹⁹ 

where $ExDFT$ is the exact KS exchange energy [the same expression as in Hartree-Fock (HF) theory, but using KS orbitals] and $EcGLn$ is the nth term in the Görling-Levy perturbation series.^18,20 The expansion for large λ of Eq. (5) has as leading term the functional $W∞DFT[ρ]$⁠, given by the minimum possible expectation value of the electron-electron repulsion in a given density ρ(r)^21–23 

while the next leading term, determined by $W∞′DFT[ρ]$⁠, corresponds to the potential energy of zero-point vibrations around the manifold corresponding to the support of the minimizing probability density in Eq. (6).¹⁹ While for the leading term there are rigorous proofs,^22,23 this next term is only a very reasonable conjecture that has been confirmed numerically in simple cases where it was possible to compute the exact integrand $WλDFT[ρ]$⁠.^24,25

Mixing KS DFT with Hartree-Fock (HF) ingredients is an approximation strategy that has a long history in chemistry, already starting with hybrids^{1–3,26–29} and double hybrids,^4–7 but also by simply inserting the HF density into a given approximate XC density functional.^30–35 Very recently, it has also been observed that rather accurate interaction energies,^36,37 particularly for non-covalent complexes,³⁸ can be obtained from models for $WλDFT[ρ]$ that interpolate between the two limits of Eq. (4)—retaining only the first term, GL2, in the GL series—and of Eq. (5), using HF densities and orbitals as input, i.e., by constructing de facto an approximate resummation of the Møller-Plesset (MP) series, a procedure that lacks so far a theoretical justification. Motivated, in particular, by these last findings, we analyze in this communication the Hartree-Fock adiabatic connection [see Eq. (7)] whose Taylor expansion around λ = 0 is the MP series [see Eq. (10)] and show that the leading term in the λ → ∞ expansion is determined by a functional of the HF density [see Eqs. (14) and (15)]. We also highlight similarities and differences with the DFT case, showing that the large λ expansion in HF theory has a structure similar to the one of Eq. (5).

We keep the notation general, as only few key properties of the HF operators are important here. We consider the adiabatic connection (see, e.g., Ref. 39)

with $V^ext$ being the (nuclear) external potential and Ĵ = Ĵ[ρ^HF] and $K^=K^[{ϕiHF}]$ being the standard HF Coulomb and exchange operators, which are fixed once for all in the initial HF calculation, and do not depend on λ, but only on the HF density ρ^HF and occupied HF orbitals ${ϕiHF}$⁠. In the ground state $ΨλHF$ of $ĤλHF$⁠, the density ρ_λ(r) changes with λ: ρ_λ=0(r) is the HF density ρ^HF(r), and ρ_λ=1(r) is the exact physical density ρ(r). Note that Teale et al.⁴⁰ have analyzed a related adiabatic connection, in which the external potential is kept fixed; in that framework in the limit λ → ∞, all the electrons but one escape to infinity. From Eq. (7), the Hellmann-Feynman theorem yields the exact formula

for the XC energy in the HF framework, with

Equation (9) has been defined to allow for a direct comparison with the DFT $WλDFT[ρ]$ of Eqs. (2) and (3), with $Wλ=0HF=ExHF$⁠, and for small λ

where $EcMPn$ the nth term in the MP series. As is well-known (see, e.g., Refs. 41 and 42), the radius of convergence of the MP series is typically smaller than 1. Here we ask the question: What happens to$ΨλHF$and$WλHF$as λ → ∞? After answering this theoretical question, we will discuss its actual relevance for constructing approximations.

When λ becomes very large, the term $λ(V^ee−Ĵ−K^)$ in Eq. (7) becomes more and more important, and we argue that the wavefunction $ΨλHF$ should end up minimizing this term alone, similarly to the DFT case²¹ of Eq. (6). The difference here is that the minimizer is not constrained to yield a fixed density and the operator to be minimized also contains $−Ĵ−K^$⁠. We further argue that the expectation value of $K^$ is subleading with respect to the one of $V^ee−Ĵ$⁠, i.e., we argue that

Before we shall support this conjecture with a variational argument, we discuss its consequences.

If Eq. (11) holds, then $ΨλHF$ for λ → ∞ ends up minimizing the even simpler operator $λ(V^ee−Ĵ)$

The “asymptotic hamiltonian” $H^∞HF=V^ee−Ĵ[ρHF]$ is completely specified by the HF density ρ^HF(r), since N = ∫drρ^HF(r) and $Ĵ[ρ]=∑i=1NvH(ri;[ρ])$⁠, with $vH(r;[ρ])≡∫ρ(r′)|r−r′|dr′$⁠. Consequently, also the minimizer in Eqs. (12) and (13) is specified solely by ρ^HF

and the minimum in Eq. (13) is a functional of ρ^HF

The minimizer in Eq. (12) could be not unique, but this does not affect the value of the minimum, which is the object of the present investigation. The functional $Eel[ρ]=minΨ⟨Ψ|V^ee−Ĵ[ρ]|Ψ⟩+U[ρ]$ has a simple classical interpretation: Since $H^∞HF=V^ee−Ĵ[ρHF]$ is a purely multiplicative operator

the square modulus $|Ψ∞HF|2$ of its minimizing wave function is a distribution in $R3N$ that is zero wherever $H^∞HF$ as a function of r₁, …, r_N does not assume its global minimum (if it were otherwise, it would not be optimal as we could always lower the energy by increasing the weight of the wave function in the global minimum of $H^∞HF$⁠). In other words,

is the minimum total electrostatic energy of N equal classical point charges (−e) in a positive background with continuous charge density (+e)ρ(r). The term U[ρ], inherited from Eq. (9), represents the background-background repulsion.

Strictly speaking, the minimizer $Ψ∞HF$ is not in the space of allowed wavefunctions so that the minimum is actually an infimum, similarly to the DFT case.^22,23,43

Equations (14) and (15) comprise a central result of this work: they show that the strong-interaction limit of the HF adiabatic connection can be determined from a functional of the HF density, providing some theoretical justification for resumming the MP series by using a DFT-like expansion at large λ with functionals of ρ^HF,^36–38 although, as we will discuss, there are still several points to be addressed.

We now analyze the functionals $Ψ∞HF[ρ]$ and E_el[ρ], comparing them with the DFT case. As λ → ∞, $ĤλDFT$ of Eq. (1) tends to $λH^∞DFT[ρ]$⁠, with²¹ 

Comparing $H^∞DFT[ρ]$ with $H^∞HF[ρ]$ of Eq. (16), we see that both hamiltonians consist of the electron-electron repulsion operator and of an attractive one-body potential. In the HF case, the attractive potential is −$vH$(r, [ρ]), which is, for typical Hartree-Fock densities, strong enough to create a classical bound crystal. To be more precise, −$vH$(r, [ρ]) is more attractive than the one-body potential $v∞$(r, [ρ]). In fact, the potential $v∞$(r, [ρ]), which has been studied in several studies,^21,44–46 is generated by a charge that integrates to N − 1^21,45

while the attractive potential −$vH$(r, [ρ]) is generated by the given density ρ(r), which integrates to N. For finite systems, the state $Ψ∞HF[ρ]$ is thus more compact than the state $Ψ∞DFT[ρ]$⁠: this is due to the density constraint in the DFT adiabatic connection, which forces $Ψ∞DFT[ρ]$ to have the given quantum mechanical density ρ(r).^21,46

We note in passing that, for given occupied HF orbitals, we have the chain of inequalities

The first one, $W∞HF≤Eel[ρHF]$⁠, is trivial since $W∞HF=Eel[ρHF]+2ExHF$ and $ExHF≤0$⁠. The second inequality holds for any density ρ(r), $Eel[ρ]≤W∞DFT[ρ]$⁠. To prove it, we introduce the bifunctional $W[ρ,v]$

for which we have

and, from the dual formulation of $W∞DFT[ρ]$⁠,^21,47,48

which clearly completes the proof.

As promised, we provide a variational argument to support the assumption of Eq. (11), sketching the main points and leaving a more detailed treatment to a subsequent paper. We start by considering the global minimum $R̲min≡{r1min,…,rNmin}$ of the function $H^∞HF$ of Eq. (16) and construct the simple trial wavefunction

where $Gα(r)=α3/4π3/4e−α2|r|2$⁠, with α being a λ-dependent variational parameter that goes to infinity for large λ, α(λ) ∼ λ^q with q > 0. By construction, when α → ∞ (i.e., when λ → ∞), we have that

where $Ψ∞HF$ was introduced in Eq. (12) (in the case of degeneracy, we can select one of the minimizers, since here we only want to obtain an upper bound to the lowest eigenvalue of $HλHF$⁠). We now analyze, for large α, the expectation value on $ΨλT$ of each term appearing in $ĤλHF$ of Eq. (7), obtaining

where t, h, and k are all positive numbers. This is obvious for t, but it is also true for k because the expectation of $−K^$ is positive for any wavefunction Ψ, as $K^$ has a negatively definite kernel. The fact that the expectation value of $K^$ on $ΨλT$ vanishes as α⁻¹ for large α is due to the non-locality of $K^$⁠, which samples the gaussians in the bra and in the ket in different points of space, and due to the regularity properties of the HF orbitals (which have no delta-function singularities). The positivity of h in Eq. (27) can be proven by expanding $H^∞HF$ around $R̲min$ up to the second order, which gives a positive definite hessian matrix.

Putting together Eqs. (26)–(29) and replacing α with λ^q, we find that, for large λ, the expectation value of $ĤλHF$ on $ΨλT$ behaves asymptotically as

With t, h, and k being positive, we see that the best variational choice to make the next leading term after O(λ) increase with the lowest possible power of λ is q = 1/2, as conjectured in Eq. (11). Although $ΨλT$ of Eq. (24) is not antisymmetric, we can always properly antisymmetrize it, which only leads to corrections O(e^−α) in the computation of the expectation values, similarly to the DFT case.²⁵ 

Thus, we have explicitly constructed a variational wavefunction that yields the minimum possible value for the leading term O(λ) in the expectation of $ĤλHF$⁠. In fact, since E_el[ρ^HF] − U[ρ^HF] is the global minimum of the multiplicative operator $V^ee−Ĵ$⁠, there is no wavefunction that can yield a lower expectation for this operator. Moreover, since $−K^$ is positive definite, the best we can do is to make its expectation zero when λ → ∞, which our wavefunction is able to do.

This variational argument also shows that the next leading term in $WλHF$ should be order λ^−1/2, similarly to the DFT case of Eq. (5). A quantitative estimate of this next leading term could be in principle obtained by using the normal modes around the minimum of $V^ee−Ĵ$⁠: a unitary transformation from the $ri−rimin$ to the normal mode coordinates ξ₁, …, ξ_3N that diagonalize the hessian of $H^∞HF$ at $R̲min$ leads to a set of uncoupled harmonic oscillators whose spring constants scale with λ

with $ωα2$ being the eigenvalues of the hessian of $H^∞HF$ at $R̲min$⁠. The ground-state of $ĤλZP$ is obtained by occupying the lowest state of each oscillator, with the product state

This wavefunction should provide the minimum possible expectation, to order λ^1/2, of $T^+λ(V^ee−Ĵ)$⁠. However, since $−λK^$ is of the same order λ^1/2, we cannot exclude at this point that the minimization of the full $T^+λ(V^ee−Ĵ−K^)$ could lead to a different set of occupied oscillator states. This investigation will be the object of future studies. From our present treatment, we have so far

with

where $WK,∞′HF$ is due to the effect of $−λK^$ at orders λ^1/2 in $ĤλHF$ and is a functional of the occupied HF orbitals. Equation (34) should be exact while Eq. (35) is for now a conjecture. We also see that both $W∞HF$ and $W∞′HF$ have a part that is a functional of the HF density only and a part that is a functional of the occupied HF orbitals. In both cases, the part that is a density functional has an origin similar to the one of the DFT functionals of Eq. (5), being, respectively, a classical electrostatic energy and the potential energy of zero-point oscillations around a classical minimum. The parts that need the knowledge of the occupied HF orbitals do not appear in the DFT case. This structure should be exact, although the detailed form of $W∞′HF$ might include a different set of occupied oscillator states.

Although the λ → ∞ limit of $WλHF$ has a structure similar to the one of DFT, there are many differences that need to be kept in mind. Both $WλDFT[ρ]$ and $WλHF$ are decreasing functions of λ

but $WλDFT[ρ]$ for λ ≥ 0 is believed to be convex or at least piecewise convex (if there are crossings of states), while $WλHF$ is for sure not always convex. In fact, the MP2 correlation energy usually underestimates (in absolute value) the total correlation energy $EcHF$⁠, implying that $WλHF$ for 0 < λ ≪ 1 must run below its tangent; thus, $WλHF$ usually starts concave for small λ and then needs to change convexity to tend to the finite asymptotic value $W∞HF$ for large λ. Moreover, while the density constraint of the DFT adiabatic connection usually mitigates the crossing of states, the HF adiabatic connection might have jumps or kinks as λ is increased. A simple example is the N = 1 case, for which $WλHF=−U[ρHF]$ for 0 ≤ λ ≤ 1, while for λ > 1, the curve starts to decrease, tending, as λ → ∞ to a well defined value, with the electrostatic energy determined by the configuration in which the electron is sitting in the minimum of −$vH$(r, [ρ]).

In conclusion, we have shown that by looking at the λ → ∞ limit of the HF adiabatic connection, we recover functionals of the HF density, revealing a new intriguing formal link between HF and DFT. However, we should also stress that the use of models for $WλHF$ taken from DFT, although somehow justified by our analysis, should at this stage still be taken with some caution. The empirical observation so far^37,38 is that these models are not accurate for total energies, but work rather well for interaction energies, with a small variance, particularly for non-covalent complexes.^37,38 This point requires further investigation, which will be the object of a paper in preparation, where $WλHF$ will be computed and analyzed for various systems, and will be compared against various models. We will also evaluate and further analyze $W∞HF$ and $W∞′HF$⁠, making a detailed numerical comparison with the corresponding DFT functionals. We can already remark that the difference between $W∞HF$ and $W∞DFT$ can be big. For example, for the He atom, we have $W∞HF≈−4.347$ Ha, while $W∞DFT≈−1.50$ Ha. We will also study whether it is possible to extract a model for the self-energy in the strong-coupling limit, to be used in the context of Green’s function approaches,^49–52 and to analyze in the same spirit adiabatic connections appearing in other theories.^39,53

This work was supported by the European Research Council under H2020/ERC Consolidator Grant corr-DFT (Grant No. 648932). S.V. acknowledges financial support from NWO through Rubicon (Grant No. 019.181EN.026).