Functional Derivative of the Zero Point Energy Functional from the Strong Interaction Limit of Density Functional Theory

TL;DR

This paper derives an explicit functional derivative for the zero point energy term in a simplified strong interaction limit of density functional theory, providing insights into strongly correlated electron systems.

Contribution

It presents a new explicit expression for the ZPE functional derivative in a specific two-electron, one-dimensional case, advancing the understanding of strong correlation in DFT.

Findings

Derived explicit ZPE functional derivative for two-electron systems

Confirmed the expression numerically and verified sum-rule compliance

Showed ZPE potential can produce a bond mid-point peak in dissociation

Abstract

We derive an explicit expression for the functional derivative of the subleading term in the strong interaction limit expansion of the generalized Levy--Lieb functional for the special case of two electrons in one dimension. The expression is derived from the zero point energy (ZPE) functional, which is valid if the quantum state reduces to strongly correlated electrons in the strong coupling limit. The explicit expression is confirmed numerically and respects the relevant sum-rule. We also show that the ZPE potential is able to generate a bond mid-point peak for homo-nuclear dissociation and is properly of purely kinetic origin. Unfortunately, the ZPE diverges for Coulomb systems, whereas the exact peaks should be finite.

Equations

$$F_{\lambda }[\rho ]$$

$$\equiv ⟨{\Psi }_{\lambda }[\rho ]|\hat{T}+\lambda {\hat{V}}_{ee}|{\Psi }_{\lambda }[\rho ]⟩.$$

$${\hat{H}}_{\lambda }[\rho ]\equiv \hat{T}+\lambda {\hat{V}}_{ee}+{\hat{V}}^{\lambda }[\rho ],$$

$$v_{ee}(x)=\frac{1}{|x|}.$$

$$⟨{\Psi }_{\lambda }|\hat{\rho }(\mathbf{r})|{\Psi }_{\lambda }⟩$$

$$E_{\lambda }[\rho ]$$

$$=\underset{\tilde{\rho }}{\min }\left(F_{\lambda }[\tilde{\rho }]+\int \mathrm{d}\mathbf{r}\tilde{\rho }(\mathbf{r})v^{\lambda }[\rho ](\mathbf{r})\right)$$

$$\frac{\delta F_{\lambda}[\tilde{\rho}]}{\delta\tilde{\rho}(\mathbf{r})}\bigg{|}_{\mathrlap{\tilde{\rho}=\rho}}=-v^{\lambda}[\rho](\mathbf{r}),$$

$$\underset{\lambda \to \infty }{\lim }\frac{{\hat{H}}_{\lambda }}{\lambda }=\underset{\lambda \to \infty }{\lim }\frac{\lambda {\hat{V}}_{ee}+{\hat{V}}^{\lambda }}{\lambda }\equiv {\hat{V}}_{ee}+{\hat{V}}^{\text{SIL}},$$

$$F_{\lambda }[\rho ]$$

$$\frac{\delta V_{ee}^{\text{SIL}}[\rho]}{\delta\rho(\mathbf{r})}\bigg{|}_{\mathrlap{\rho=\rho_{0}}}=-v^{\text{SIL}}[\rho_{0}](\mathbf{r}).$$

$$\begin{gathered} E_{\text{pot}}({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N) \\ \equiv \sum _{i>j}{v}_{ee}(|{\mathbf{r}}_i-{\mathbf{r}}_j|)+\sum _{i=1}^N{v}^{\text{SIL}}({\mathbf{r}}_i). \end{gathered}$$

$$|{\Psi }_{\text{SCE}}({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N){|}^2\equiv |{\Psi }_{\lambda \to \infty }({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N){|}^2$$

$$=\frac{1}{N!}\sum _{\wp }\int \mathrm{d}\mathbf{s}\prod _{i=1}^N\frac{\rho (\mathbf{s})}{N}\delta ({\mathbf{r}}_i-{\mathbf{f}}_{\wp (i)}(\mathbf{s}))$$

$${\Omega }_0(\mathbf{s})$$

$${\mathbf{f}}_1(\mathbf{r})$$

$${\mathbf{f}}_2(\mathbf{r})$$

$${\mathbf{f}}_3(\mathbf{r})$$

$$\text{mathmakebox}[\text{widthof}=][c]⋮$$

$${\mathbf{f}}_N(\mathbf{r})$$

$$\displaystyle\underbrace{\mathbf{f}\bigl{(}\mathbf{f}(\ldots\mathbf{f}(\mathbf{r})\ldots)\bigr{)}}_{N\text{ times}}=\mathbf{r}.$$

$$\rho (\mathbf{r})\mathrm{d}\mathbf{r}$$

$$\displaystyle E_{\text{pot}}\bigl{(}\Omega_{0}(\mathbf{s})\bigr{)}$$

$$V_{ee}^{\text{SCE}}[\rho ]$$

$$=\frac{1}{N}\sum _{i=1}^{N-1}\sum _{j=i+1}^N{\int }_{{\mathbb{R}}^D}^{}\mathrm{d}\mathbf{r}\rho (\mathbf{r})v_{ee}(|{\mathbf{f}}_i(\mathbf{r})-{\mathbf{f}}_j(\mathbf{r})|)$$

$$=\frac{1}{2}\sum _{i=1}^{N-1}{\int }_{{\mathbb{R}}^D}^{}\mathrm{d}\mathbf{r}\rho (\mathbf{r})v_{ee}(|\mathbf{r}-{\mathbf{f}}_i(\mathbf{r})|).$$

$$\mathrm{Tr}\left(\sqrt{\mathbb{H}(\mathbf{s})}\right)\equiv \sum _{\mu =D+1}^{DN}{\omega }_{\mu }(\mathbf{s}).$$

$$v^{\lambda }(\mathbf{r})$$

$${\hat{H}}_{\lambda }$$

$${\hat{H}}^{\text{ZPE}}=$$

$$\displaystyle{}+\hat{V}^{\text{ZPE}}\bigl{(}\mathbf{s},\mathbf{f}_{2}(\mathbf{s}),\dotsc,\mathbf{f}_{N}(\mathbf{s})\bigr{)}.$$

$$\displaystyle\hat{V}^{\text{ZPE}}\bigl{(}\mathbf{s},\mathbf{f}_{2}(\mathbf{s}),\dotsc,\mathbf{f}_{N}(\mathbf{s})\bigr{)}$$

$$=-\sum _{\mu =D+1}^{DN}\frac{{\omega }_{\mu }(\mathbf{s})}{2}+\text{const}.$$

$$F_{\lambda }[\rho ]$$

$$F^{\text{ZPE}}[\rho ]$$

$$=\frac{1}{2}{\int }_{{\mathbb{R}}^D}^{}\mathrm{d}\mathbf{s}\frac{\rho (\mathbf{s})}{N}\mathrm{Tr}\left(\sqrt{\mathbb{H}(\mathbf{s})}\right),$$

$$f[\rho](s)=\begin{cases}N_{e}^{-1}\bigl{(}N_{e}(s)+1\bigr{)}&s<N_{e}^{-1}(1)\\ N_{e}^{-1}\bigl{(}N_{e}(s)-1\bigr{)}&s>N_{e}^{-1}(1),\end{cases}$$

Full text

I Introduction

Kohn-Sham (KS) density functional theory (DFT) is the workhorse of electronic structure calculations in physics and chemistry, thanks to its good compromise between accuracy and computational cost. Although exact in principle, KS DFT must rely in practice on approximations for the exchange-correlation (XC) functional, which, despite their many successes, have still problems in describing strongly correlated systems, whose physics is very different than the one of the non-interacting KS reference system (Mardirossian and Head-Gordon, 2017; Cohen et al., 2012; Burke, 2012; Becke, 2014).

The strong interaction limit (Seidl, 1999; Seidl et al., 2007; Gori-Giorgi and Seidl, 2010) (SIL) of the universal part of the ground-state energy density functional (Levy, 1979, 1982; Lieb, 1983) is a semi-classical limit in which the electron-electron energy dominates over the kinetic energy, and it is the first term of an expansion of the generalized Levy–Lieb functional in the form of an asymptotic series for the electronic coupling constant $\lambda\rightarrow\infty$. This same expansion also determines the asymptotic behaviour of the exact XC functional of KS DFT (Malet and Gori-Giorgi, 2012; Malet et al., 2013; Mendl et al., 2014) at strong coupling. In specific cases, the SIL solution reduces to a particular simple form of strictly correlated electrons (SCE), in which the position of one electron dictates the position of all other electrons (Seidl et al., 2007; Seidl, 1999; Gori-Giorgi and Seidl, 2010). In more general cases, the search over SCE-type solutions yields the exact SIL as an infimum (Colombo and Di Marino, 2013). The SCE solution unveils how the exact XC functional mathematically transforms the information on the electronic density into an expectation value of the electron-electron repulsion, even if only in the case of its $\lambda\to\infty$ asymptotic expansion. Its investigation lead to the construction of new non-local density functionals, based on particular integrals of the density (Wagner and Gori-Giorgi, 2014; Zhou et al., 2015; Vuckovic and Gori-Giorgi, 2017) rather than on the traditional ingredients of standard approximations (local density and gradients, occupied and unoccupied KS orbitals).

The first subleading term for SCE-type solutions introduces kinetic energy in the form of zero-point oscillations (ZPE). It has been first evaluated in 2009 (Gori-Giorgi et al., 2009) and received numerical confirmations only recently (Cort et al., 2017; Grossi et al., 2017). Little is known yet on the third leading term, for which scaling arguments suggest it to be of purely kinetic nature (Gori-Giorgi et al., 2009; Ying et al., 2016). This third term should incorporate exact pieces of information on the ionization energy of the system under exam (Giarrusso et al., 2018a).

Besides the XC functional itself, another quantity that plays an important role in KS DFT is its functional derivative with respect to the density, which determines the XC potential entering in the KS equations. The exact (or very accurate) XC potential has been studied for small systems in several works, using various reverse-engineering procedures (Ospadov et al., 2017; Cuevas-Saavedra et al., 2015; Cuevas-Saavedra and Staroverov, 2016): these works have shown that for strongly-correlated systems the XC potential must display very peculiar features, such as “peaks” and “steps” (Baerends and Gritsenko, 1997; Hodgson et al., 2016; van Leeuwen et al., 1995). While the functional derivative of the SCE leading term has been evaluated and used as an approximation for the XC potential in the self-consistent KS equations in various works (Malet and Gori-Giorgi, 2012; Malet et al., 2013; Mendl et al., 2014; Malet et al., 2014), the potential associated to the next leading term has never been investigated in an exact manner (only very recently, a semi-local approximation for the ZPE has been used to look at KS potentials coming from functionals that interpolate between the weak- and strong-coupling limits of the XC functional (Fabiano et al., 2019)). It is the purpose of this paper to fill this gap, by starting an investigation of the exact ZPE functional derivative. The SIL functionals have a density dependence that is rather complicated and unusual, making it actually difficult to evaluate functional derivatives. The reason why the functional derivative of the leading SIL term (the SCE term) could be easily computed is that it can be obtained from an exact shortcut (Malet and Gori-Giorgi, 2012; Malet et al., 2013), which seems to be missing at the next leading order. For this reason, our investigation starts from a simple, yet non trivial, case: two electrons confined in one dimension (1D). Similar 1D models have been widely used to investigate features in exact KS DFT, proving to provide a good qualitative description of the relevant features of their 3D counterparts (Mori-Sánchez and Cohen, 2018; Tempel et al., 2009; Benítez and Proetto, 2016; Magyar and Burke, 2004; Hodgson et al., 2016; Helbig et al., 2009).

Besides its interest as an XC potential at strong coupling, the ZPE functional derivative that we compute here is also a crucial ingredient to analyze the third term in the large-$\lambda$ expansion of the exact Levy–Lieb functional. This next term, in fact, requires solving a hierarchy of Schrödinger equations for which knowledge of the asymptotic expansion at strong coupling of $v^{\lambda}$ (the 1-body potential that keeps the density fixed at each $\lambda$) is needed; the potential $v^{\lambda}$ at orders $\lambda^{1/2}$ should be given by minus the ZPE functional derivative (Gori-Giorgi et al., 2009).

The paper is organized as follows: we first briefly cover the key concepts of SCE and ZPE formalism in section II. The core of the paper, section III, hosts an analytical expression of the functional derivative of the ZPE functional (35). Its features are discussed and numerical calculation is provided to verify the consistency of our results. Last, in section V we draw our conclusions and outline future steps.

II Theoretical Overview

Let us consider the universal functional $F_{\lambda}[\rho]$, defined in the Levy constrained search formulation for any $\lambda\in\mathbb{R}$ as

[TABLE]

Under the assumption of a ground state $v$-representable density, the minimizing wavefunction $\Psi_{\lambda}[\rho]$ in (II) is also a ground state (Levy and Perdew, 1985; Levy, 1982) of the $\lambda$-dependent Hamiltonian

[TABLE]

where $\hat{T}$ is the familiar kinetic energy operator, and $\hat{V}_{ee}=\frac{1}{2}\sum_{i\neq j}^{N}v_{ee}(\lvert\mathbf{r}_{i}-\mathbf{r}_{j}\rvert)$ is the electron-electron interaction operator. For realistic electrons in 3D space,

[TABLE]

For the 1D case, see (38) below. Generally, we choose piecewise convex functions $v_{ee}(x)$. The local one body operator $\hat{V}^{\lambda}[\rho]=\sum_{i=1}^{N}v^{\lambda}[\rho](\mathbf{r}_{i})$ is the Lagrange multiplier that enforces the constraint

[TABLE]

The $\lambda$-dependent energy

[TABLE]

connects (II) and (2). The minimization over $\tilde{\rho}$ implies that (Dreizler and Gross, 1990)

[TABLE]

modulo a constant. In what follows, we recall the basic ideas needed to apply these concepts to the regime $\lambda\gg 1$.

II.1 Strictly Correlated Electrons

From physical arguments, one suspects that $\lim_{\lambda\to\infty}F_{\lambda}[\rho]/\lambda=\langle{\Psi_{\lambda\to\infty}[\rho]}|{\hat{V}_{ee}}|{\Psi_{\lambda\to\infty}[\rho]}\rangle$ (Seidl, 1999; Seidl et al., 2007; Gori-Giorgi et al., 2009); this result was proved rigorously only recently (Cotar et al., 2018; Lewin, 2018). As a consequence, to satisfy the density constraint (4) we must have to leading order that in the limit $\lambda\to\infty$ the force exerted by the external potential is of the same order in $\lambda$ as the electron-electron repulsion. In the SIL regime we hence define the local one-body operator $\hat{V}^{\text{SIL}}=\sum_{i=1}^{N}v^{\text{SIL}}[\rho](\mathbf{r}_{i})$

[TABLE]

and the functional $V_{ee}^{\text{SIL}}[\rho]$ as

[TABLE]

These two quantities are connected by (6), i.e.

[TABLE]

Equation (7) defines a function in configuration space:

[TABLE]

The minimization problem in (8) can be regarded as an optimal transport problem with repulsive cost (Buttazzo et al., 2012).

A candidate solution to this problem was first the so-called strictly correlated electrons (SCE) ansatz and satisfies $V_{ee}^{\text{SCE}}[\rho]\geq V_{ee}^{\text{SIL}}[\rho]$. The SCE ansatz was suggested on physical grounds by Seidl and co-workers (Seidl, 1999; Seidl et al., 2007), and has been proved rigorously to be exact for $D=1$ or $N=2$ in $D>1$, provided the interaction $v_{ee}(x)$ is convex and bounded from below (Colombo et al., 2015). In the following, we assume that the proposed SCE solution is the exact SIL solution, so we replace SIL by SCE.

The underlying idea of SCE is that the positions of the electrons become strictly correlated, i.e. the position of one electron dictates the whereabouts of all other electrons. This means that the minimizer of (8) is a distribution that is zero in the whole configuration space except for a subset $\Omega_{0}$ of dimension $D$

[TABLE]

where $\wp$ denotes any permutation of $N$ elements, and

[TABLE]

The optimal maps or co-motion functions $\mathbf{f}_{i}[\rho]$ are non-local functionals of the density and their physical meaning is to provide the position of $N-1$ electrons as a function of the position of the first electron. Indistinguishability can be guaranteed by requiring the following group properties (Seidl et al., 2007; Malet et al., 2014)

[TABLE]

Furthermore, the density constraint implies the differential equation

[TABLE]

The minimum of (7) must be degenerate in $\Omega_{0}(\mathbf{s})$: a hypothetical minimum in a specific point $\mathbf{s^{*}}$ would collapse the system into a frozen configuration of positions $\{\mathbf{f}_{1}(\mathbf{s^{*}}),\mathbf{f}_{2}(\mathbf{s^{*}}),\dotsc,\mathbf{f}_{N}(\mathbf{s^{*}})\}$, in violation of the smooth density constraint (4). Hence we must have

[TABLE]

Finally, with (II.1) $V_{ee}^{\text{SCE}}[\rho]$ reads

[TABLE]

II.2 Zero Point Energy

As anticipated, at finite large $\lambda$ the characterization of the ground state of Hamiltonian (2) departs from the semi-classical picture, as the kinetic energy starts to play a relevant role in the description of the underlying physics in the form of zero-point oscillations performed near $\Omega_{0}$.

Consider $\mathbb{H}(\mathbf{s})$, the Hessian of $E_{\mathrm{pot}}(\mathbf{r}_{1},\dotsc,\mathbf{r}_{N})$ evaluated in $\Omega_{0}(\mathbf{s})$. This matrix has $D$ zero eigenvalues and $DN-D$ positive $\mathbf{s}$-dependent eigenvalues, $\omega_{\mu}(\mathbf{s})^{2}$,

[TABLE]

The corresponding eigenvectors induce a set of curvilinear coordinates $u_{\mu}$ in terms of which $\hat{H}^{\lambda}$ can be expanded (Gori-Giorgi et al., 2009; Grossi et al., 2017).

Retaining the leading order in the expansion of the Laplace–Beltrami operator for the kinetic energy, we have argued (Gori-Giorgi et al., 2009; Grossi et al., 2017)

[TABLE]

This allows one to write

[TABLE]

where the operator $\hat{H}^{\text{ZPE}}$ reads

[TABLE]

For each fixed $\mathbf{s}$, $\hat{H}^{\text{ZPE}}$ has the structure of a set of harmonic oscillators in the coordinates $u_{\mu}$. The term denoted $\hat{V}^{\text{ZPE}}$, depending only on $\mathbf{s}$, does not affect the harmonic nature of its solution and, by correcting the external potential computed in (8), keeps the degeneracy of the energy with respect to $\mathbf{s}$ through order $\sqrt{\lambda}$, provided the following constraint (Gori-Giorgi et al., 2009) is imposed

[TABLE]

This allows us to give an explicit expression for the subleading term of the generalized universal functional in the strong interaction limit

[TABLE]

with

[TABLE]

and $|\Psi_{\text{ZPE}}[\rho]\rangle$ denotes the ground state of (II.2). Notice that in previous works (Seidl, 1999; Gori-Giorgi et al., 2009; Grossi et al., 2017) $F^{\text{ZPE}}[\rho]$ was denoted as $2W_{\infty}^{\prime}[\rho]$, in analogy with the linear coefficient in the expansion at small $\lambda$ of $F_{\lambda}[\rho]$ (see also (48) below).

III Functional derivative of $F^{\text{ZPE}}[\rho]$ for $N=2$, $D=1$

III.1 SCE + ZPE for $N=2$ electrons in 1D

This brief paragraph is devoted to provide the quantities described so far in the case of 2 electrons in $D=1$ , as well as set of useful relations that help to considerably simplify the calculation outlined in the next sections.

Defining $f_{1}(s)\equiv s$, $f_{2}(s)\equiv f(s)$, we have (Malet et al., 2014)

[TABLE]

where

[TABLE]

The comotion function is such that the integral of the density between $x$ and $f(x)$ always integrates to 1 independently of $x$. Therefore, when $x<0$, for a symmetric density $f(x)$ must necessarily be positive, and vice versa. As the reference electron approaches [math] from the left, the second electron is pushed towards $+\infty$. When the reference electron crosses the origin, the second electron must ”jump” to $-\infty$.

The only non-zero frequency (eigenvalue of the $2\times 2$ matrix $\mathbb{H}(s)$) is given by (Malet et al., 2014)

[TABLE]

Notice that $v_{ee}(x)$ is convex, $v^{\prime\prime}_{ee}(x)>0$, and that $f^{\prime}(x)>0$, see (28a) below. Equation (II.2) reads explicitly

[TABLE]

Moreover, equations (13) and (14) read

[TABLE]

implying \omega\bigl{(}f(s)\bigr{)}=\omega(s).

III.2 Explicit expression

Inserting (18) and (22) in (6) and comparing the terms proportional to $\sqrt{\lambda}$, we have

[TABLE]

The derivation of an explicit form for $\delta F^{\text{ZPE}}/\delta\rho(x)$ starts from noticing that

[TABLE]

The frequency function $\omega(x)$ is an implicit functional of the density, via the co-motion function and its derivative. Even for 2 electrons in $D=1$, computing the functional derivatives of $f(x)$ can be delicate, as it changes sign when $N_{e}(s)=1$: perturbing the density in this point implies taking into account a step function, for which the chain rule does not apply (see Appendix in (Lani et al., 2016) for further details). Step functions are also expected whenever there is a step in $\rho(x)$ or a difference in the values of the density at the boundaries in a compact support. This is not our case however, since we assume $\rho(x)$ to be a continuous integrable function defined on the whole real axis. As a consequence, $\lim_{x\rightarrow N_{e}^{-1}(1)\uparrow}\omega(x)=\lim_{x\rightarrow N_{e}^{-1}(1)\downarrow}\omega(x)$, there is no step to be taken into account. Hence, we can simply apply the chain rule and write

[TABLE]

which reads

[TABLE]

For the chain rule, only the regular part of the functional derivative of $f(x)$, which can be found in (Lani et al., 2016), is relevant, and reads in 1D

[TABLE]

$\Theta(x)$ being the Heaviside step function.

For the functional derivative of $f^{\prime}[\rho](x)$, we make use of (28b),

[TABLE]

In the appendix, we show that, using (33) and (III.2) in (31) and inserting the result in (III.2), $\delta F^{\text{ZPE}}/\delta\rho(x)$ can be expressed as (see Appendix for details)

[TABLE]

where $\Lambda(y)$ is an odd, well behaved function (see also Fig. 1) and reads explicitly

[TABLE]

Equation (II.2) implies a sum rule on $\delta F^{\text{ZPE}}/\delta\rho(x)$. Inserting (29) in (II.2), and remembering that \omega(s)=\omega\bigl{(}f(s)\bigr{)}, we see that we must have

[TABLE]

Since I\bigl{(}f(x)\bigr{)}=-I(x), this is consistent with our result (35).

III.3 Numerical results for selected densities

In this section, we are going to verify (35) numerically, using the effective convex Coulomb interaction renormalized at the origin

[TABLE]

(See (Grossi et al., 2017) for a brief discussion on the importance of convexity of the interaction in SCE-DFT.) We pick 3 test densities, peaked at $x=0$

[TABLE]

All the respective co-motion functions can be evaluated analytically since the inverse function of (25) can be written explicitly. In Fig. 2, we provide the profile of $\delta F^{\text{ZPE}}/\delta\rho(x)$ for the test densities (39). The plots show that the shape of the curve can vary drastically depending on the density chosen. In particular, in all the densities we chose (excluding $\rho_{3}$) the functional derivative shows divergences both in the origin and in the large $x$ limit. The nature of this divergences shall be investigated deeper in Sec. III.4.

Since the derivation of (35) was cumbersome, we decided to verify it numerically to exclude any possible error. We thus simply use the definition of functional derivative

[TABLE]

If our expression for the functional derivative is correct, we should have that the slope of the l.h.s. of (40) at $\epsilon=0$ coincides with the straight line on the r.h.s. of (40). For the numerical verification, we consider the following perturbations

[TABLE]

The shape of these functions has been chosen arbitrarily, though they are symmetric, integrate to 0 (thus not changing the number of particles) and, thanks to their fast decay at large $x$, are such that $\rho_{i}(x)+\epsilon\phi(x)>0\quad\forall x\in\mathbb{R}$, for at least $\epsilon\in[-0.5,0.5]$ for the chosen densities. In Fig. 3 we show the l.h.s. of (40) as a function of $\epsilon$ and the corresponding r.h.s., linear in $\epsilon$. In all cases the tangent of the l.h.s. of (40) shows an excellent agreement with (35).

III.4 Divergencies of $\delta F^{\text{ZPE}}/\delta\rho(x)$ in 1D

In what follows, we study the behaviour of the functional derivative at large $x$. The same behaviour can be deduced for small $x$, due to the fact that \omega(x)=\omega\bigl{(}f(x)\bigr{)} and that $\lim_{x\to 0^{\pm}}f(x)=\mp\infty$ (see text after (25)). Keeping in mind that $\lim_{x\to\infty}I(x)=\mathrm{const}$, it is clear from (35) that for $x\gg 1$

[TABLE]

The behaviour of $\delta F^{\text{ZPE}}/\delta\rho(x)$ at large $x$ is dominated by $\omega(x)$, which in turn is determined by the interplay between the electron-electron interaction and the density decay at large $x$, cf. (26). With interaction (38) $v^{\prime\prime}_{ee}(x)\sim x^{-3}$ at large $x$, hence the frequency will diverge whenever $\rho(x)=o(x^{-3})$ for $x\gg 1$. This is the case for densities $\rho_{1,2}$ in (39) which both decay exponentially (or faster). Such a divergence of $\omega(x)$ makes the interpretation of the expansion of $v^{\lambda}$ less straightforward: for what just stated in (42), at large distances, its asymptotic expansion reads

[TABLE]

At first glance, it seems that the expansion at large $\lambda$ for $v^{\lambda}$ is not consistent with the requirement $v^{\lambda}\in L^{3/2}+L^{\infty}$: if $\omega(x)$ diverges to $+\infty$ then, for every fixed $\lambda$, there is a point $x$ after which the second term in (43) becomes dominant and the minimum of $v^{\lambda}(x)$ is at $x=\pm\infty$ (since $v^{\text{SCE}}(x)\sim-(N-1)/\lvert x\rvert$ for large $x$ for the chosen interaction). To make sense of (43), one has to be careful in taking the correct order of limits: what we mean here is that for each fixed $x$ the expansion of $v^{\lambda}$ as a function of $\lambda$ follows (43).

On the other hand, 1D models often assume an effective electron-electron interaction depending on the physics they aim to describe, often leading to short range interactions. From the preceding discussion, it is clear that a short-range interaction should lead to a better behaviour of $\omega(x)$ and hence, the behaviour of $v^{\text{ZPE}}$. We have tried two different short range interactions, namely a modified Yukawa potential

[TABLE]

As an example, in Fig. 4 we plot how the profile of $\omega[\rho_{2}](x)$ varies as we pick different interactions. If we pick a sufficiently high $\alpha$, $\omega(x)$ is damped (and consequently the functional derivative $\delta F^{\text{ZPE}}/\delta\rho(x)$). We choose $\alpha=2$, since $\alpha\geq 1$ leads to a convergent frequency ($\omega[\rho_{2}](x)\sim\frac{\sqrt{\pi}\alpha}{2}x^{-\frac{1}{2}}\mathrm{e}^{(1-\alpha)x}$).

Notice that neither (44a) nor (44b) would provide a finite $\omega(x)$ when using density $\rho_{1}$, as the Gaussian decay would prevail on both interactions with any choice of parameters. A faster decaying interaction would be needed, e.g. $\sim\mathrm{e}^{-x^{2}}$.

IV Exchange-correlation potential for a 1D dimer

It is known (Buijse et al., 1989; Gritsenko and Baerends, 1996; Baerends and Gritsenko, 1997; van Leeuwen and Baerends, 1994; Helbig et al., 2009; Ying et al., 2016) that the exact exchange-correlation (XC) potential of a homo-nuclear dimer builds a peak in the mid-bond region that, in the dissociating limit, must be proportional in height to the ionization potential of each fragment. Although some GGA functionals build peak-like features in the bond mid-point (van Leeuwen and Baerends, 1994), they miss its peculiar scaling properties (Ying et al., 2016) which in general are not recovered by local, semilocal or hybrid functionals (Ying et al., 2016). Using only $-v^{\text{SCE}}$ as an approximation to the true XC correlation potential does not allow to recover exactly this feature, which is of purely kinetic nature (Ying et al., 2016; Malet et al., 2014; Giarrusso et al., 2018a). It is the purpose of this section to investigate whether the expression obtained so far can help in reproducing, at least qualitatively, this characteristic.

Consider the density $\rho_{D}$

[TABLE]

Having two equal maxima located at $\pm\frac{R}{2}$, $\rho_{D}$ can be considered as a 1D model for a homo-nuclear dimer whose density profile is parametrically dependent on the internuclear distance $R$. This model has been used several times (Tempel et al., 2009; Baker et al., 2015; Helbig et al., 2009; Hodgson et al., 2016; Benítez and Proetto, 2016) since it has been proved to mimic many exact features of the exact KS potential for real molecules; in particular, it gives us the opportunity to model the bond stretching and analyse the kinetic contributions to the XC potential.

To write an expression for the XC potential, we start from the adiabatic connection formalism (Harris, 1984). The XC energy can be written exactly in terms of an integral over the coupling $\lambda$

[TABLE]

where

[TABLE]

$U_{\text{H}}[\rho]$ being the Hartree functional. Using the large $\lambda$ expansion of the adiabatic connection integrand (Gori-Giorgi and Seidl, 2010)

[TABLE]

we obtain

[TABLE]

and

[TABLE]

In Fig. 5 we show the potential in (50) for $R=15$. Via (42), $\delta F^{\text{ZPE}}/\delta\rho(x)$ indeed introduces a correction in the mid-bond region. In fact, since we have (Lani et al., 2016)

[TABLE]

from our treatment in section III.4, the divergence in the mid-bond region can be readily evaluated inserting (51) in (26)

[TABLE]

For any fixed $x\neq 0$, we find $\delta F^{\text{ZPE}}/\delta\rho_{D}(R;x)\to 0$ as $R\to\infty$ while similarly, due to the fact that \omega(x)=\omega\bigl{(}f(x)\bigr{)}, a divergence of the XC potential appears also at large $x$. Thus, the kinetic correlation energy introduced by the ZPE creates a divergence instead of a finite peak in the bond mid-point, and this divergence occurs on a region that shrinks when $R\to\infty$. This divergence is due to the extreme correlation between the two electrons: when, say, electron 1 oscillates around the origin, electron 2 jumps from plus to minus infinity. In the exact wavefunction, when one electron crosses the bond mid-point, the conditional position of the other electron also “jumps” from one atom to the other (which is the origin of the peak (Buijse et al., 1989; Helbig et al., 2009; Ying et al., 2016)), but it is distributed according to the one-electron density on each atom.

The ZPE correction to the SCE approximation of $v_{\text{XC}}$ in (50) includes a positive contribution from the region $\lambda\sim 0$ that, although integrable, is too large to provide a reasonable estimate of the XC energy $E_{\text{XC}}[\rho]$. This is due to the fact that we are using only pieces of information from the high coupling limit to approximate $W_{\lambda}[\rho]$. A way to improve this approximation is to include also exact ingredients from the $\lambda\rightarrow 0$ limit (Seidl et al., 1999, 2000; Malet et al., 2014; Giarrusso et al., 2018b), by writing an expression that reproduces the correct behaviour of $W_{\lambda}$ at small and strong couplings; among these, the interaction strength interpolation (ISI) (Seidl et al., 1999) has been object of study in recent years (Giarrusso et al., 2018b; Vuckovic et al., 2018; Fabiano et al., 2019). In this final paragraph, we investigate the effect of one of a simplified ISI as proposed in (Malet et al., 2014), which is size consistent for the dissociation of a system into two equal fragments. Hence we approximate $W_{\lambda}[\rho]$ to

[TABLE]

with

[TABLE]

The XC energy reads then

[TABLE]

While the XC potential is changed considerably at small $R$ (see Fig. 6), for large internuclear distances the effect of the ISI becomes negligible: already at $R=5$ (Fig. 7) we see that the effect is small and at $R=15$ the two curves becomes indistinguishable.

V Conclusions

In this work we worked out an explicit expression for the functional derivative of the subleading term of the generalized universal functional $F^{\text{ZPE}}[\rho]$ in the strong coupling limit of DFT for 2 electrons in 1D. Our expression respects the sum rules deduced first in (Gori-Giorgi et al., 2009) on physical grounds, and has been verified numerically.

We found that the asymptotic behaviour of $\delta F^{\text{ZPE}}/\delta\rho(x)$ for $x\to\infty$ is dictated by the asymptotic behaviour of the ZPE frequency $\omega(x)$. The asymptotic behaviour of $\omega(x)$ is dominated by the ratio $v^{\prime\prime}_{ee}(x)/\rho(x)$ for large $x$, so typically depends on the relative decay of the density compared to the interaction. For relatively fast decaying densities, $\omega(x)$ and hence $v^{\text{ZPE}}$ diverges for $x\to\infty$ and $x\to N_{e}(1)$. We expect similar features to be present in more general cases (higher dimensions and more particles). Though we do not have an explicit expression of $F^{\text{ZPE}}$ to directly evaluate its functional derivative, the sum rule (II.2) is generally valid and indicates that $v^{\text{ZPE}}$ should have at least the same divergencies as the ZPE frequencies $\omega_{\mu}(x)$. So in the general 1D case, we expect divergencies of the ZPE potential at the points where the density integrates to an integer particle number.

By studying the dissociation of a symmetric dimer, we have demonstrated that the ZPE correctly generates a peak in the mid-point region, properly purely built by the kinetic energy. Unfortunately, the diverging features of $\omega(x)$ also make the peak diverging for Coulomb systems, instead of reaching a finite value as in the exact case.

In the future, we aim to investigate the next leading term of the generalized universal functional. This should include exact pieces of information on the ionization energy, hence “curing” the divergencies appearing at the ZPE order (Giarrusso et al., 2018a). Another promising research line is the calculation of the kernel of $F^{\text{ZPE}}[\rho]$, i.e. its second functional derivative, which can be used as an adiabatic but spatially non-local XC-kernel in the response formulation of TD-DFT.

Acknowledgements.

Financial support was provided by the European Research Council under H2020/ERC Consolidator Grant “corr-DFT” [grant number 648932]. K.J.H.G. also acknowledges funding by Stichting voor Fundamenteel Onderzoek der Materie FOM Projectruimte [project 15PR3232].

Appendix A Calculation details for $\delta F^{\text{ZPE}}/\delta\rho(x)$

We have from (27)

[TABLE]

Using the chain rule in (31), the integral in the last equation can be written as

[TABLE]

With the substitution $u=f(y)$, the second delta function and step function can be combined with the first ones to yield

[TABLE]

The integrand of last integral is not well behaved due to the presence of $\omega(y)$, and is prone to numerical instabilities when evaluated. In our investigation we found that both integrals have opposite divergences, which can be eliminated by combining them. In order to do so, we proceed along two lines: first, we integrate the Dirac deltas in the first term and then rewrite the result as an integral, effectively performing an integration by parts of the Dirac delta. Secondly, remembering that the functional derivative is only defined modulo a constant, we can shift the region of integration, as this only gives a constant contribution and write

[TABLE]

where we defined $b_{+}>0$ as the point where $b_{+}=-f(b_{+})$. As outlined, we can now use the fundamental theorem of calculus to rewrite the first term as

[TABLE]

We make use of

[TABLE]

to write

[TABLE]

Combining these results we can write the integral in (56) as

[TABLE]

It is not transparent from this expression that (63) is odd under the exchange $x\to f(x)$. Moreover, the term $\sim v^{\prime\prime\prime}_{ee}f/\omega$ might not be bounded. To make it more clear, we apply again the transformation $u=f(y)$ to the first two terms in the integrand above and rewrite them as

[TABLE]

Now the integrands can be summed to yield

[TABLE]

and adding the integration between $-b_{+}$ and $b_{+}$, which amounts to adding only an immaterial constant to the functional derivative, yields (35).