Strictly correlated electrons approach to excitation energies of dissociating molecules

TL;DR

This paper investigates the application of the strictly correlated electrons (SCE) approach within density functional theory to accurately compute excitation energies in dissociating molecules, highlighting its strengths and limitations.

Contribution

It provides an analytic and numerical analysis of the SCE-based approximation for excitation energies, revealing its potential and challenges in modeling dissociation.

Findings

ASCE kernel exhibits spatial divergence at large bond distances.

ASCE yields non-zero excitation energies in dissociation regime, unlike ALDA.

Exact excitation energies are recovered only at very strong interactions.

Abstract

In this work we consider a numerically solvable model of a two-electron diatomic molecule to study a recently proposed approximation based on the density-functional theory of so-called strictly correlated electrons (SCE). We map out the full two-particle wave function for a wide range of bond distances and interaction strengths and obtain analytic results for the two-particle states and eigenenergies in various limits of strong and weak interactions, and in the limit of large bond distance. We then study the so-called Hartree-exchange-correlation (Hxc) kernel of time-dependent density functional theory which is a key ingredient in calculating excitation energies. We study an approximation based on adiabatic SCE (ASCE) theory which was shown to display a particular feature of the exact Hxc-kernel, namely a spatial divergence as function of the bond distance. This makes the ASCE kernel a candidate for correcting a notorious failure of the commonly used adiabatic local density approximation (ALDA) in the calculation of excitation energies of dissociating molecules. Unlike the ALDA, we obtain non-zero excitation energies from the ASCE kernel in the dissociation regime but they do not correspond to those of the true spectrum unless the interaction strength is taken to be very large such that the SCE theory has the right regime of validity, in which case the excitation energies become exact and represent the so-called zero point oscillations of the strictly correlated electrons. The commonly studied physical dissociation regime, namely large molecular separation at intermediate interaction strength, therefore remains a challenge for density functional approximations based on SCE theory.

Figures

Equations

$${\hat{H}}_{\lambda }=\hat{T}+{\hat{V}}_{\lambda }+\lambda \hat{W}$$

$$\left(-\frac{1}{2}{\nabla }^2+v_s[n](\mathbf{r})\right){\phi }_i(\mathbf{r},\sigma )={ϵ}_i{\phi }_i(\mathbf{r},\sigma )$$

$$v_{\text{Hxc}}^{\lambda }[n](\mathbf{r})=v_s[n](\mathbf{r})-v_{\lambda }[n](\mathbf{r}).$$

$$E_{\text{Hxc}}^{\lambda }[n]={\int }_0^{\lambda }d{\lambda }^{\prime }{W}_{{\lambda }^{\prime }}[n]$$

$$W_{\lambda }[n]=⟨{\Psi }_{\lambda }[n]|\hat{W}|{\Psi }_{\lambda }[n]⟩.$$

$$W_{\lambda }[n]=V_{\text{SCE}}[n]+\frac{V_{\text{ZPE}}[n]}{\sqrt{\lambda }}+O({\lambda }^{-3/2})$$

$$E_{\text{Hxc}}^{\lambda }[n]=\lambda V_{\text{SCE}}[n]+2\sqrt{\lambda }{V}_{\text{ZPE}}[n]+E_2[n]+O({\lambda }^{-1/2}).$$

$$v_{\text{Hxc}}^{\lambda }(\mathbf{r})=\lambda v_{\text{SCE}}(\mathbf{r})+\sqrt{\lambda }{v}_{\text{ZPE}}(\mathbf{r})+v_2(\mathbf{r})+O\left({\lambda }^{-1/2}\right)$$

$$f_{\text{Hxc}}^{\lambda }(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{\delta v_{\text{Hxc}}^{\lambda }(\mathbf{r})}{\delta n({\mathbf{r}}^{\prime })}$$

$$f_{\text{Hxc}}^{\lambda }(\mathbf{r},{\mathbf{r}}^{\prime })=\lambda \frac{\delta v_{\text{SCE}}(\mathbf{r})}{\delta n({\mathbf{r}}^{\prime })}+O(\sqrt{\lambda })$$

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

$$+\lambda {\cos }^2\left[\frac{\pi }{L}(x_1-x_2)\right]$$

$${\hat{H}}_s=-\frac{1}{2}\left({\partial }_{x_1}^2+{\partial }_{x_2}^2\right)+v_s(x_1)+v_s(x_2).$$

$$v_s(x)=V_0\left[1+\cos \left(\frac{4\pi x}{L}\right)\right]$$

$$v_s^{\prime \prime }(x_0)={\left(\frac{4\pi }{L}\right)}^2{V}_0=\alpha$$

$$\left[-\frac{1}{2}{\partial }_x^2+v_s(x)\right]{\phi }_l^{\pm }(x)={\epsilon }_l^{\pm }{\phi }_l^{\pm }(x)$$

$$\left[-{\partial }_z^2+2\nu \cos (2z)\right]M_l^{\pm }(z)=a_l^{\pm }(\nu )M_l^{\pm }(z)$$

$$a_l^{\pm }(\nu )$$

$$\nu (L)$$

$${\int }_{-\pi }^{\pi }dz(M_l^{\pm }(z){)}^2=\pi .$$

$${\phi }_l^{+}(x)$$

$${\phi }_l^{-}(x)$$

$${\epsilon }_1^{-}-{\epsilon }_0^{+}$$

$$=\frac{32}{L^2}(2\pi {)}^{3/2}{\nu }^{3/4}{e}^{-4\sqrt{\nu }}\underset{L\to \infty }{\to }0$$

$$n(0)=\frac{16}{L}{\left(2\pi \right)}^{\frac{1}{2}}{\nu }^{1/4}{e}^{-4\sqrt{\nu }}(L\to \infty )$$

$$n(x)=2|{\phi }_0^{+}(x){|}^2=\frac{4}{L}{C}_0^2(\frac{2\pi x}{L};\nu )$$

$$\Psi (x_1{\sigma }_1,x_2{\sigma }_2)=\psi (x_1,x_2)\frac{1}{\sqrt{2}}\left({\delta }_{{\sigma }_1\uparrow }{\delta }_{{\sigma }_2\downarrow }-{\delta }_{{\sigma }_1\downarrow }{\delta }_{{\sigma }_2\uparrow }\right)$$

$${\Psi }_{\lambda }^{\text{GS}}(x_1,x_2)$$

$$=\frac{1}{\sqrt{2}}\left[{\phi }_0^{+}(x_1){\phi }_0^{+}(x_2)-{\phi }_1^{-}(x_1){\phi }_1^{-}(x_2)\right]$$

$$E_{\lambda }^{\text{GS}}=2{\epsilon }_0^{+}=\sqrt{\alpha }-\frac{{\pi }^2}{L^2}+O(L^{-4})(L\to \infty )$$

$${\Psi }_{\lambda }^{(1)}(x_1,x_2)=$$

$$-{\phi }_0^{+}(x_1){\phi }_2^{-}(x_2)-{\phi }_2^{-}(x_1){\phi }_0^{+}(x_2)].$$

$$E_{\lambda }^{(1)}={\epsilon }_0^{+}+{\epsilon }_1^{+}=2\sqrt{\alpha }-\frac{3{\pi }^2}{L^2}+O(L^{-4})(L\to \infty )$$

Full text

I Introduction

Density-functional theory (DFT) is a commonly used electronic structure method. Its ground state version is mainly used to calculate energies and structures of electronic systems Barth:review , while its time-dependent (TD) counterpart TDDFT also allows for the calculation of dynamic properties and excitation energies Ullrich:book . Virtually all density-functional calculations are based on the Kohn-Sham (KS) system, a non-interacting system that produces the same electronic density as the true system of interest. The KS system provides a considerable simplification of the many-body problem which is advantageous for numerical implementations. However, all the complications of the true many-body system are hidden in the effective potential of the KS-system. This KS potential is typically expressed as a sum of the external potential of the interacting system of interest and the Hartree-exchange-correlation (Hxc) potential containing implicitly the many-body effects of the interacting system. The KS formalism is equally applicable in ground state and time-dependent DFT but in this work we will focus on the calculation of excitation energies which are obtained in TDDFT using a linear response formalism. For this purpose, it is enough to know the functional derivative of the Hxc potential with respect to the density which yields a quantity known as the Hxc-kernel. The simplest possible approximation for the Hxc kernel is the adiabatic local-density approximation (ALDA), for which the kernel is local in space and time. Although this approximation has been used successfully Ullrich:book it has a number of important deficiencies, such as the inability to reproduce Born-Oppenheimer surfaces of excited states in dissociating molecules GritsenkoBaerends ; GritsenkoBaerends2 .

When a molecule separates into fragments its excitation energies should approach those of the separate fragments. This behavior is not reproduced by the ALDA since upon dissociation the gap between the bonding and anti-bonding KS eigenvalues decreases exponentially fast with the bond distance, and the ALDA kernel is unable to correct for this thereby rendering many of the excitation energies to become zero in the dissociation limit. To correct for this, asymptotic corrections have been devised GritsenkoBaerends ; GritsenkoBaerends2 that introduce exponentially growing terms in the kernel that compensate for the closing of the bonding-antibonding gap. Although such corrections can reproduce the main features of the exact bonding curve for the lowest excited state GritsenkoBaerends ; GritsenkoBaerends2 , there is no systematic way to construct such functionals. Other more systematic approximations often rely on perturbative expansions, which makes them questionable in the multi-configuration regime required to describe molecular dissociation.

In recent work Lani an approximate kernel was derived within the framework of so-called strictly correlate electrons (SCE). This is a ground state DFT formalism that becomes exact in the limit of very large two-body interactions. When the simplest approximation within this formalism is applied within the adiabatic approximation an approximate Hxc kernel can be derived . This so-called adiabatic SCE (ASCE) kernel was shown to have a number of desirable features. It was shown to obey the so-called zero-force theorem Ullrich:book ; Mundt:PRA and it was shown that in the case of molecular dissociation it exhibits an exponential growth with the bond distance Lani . The kernel therefore displays a very non-local spatial behavior that has the potential to cure the deficiency of the ALDA kernel for molecular dissociation. We recently investigated the ASCE kernel Cort in a model system for which the exchange-correlation kernel can be obtained exactly for various two-body interaction strengths. It was found that the leading order and the next to leading order of the asymptotic expansion for the exact Hxc kernel in terms of the interaction strength agreed with that one predicted by the adiabatic SCE formalism. This result shows that the SCE formalism is a promising method for describing the linear response properties in the strong interaction limit. Moreover, these terms were also shown to be frequency independent in the exact theory such that the adiabatic approximation in this limit is in fact exact. In view of these favorable properties of the ASCE kernel the natural question arises whether this kernel can be used to correctly predict the excitation energies of dissociating molecules. Answering this question is the main aim of the present work.

To attack this problem, we developed a simplified one-dimensional model of a diatomic molecule having the main physical characteristics of a real three-dimensional hydrogen molecule and for which we can perform analytical and numerical calculations for arbitrary bond distance and interaction strength. In particular the KS orbitals and eigenvalues are known analytically, a feature that is very desirable as it provides an analytic expression for the KS gap upon dissociation. The model is used to benchmark the performance of the ASCE kernel as well as to discuss many features of the SCE formalism in the limit of large interactions.

The paper is organized as follows: In Sec. II we give a brief introduction to the main elements of SCE theory that we need. In Sec. III we introduce the model system and discuss its properties. In Sec. IV we discuss the ASCE kernel for our model density and obtain the excitation energies. In Sec. V we present our conclusions.

II Density-functional theory in the large interaction limit

The main motivation of this work is to benchmark the recently proposed approximations for the exchange-correlation (xc) potential and xc-kernel based on the so-called theory of strictly correlated electrons Lani ; Cort . To provide a self-contained minimal background for the reader we briefly review some basic aspects of DFT. Our starting point is the time-independent N-body Hamiltonian of a system which we write as Barth:review :

[TABLE]

where $\hat{T}$ is the kinetic energy and $\hat{W}$ the two-body interaction, the strength of which is regulated by a real parameter $\lambda$. Finally, $\hat{V}_{\lambda}$ represents the external potential and is the sum of one-body potentials $v_{\lambda}(\mathbf{r})$. The latter potential depends on the interaction strength $\lambda$ via the requirement that for each value of $\lambda$ the same electronic density $n(\mathbf{r})$ is obtained from the ground state of Eq.(1). This makes $v_{\lambda}$ a functional of the density via the Hohenberg-Kohn theorem HK and we will therefore sometimes write $v_{\lambda}[n]$ to stress this fact when necessary.

Typically the Hamiltonian is given at $\lambda=1$ with a known external potential and the key many-body problem is to solve for its eigenstates. However, consideration of the full $\lambda$-dependence is useful in formal derivations in DFT and is particularly relevant for the present work. An important limit is obtained by taking $\lambda=0$, in which case the system becomes non-interacting while retaining the density of the interacting system. This system is denoted as the Kohn-Sham (KS) system and its external potential is commonly denoted by $v_{s}(\mathbf{r})$. The ground state of the KS system is a Slater determinant consisting of KS orbitals $\varphi_{i}$ satisfying

[TABLE]

where $\sigma$ is a spin index. The KS equations are a device for obtaining the density of the interacting system by solving one-particle equations. However, to make the procedure useful we need to make a connection to the interacting system which we will take at a general interaction strength $\lambda$. To do this we define the Hxc potential as

[TABLE]

A given approximation for this quantity allows us to obtain the density of the interacting system by using the potential $v_{\textrm{KS}}[n,v_{\lambda}]=v_{\lambda}+v_{\textrm{Hxc}}^{\lambda}[n]$ in Eq.(2) instead of $v_{s}[n]$ where $v_{\lambda}$ is a given and known potential of the interacting system at interaction strenght $\lambda$ (which is commonly taken to be $\lambda=1$ but we would like here to use a general interaction strength for the discussion below) Ruggenthaler:2013 . The central object of DFT is therefore the Hxc potential. This quantity in turn is given by the functional derivative of the Hxc-energy with respect to the density $v_{\textrm{Hxc}}^{\lambda}(\mathbf{r})=\delta E_{\textrm{Hxc}}^{\lambda}/\delta n(\mathbf{r})$. The Hxc-energy can be obtained from

[TABLE]

where we defined

[TABLE]

where $\Psi_{\lambda}[n]$ is the ground-state of Hamiltonian Eq.(1). The quantity $W_{\lambda}$ has been studied in limiting cases. For small values of $\lambda$ it is accessible via perturbation theory while in the limit of large values of $\lambda$ there is an asymptotic expansion that is derived from SCE theory. This expansion has the form GoriGiorgi:2009

[TABLE]

where the leading term is the interaction energy of the strictly correlated electrons and the next term arises from their zero-point energy (ZPE) in vibrations around their equilibrium positions. Correspondingly the asymptotic expansion of the Hxc energy for large $\lambda$ is given by:

[TABLE]

as can be checked by differentiation with respect to $\lambda$ and comparison to Eq.(6). This expression further introduces a density functional $E_{2}[n]$ the relevance of which will become clear later. The functional derivative with respect to the density gives an expansion of the Hxc-potential in powers of $\sqrt{\lambda}$

[TABLE]

which is valid for large value of $\lambda$. A very interesting point is that, at least for one-dimensional systems many-electron systems, the two leading terms are explicitly known functionals of the density and can be calculated explicitly in a rather simple way from so-called co-motions functions Seidl1 . Before we discuss the applicability of this expansion let us further define the adiabatic Hxc kernel by

[TABLE]

which according to Eq.(8) has the expansion

[TABLE]

