Excited states from range-separated density-functional perturbation theory

We explore the possibility of calculating electronic excited states by using perturbation theory along a range-separated adiabatic connection. Starting from the energies of a partially interacting Hamiltonian, a first-order correction is defined with two variants of perturbation theory: a straight-forward perturbation theory, and an extension of the G{\"o}rling--Levy one that has the advantage of keeping the ground-state density constant at each order in the perturbation. Only the first, simpler, variant is tested here on the helium and beryllium atoms and on the dihydrogene molecule. The first-order correction within this perturbation theory improves significantly the total ground-and excited-state energies of the different systems. However, the excitation energies are mostly deterio-rated with respect to the zeroth-order ones, which may be explained by the fact that the ionization energy is no longer correct for all interaction strengths. The second variant of the perturbation theory should improve these results but has not been tested yet along the range-separated adiabatic connection.


Introduction
In density-functional theory (DFT) of quantum electronic systems, the most widely used approach for calculating excitation energies is nowadays linear-response timedependent density-functional theory (TDDFT) (see, e.g., Refs. [1,2]). However, in spite of many successes, when applied with the usual adiabatic semilocal approximations, linear-response TDDFT has serious limitations for describing systems with static (or strong) correlation [3], double or multiple excitations [4], and Rydberg and charge-transfer excitations [5,6]. Besides, the Hohenberg-Kohn theorem [7] states that the time-independent ground-state density contains all the information about the system implying that time-dependence is in principle not required to describe excited states.
Several time-independent DFT approaches for calculating excitation energies exist and are still being developed. A first strategy consists of simultaneously optimising an ensemble of states. Such an ensemble DFT was pioneered by Theophilou [8] and by Gross, Oliveira and Kohn [9] and is still a subject of research [10][11][12][13], but it is hampered by the absence of appropriate approximate ensemble functionals. A second strategy is to apply the self-consistent field (SCF) method to directly optimise a single excited state. This approach was started by Gunnarsson and Lundqvist [14], who extended ground-state DFT to the lowest-energy state in each symmetry class, and developed into the pragmatic (multiplet-sum) SCF method [15,16] (still in use today [17]) and related methods [18][19][20]. Great efforts have been made by Nagy, Görling, Levy, Ayers and others to formulate a rigorous self-consistent DFT of an arbitrary individual excited state [21][22][23][24][25][26][27][28][29][30][31][32][33] but a major difficulty is the need to develop approximate functionals for a specific excited state (see Ref. [34] for a proposal of such excited-state functionals). A third strategy, first proposed by Grimme, consists of using configuration-interaction (CI) schemes in which modified Hamiltonian matrix elements include information from DFT [35][36][37][38].
Finally, a fourth possible approach, proposed by Görling [39], is to calculate the excitation energies from Görling-Levy (GL) perturbation theory [40,41] along the adiabatic connection using the non-interacting Kohn-Sham (KS) Hamiltonian as the zeroth-order Hamiltonian. In this approach, the zeroth-order approximation to the exact excitation energies is provided by KS orbital energy differences (which, for accurate KS potentials, is known to be already a fairly good approximation [42][43][44]). It can be improved upon by perturbation theory at a given order in the coupling constant of the adiabatic connection. Filippi, Umrigar, and Gonze [45] showed that the GL first-order corrections provide a factor of two improvement on the KS zeroth-order excitation energies for the He, Li + and Be atoms when using accurate KS potentials. For (nearly) degenerate states, Zhang and Burke [46] proposed to use degenerate second-order GL perturbation theory, showing that it works well for a simple one-dimensional model. This approach is conceptually simple as it uses the standard adiabatic connection along which the ground-state density is kept constant (in contrast to approaches introducing generalised adiabatic connections keeping an excited-state density constant [21,22,24,29]). In spite of promising early results, this approach has not been pursued further, perhaps because it can be considered an approximation to TDDFT [47].
In this work, we explore further this density-functional perturbation-theory approach for calculating excitation energies, introducing one key modification in comparison to the earlier work of Refs. [39,45]: as a zeroth-order Hamiltonian, instead of using the non-interacting KS Hamiltonian, we use a partially interacting Hamiltonian incorporating the long-range part only of the Coulomb electron-electron interaction, corresponding to an intermediate point along a range-separated adiabatic connection [48][49][50][51][52][53]. The partially interacting zeroth-order Hamiltonian is of course closer to the exact Hamiltonian than is the non-interacting KS Hamiltonian, thereby putting less demand on the perturbation theory. In fact, the zeroth-order Hamiltonian can already incorporate some static correlation.
The downside of this approach is that a many-body method such as CI theory is required to generate the eigenstates and eigenvalues of the zeroth-order Hamiltonian. However, if the partial electron-electron interaction is only a relatively weak long-range interaction, we expect a faster convergence of the eigenstates and eigenvalues with respect to the one-and many-electron CI expansion than for the full Coulomb interaction [52,54], so that a small CI or multi-configuration self-consistent field (MCSCF) description would be sufficiently accurate.
When using a semi-local density-functional approximation for the effective potential of the range-separated adiabatic connection, the presence of an explicit long-range electron-electron interaction in the zeroth-order Hamiltonian has the additional advantage of preventing the collapse of the high-lying Rydberg excitation energies [48,55,56]. In contrast to adiabatic TDDFT, double and multiple excitations can be described with this density-functional perturbation-theory approach, although this possibility was not explored in Refs. [39,45]. Finally, approximate excitedstate wave functions are obtained in the course of the calculations, which is useful for interpretative analysis and for the calculation of properties. We envisage using this densityfunctional perturbation theory to calculate excited states after a range-separated ground-state calculation combining a long-range CI [57,58] or long-range MCSCF [59,60] treatment with a short-range density functional. This would be a simpler alternative to linear-response range-separated MC-SCF theory [61,62] for calculations of excitation energies.
In this work, we study in detail two variants of rangeseparated density-functional perturbation theory based either on the Rayleigh-Schrödinger (RS) or GL perturbation theories and test the first, simpler variant on the He and Be atoms and the H 2 molecule, performing accurate calculations along a range-separated adiabatic connection without introducing density-functional approximations.
The two variants of the range-separated perturbation theory are presented in Section 2. Except for the finite basis approximation, no other approximation is introduced and the computational details can be found in Section 3. Finally, the results obtained for the He and Be atoms and for the H 2 molecule are discussed in Section 4. Section 5 contains our conclusions.

Range-separated ground-state density-functional theory
In range-separated DFT (see, e.g., Ref. [52]), the exact ground-state energy of an N-electron system is obtained by the following minimisation over normalised multideterminantal wave functions where we have introduced the kinetic-energy operatorT , the nuclear attraction operatorV ne = v ne (r)n(r)dr written in terms of the density operatorn(r), a long-range (lr) electron-electron interaction W lr,μ ee = 1 2 w lr,μ ee (r 12 )n 2 (r 1 , r 2 )dr 1 dr 2 , written in terms of the error-function interaction w lr,μ ee (r) = erf (μr)/r and the pair-density operatorn 2 (r 1 , r 2 ) and finally the corresponding complement short-range (sr) Hartree-exchange-correlation density functionalĒ sr,μ Hxc [n ] evaluated at the density of . The density and pair density are obtained as expectation values of the density and pair density operators, respectively, The parameter μ in the error function controls the separation range, with 1/μ acting as a smooth cut-off radius. The Euler-Lagrange equation for the minimisation of Equation (1) leads to the (self-consistent) eigenvalue equation,Ĥ lr,μ μ where μ 0 and E μ 0 are taken as the ground-state wave function and associated energy of the partially interacting Hamiltonian (with an explicit long-range electron-electron interaction)Ĥ lr,μ =T +V ne +Ŵ lr,μ ee +ˆV sr,μ Hxc , which contains the short-range Hartree-exchangecorrelation potential operator,V sr,μ For μ = 0, the HamiltonianĤ lr,μ of Equation (6) reduces to the standard non-interacting KS Hamiltonian, H lr,μ=0 =Ĥ KS , whereas, for μ → ∞, it reduces to the physical HamiltonianĤ lr,μ→∞ =Ĥ . Therefore, when varying the parameter μ between these two limits, the HamiltonianĤ lr,μ defines a range-separated adiabatic connection, linking the non-interacting KS system to the physical system with the ground-state density kept constant (assuming that the exact short-range Hartree-exchangecorrelation potentialv sr,μ Hxc (r) is used).

Excited states from perturbation theory
Excitation energies in range-separated DFT can be obtained by linear-response theory starting from the (adiabatic) timedependent generalisation of Equation (1) [63], where the excited states and their associated energies are obtained from time-independent many-body perturbation theory. In standard KS theory, the single-determinant eigenstates and associated energies of the non-interacting KS Hamiltonian, whereĤ KS =T +V ne +V Hxc , give a first approximation to the eigenstates and associated energies of the physical Hamiltonian. To calculate excitation energies, two variants of perturbation theory using the KS Hamiltonian as zerothorder Hamiltonian have been proposed [39,45]. We here extend these two variants of perturbation theory to rangeseparated DFT. As a first approximation, it is natural to use the excited-state wave functions and energies of the long-range interacting Hamiltonian whereĤ lr,μ is the Hamiltonian of Equation (6) with the short-range Hartree-exchange-correlation potentialˆV sr,μ Hxc evaluated at the ground-state density n 0 . These excitedstate wave functions and energies can then be improved upon by defining perturbation theories in which the Hamil-tonianĤ lr,μ is used as the zeroth-order Hamiltonian.

RS-based variant of perturbation theory
The simplest way of defining such a perturbation theory is to introduce the following Hamiltonian dependent on the coupling constant λ, where the short-range perturbation operator iŝ with the short-range electron-electron interaction W sr,μ ee = (1/2) w sr,μ ee (r 12 )n 2 (r 1 , r 2 )dr 1 dr 2 (12) defined with the complementary error-function interaction w sr,μ ee (r) = erfc(μr)/r. When varying λ, Equation (10) sets up an adiabatic connection linking the long-range interacting Hamiltonian atĤ μ,λ=0 =Ĥ lr,μ , to the physical Hamil-tonianĤ μ,λ=1 =Ĥ , for all μ. Importantly, the ground-state density is not kept constant along this adiabatic connection.
The exact eigenstates and associated eigenvalues of the physical Hamiltonian can be obtained by standard RS perturbation theory -that is by Taylor expanding the eigenstates and eigenvalues of the HamiltonianĤ μ,λ in λ and setting λ = 1: act as zeroth-order eigenstates and energies. Using orthonormalised zerothorder eigenstates μ k | μ l = δ kl and assuming nondegenerate zeroth-order eigenstates, the first-order energy correction for the state k becomes As usual, the zeroth + first-order energy is simply the expectation value of the physical Hamiltonian over the zerothorder eigenstate This expression is a multi-determinantal extension of the exact-exchange KS energy expression for the state k, proposed and studied for the ground state in Refs. [64][65][66].
The second-order energy correction is given by where the first-order wave-function correction is given by For μ = 0, this perturbation theory reduces to the first variant of the KS perturbation theory studied by Filippi et al., see Equation (5) of Ref. [45].
To understand the numerical results in Section 4, we now consider how the zeroth + first-order energies behave with respect to μ near the KS system (μ = 0) and near the physical system (μ → ∞). The total energies up to the first order in perturbation theory are given by the expectation value of the full Hamiltonian over the zeroth-order wave functions in Equation (14). Using the Taylor expansion of the wave function around the KS wave function [53], the zeroth + first-order energies are thus given by where (3) k is the contribution entering at the third power of μ in the zeroth-order wave function.
From the asymptotic expansion of the wave function , which is valid almost everywhere when μ → ∞ (the electron-electron coalescence needs to be treated carefully) [53], the first correction to the zeroth + first-order energies are seen to enter at the fourth power in 1/μ where E (0+1,−4) k is the contribution entering at the fourth power of 1/μ.

GL-based variant of perturbation theory
A second possibility is to define a perturbation theory based on a slightly more complicated adiabatic connection, in which the ground-state density is kept constant between the long-range interacting Hamiltonian and the physical Hamiltonian, see Appendix 1. The Hamiltonian of Equation (10) is then replaced bŷ whereŴ sr,μ is now defined aŝ that depends non-linearly on λ so that the ground-state density n 0 is kept constant for all μ and λ. The density functionals E sr,μ Hx,md [n] and E sr,μ,λ c,md [n] are defined in Appendix 1.
One can show that, for non-degenerate ground-state wave functions μ 0 , the expansion ofV sr,μ,λ c,md in λ for λ → 0 starts at second order V sr,μ,λ c,md = λ 2V sr,μ, (2) c,md Hence, the Hamiltonian of Equation (20) properly reduces to the long-range Hamiltonian at λ = 0,Ĥ μ,λ=0 =Ĥ lr,μ , whereas, at λ = 1, it correctly reduces to the physical Hamiltonian,Ĥ μ,λ=1 =Ĥ . This is so because the short-range Hartree-exchange-correlation potential in the Hamiltonian H lr,μ can be decomposed aŝV sr,μ whereˆV sr,μ c,md =V sr,μ,λ=1 c,md is cancelled by the perturbation terms for λ = 1. Equation (25) corresponds to an alternative decomposition of the short-range Hartree-exchangecorrelation energy into Hartree-exchange and correlation contributions based on the multi-determinantal wave function μ 0 instead of the single-determinant KS wave function KS 0 [64][65][66], which is more natural in range-separated DFT. This decomposition is especially relevant here since it separates the perturbation into a Hartree-exchange contribution that is linear in λ and a correlation contribution containing all the higher order terms in λ.
As before, the first-order energy correction is given by Equation (14) but with the perturbation operator of Equation (21), yielding the following energy up to first order: The second-order energy correction of Equation (16) whereas the expression of the first-order wave function correction is still given by Equation (17) but with the perturbation operator of Equation (21).
For μ = 0, this density-fixed perturbation theory reduces to the second variant of the KS perturbation theory proposed by Görling [39] and studied by Filippi et al. [Equation (6) of Ref. [45]], which is simply the application of GL perturbation theory [40,41] to excited states. In Ref. [45], it was found that the first-order energy corrections in densityfixed KS perturbation theory provided on average a factor of two improvement on the KS zeroth-order excitation energies for the He, Li + and Be atoms when using accurate KS potentials. By contrast, the first-order energy corrections in the first variant of KS perturbation theory, without a fixed density, deteriorated on average the KS excitation energies.
The good results obtained with the second variant of KS perturbation theory may be understood from the fact that, in GL perturbation theory, the ionisation potential remains exact to all orders in λ. In fact, this nice feature of GL theory holds also with range separation, so that the GLbased variant of range-separated perturbation theory should in principle be preferred. However, it requires the separation of the short-range Hartree-exchange-correlation potential into the multi-determinantal Hartree-exchange and multi-determinantal correlation contributions (according to Equation (25)), which we have not done for accurate potentials or calculations along the double adiabatic connection with a partial interaction defined byŴ lr,μ ee + λŴ sr,μ ee (cf. Appendix 1). We, therefore, consider only the RS-based variant of range-separated perturbation theory here but note that the GL-based variant can be straightforwardly applied with density-functional approximations -using, for example, the local-density approximation that has been constructed for the 'multi-determinantal correlation' functional [64,67].

Computational details
Calculations were performed for the He and Be atoms and the H 2 molecule with a development version of the DALTON program [68], see Refs. [69][70][71]. Following the same settings as in Ref. [53], a full CI (FCI) calculation was first carried out to get the exact ground-state density within the basis set considered. Next, a Lieb optimisation of the short-range potential v sr,μ (r) was performed to reproduce the FCI density with the long-range electron-electron interaction w lr,μ ee (r 12 ). Then, an FCI calculation was done with the partially interacting Hamiltonian constructed from w lr,μ ee (r 12 ) and v sr,μ (r) to obtain the zeroth-order energies and wave functions according to Equation (9). Finally, the zeroth + first-order energies were calculated according to Equation (15). The second-order correction of Equation (16) is not calculated in this work. The basis sets used were: uncontracted t-aug-cc-pV5Z for He, uncontracted daug-cc-pVDZ for Be and uncontracted d-aug-cc-pVTZ for H 2 .

Results and discussion
All the zeroth-order curves shown hereinafter correspond to the curves of Ref. [53] as the partially interacting Hamiltonian acts as starting point for the perturbation theory.

Helium atom
The ground-and excited-state total energies to first order along the range-separated adiabatic connection of helium are shown in Figure 1. In the KS limit, when μ = 0, the total energies are significantly improved with respect to the zeroth-order ones. In fact, as shown for the groundstate energy, the zeroth-order total energies are off by approximately 1.2 hartree with respect to the energies of the physical system. When the first-order correction is added, the error becomes smaller than 0.06 hartree for all states. Moreover, the singlet and triplet excited-state energies are no longer degenerate. With increasing range-separation parameter μ, a faster convergence towards the total energies of the physical system is also observed at first order for all states.
The description of the total energies is, therefore, much improved with the addition of the first-order correction. The linear term in μ present in the zeroth-order total energies [53] vanishes for the zeroth + first-order total energies, which instead depend on the third power of μ for small μ (cf. Equation (18)). At large μ, the error relative to the physical energies enters as 1/μ 4 rather than as 1/μ 2 in the zerothorder case, explaining the observed faster convergence of the first-order energies.
The excitation energies of the helium atom correct to zeroth and first orders are plotted in Figure 2. As previously noted, at μ = 0, the degeneracy of the zerothorder singlet and triplet excitation energies is lifted by the first-order correction. However, the excitation energies correct to first order overestimate the physical excitation energies by 0.1-0.2 hartree so that the error is actually larger than at the zeroth order. For the 1 1 S → 1 3 P excitation energy, the correction is even going in the wrong direction and the singlet-triplet splitting is too large by about a factor 1.5.
When the extreme long-range part of the Coulombic interaction is switched on with positive μ close to 0, this initial overestimation is corrected. In fact, for small μ, all excitation energies decrease in the third power of μ, in agreement with Equation (18). When μ 0.5-1, this correction becomes too large and the excitation energies of the partially interacting system become lower than their fully interacting limits. As μ increases further so that more interaction is included, the excitation energies converge toward their fully interacting values from below. The zeroth-order excitation energies, which do not oscillate for small μ, converge monotonically toward their physical limit and are on average more accurate than the zeroth + first order excitation energies. In short, the first-order correction does not improve excitation energies, although total energies are improved.
The inability of the first-order correction to improve excitation energies should be connected to the fact that, since the ground-state density is not kept constant at each order in the perturbation, the ionisation potential is no longer constant to first order along the adiabatic connection. This behaviour results in an unbalanced treatment of the ground and excited states. Moreover, high-energy Rydberg excitation energies should be even more sensitive to this effect, as observed for transitions to the P state. The second GL-based variant of perturbation theory should correct this behaviour Excitation energies in hartree 1 1 S → 1 3 P 0th-order 1 1 S → 1 3 P 0th+1st order 1 1 S → 1 1 P 0th-order 1 1 S → 1 1 P 0th+1st order Excitation energies in hartree 1 1 S → 1 3 P 0th order 1 1 S → 1 3 P 0th+1st order 1 1 S → 1 1 P 0th order 1 1 S → 1 1 P 0th+1st order by keeping the density constant at each order, as shown in the KS case [41,45].

Beryllium atom
When the first-order perturbation correction is applied to the ground-state and valence-excited states of beryllium, total energies are again improved (not illustrated here). In Figure 3, we have plotted the zeroth-and first-order valence excitation energies of beryllium against the rangeseparation parameter μ.
Since valence excitation energies should be less sensitive to a poor description of the ionisation energy than Rydberg excitation energies, the first-order correction should work better for the beryllium valence excitations than for the helium Rydberg excitations. However, although the singlet excitation energy of beryllium is improved at μ = 0 at first order, the corresponding triplet excitation energy is not. In fact, whereas the triplet excitation energy is overestimated at zeroth order, it is underestimated by about the same amount at first order.
As the interaction is switched on, a bump is observed for small μ for the singlet excitation energy but not the triplet excitation energy, which converges monotonically to its physical limit. The convergence of the excitation energies with μ is improved by the first-order excitation energies, especially in the singlet case.

Hydrogen molecule
In Figure 4, we have plotted the excitation energies of H 2 as a function of μ at the equilibrium distance R eq and at 3R eq . At the equilibrium geometry, the first-order correction works well. At μ = 0, the correction is in the right direction (singlet and triplet excitation energies being raised and lowered, respectively); for nearly all μ > 0, the error is smaller than for the zeroth-order excitation energies.
Unfortunately, when the bond is stretched, this is no longer the case. At the stretched geometry, the first excitation energy 1 1 + g → 1 3 + u becomes negative for small values of μ and the error with respect to the physical excitation energy is higher than in the zeroth-order case. Moreover, the ordering of the two singlet excitation energies is incorrect at small μ and they exhibit strong oscillations when   the interaction is switched on. In this case, therefore, the zeroth-order excitation energies are better approximations to the physical excitation energies.

Conclusion
We have considered two variants of perturbation theory along a range-separated adiabatic connection. The first and simpler variant, based on the usual RS perturbation theory, was tested on the helium and beryllium atoms and on the hydrogen molecule at equilibrium and stretched geometries. Although total energies are improved to first order in the perturbation, excitation energies are not improved since the theory does not keep the density constant along the adiabatic connection at each order of perturbation. It would be interesting to examine the evolution of the ionisation potential to understand better the effect of this variant of the perturbation theory on our systems of interest.
The second variant of the perturbation theory, based on GL theory, should improve the results significantly by keeping the ground-state density constant at each order in the perturbation [41], as already observed on the KS system [45]. However, this more complicated theory has not yet been implemented for μ > 0.
An interesting alternative to perturbation theory is provided by extrapolation, which make use of the behaviour of the energies with respect to μ near the physical system to estimate the exact energies from the energy of the partially interacting system at a given μ and its first-order or higher order derivatives with respect to μ [72,73]. Work using this approach will be presented elsewhere.

Disclosure statement
No potential conflict of interest was reported by the authors. We here present a double adiabatic connection, depending on two parameters, that keeps the ground-state density constant. It is the basis for the perturbation theory presented in Section 2.2.2. A different density-fixed double adiabatic connection was considered in Refs. [74,75].

Funding
The Levy-Lieb universal density functional for the Coulomb electron-electron interactionŴ ee is given by [76][77][78]  We here generalise it to the interactionŴ lr,μ ee + λŴ sr,μ ee , wherê W lr,μ ee andŴ sr,μ ee are long-range and short-range electron-electron interactions, respectively, that depend on both a range-separation parameter μ and on a linear parameter λ: