Excitation energies from G{\"o}rling-Levy perturbation theory along the range-separated adiabatic connection

A G{\"o}rling-Levy (GL)-based perturbation theory along the range-separated adiabatic connection is assessed for the calculation of electronic excitation energies. In comparison with the Rayleigh-Schr{\"o}dinger (RS)-based perturbation theory introduced in a previous work [E. Rebolini, J. Toulouse, A. M. Teale, T. Helgaker, A. Savin, Mol. Phys. 113, 1740 (2015)], this GL-based perturbation theory keeps the ground-state density constant at each order and thus gives the correct ionization energy at each order. Excitation energies up to first order in the perturbation have been calculated numerically for the helium and beryllium atoms and the hydrogen molecule without introducing any density-functional approximations. In comparison with the RS-based perturbation theory, the present GL-based perturbation theory gives much more accurate excitation energies for Rydberg states but similar excitation energies for valence states.


Introduction
Within the framework of density-functional theory (DFT), the calculation of molecular excitation energies is nowadays mostly performed using linear-response timedependent density-functional theory (TDDFT) (see, e.g. Refs. [1,2]) within the adiabatic local or semilocal approximations. Despite its many successes, linearresponse TDDFT within these approximations suffers CONTACT Julien Toulouse julien.toulouse@upmc.fr from serious limitations, especially for describing systems with static (or strong) correlation [3], double or multiple excitations [4], and Rydberg and charge-transfer excitations [5,6]. These deficiencies have been attributed to the locality of the approximated exchange-correlation potential and kernel either in space (local and semi-local approximations) or in time (adiabatic approximation). While the former is directly linked to functional development in time-independent DFT, the latter is a problem specific to the time-dependent formulation. However, time dependence is in principle not required to describe excited states since by the Hohenberg-Kohn theorem [7], the time-independent ground-state density contains all the information about the system including information about its excited states. Over the years, several time-independent DFT approaches for calculating excitation energies have emerged and are still being developed: ensemble DFT [8][9][10][11][12][13][14][15][16][17][18], state-specific self-consistent DFT and related methods , hybrid DFT/configuration interaction (CI) methods [41][42][43][44][45] and perturbation theory starting from the non-interacting Kohn-Sham (KS) Hamiltonian [46][47][48][49]. In a previous work [50], we have explored further this density-functional perturbation-theory approach with one key modification: 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 [51][52][53][54][55][56][57]. 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 this previous work [50], a Rayleigh-Schrödinger (RS)-based perturbation theory was tested on a few atoms and molecules and it was found that the first-order excitation energies are not overall improved in comparison with the zeroth-order excitation energies. This finding was rationalised by the fact that this perturbation theory does not keep the ground-state density constant at each order. In the present work, we explore an alternative approach, based on Görling-Levy (GL) perturbation theory along a range-separated adiabatic connection, which keeps the ground-state density constant at each order and is expected to give more accurate excitation energies.
The paper is organised as follows. The main equations of our GL-based range-separated perturbation theory are given in Section 2. The computational details for the calculations carried out, involving no other approximations than the use of a finite basis, can be found in Section 3. The results obtained for the He and Be atoms and for the H 2 molecule are discussed in Section 4. Finally, Section 5 contains our conclusions.

Theoretical background
This section consists of two parts. We first review rangeseparated DFT for ground states in Section 2.1 and then GL-based perturbation theory for excited states in Section 2.2.

Range-separated ground-state density-functional theory
In range-separated DFT (see, e.g. Ref. [55]), the electronelectron interaction is partitioned into long-range and short-range contributions by means of the error function and a range-separation parameter μ which controls the range of the separation. The long-range (lr) interaction operator is defined aŝ and is written in terms of the error-function interaction w lr,μ ee (r) = erf (μr)/r and the pair-density operator n 2 (r 1 , r 2 ), where r refers to the electron coordinates. The exact ground-state energy of an N-electron system is then expressed as where the minimisation is performed over normalised multi-determinantal wave functions . In Equation (2), we have introduced the kinetic-energy opera-torT , the nuclear-electron attraction operatorV ne = v ne (r)n(r)dr written in terms of the density operator n(r) and the complement short-range (sr) Hartreeexchange-correlation density functionalĒ sr,μ Hxc [n ] evaluated at the density of , n (r) = |n(r)| . The minimising wave function μ 0 in Equation (2) is the ground-state wave function of the following eigenvalue problem:Ĥ associated with the energy E μ 0 . In Equation (3), we have introduced the partially interacting Hamiltonian H lr,μ =T +V ne +Ŵ lr,μ ee +V sr,μ Hxc , which contains the operator V sr,μ Hxc = v sr,μ Hxc (r)n(r)dr, with the short-range Hartree-exchange-correlation potentialv sr,μ Hxc (r) = δĒ sr,μ Hxc [n 0 ]/δn(r) evaluated at the exact ground-state density n 0 (r). This is the potential that keeps the ground-state density constant for all μ, i.e. μ 0 |n(r)| μ 0 = n 0 (r). In this paper, contrary to what was sometimes done in previous papers on rangeseparated DFT, we use an overline in the notation for the short-range Hartree-exchange-correlation energȳ E sr,μ Hxc [n] and its associated potential operatorV sr,μ Hxc to emphasise that these quantities are defined as complements to their long-range analogues, i.e. they include a mixed long-range/short-range correlation contribution [58,59]. At μ = 0, the HamiltonianĤ lr,μ in Equation (4) reduces to the standard non-interacting KS Hamiltonian, whereas for μ → Ý it reduces to the physical Hamiltonian. The HamiltonianĤ lr,μ therefore defines a rangeseparated adiabatic connection, linking the KS and the physical systems when varying μ, keeping the groundstate density constant.

Excited states from GL-based perturbation theory
For a given value of μ, the excited-state wave functions and energies of the long-range interacting Hamiltonian can be used as zeroth-order approximations to the physical excited-state wave functions and energies. They can then be improved upon by setting up perturbation theories withĤ lr,μ as the zeroth-order Hamiltonian. As shown in Ref. [50], an RS-based perturbation theory in which the ground-state density is not kept constant gives first-order excitation energies that overall do not improve upon the zeroth-order excitation energies. Here, we instead use a GL-based perturbation theory, with the ground-state density kept constant.
To formulate such a GL-based perturbation theory, we define the following Hamiltonian depending on a coupling constant λ: which contains the operator wherev sr,μ,λ Hxc (r) is the λ-dependent short-range Hartreeexchange-correlation potential that keeps the groundstate density constant for all μ and all λ -that is, is the groundstate wave function of the HamiltonianĤ μ,λ in Equation (7). The HamiltonianĤ μ,λ thus sets up a double adiabatic connection with a constant ground-state density. To clearly separate the linear and non-linear dependence in λ, we then rewrite Equation (7) and we define the potential operator which we choose to denote without an overline because it is a 'double complement' in the sense that it is the complement to the complement potentialV sr,μ,λ Hxc with respect to the potentialV sr,μ Hxc . This potential can be decomposed into a linear contribution with respect to λ and a term containing all higher-order terms (see Ref. [50]) where the potentialV sr,μ Hx,md is the short-range 'multideterminantal (md) Hartree-exchange' potential operator, whileV sr,μ,λ c,md is the λ-dependent short-range 'multi-determinantal correlation' potential operator (see Ref. [60]). For non-degenerate ground states of the long-range interacting Hamiltonian in Equation (4), the expansion ofV sr,μ,λ c,md in λ around λ = 0 begins at second-order: The partially interacting Hamiltonian can then be rewritten aŝ and the corresponding zeroth+first-order energy is (17) whereV sr,μ c,md =V sr,μ Hxc −V sr,μ Hx,md =V sr,μ,λ=1 c,md is the short-range multi-determinantal correlation potential operator. This operator can be expressed asV sr,μ c,md = v sr,μ c,md (r)n(r)dr wherev sr,μ c,md (r) = δĒ sr,μ c,md [n 0 ]/δn(r) is the functional derivative of the short-range multideterminantal correlation energy functionalĒ sr,μ c,md [n] introduced in Refs. [58,60] (see also Refs. [61][62][63]), evaluated at the exact ground-state density n 0 . Equation (17) contains two contributions; the zeroth+first-order energy in standard RS perturbation theory, and an additional term, which is only present in GL-based perturbation theory. For μ = 0, the zeroth-order wave functions are KS single determinants, μ=0 k = k , and the short-range multi-determinantal correlation potential reduces to the standard KS potential,v sr,μ=0 c,md (r) = v c (r), and we thus recover the standard GL perturbation theory for excitation energies [48,49].
A local-density approximation has been developed for the short-range multi-determinant correlation energy functionalĒ sr,μ c,md [n] (and thus for its associated potential v sr,μ c,md (r)) [60,64]. In the present study, we will test the GL-based first-order perturbation theory without introducing any density-functional approximations, providing benchmark data for future approximations.

Computational details
Calculations were performed for the He and Be atoms and the H 2 molecule with a development version of the DAL-TON program [65,66]; see Refs. [56,67,68]. Following the same procedure as in Ref. [57], a full CI (FCI) calculation was first carried out to obtain the exact ground-state density in the chosen basis set. Next, for several values of μ and λ, a Lieb optimisation [67,69,70] was carried out to obtain the short-range potentialv sr,μ,λ (r) = v ne (r) + v sr,μ,λ Hxc (r) needed in Equation (7) to reproduce this FCI ground-state density with the partial electron-electron interaction w lr,μ ee (r 12 ) + λw sr,μ ee (r 12 ). Then, an FCI calculation was performed with the partially interacting Hamiltonian constructed from w lr,μ ee (r 12 ) + λw sr,μ ee (r 12 ) andv sr,μ,λ (r) to obtain, for a few states, the energies E μ,λ k and wave functions μ,λ k , for each values of μ and λ. For each system, each excited state and each value of μ, a third-degree polynomial fit in λ was performed on the excitation energy E μ,λ 0 for values of λ ranging from 0 to 0.3 and the first-order correction in the GL-based perturbation theory was obtained as the firstorder derivative with respect to λ. The zeroth-order excitation energies and the first-order excitation energies in the RS-based perturbation theory were already available from Refs. [57] and [50], respectively. Because of issues with loss of numerical precision, we refrain from extracting second-and higher-order derivatives with respect to λ, which could be used to test higher-order perturbation theories or extrapolation schemes [71].
The basis sets used were uncontracted t-aug-cc-pV5Z for He, uncontracted d-aug-cc-pVDZ for Be and uncontracted d-aug-cc-pVTZ for H 2 .

First Rydberg excitation energies of the helium atom
The first singlet and triplet excitation energies of the helium atom correct to zeroth and first-orders along the range-separated adiabatic connection in the GL-and RSbased perturbation theories are shown in Figure 1.
In the KS limit, at μ = 0, the zeroth-order singlet and triplet excitation energies are degenerate. When increasing μ, this degeneracy is lifted and the zeroth-order singlet and triplet excitation energies eventually converge to their physical excitation energies for μ → Ý. With the introduction of the perturbation, the singlet/triplet degeneracy is lifted already at μ = 0. As found in Ref. [50], the RS-based first-order perturbative correction globally deteriorates the zeroth-order excitation energies, leading to large errors that do not decrease monotonically with μ.
Interestingly, the zeroth+first-order excitation energies obtained from the GL-based perturbation theory have smaller errors than those obtained by the RS-based perturbation theory. At μ = 0, the contributions to the excitation energies in the GL-based first-order perturbation theory coming from the additional term in Equation (19) are equal to v c (r)[n 2s (r) − n 1s (r)]dr for the 1 1 S → 2 1 S and 1 1 S → 2 3 S transitions and v c (r)[n 2p (r) − n 1s (r)]dr for the 1 1 S → 1 1 P and 1 1 S → 1 3 P transitions, where n k (r) is the density of the KS orbital k. At the scale of the plot, these corrections to the RS-based first-order perturbative excitation energies are significant (more than −0.01 eV) and strongly reduce the errors. As μ increases, the corrections coming from the additional term in Equation (19) change sign but are still efficient to reduce the errors. For sufficiently large μ, in comparison to the zeroth-order excitation energies and to the zeroth+first-order RS excitation energies, the zeroth+first-order GL excitation energies systematically converge faster with respect to the rangeseparation parameter to the physical excitation energies: a 1 millihartree accuracy is reached for μ larger than 1.2 bohr −1 , while values of 3-4 bohr −1 are necessary for the zeroth-order curves. However, near μ = 0, the first-order GL correction does not always improve the zeroth-order excitation energy (see the 1 1 S → 1 1 P transition).
These comparisons show that the GL-based perturbation theory (which keeps the ground-state density constant, ensuring a correct ionisation energy) is a much better strategy for calculating Rydberg excitation energies than the RS-based perturbation theory. This was expected since these Rydberg excitation energies are close to the ionisation threshold.

Valence excitation energies of the beryllium atom
The valence excitation energies of the beryllium atom correct to zeroth and first-orders in the GL-based and RS-based perturbation theories are plotted in Figure 2.
For these excitations, the two approaches give very similar first-order corrections. This behaviour can be rationalised from the fact that these valence excited states are far from the ionisation threshold so that imposing the correct ionisation energy (by keeping the ground-state density constant) has much less impact than for the Rydberg excitation energies. This can also be understood by looking at the expression of the difference between the zeroth+first-order GL and RS excitation energies coming from the additional term in Equation (19) which is, for μ = 0, equal to v c (r)[n 2p (r) − n 2s (r)]dr for the 1 1 S → 1 1 P and 1 1 S → 1 3 P transitions. This quantity is necessarily small since the 2s and 2p orbitals are localised in the same region of space.
We note that the first-order correction systematically improves the singlet excitation energy but not the triplet excitation energy, which is overestimated at zerothorder and underestimated by about the same amount at zeroth+first-order for small μ.

First valence excitation energies of the hydrogen molecule
In Figure 3, we have plotted the first valence excitation energies of the hydrogen molecule at its equilibrium internuclear distance R eq and at 3R eq .
As for the valence excitation energies in beryllium, the first-order RS-and GL-based corrections are overall quite similar, but with a discernible faster μ-convergence of the zeroth+first-order GL excitation energies to the physical excitation energies at the equilibrium distance. Again, the fact that the difference between the zeroth+first-order GL and RS excitation energies is small can be understood from its expression at μ = 0 which is v c (r)[n 1σ u (r) − n 1σ g (r)]dr for the 1 1 + g → 1 1 + u and 1 1 + g → 1 3 + u transitions, and which involves the difference between two similar orbital densities. At the equilibrium distance, the GL-based first-order perturbation theory overshoots the correction to both the singlet and triplet zeroth-order excitation energies for small values of μ, but nevertheless improves upon the zeroth-order correction for all μ values.
When the bond is stretched, the first-order correction no longer systematically improves on the zerothorder excitation energies for small μ. The zeroth+firstorder excitation energy of the first transition 1 1 + g → 1 3 + u becomes negative for small μ and the error with respect to the physical excitation energy is higher than in the zeroth-order case. The zeroth+first-order excitation energies for the two singlet states, 1 1 + u and 2 1 + g , are incorrectly ordered at small μ and do not monotonically converge to the physical energies, passing through a maximum at around μ 0.5 bohr −1 . At μ = 0, in comparison to the RS-based first-order perturbation theory, the GL-based first-order perturbation theory gives a better estimate of the excitation energy for the 2 1 + g state of double-excitation character, the additional term in the first-order GL-based correction being 2 v c (r)[n 1σ u (r) − n 1σ g (r)]dr. However, in the dissociation limit, it is easy to show that both RS-based and GLbased first-order perturbation theories will incorrectly lead to a vanishing excitation energy at μ = 0 for this state. Only for values of μ greater than about 1 bohr −1 , both RS-and GL-based first-order perturbation theories provide accurate estimates of the excitation energies for all the states considered.

Conclusion
We have applied a GL-based perturbation theory along a range-separated adiabatic connection for the calculation of electronic excitation energies. Unlike the RS-based perturbation theory that we explored in a previous work, the GL-based perturbation theory keeps the ground-state density (and thus the ionisation energy) constant at each order. Excitation energies up to first-order in the perturbation have been calculated numerically for the helium and beryllium atoms and the hydrogen molecule without introducing any density-functional approximations.
In comparison with the RS-based perturbation theory, the GL-based perturbation theory gives much more accurate excitation energies for the Rydberg states of the helium atom but similar excitation energies for the valence states of the beryllium atom and of the hydrogen molecule. This can be rationalised by observing that the Rydberg states are close to the ionisation threshold and therefore sensitive to having the correct ionisation energy.
This first-order GL-based perturbation theory works reasonably well for calculating the first valence excitation energies of the hydrogen molecule at its equilibrium distance. However, results are less satisfactory for the valence excitation energies of the beryllium atom and stretched hydrogen molecule at small range-separation parameter. For such systems, with small gaps, it may be necessary to go beyond single-reference first-order perturbation theory for small range-separation parameters.
One possible extension of this work would be to test density-functional approximations for short-range exchange-correlation potential and approximations for the wave-function part of the calculation. In particular, the present approach could be applied in practice by first performing a range-separated DFT groundstate calculation using a long-range multiconfiguration self-consistent-field (MCSCF) wave function and a short-range (semi-)local density-functional approximation in Equation (2), as developed in Refs. [72][73][74][75], then doing a CI calculation for the excited states in Equation (6) and evaluating the excitation energies via Equation (17) using a density-functional approximation forv sr,μ c,md (r). We would then obtain a timeindependent range-separated DFT method for calculating excitation energies, as an alternative to the more usual linear-response time-dependent range-separated DFT approaches [76][77][78][79]. Finally, for nearly degenerate systems, it would also be interesting to explore the extension of this approach to range-separated ensemble DFT [14][15][16]18] which would involve replacing the initial ground-state MCSCF calculation by a state-average MCSCF.