The first term on the right hand side represents the so-called adiabatic SCE kernel $\lambda f^{\textrm{ASCE}}_{\textrm{Hxc}}$ which has been studied in detail in Refs.Lani ; Cort which we refer to for more details. So far our discussion has been very general and, apart from the adiabatic approximation to the time-dependent kernel of TDDFT in Eq.(10), no approximations have been used. The main question is, however, how reliable the asymptotic expansions in Eqs.(7) and (8) are for values close to the physically relevant interaction strenght $\lambda=1$. Since the expansion is asymptotic, retaining higher order terms typically worsen the approximation unless we increase the value of $\lambda$. This means that for values of $\lambda$ close to one the best approximation may be obtained by only retaining the term $v_{\textrm{SCE}}$. Indeed, it was pointed out in Ref.Malet that in this interaction regime adding the ZPE contribution generally will give a worsening of the result. It was found that at the lowest SCE level for a model one-dimensional diatomic molecule the bonding curve is correct at large separation but inaccurate at equilibrium separation, while adding the ZPE contribution gives an overall worse result for the bonding curve. The asymptotic expansion can therefore not been applied as such and consequently Ref.Malet considers various amendments. A similar conclusion was obtained from our previous work on the model system of a quantum ring Cort where we found the ZPE contribution to worsen the results at smaller interaction strengths. This work was done for a homogeneous system in which we mainly studied the properties of the kernel itself. In the present work we extend that work to an inhomogeneous model system in which again the kernel will be at the focus of attention. The equations derived in the present section will be referenced in later sections.

III The molecular model

III.1 Definition of the model

For our description of the simplified molecular model we consider two electrons with spatial coordinates $x_{1}$ and $x_{2}$ both in the domain $\left[-\frac{L}{2},\frac{L}{2}\right]$ on a ring of length $L$. The Hamiltonian of our system is given by

[TABLE]

where the first two terms are the kinetic energy of each electron, $v_{\lambda}$ is the one body external potential and $w(x)=\lambda\cos^{2}(\pi x/L)$ is the electron-electron repulsion. We impose periodic boundary conditions such that the particles effectively move on a ring which is commonly referred to as a quantum ring (QR) system Ruggenthaler:2013 . The strength of the interaction $\lambda$ is a parameter which we will take to be positive. The interaction tends to keep particles on opposing parts of ring and has a convenient form for numerical considerations. In accordance with Eq.(1) the potential $v_{\lambda}$ is chosen in such a way that for each value of $\lambda$ the same ground state density is produced. For our model it turns out to be useful to specify the external potential at $\lambda=0$ which corresponds to the KS-potential. In this way we can choose the potential in such a way that we obtain an analytic solution for the KS orbitals. The potential at all other interaction strengths, including the physically relevant case $\lambda=1$, is subsequently determined by the constraint that the density is the same for all values of $\lambda$ as we will discuss in more detail later.

III.2 The Kohn-Sham system

The KS system is obtained from Eq.(11) by taking $\lambda=0$ and we adopt the common notation of denoting the KS -potential by $v_{s}$, i.e. $v_{s}=v_{\lambda=0}$. In this limit the Hamiltonian of Eq.(11), which we now denote by $\hat{H}_{s}$, attains the form

[TABLE]

We now specify an explicit choice for $v_{s}$ which we take to be

[TABLE]

where $V_{0}$ is a constant with units of energy. This potential has two minima located at $x_{0}=\pm L/4$ where $v_{s}(x_{0})=0$ and is positive everywhere else. The ground-state density has two maxima at the potential minima and therefore represents a simple model of a diatomic molecule in which the atoms are separated by a bond distance $L/2$.

We want to use this model to describe molecular dissociation and therefore vary the bond length $L$. While doing this we want to guarantee that the width of each atomic density remains fixed upon separation, which can be achieved by requiring that the curvature of the potential at $x_{0}=\pm L/4$ is independent of $L$. This condition reads

[TABLE]

where $\alpha$ length independent which gives $V_{0}=\alpha\left(L/(4\pi)\right)^{2}$ for an arbitrary $\alpha$ (in this paper we will always take $\alpha=1$). The KS orbitals of our system satisfy the eigenvalue equation

[TABLE]

where we added a symmetry label $\pm$ for orbitals that are even or odd with respect to reflection in the origin, i.e. $\varphi_{l}^{\pm}(-x)=\pm\varphi_{l}^{\pm}(x)$. These equations must be solved together with the boundary conditions $\varphi_{l}(-L/2)=\varphi_{l}(L/2)$ and the same for their derivatives. It is convenient to define the dimensionless coordinate $z=2\pi x/L$ and use the explicit form of the potential to rewrite Eq.(15) as

[TABLE]

where we have defined the following constants

[TABLE]

We recover the KS orbitals from $\varphi_{l}^{\pm}(x)=M_{l}^{\pm}(2\pi x/L)$. Equation (16) is the well-known Mathieu equation and its eigenfunctions and eigenvalues have been intensively studied NIST . The functions $M_{l}^{+}$ and $M_{l}^{-}$ are commonly denoted as the Mathieu-cosine $C_{l}$ and the Mathieu-sine $S_{l}$ functions respectively, while the values $a_{l}^{\pm}$ are called the Mathieu characteristic values. The convention is that the label of the even states start at $l=0$ whereas the labels of the odd states start at $l=1$. The Mathieu functions satisfy $M_{l}^{\pm}(z+\pi)=(-1)^{l}M_{l}^{\pm}(z)$ and are therefore $2\pi$-periodic. They are commonly normalized as follows

[TABLE]

Correspondingly the normalized (to one) KS orbitals are expressed in terms of Mathieu functions as

[TABLE]

while the Kohn-Sham eigenenergies can be recovered from the Mathieu characteristic values by means of Eq.(17). In FIG. 1 we plot the KS potential and the ground state density for different bond distances to illustrate the main features that we mentioned, in particular the fact that the width of the maxima becomes independent of the bond distance for large $L$. Although we are not particularly interested in the case of very short bond distances we note that in the limit $L\rightarrow 0$ the parameter $\nu$ becomes equal to zero and the ground state KS orbital is given by the constant function $\varphi_{0}(x)=1/\sqrt{L}$ representing a system of constant density. We will not investigate this limit in detail; a homogenous QR at various interaction strength has been studied in detail in Ref.Cort .

In Fig.2 we plot the ground state and the first few excited state KS orbitals. Of particular interest for our later discussion of the Hxc kernel is the lowest pair of bonding and anti-bonding states represented by the pair of Mathieu functions $C_{0}$ and $S_{1}$.

The corresponding energy gap between the KS eigenvalues closes exponentially fast with increasing bond distance:

[TABLE]

where we used the asymptotic expansion for the Mathieu characteristic value given in (A). We remind the reader that $\nu$ is an increasing function of $L$ given by Eq.(18). The density in the bond midpoint has a similar exponential decay (see Eq.(71) in Appendix A) given by

[TABLE]

The knowledge of this precise behavior of the KS gap as well as the density in the bond midpoint will facilitate considerably the calculation of the excitation energy from the ASCE kernel in Sec.IV.

III.3 Exact solution of the model

After having considered the model in the KS limit we will now consider the case of finite interaction strength $\lambda$. The potential $v_{\lambda}$ in Eq.(11) can not be obtained analytically except in some limiting cases that we will discuss below. We therefore obtain $v_{\lambda}$ directly from the constraint that the density is independent of $\lambda$ using the numerical algorithm outlined in Ref.Nielsen:EPJB2018 . In our case the density is given by the ground state KS orbital from Eq.(20) to be

[TABLE]

for all $\lambda$ where we remind the reader that $\nu$ depends on $L$ via Eq.(18). Therefore for a given value of $L$ we have the numerical task to find $v_{\lambda}$ for a range of interaction strengths of interest. The ground state is a spin singlet state and consequently we will mostly be interested in the singlet excited states. The singlet wave function has the structure

[TABLE]

where $\sigma_{i}$ for $i=1,2$ are spin variables and where the spatial part of the wave function is symmetric $\psi(x_{1},x_{2})=\psi(x_{2},x_{1})$ to ensure anti-symmetry of the full space-spin wave function. To obtain deeper insight in the results we will also derive analytic results in the regime of large bond distance $L/2$ for fixed interaction strength $\lambda$ which is the common molecular dissociation regime and the complementary regime of large interaction strength $\lambda$ for fixed bond distance $L/2$ which is the SCE regime. We will start in the next subsection with the first regime.

III.3.1 Large bond distance for fixed interaction strength

We first consider the regime of large bond distance $L/2$ at fixed values of $\lambda$. In this regime the molecule is typically dissociated in two one-electron atoms (unless the interacting strength $\lambda$ is very small such that there are contributions from the ionic states with two or zero electrons on each atom). For a one-electron atom the KS potential is equal to the true external potential and therefore we have $v_{\lambda}(x)=v_{s}(x)$ for $x$ in the neighbourhood of each atom at large separation. The ground state atomic orbitals $A(x)$ and $B(x)$ on atoms $A$ and $B$ are localized around $x=\pm L/4$ and can be expressed in terms of the first bonding and anti-bonding molecular KS orbitals as $A(x)=(\varphi_{0}^{+}(x)+\varphi_{1}^{-}(x))/\sqrt{2}$ and $B(x)=(\varphi_{0}^{+}(x)-\varphi_{1}^{-}(x))/\sqrt{2}$ (see for example Fig.2). The exact ground-state (GS) wave function for the large bond distance limit is the well known Heitler-London (HL) wave function

[TABLE]

The ground state energy is given by

[TABLE]

where we used that $\varepsilon_{0}^{+}=\varepsilon_{1}^{-}$ in the large $L$ limit and the asymptotic expansion of the Mathieu characteristic values in Appendix A. This result is easy to understand. Since at the atomic positions $x_{0}=\pm L/4$ we have that $v_{s}^{{}^{\prime\prime}}(x_{0})=\alpha$ the potential around each atom is given by $v_{s}(x)=\alpha(x-x_{0})^{2}/2$ which corresponds to a harmonic well with harmonic frequency $\sqrt{\alpha}$. Each atomic oscillator has ground state energy $\sqrt{\alpha}/2$ thereby adding up to the molecular ground state energy $\sqrt{\alpha}$.

Let us now consider the first excited state which in the large $L$ limit is given by

[TABLE]

The orbitals used in this expression are displayed in Fig.2. For large $L$ the states $\varphi_{0}^{+}$ and $\varphi_{1}^{-}$ become degenerate and the same is true for the states $\varphi_{1}^{+}$ and $\varphi_{2}^{-}$. These orbitals can be used to construct localized ground and excited state atomic orbitals from the combinations $\varphi_{0}^{+}\pm\varphi_{1}^{-}$ and $\varphi_{1}^{+}\pm\varphi_{2}^{-}$ if desired. The energy of the two-particle state of Eq.(28) is given by

[TABLE]

Again it is straightforward to interpret the energy. The system is a superposition of two states in which one atom is a ground state oscillator with energy $\sqrt{\alpha}/2$ and the other one a first excited oscillator with energy $3\sqrt{\alpha}/2$ giving a total molecular energy of $2\sqrt{\alpha}$.

To judge the accuracy of these limiting wave functions we plot the exact $\psi_{\lambda}$ for $\lambda=1$ and $L=9$ and $21$ (corresponding to bond lengths $4.5$ and $10.5$) in Fig. 3. We see that for $L=21$ the wave functions Eq.( 26) and (28) are a good approximation to the true wave functions (as we also checked numerically). At $L=9$ the system still has a considerable density at the bond midpoint and the HL-type wave functions are a less good approximation.

Finally we compare in Fig.4 the exact external potential $v_{\lambda}$ to $v_{s}$. We see that around the atoms both potentials agree but that around the bond midpoint there is a considerable deviation. This amounts to a peak in the Hxc-potential $v_{\textrm{Hxc}}^{\lambda}=v_{s}-v_{\lambda}$ at the bond midpoint. This is a well-known feature of the Hxc-potential BuijseBaerends and is related to the so-called left-right correlation in the system. We refer to the cited reference for a more in-depth discussion.

III.3.2 Large interaction strength at fixed bond distance

We now turn our attention to the complementary regime of larger interaction strength $\lambda$ for fixed bond distance. This is the regime in which SCE become exact. From our numerical work we find that in this limit the two-particle wave function localizes in a region where $|x_{1}-x_{2}|\approx L/2$ as displayed in FIG. 5.

This is in accordance with SCE theory which tells that in the very strong interaction limit the position of a single electron determines the positions of the remaining electrons uniquely. For this reason it is convenient to introduce the center of mass $R=(x_{1}+x_{2})/2$ and relative coordinate $r=x_{1}-x_{2}$, where $R\in[-L/2,L/2]$ and $r\in[-L,L]$. The Hamiltonian (11) in the new coordinates attains the form

[TABLE]

We want to give an explicit approximate expression of the hamiltonian (11) for the limit $\lambda\to\infty$ for any fixed bond distance $L/2$. Since the wave function is localized around the lines $r=\pm L/2$ it is natural to expand the external potential $v_{\lambda}$ around these values. For example, for $r=L/2$ we have to second order

[TABLE]

where we defined

[TABLE]

with an essentially identical result for the expansion around $r=-L/2$, and where we used the property $v_{\lambda}(x)=v_{\lambda}(x+L/2)$ in the definitions of $\bar{v}_{\lambda}$ and $\beta_{\lambda}$ and in the cancellation of the linear term. With the expansion of Eq.(31) the Hamiltonian becomes

[TABLE]

with a similar expansion around $r=-L/2$. We see that this Hamiltonian becomes separable when we neglect the term $\beta_{\lambda}$. However, the two-body interaction has form $w(r)=\lambda(\pi/L)^{2}(r-L/2)^{2}$ around $r=L/2$ and the question is therefore whether we can neglect $\beta_{\lambda}$ compared to $\lambda(\pi/L)^{2}$. From our calculation we find that $v_{\lambda}$ and therefore also $\beta_{\lambda}$ converges to a finite value for large $\lambda$. Therefore for fixed $L$ and large enough $\lambda$ we can neglect $\beta_{\lambda}$ and the system becomes approximately separable. If we write the wave function in this limit as $\Psi_{\lambda}(r,R)=\chi_{\lambda}(r)\varphi_{\lambda}(R)$ then its factors are determined from the equations

[TABLE]

These equations determine all the eigenstates in the large $\lambda$ limit. Let us, however, focus on the ground state and take $\chi_{\lambda}$ and $\varphi_{\lambda}$ to be ground states of their corresponding Hamiltonians. The ground state density is then obtained from

[TABLE]

The function $|\chi_{\lambda}(r)|^{2}$ becomes very narrowly peaked around $r=\pm L/2$ as $\lambda$ becomes very large. We can therefore normalize it such that for the limit that $\lambda\rightarrow\infty$

[TABLE]

from which we obtain, using Eq.(37), that for large interaction strength

[TABLE]

The ground state density is also given by $n(x)=2|\varphi_{0}^{+}(x)|^{2}$ in which $\varphi_{0}^{+}(x)$ solves Eq.(15). Comparison of this equation to Eq.(35) then immediately yields that

[TABLE]

and $\varphi_{0}^{+}(x)=\sqrt{2}\,\varphi_{\lambda}(x-L/4)$. From our derivation we therefore deduce that in our system

[TABLE]

A comparison with the general Eq.(8) from SCE theory shows that in our case $v_{\textrm{SCE}}$ and $v_{\textrm{ZPE}}$ are zero and that $v_{2}(x)=v_{s}-v_{\lambda}=3v_{s}(x)/4$. The fact that $v_{\textrm{SCE}}$ and $v_{\textrm{ZPE}}$ vanish can also be directly derived from SCE theory and is a consequence of the symmetry of our system. In Fig. 6 we compare $v_{\lambda}$ to $v_{s}/4$ for various large values of $\lambda$ and note a good agreement between them with the exception of some deviations around the bond midpoint. This discrepancy becomes smaller for higher values of $\lambda$.

Let us now consider the energies of the system. The eigenenergies of the two-particle state are given by $E=\epsilon+\tilde{\epsilon}$ where $\epsilon$ and $\tilde{\epsilon}$ are the eigenvalues of Hamiltonians in Eq.(35) and Eq.(36). From the fact that $v_{\lambda}=v_{s}/4$ in Eq.(35) we see that the eigenvalues $\epsilon$ are half of the KS eigenvalues of Eq.(15). These eigenvalues correspond to an excitation which only involves a change of the center-of-mass wave function without changing the relative wave function. The eigenvalues $\tilde{\epsilon}$ are calculated from Eq.(36). The transformation $z=\pi r/L$ transforms this Hamiltonian to

[TABLE]

where

[TABLE]

Eq.(42) is again the Mathieu equation with this time a parameter $q$ that depends on the interaction strength. The eigenvalues in the limit of large interactions have the form

[TABLE]

which is a harmonic spectrum with harmonic frequency $\omega_{\lambda}=2\pi\sqrt{\lambda}/L$. These excitations of involve a change of the relative wave function and represent the zero point vibrations of the strictly correlated electrons of SCE theory. The lowest excitation energy for this mode is therefore $\omega_{\lambda}$. This will be relevant of our discussion of the excitation energy obtained from the ASCE kernel.

IV The adiabatic SCE kernel

IV.1 Definition and properties

We have in studied in detail the excitation properties of our model system in two different regimes. We will now investigate the adiabatic SCE kernel. As was discussed below Eq.(10) the ASCE kernel is defined as

[TABLE]

The SCE potential vanishes for our system, but its functional derivative does not. As was discussed in detail in Refs.Lani ; Cort it is explicitly given by the expression

[TABLE]

where $\theta$ is the usual Heaviside function and $w(x)$ the two body interaction. The function $f(x)$ is the so-called co-motion function which specifies the position of another electron given the position of a reference electron. For our system the co-motion function attains the simple form

[TABLE]

If we define the function $\mathcal{P}(x)$ to be

[TABLE]

then the integrand contains $d\mathcal{P}/dx$ and we can obtain $f^{\textrm{ASCE}}$ by partial integration while usefully manipulating the results using the fact that $\mathcal{P}(x)-\mathcal{P}(0)$ is an odd function. In the quadrant $x,x^{\prime}>0$ we obtain

[TABLE]

while in the quadrant $x<0,x^{\prime}>0$ we have

[TABLE]

The function in the remaining quadrants is determined from the symmetry $f^{\textrm{ASCE}}(x,x^{\prime})=f^{\textrm{ASCE}}(-x,-x^{\prime})$. For our system the function $\mathcal{P}(x)$ can be written more explicitly as:

[TABLE]

where for our two-body potential $w^{\prime\prime}(L/2)=2\pi^{2}/L^{2}$. In the Appendix B we show that

[TABLE]

for $x\neq 0$. This equation implies that for large values of $L$ the kernel assumes the form

[TABLE]

for $x,x^{\prime}\neq 0$. The function exhibits plateaux of height $\mathcal{P}(0)$ in the quadrants in which both coordinates have the same sign and is zero otherwise. In the Appendix B we show that this height grows exponentially fast with $L$ according to

[TABLE]

(we remind the reader that $\nu$ depends on $L$ according to Eq.(18)). With these results we are ready to calculate excitation energies from the ASCE kernel.

IV.1.1 Lowest excitation energy

We now address the issue of calculation the excitation energy of the system. To make our point it is sufficient to restrict ourselves to the so-called small matrix approximation Ullrich:book in which the singlet excitation energy $\Omega$ from an occupied state $i$ to an unoccupied state $a$ is given by

[TABLE]

where $\omega_{ia}=\epsilon_{a}-\epsilon_{i}$ is the difference in KS energies. and

[TABLE]

where $\Phi_{ia}(x)=\varphi_{i}(x)\varphi_{a}(x)$ is an excitation function (in which we take the orbitals to be real for simplicity) and $f_{\textrm{Hxc}}$ the Hxc kernel which we took in an adiabatic approximation relevant to the discussion below. In our particular case we consider the excitation from the lowest KS orbital $\varphi_{0}^{+}$ to $\varphi_{1}^{-}$. For easy of notation and to be in accordance with adopted language we denote the orbitals by the gerade and ungerade sigma orbitals $\sigma_{g}(x)$ and $\sigma_{u}(x)$ and their eigenvalues by $\epsilon_{g}$ and $\epsilon_{u}$. We know that in the dissociation limit the KS gap $\omega_{gu}$ vanishes . The excitation energy is therefore given by

[TABLE]

In the ALDA this expression vanishes as the kernel can not compensate for the decay of the KS gap. However, as we will show now, the ASCE kernel (we remind the reader of Eq.(10) ) will lead to a finite contribution. The matrix element in the large separation limit is readily calculated from Eq.(57) to be

[TABLE]

where we used the symmetry and normalization of the KS orbitals. If we use this in Eq.(61) we find that in the large $L$ limit

[TABLE]

For our system we have $w^{\prime\prime}(L/2)=2\pi^{2}/L^{2}$ and we obtain $\Omega=2\pi\sqrt{\lambda}/L$ which is exactly the harmonic frequency of the zero point oscillation of Eq.(45). We therefore deduce that the excitations that we recover from the ASCE kernel are exactly the ones that correspond to the zero point oscillations. With hindsight this may not be surprising as, after all, the zero point oscillations represent an always present set of excitations in SCE theory. Note that in the derivation of Eq.(63) it is important to consider a fixed but arbitrary large $L$ and then take the limit $\lambda\rightarrow\infty$, i.e. the standard SCE regime, and not the other way around otherwise $\Omega=2\pi\sqrt{\lambda}/L\rightarrow 0$.

IV.2 The ASCE kernel in the conventional molecular dissociation regime

In the previous subsection we found that in the limit that the interaction strength $\lambda$ becomes very large at fixed bond distance $L/2$ the lowest excitation energy is that of the lowest zero point oscillation of the strictly correlated electrons, and in that regime the ASCE kernel gives an exact result. Let us now see how the ASCE kernel performs in the opposite regime in which the bond distance becomes large at fixed interaction strength, in particular for the chemically relevant case of interaction strength $\lambda=1$. This is the conventional dissociation regime as commonly studied in bond breaking in chemistry. Note that we now apply the ASCE kernel outside its formal range of applicability and therefore the approximation becomes uncontrolled. The consideration is nevertheless illuminating as it illustrates the reasons for the breakdown of the approximation. For $\lambda=1$ the matrix element Eq. (62) of the ASCE kernel is given by $\mathcal{P}(0)/4$ and we have for the lowest excitation energy

[TABLE]

Let us compare this to the exact excitation energy

[TABLE]

as follows directly from Eqs.(27) and (29). We remind the reader that the parameter $\alpha$ (see Eq.(14)) is given by the curvature of the external potential at its minima (as $v_{s}$ becomes the true external potential around the atoms in the dissociation limit). Since upon dissociation the separate atoms become independent single particle oscillators, Eq.(65) is a natural result. If we consider the ASCE approximation, on the other hand, we see that according to Eq.(64) the lowest excitation energy is determined solely by the curvature of the interaction potential $w^{\prime\prime}(L/2)$. This is because, by using the ASCE kernel, we pretend that the separated atoms still behave as strictly correlated electrons with an excitation energy determined by the zero point oscillations. This is the wrong physical picture in this regime and therefore the ASCE approximation fails to describe the right physics. In fact, in our system $w^{\prime\prime}(L/2)=2\pi^{2}/L^{2}\rightarrow 0$ for $L\rightarrow\infty$ and therefore the ASCE excitation energy becomes zero in the dissociation limit. For other forms of the two-body interaction this may not be the case but this does not change our conclusion regarding the physical picture. The ASCE approximation is therefore not an improvement over the ALDA in the dissociation regime. Both approximations attain the wrong dissociation limit; in the case of the ALDA the excitation energy becomes zero whereas in the case of the ASCE approximation the excitation energy is determined by the two-body interaction potential rather than by the external potential of the separated atoms. This result is not surprising as we have used the ASCE kernel outside its regime of applicability. The ASCE kernel is therefore not of use if one is interested in regime of large bond length at intermediate interaction strength which is the relevant case for bond breaking in most common chemical applications. To correct these problems within the present formalism a natural way to proceed would be include ZPE and higher order kernels in the expansion of the Hxc kernel as was done in Ref.Cort . However, that work showed that the extra terms lead to worse approximation than just the ASCE approximation for low interaction strengths, as is typical for an asymptotic expansion. The description of the conventional dissociation regime using density-functional methods therefore remains a challenging task.

V Conclusions

In this work we studied the properties of an approximate adiabatic Hxc kernel based on the theory of strictly correlated electrons. To benchmark this approximation we studied a numerically and analytically solvable system which is able to simulate the main features of a dissociating molecule. We studied in detail the two-particle eigenstates in various limits and calculated the excitation spectrum in the limit of large interaction strength. The ASCE kernel was shown to reproduce the so-called zero-point oscillation part of the spectrum. The attainment of this exact result shows that the ASCE kernel becomes exact in the this regime as we also concluded from earlier work Cort . However, most current interest in molecular dissociation in chemistry is devoted to the complementary regime of large bond distance at intermediate interaction strength. In this regime the ASCE kernel is not suitable for obtaining the excitation spectrum. We conclude that the description of molecular dissociation based on functionals founded on SCE theory remains a challenge for the future.

Appendices

Appendix A Properties of Mathieu functions

In this Appendix we describe a few useful properties the Mathieu functions and their characteristic values that we use in the main text. Many properties of these functions can be found in NIST . The Mathieu characteristic values have the following expansion for large $q$ (where $q$ is the parameter in the Mathieu equation)

[TABLE]

The difference $a_{l+1}^{-}(q)-a_{l}^{+}(q)$ is exponentially small in the large $q$ limit NIST

[TABLE]

We note that in our previous work Cort we denoted $a_{l+1}^{-}$ by $a_{l}^{-}$ in the asymptotic formula Eq.(66) which amounts to a different labeling convention for the characteristic values. Here we stick to a more common convention.

For this work we need an accurate representation of $C_{0}(z;q)$ for small values of $z$. A representation that is valid for large $q$ in the interval $|z|<\pi/2$ is given by

[TABLE]

To the determine this function we also need to know its prefactor $C_{0}(0;q)$ which is given by Gertrud

[TABLE]

This equation involves yet another prefactor which is obtainable from Sips’ expansion Cort and given in leading order in $q$ to be

[TABLE]

In particular we find that

[TABLE]

from which we obtain the density in the bond midpoint of Eq.(24).

Appendix B Analysis of the function $\mathcal{P}(x)$

We study here the properties of the function $\mathcal{P}(x)$ defined in Eq.(55) we rewrite here as

[TABLE]

where we used the explicit form of the density and we defined

[TABLE]

It will be convenient to further introduce the functions

[TABLE]

and $\mathcal{J}(\nu)=\mathcal{I}(\pi/2,\nu)$ such that we can write

[TABLE]

where we used the symmetry of the integrand. Using then the asymptotic expansion of Mathieu functions functions Eq.(68), the $f(t;\nu)$ reads

[TABLE]

For $\nu$ very large this function has its main contributions from $z=0$ and we can approximate

[TABLE]

which inserted into Eq.(74) gives

[TABLE]

and consequently

[TABLE]

From this we can evaluate $\mathcal{P}(0)$. Using Eq.(71) and (73) we find

[TABLE]

which yields Eq.(58). Finally we consider the quantity

[TABLE]

and therefore

[TABLE]

This behavior of the function $\mathcal{P}(x)$ is illustrated in Fig.8 where we plotted $\mathcal{P}(x)/\mathcal{P}(0)$. In this figure we clearly see the step appearing with increasing $L$.