Linear-response range-separated density-functional theory for atomic photoexcitation and photoionization spectra
LLinear-response range-separated density-functional theoryfor atomic photoexcitation and photoionization spectra
Felipe Zapata, ∗ Eleonora Luppi, † and Julien Toulouse ‡ Laboratoire de Chimie Th´eorique (LCT), Sorbonne Universit´e and CNRS, F-75005 Paris, France (Dated: March 14, 2019)We investigate the performance of the range-separated hybrid (RSH) scheme, which combineslong-range Hartree-Fock (HF) and a short-range density-functional approximation (DFA), for cal-culating photoexcitation/photoionization spectra of the H and He atoms, using a B-spline basisset in order to correctly describe the continuum part of the spectra. The study of these simplesystems allows us to quantify the influence on the spectra of the errors coming from the short-rangeexchange-correlation DFA and from the missing long-range correlation in the RSH scheme. We studythe differences between using the long-range HF exchange (nonlocal) potential and the long-rangeexact exchange (local) potential. Contrary to the former, the latter supports a series of Rydbergstates and gives reasonable photoexcitation/photoionization spectra, even without applying linear-response theory. The most accurate spectra are obtained with the linear-response time-dependentrange-separated hybrid (TDRSH) scheme. In particular, for the He atom at the optimal value of therange-separation parameter, TDRSH gives slightly more accurate photoexcitation and photoioniza-tion spectra than standard linear-response time-dependent HF. More generally, the present workshows the potential of range-separated density-functional theory for calculating linear and nonlinearoptical properties involving continuum states.
I. INTRODUCTION
Nowadays, time-dependent density-functional theory(TDDFT) [1], applied within the linear-response formal-ism [2–4], is a widely used approach for calculating pho-toexcitation spectra (transitions from bound to boundstates) of electronic systems. In spite of many successes,it is however well known that usual (semi-)local density-functional approximations (DFAs), i.e. the local-densityapproximation (LDA) and generalized-gradient approx-imations (GGAs), for the exchange-correlation poten-tial and its associated exchange-correlation kernel do notcorrectly describe long-range electronic transitions, suchas those to Rydberg [5] and charge-transfer [6] statesin atomic and molecular systems. A better descriptionof Rydberg excitations can be obtained with exchange-correlation potential approximations having the correct − /r long-range asymptotic decay [7–10], even thoughit has been shown that accurate Rydberg excitation en-ergies and oscillator strengths can in fact be extractedfrom LDA calculations in small atoms [11, 12]. A moregeneral solution for correcting both Rydberg and charge-transfer excitations is given by range-separated TDDFTapproaches [13–19] which express the long-range partof the exchange potential and kernel at the Hartree-Fock (HF) level. These range-separated approaches alsogive reasonably accurate values for the ionization energythreshold [14, 20, 21].Linear-response TDDFT has also been used for calcu-lating photoionization spectra (transitions from boundto continuum states) of atoms and molecules [22–33]. ∗ [email protected] † [email protected] ‡ [email protected] These calculations are less standard in quantum chem-istry since they involve spatial grid methods or B-splinebasis sets for a proper description of the continuumstates. In this case as well, usual (semi-)local DFAsprovide a limited accuracy and asymptotically correctedexchange-correlation potential approximations give moresatisfactory results. More accurate still, but less com-mon, are photoionization spectra calculated with theexact-exchange (EXX) potential [28] or the localized HFexchange potential and its associated kernel [33]. Re-cently, range-separated approximations have been suc-cessfully used for calculating photoexcitation and pho-toionization spectra of molecular systems using time-propagation TDDFT with Gaussian basis sets togetherwith an effective lifetime model compensating for themissing continuum states [34–36]. However, to the bestof our knowledge, range-separated approximations havenot yet been used in frequency-domain linear-responseTDDFT calculations of photoionization spectra.In this work, we explore the performance of thelinear-response time-dependent range-separated hybrid(TDRSH) scheme [19, 37] for calculating photoexcitationand photoionization spectra of the H and He atoms us-ing a B-spline basis set to accurately describe the con-tinuum part of the spectra. The TDRSH scheme al-lows us to treat long-range exchange effects at the HFlevel and short-range exchange-correlation effects within(semi-)local DFAs. First, the dependence of the range-separated hybrid (RSH) orbital energies on the range-separation parameter is investigated, as well as the ef-fect of replacing the long-range HF exchange nonlocalpotential by the long-range EXX local potential (result-ing in a scheme that we refer to as RSH-EXX). Sec-ond, oscillator strengths directly computed with the RSHand the RSH-EXX orbitals are compared with oscilla-tor strengths obtained with the linear-response TDRSH a r X i v : . [ phy s i c s . c h e m - ph ] M a r scheme. The study of the H atom allows us to quantifythe residual self-interaction error coming from the short-range exchange-correlation DFA, and the study of the Heatom permits to quantify the effect of the missing long-range correlation in the RSH scheme. This work consti-tutes a first step for applying range-separated TDDFTto strong-field phenomena, such as high-harmonic gen-eration or above-threshold ionization, where long-rangeeffects and continuum states play an important role.The outline of the paper is as follows. In Sec. II, firstly,we briefly review the RSH scheme and introduce theRSH-EXX variant, and, secondly, we review the linear-response TDRSH method. In Sec. III, the basis set ofB-spline functions is defined, and we indicate how therange-separated two-electron integrals are computed us-ing an exact spherical harmonic expansion for the range-separated interaction. In Sec. IV results are presentedand discussed. Firstly, we show the performance of theB-spline basis set for describing the density of contin-uum states of the H atom within the different meth-ods. Secondly, the dependence of the orbital energies ofthe H and He atoms on the range-separation parameteris analyzed. Thirdly, different calculated photoexcita-tion/photoionization spectra for the H and He atoms arediscussed and compared with exact results. In Sec. V,conclusions and perspectives are given. Unless otherwiseindicated, Hartree atomic units are used throughout thepaper. II. RANGE-SEPARATEDDENSITY-FUNCTIONAL THEORYA. Range-separated hybrid scheme
Range-separated density-functional theory (see, e.g.,Refs. 38 and 39) is based on the splitting of the Coulombelectron-electron interaction w ee ( r ) = 1 /r into long-range (lr) and short-range (sr) contributions w ee ( r ) = w lree ( r ) + w sree ( r ) , (1)and the most common forms for the long-range and short-range interactions are w lree ( r ) = erf( µr ) r , (2)and w sree ( r ) = erfc( µr ) r . (3)where erf and erfc are the error function and the comple-mentary error function, respectively, and µ is a tunablerange-separation parameter controlling the range of theseparation. Using this decomposition, it is possible torigorously combine a long-range wave-function approachwith a complementary short-range DFA. The simplest approach in range-separated density-functional theory consists in using a single-determinantwave function for the long-range interaction. This leadsto the RSH scheme [40] which spin orbitals { ϕ p ( x ) } (where x = ( r , σ ) are space-spin coordinates) and or-bital energies ε p can be determined for a given system bythe following eigenvalue problem, (cid:18) − ∇ + v ne ( r ) + v H ( r ) + v srxc ( x ) (cid:19) ϕ p ( x )+ (cid:90) v lr , HFx ( x , x (cid:48) ) ϕ p ( x (cid:48) )d x (cid:48) = ε p ϕ p ( x ) , (4)where v ne ( r ) is the nuclei-electron potential, v H ( r ) is theHartree potential for the Coulomb electron-electron in-teraction, v H ( r ) = (cid:90) n ( x (cid:48) ) w ee ( | r − r (cid:48) | )d x (cid:48) , (5)where n ( x ) = (cid:80) occ i | ϕ i ( x ) | are the spin densities ( i refersto occupied spin orbitals), v lr , HFx ( x , x (cid:48) ) is the nonlocal HFexchange potential for the long-range electron-electroninteraction, v lr , HFx ( x , x (cid:48) ) = − occ (cid:88) i ϕ ∗ i ( x (cid:48) ) ϕ i ( x ) w lree ( | r − r (cid:48) | ) , (6)and v srxc ( x ) is the short-range exchange-correlation poten-tial v srxc ( x ) = δ ¯ E srxc δn ( x ) , (7)where ¯ E srxc is the complement short-range exchange-correlation density functional. In this work, we use theshort-range spin-dependent LDA exchange-correlationfunctional of Ref. 41 for ¯ E srxc . The long-range and short-range potentials, v lr , HFx ( x , x (cid:48) ) and v srxc ( x ), explicitly de-pend on the range-separation parameter µ , and con-sequently the spin orbitals, the orbital energies, andthe density also implicitly depend on it. For µ = 0, v lr , HFx ( x , x (cid:48) ) vanishes and v srxc ( x ) becomes the usual full-range LDA exchange-correlation potential, and thus theRSH scheme reduces to standard Kohn-Sham LDA. For µ → ∞ , v lr , HFx ( x , x (cid:48) ) becomes the usual full-range HF ex-change potential and v srxc ( x ) vanishes, and thus the RSHscheme reduces to standard HF.In the present paper, we also consider the followingvariant of the RSH scheme, (cid:18) − ∇ + v ne ( r ) + v H ( r ) + v srxc ( x ) + v lr , EXXx ( x ) (cid:19) ϕ p ( x )= ε p ϕ p ( x ) , (8)in which the long-range nonlocal HF exchange potentialhas been replaced by the long-range local EXX [42–44]potential v lr , EXXx ( x ) = δE lrx δn ( x ) , (9)where E lrx is the long-range exchange density func-tional [45, 46]. We will refer to this scheme as RSH-EXX.The calculation of the EXX potential is involved [47–49],with the exception of one- and two-electron systems. In-deed, for one-electron systems, the long-range EXX po-tential is simply v lr , EXXx ( x ) = − v lrH ( r ) , (10)and for systems of two electrons in a single spatial orbital,it is v lr , EXXx ( x ) = − v lrH ( r ) , (11)where v lrH ( r ) = (cid:82) n ( x (cid:48) ) w lree ( | r − r (cid:48) | )d x (cid:48) is the long-rangeHartree potential. For these one- and two-electron cases,it can be shown that Eqs. (4) and (8) give identical oc-cupied orbitals but different unoccupied orbitals. Moregenerally, for systems with more than two electrons, theHF and EXX exchange potentials give similar occupiedorbitals but very different unoccupied orbitals.Once orbitals and orbital energies are obtained fromEqs. (4) and (8), the bare oscillator strengths can becalculated. They are defined as f ia = 23 ω ia (cid:88) ν = x,y,z | d ν,ia | , (12)where i and a refer to occupied and unoccupied spin or-bitals, respectively, ω ia = ε a − ε i are the bare excitationenergies and d ν,ia = (cid:82) ϕ ∗ i ( x ) r ν ϕ a ( x )d x are the dipole-moment transition integrals. We will consider these bareexcitation energies ω ia and oscillator strengths f ia fora first approximation to photoexcitation/photoionizationspectra. B. Linear-response time-dependentrange-separated hybrid scheme
In the time-dependent extension of the RSH schemewithin linear response (referred to as TDRSH) [18, 19,37], one has to solve the following pseudo-Hermitianeigenvalue equation (cid:18)
A B − B ∗ − A ∗ (cid:19) (cid:18) X n Y n (cid:19) = ω n (cid:18) X n Y n (cid:19) , (13)whose solutions come in pairs: excitation energies ω n > X n , Y n ), and de-excitation energies ω n < Y ∗ n , X ∗ n ). The elements of thematrices A and B are A ia,jb = ( ε a − ε i ) δ ij δ ab + K ia,jb , (14) B ia,jb = K ia,bj , (15)where i, j and a, b refer to occupied and unoccupiedRSH spin orbitals, respectively, and the coupling ma-trix K contains the contributions from the Hartree kernel f H ( r , r ) = w ee ( | r − r | ), the long-range HF exchangekernel f lr , HFx ( x , x ; x (cid:48) , x (cid:48) ) = − w lree ( | r − r | ) δ ( x − x (cid:48) ) δ ( x (cid:48) − x ), and the adiabatic short-range exchange-correlation kernel f srxc ( x , x ) = δv srxc ( x ) /δn ( x ) K ia,jb = (cid:104) aj | f H | ib (cid:105) + (cid:104) aj | f lr , HFx | ib (cid:105) + (cid:104) aj | f srxc | ib (cid:105) = (cid:104) aj | w ee | ib (cid:105) − (cid:104) aj | w lree | bi (cid:105) + (cid:104) aj | f srxc | ib (cid:105) , (16)where (cid:104) aj | w ee | ib (cid:105) and (cid:104) aj | w lree | bi (cid:105) are the two-electronintegrals associated with the Coulomb and long-range interactions, respectively, and (cid:104) aj | f srxc | ib (cid:105) = (cid:82)(cid:82) ϕ ∗ a ( x ) ϕ ∗ j ( x ) f srxc ( x , x ) ϕ i ( x ) ϕ b ( x )d x d x . Sincewe use the short-range LDA exchange-correlation den-sity functional, for µ = 0 the TDRSH scheme re-duces to the usual linear-response time-dependent local-density approximation (TDLDA). For µ → ∞ , theTDRSH scheme reduces to standard linear-responsetime-dependent Hartree-Fock (TDHF).The time-dependent extension of the RSH-EXXvariant within linear response (referred to asTDRSH-EXX) leads to identical equations withthe exception that the long-range HF exchangekernel f lr , HFx ( x , x ; x (cid:48) , x (cid:48) ) is replaced by the long-range frequency-dependent EXX kernel [50, 51] f lr , EXXx ( x , x ; ω ) = δv lr , EXXx ( x , ω ) /δn ( x , ω ). Forone-electron systems, the long-range EXX kernel issimply f lr , EXXx ( x , x ; ω ) = − f lrH ( r , r ) , (17)and, for systems with two electrons in a single spatialorbital, it is f lr , EXXx ( x , x ; ω ) = − f lrH ( r , r ) , (18)where f lrH ( r , r ) = w lree ( | r − r | ) is the long-rangeHartree kernel. For these one- and two-electron cases,TDRSH and TDRSH-EXX give rise to identical excita-tion energies and oscillator strengths.Finally, we can calculate the corresponding TDRSH(or TDRSH-EXX) oscillator strengths as f n = 23 ω n (cid:88) ν = x,y,z | d ν,ia ( X n,ia + Y n,ia ) | . (19)In the limit of a complete basis set, the linear-responseoscillator strengths in Eq. (19) always fulfill the Thomas-Reiche-Kuhn (TRK) sum rule, (cid:80) n f n = N where N isthe electron number. The bare oscillator strengths ofEq. (12) fulfill the TRK sum rule only in the case wherethe orbitals have been obtained from an effective localpotential, i.e. for LDA and RSH-EXX but not for HFand RSH (see Ref. 37). III. IMPLEMENTATION INA B-SPLINE BASIS SET
In practice, each spin orbital is decomposed into aproduct of a spatial orbital and a spin function, ϕ p ( x ) = ϕ p ( r ) δ σ p ,σ where σ p is the spin of the spin orbital p , andwe use spin-adapted equations. As we investigate atomicsystems, the spatial orbitals are written in spherical co-ordinates, ϕ p ( r ) = R n p l p ( r ) Y m p l p (Ω) , (20)where Y m p l p (Ω) are the spherical harmonics (Ω stands forthe angles θ, φ ) and the radial functions R n p l p ( r ) are ex-pressed as linear combinations of B-spline functions oforder k s , R n p l p ( r ) = N s (cid:88) α =1 c n p l p α B k s α ( r ) r , (21)where N s is the dimension of the basis. To completelydefine a basis of B-spline functions, a non-decreasing se-quence of N s + k s knot points (some knot points are pos-sibly coincident) must be given [52]. The B-spline func-tion B k s α ( r ) is non zero only on the supporting interval[ r α , r α + k s ] (containing k s + 1 consecutive knot points)and is a piecewise function composed of polynomialsof degree k s − k s − m deriva-tives across each knot of multiplicity m . We have cho-sen the first and the last knots to be k s -fold degener-ate, i.e. r = r = · · · = r k s = R min and r N s +1 = r N s +2 = · · · = r N s + k s = R max , while the multiplicity ofthe other knots is unity. The spatial grid spacing waschosen to be constant in the whole radial space betweentwo consecutive non-coincident points and is given by∆ r = R max / ( N s − k s + 1). In the present work, thefirst and the last B-spline functions were removed fromthe calculation to ensure zero boundary conditions at r = R min and r = R max . The results presented in thispaper have been obtained using the following parame-ters: k s = 8, N s = 200, R min = 0, and R max = 100bohr. Moreover, we need to use only s and p z sphericalharmonics.Working with such a B-spline representation, one mustcompute matrix elements involving integrals over B-spline functions. The principle of the calculation of one-electron and two-electron integrals over B-spline func-tions are well described by Bachau et al. in Ref. 53. Wewill now briefly review the computation of the standardCoulomb two-electron integrals over B-spline functions,and then we will present the calculation of the long-rangeor short-range two-electron integrals over B-spline func-tions, the latter being original to the present work. A. Coulomb two-electron integrals
The Coulomb electron-electron interaction is given by w ee ( | r − r (cid:48) | ) = 1( | r | + | r (cid:48) | − | r || r (cid:48) | cos γ ) / , (22)where r and r (cid:48) are electron vector positions and γ is theangle between them. The multipolar expansion for this interaction is w ee ( | r − r (cid:48) | ) = ∞ (cid:88) k =0 (cid:20) r k< r k +1 > (cid:21) k (cid:88) m k = − k ( − m k C k − m k (Ω) C km k (Ω (cid:48) ) , (23)where r < = min( | r | , | r (cid:48) | ) and r > = max( | r | , | r (cid:48) | ) and C km k (Ω) = (4 π/ (2 k + 1)) / Y m k k (Ω) are the renormal-ized spherical harmonics. The Coulomb two-electron in-tegrals, in the spatial orbital basis, can then be expressedas the sum of products of radial integrals and angular fac-tors (cid:104) pq | w ee | tu (cid:105) = ∞ (cid:88) k =0 R k ( p, q ; t, u ) k (cid:88) m k = − k δ m k ,m p − m t δ m k ,m q − m u × ( − m k c k ( l p , m p , l t , m t ) c k ( l q , m q , l u , m u ) , (24)where R k ( p, q ; t, u ) are the two-dimensional radial Slaterintegrals and the angular coefficients c k ( l p , m p , l t , m t )and c k ( l q , m q , l u , m u ) are obtained from the Gaunt co-efficients [54, 55]. The coefficient c k ( l, m, l (cid:48) , m (cid:48) ) is nonzero only if | l − l (cid:48) | ≤ k ≤ l + l (cid:48) and if l + l (cid:48) + k is an eveninteger, which makes the sum over k in Eq. (24) exactlyterminate. The Slater integrals are defined as R k ( p, q ; t, u ) = N s (cid:88) α =1 N s (cid:88) λ =1 N s (cid:88) β =1 N s (cid:88) ν =1 c n p l p α c n q l q λ c n t l t β c n u l u ν × R k ( α, λ ; β, ν ) , (25)where R k ( α, λ ; β, ν ) are the Slater matrix elements givenby R k ( α, λ ; β, ν ) = (cid:90) ∞ (cid:90) ∞ B k s α ( r ) B k s λ ( r (cid:48) ) (cid:20) r k< r k +1 > (cid:21) × B k s β ( r ) B k s ν ( r (cid:48) )d r d r (cid:48) . (26)In order to compute the Slater matrix elements R k ( α, λ ; β, ν ), we have implemented the integration-cellalgorithm developed by Qiu and Froese Fischer [56]. Thisalgorithm exploits all possible symmetries and B-splineproperties to evaluate efficiently the integrals in each two-dimensional radial region on which the integrals are de-fined. Gaussian quadrature is used to compute the inte-grals in each cell. B. Long-range and short-range two-electronintegrals
A closed form of the multipolar expansion of the short-range electron-electron interaction defined in Eq. (3) wasdetermined by ´Angy´an et al. [57], following a previouswork of Marshall [58] who applied the Gegenbauer addi-tion theorem to the Laplace transform of Eq. (3). Thisexact expansion is w sree ( | r − r (cid:48) | ) = ∞ (cid:88) k =0 S k ( r > , r < ; µ ) × k (cid:88) m k = − k ( − m k C k − m k (Ω) C km k (Ω (cid:48) ) , (27)where the µ -dependent radial function is written in termsof the scaled radial coordinates Ξ = µ r > and ξ = µ r < as S k ( r > , r < ; µ ) = µ Φ k (Ξ , ξ ) , (28)with Φ k (Ξ , ξ ) = H k (Ξ , ξ ) + F k (Ξ , ξ )+ k (cid:88) m =1 F k − m (Ξ , ξ ) Ξ m + ξ m ( ξ Ξ) m , (29)and the introduced auxiliary functions H k (Ξ , ξ ) = 12( ξ Ξ) k +1 (cid:2)(cid:0) Ξ k +1 + ξ k +1 (cid:1) erfc(Ξ + ξ ) − (cid:0) Ξ k +1 − ξ k +1 (cid:1) erfc(Ξ − ξ ) (cid:3) , (30)and F k (Ξ , ξ ) = 2 π / k (cid:88) p =0 (cid:18) − ξ Ξ) (cid:19) p +1 ( k + p )! p !( k − p )! × (cid:104) ( − k − p e − (Ξ+ ξ ) − e − (Ξ − ξ ) (cid:105) . (31)In order to arrive at a separable expression in Ξ and ξ , ´Angy´an et al. [57] also introduced a power series ex-pansion of the radial function Φ k (Ξ , ξ ) in the smallerreduced variable ξ . However, the range of validity ofthis expansion truncated to the first few terms is lim-ited to small values of ξ , i.e. ξ (cid:46) .
5, and higher-orderexpansions show spurious oscillations. After some tests,we decided to use the exact short-range radial functionΦ k (Ξ , ξ ) without expansion in our work.The expression of the short-range two-electron inte-grals (cid:104) pq | w sree | tu (cid:105) is then identical to the one in Eq. (24)with the simple difference that the radial term is notgiven by the standard Slater matrix elements. Now, theradial kernel in Eq. (26) is changed to that of Eq. (28).Due to the fact that the radial kernel is not multi-plicatively separable in the variables r > and r < , theintegration-cell algorithm is modified in order to calculateall integrals as non-separable two-dimensional integrals.In a second step, the long-range two-electron integralscan be simply obtained by difference (cid:104) pq | w lree | tu (cid:105) = (cid:104) pq | w ee | tu (cid:105) − (cid:104) pq | w sree | tu (cid:105) . (32) IV. RESULTS AND DISCUSSION
In this section, photoexcitation and photoionizationspectra for the H and He atoms are presented. Pho-toexcitation and photoionization processes imply transi-tions from bound to bound and from bound to continuumstates, respectively. For this reason, we first check thedensity of continuum states obtained with our B-splinebasis set. After that, we show how orbital energies for theH and He atoms are influenced by the range-separationparameter µ . Finally, having in mind these aspects, wediscuss the different calculated spectra. All the stud-ied transitions correspond to dipole-allowed spin-singlettransitions from the Lyman series, i.e. 1s → n p. A. Density of continuum states
In Fig. 1, the radial density of states (DOS) of a freeparticle in a spherical box is compared with the radialDOS of the continuum p orbitals of the H atom computedwith the exact Hamiltonian or with the HF or LDA effec-tive Hamiltonian using the B-spline basis set. The radialDOS of a free particle is given by [53] ρ ( ε ) = R max /π √ ε where R max is the radial size of the box, while for the dif-ferent Hamiltonians using the B-spline basis set (with thesame R max ) the radial DOS is calculated by finite differ-ences as ρ ( ε p ) = 2 / ( ε p +1 − ε p − ) where ε p are positiveorbital energies.As one can observe, the radial DOS computed with theLDA or the HF Hamiltonian is essentially identical to theDOS of the free particle. This can be explained by thefact that since the unoccupied LDA and HF orbitals donot see a − /r attractive potential they are all unboundand they all contribute to the continuum, similarly to thefree-particle case. By contrast, for the exact Hamiltonianwith the same B-spline basis set, one obtains a slightlysmaller DOS in the low-energy region. This is due to thepresence of the − /r attractive Coulomb potential whichsupports a series of bound Rydberg states, necessarilyimplying less unoccupied orbitals in the continuum for agiven basis.We have checked that, by increasing the size of the sim-ulation box, together with the number of B-spline func-tions in the basis so as to keep constant the density ofB-spline functions, the DOS of the exact Hamiltonianconverges, albeit slowly, to the free-particle DOS. Thismust be the case since, for potentials vanishing at infin-ity, the global density of unbound states is independent ofthe potential for an infinite simulation box (only the localDOS depends on the potential, see e.g. Ref. 59). Froma numerical point of view, the computation of the DOScan be seen as a convergence test. With the present basisset, a huge energy range of the continuum spectrum is de-scribed correctly, and the difference between the DOS ofthe exact Hamiltonian and the free-particle DOS at lowenergies (0 . − . − Ha − . Thisdifference is small enough to fairly compare the different orbital energy ε p (Ha) D O S ( H a − ) R max /π q ε p ExactHFLDA ε p (Ha) FIG. 1. Radial density of states (DOS) for a free particle, ρ ( ε p ) = R max /π √ ε p , in a spherical box of size R max = 100bohr, and for the continuum p orbitals of the H atom com-puted with the exact Hamiltonian, or with the HF or LDAeffective Hamiltonian using the B-spline basis set with thesame R max . methods considered in this paper.The calculation of the DOS is also important in or-der to compute proper oscillator strengths involving con-tinuum states. Because of the use of a finite simula-tion box, the calculated positive-energy orbitals form, ofcourse, a discrete set and not strictly a continuum. Thesepositive-energy orbitals are thus not energy normalizedas the exact continuum states should be. To better ap-proximate pointwise the exact continuum wave functions,the obtained positive-energy orbitals should be renormal-ized. Following Mac´ıas et al. [60], we renormalize thepositive-energy orbitals by the square root of the DOSas ˜ ϕ p ( r ) = (cid:112) ρ ( ε p ) ϕ p ( r ). B. Range-separated orbital energies
In Fig. 2 we show the 1s and the low-lying p orbital en-ergies for the H atom calculated with both the RSH andRSH-EXX methods as a function of the range-separationparameter µ .As one observes in Fig. 2a, with the RSH method onlythe 1s ground state is bound, and the energy of this stateis strongly dependent on µ . At µ = 0, the self-interactionerror introduced by the LDA exchange-correlation po-tential is maximal. But, when µ increases, the long-range HF exchange potential progressively replaces thelong-range part of the LDA exchange-correlation poten-tial and the self-interaction error is gradually eliminateduntil reaching the HF limit for µ → ∞ , where one ob-tains the exact 1s orbital energy. The p orbitals (and allthe other unoccupied orbitals) are always unbound andtheir (positive) energies are insensible to the value of µ .One also observes that the approximate continuum of porbitals has a DOS correctly decreasing as the energyincreases, as previously seen in Fig. 1. In Fig. 2b, one sees that the 1s orbital energy com-puted with the RSH-EXX method is identical to the 1sorbital energy obtained by the RSH scheme, as expected.However, a very different behavior is observed for theunoccupied p orbitals. Starting from the LDA limit at µ = 0 where all unoccupied orbitals are unbound, whenthe value of µ increases one sees the emergence of a seriesof bound Rydberg states coming down from the contin-uum. This is due to the introduction of an attractive − /r term in the long-range EXX potential, which sup-ports a Rydberg series. For µ → ∞ , we obtain the spec-trum of the exact hydrogen Hamiltonian calculated withthe B-spline basis set. Necessarily, with the finite ba-sis used, the appearance of the discrete bound states isaccompanied by a small reduction of the density of con-tinuum states, as we already observed in Fig. 1 with theexact Hamiltonian.Another interesting aspect that can be observed in Fig.2b is the fact that the different bound-state energies reachtheir exact µ → ∞ values at different values of µ . Thus,for a fixed small value of µ , each bound-state energy isaffected differently by the self-interaction error. For thecompact 1s orbital, the self-interaction error is eliminatedfor µ (cid:38) − . For the more diffuse 2p Rydbergstate, the self-interaction error is essentially eliminatedwith µ (cid:38) . − . When we continue to climb inthe Rydberg series, the orbitals become more and morediffuse and the self-interaction error is eliminated fromsmaller and smaller values of µ .In Fig. 3, the 1s and low-lying p orbital energies for theHe atom are shown. Again, for the RSH method, one seesin Fig. 3a that only the occupied 1s orbital is bound andall the unoccupied p orbitals are in the continuum. Sim-ilarly to the case of the H atom, at µ = 0 the 1s orbitalenergy is too high, which can essentially be attributed tothe self-interaction error in the LDA exchange-correlationpotential. This error decreases when µ increases and the1s orbital energy converges to its HF value for µ → ∞ .However, contrary to the case of the H atom, for this two-electron system, the 1s HF orbital energy is not equal tothe opposite of the exact ionization energy but is slightlytoo low due to missing correlation effects. In the spirit ofthe optimally tuned range-separated hybrids [16, 62, 63],the range-separation parameter µ can be chosen so thatthe HOMO orbital energy is equal to the opposite of theexact ionization energy, which gives µ = 1 .
115 bohr − for the He atom.As regards the RSH-EXX method, one sees again inFig. 3b that, for this two-electron system, the 1s RSH-EXX orbital energy is identical to the 1s RSH orbitalenergy. As in the case of the H atom, the introduction ofthe long-range EXX potential generates a series of boundRydberg states, whose energies converge to the Kohn-Sham EXX orbital energies for µ → ∞ . For the Rydbergstates of the He atom, it turns out that the Kohn-ShamEXX orbital energies are practically identical to the ex-act Kohn-Sham orbital energies [61], implying that theKohn-Sham correlation potential has essentially no effect µ (bohr − ) -0.50-0.250.000.20 o r b i t a l e n e r g y ε p ( H a ) (a) µ (bohr − ) -0.50-0.250.000.20 o r b i t a l e n e r g y ε p ( H a ) (b) FIG. 2. Orbital energies obtained with the RSH (a) and with the RSH-EXX (b) methods as a function of range-separationparameter µ for the H atom. The occupied 1s orbital energy is plotted in red and the unoccupied p orbital energies are plottedin blue. Horizontal dotted lines indicate the exact 1s orbital energy ( − . µ (bohr − ) -1.0-0.50.00.2 o r b i t a l e n e r g y ε p ( H a ) (a) µ (bohr − ) -1.0-0.50.00.2 o r b i t a l e n e r g y ε p ( H a ) (b) FIG. 3. Orbital energies obtained with the RSH (a) and with the RSH-EXX (b) methods as a function of range-separationparameter µ for the He atom. The occupied 1s orbital energy is plotted in red and the unoccupied p orbital energies are plottedin blue. Horizontal dotted lines indicate exact Kohn-Sham orbital energies [61], including the opposite of the exact ionizationenergy ( − . on these Rydberg states. As we will see, contrary to theRSH case, the set of unoccupied RSH-EXX orbitals canbe considered as a reasonably good first approximationfor the computation of photoexcitation and photoioniza-tion spectra, even before applying linear-response theory. C. Photoexcitation and photoionizationspectra for the hydrogen atom
In Fig. 4, photoexcitation/photoionization spectra forthe H atom calculated with different methods are shown.For the calculation using the exact Hamiltonian, thespectrum is correctly divided into a discrete and a con-tinuum part, corresponding to the photoexcitation andphotoionization processes, respectively. As already dis-cussed in Sec. IV A, for all calculations, the continuumstates have been renormalized, or equivalently the oscil-lator strengths of the continuum part of the spectrum have been renormalized as ˜ f → n p = ρ ( ε n p ) f → n p where ρ ( ε n p ) is the DOS at the corresponding positive orbitalenergy ε n p . Moreover, for better readability of the spec-tra, following Refs. 11, 64, and 65, we have also renor-malized the oscillator strengths of the discrete part of thespectrum as ˜ f → n p = n f → n p where n is the principalquantum number of the excited p orbital. This makes thetransition between the discrete and the continuum partof the spectrum smooth. Another thing is, since we areworking with a finite B-spline basis set principally target-ing a good continuum, we obtain only a limited numberof Rydberg states and the last Rydberg states near theionization threshold are not accurately described. In par-ticular, the corresponding oscillator strengths are overes-timated (not shown). To fix this problem, we could forexample use quantum defect theory in order to accuratelyextract the series of Rydberg states [64, 66–68]. However,for the propose of the present work, we did not find neces-sary to do that, and instead we have simply corrected the excitation energy ω n (Ha) o s c ill a t o r s t r e n g t h ˜ f n (a) ExactHFLDATDLDA excitation energy ω n (Ha) o s c ill a t o r s t r e n g t h ˜ f n (b) ExactRSH ( µ =0 . RSH-EXX ( µ =0 . TDRSH ( µ =0 . FIG. 4. Photoexcitation/photoionization spectra calculatedwith different methods for the H atom. In (a) comparisonof the HF, LDA, and TDLDA methods with respect to thecalculation with the exact Hamiltonian. In (b) comparisonof the RSH, RSH-EXX, and TDRSH methods (all of themwith a range-separation parameter of µ = 0 . − ) withrespect to the calculation with the exact Hamiltonian. r (bohr) -0.250.000.50 r a d i a l a m p li t ud e ExactHFLDARSH ( µ =0 . RSH-EXX ( µ =0 . FIG. 5. Comparison of the renormalized radial amplitude˜ R ( r ) = (cid:112) ρ ( ε ) R ( r ) of the continuum p orbital involved in thetransition energy ω n = ε − ε = 0 . µ = 0 . − ) with respect to the exact calculationfor the H atom. oscillator strengths of the last Rydberg states by inter- polating between the oscillator strengths of the first fiveRydberg states and the oscillator strength of the firstcontinuum state using a second-order polynomial func-tion of the type ˜ f n = c + c ω n + c ω n . This procedurewas applied for all spectra having a discrete part.Let us first discuss the spectra in Fig. 4a. The LDAspectrum, calculated using the bare oscillator strengthsof Eq. (12), does not possess a discrete photoexcitationpart, which was of course expected since the LDA po-tential does not support bound Rydberg states, as seenin the µ = 0 limit of Fig. 2. The ionization thresh-old energy, giving the onset of the continuum spec-trum, is much lower than the exact value (0.5 Ha) dueto the self-interaction error in the ground-state orbitalenergy. At the ionization threshold, the LDA oscilla-tor strengths are zero, in agreement with the Wigner-threshold law [69, 70] for potentials lacking a long-rangeattractive − /r Coulomb tail. Close above the ionizationthreshold, the LDA spectrum has an unphysical largepeak, which corresponds to continuum states with animportant local character. However, as noted in Ref. 11,at the exact Rydberg transition energies, the LDA con-tinuum oscillator strengths are actually reasonably goodapproximations to the exact discrete oscillator strengths,which was explained by the fact that the LDA potential isapproximately the exact Kohn-Sham potential shifted bya constant. Moreover, above the exact ionization energy,LDA reproduces relatively well the exact photoionizationspectrum and becomes essentially asymptotically exactin the high-energy limit. This is consistent with the factthat, at a sufficiently high transition energy, the LDAcontinuum orbitals are very similar to the exact ones, atleast in the spatial region relevant for the calculation ofthe oscillation strengths, as shown in Fig. 5.The TDLDA spectrum differs notably from the LDAspectrum only in that the unphysical peak at around0 . µ → ∞ limit of Fig. 2a),but does not even look as a photoionization spectrum.The HF unoccupied orbitals actually represent approx-imations to the continuum states of the H − anion, andare thus much more diffuse than the exact continuumstates of the H atom, as shown in Fig. 5. Consequently,the HF spectrum has in fact the characteristic shape ofthe photodetachment spectrum of the H − anion [71, 72](with the caveat that the initial state is the 1s orbital ofthe H atom instead of the 1s orbital of the H − anion).Finally, note that, for the H atom, linear-response TDHFgives of course the exact photoexcitation/photoionizationspectrum.Let us now discuss the spectra obtained with the range-separated methods in Fig. 4b. The common value ofthe range-separation parameter µ = 0 . − has beenused [20]. The RSH spectrum looks like the photode-tachment spectrum of the H − anion. This is not sur-prising since the RSH effective Hamiltonian contains along-range HF exchange potential. The RSH continuumorbitals are similarly diffuse as the HF continuum or-bitals, as shown in Fig. 5. The RSH ionization thresh-old energy is slightly smaller than the exact value (0.5Ha) due to the remaining self-interaction error in the1s orbital energy stemming from the short-range LDAexchange-correlation potential at this value of µ . TheRSH-EXX ionization threshold is identical to the RSHone, but, contrary to the RSH spectrum, the RSH-EXXspectrum correctly shows a discrete photoexcitation partand a continuum photoionization part. Beside the smallredshift of the spectrum, the self-interaction error at thisvalue of µ manifests itself in slightly too small RSH-EXXoscillator strengths. The RSH-EXX continuum orbitalsare very similar to the exact continuum orbitals, as shownin Fig. 5. Finally, at this value of µ , TDRSH gives a pho-toexcitation/photoionization spectrum essentially identi-cal to the RSH-EXX spectrum. D. Photoexcitation and photoionizationspectra for the helium atom
In Fig. 6, different photoexcitation/photoionizationspectra for the He atom are shown. As in the H atomcase, the oscillator strengths of the discrete part of theTDHF, RSH-EXX, and TDRSH spectra have been in-terpolated (using again the oscillator strengths of firstfive Rydberg states and of the first continuum state)to correct the overestimation of the oscillator strengthsfor the last Rydberg transitions. The excitation energiesand the (non-interpolated) oscillator strengths of the firstfive discrete transitions are reported in Table I and com-pared with exact results. The photoionization part ofsome of the calculated spectra are compared with fullconfiguration-interaction (FCI) calculations and experi-mental results in Fig. 7.In Fig. 6a, one sees that the HF spectrum looks againlike a photodetachment spectrum, corresponding in thiscase to the He − anion. By contrast, TDHF gives a rea-sonable photoexcitation/photoionization spectrum. Inparticular, for the first discrete transitions listed in Ta-ble I, TDHF gives slightly too large excitation energiesby at most about 0.02 Ha (or 0.5 eV) and slightly toosmall oscillator strengths by at most about 0.025. Theionization energy is also slightly too large by about 0.015Ha, as already seen from the HF 1s orbital energy in the µ → ∞ limit of Fig. 3. As regards the photoionizationpart of the spectrum, one sees in Fig. 7 that TDHF gives excitation energy ω n (Ha) o s c ill a t o r s t r e n g t h ˜ f n (a) HFTDHFLDATDLDA ω n (Ha) o s c ill a t o r s t r e n g t h (b) TDHFRSH ( µ =0 . TDRSH ( µ =0 . RSH-EXX ( µ =0 . FIG. 6. Photoexcitation and photoionization spectra calcu-lated with different methods for the He atom. In (a) compar-ison of HF, TDHF, LDA, and TDLDA methods. In (b) com-parison of RSH, RSH-EXX, and TDRSH methods (all of themwith a range-separation parameter of µ = 1 .
115 bohr − ).
25 30 35 40 45 50 55 excitation energy ω n (eV) c r o ss s e c t i o n ( M b ) Exp. FCITDLDATDHFTDRSH ( µ =1 . ) FIG. 7. Photoionization cross-section profile for the He atom.Normalized cross sections are given (in Hartree atomic units)by σ n = (2 π /c ) ˜ f n where ˜ f n are the renormalized oscillatorstrengths and c is the speed of light. Conversion factors 1 Ha= 27.207696 eV and 1 bohr = 28 . slightly too large photoionization cross sections.The LDA spectrum in Fig. 6a is also similar to the0 TABLE I. Excitation energies ( ω n in Ha) and oscillator strengths ( f n ) of the first discrete transitions calculated with differentmethods for the He atom. The ionization energy is also given.Transition Exact a TDHF RSH-EXX ( µ = 1 . µ = 1 . ω n f n ω n f n ω n f n ω n f n S → P 0.7799 0.2762 0.7970 0.2518 0.7766 0.3303 0.7827 0.25471 S → P 0.8486 0.0734 0.8636 0.0704 0.8474 0.0857 0.8493 0.07081 S → P 0.8727 0.0299 0.8872 0.0291 0.8721 0.0344 0.8729 0.02921 S → P 0.8838 0.0150 0.8982 0.0148 0.8835 0.0172 0.8839 0.01481 S → P 0.8899 0.0086 0.9042 0.0087 0.8897 0.0100 0.8899 0.0087Ionization energy 0.9036 0.9180 0.9036 0.9036 a From Ref. 74.
LDA spectrum for the H atom. The ionization thresholdenergy is much too low, and the spectrum lacks a dis-crete part and has an unphysical maximum close abovethe ionization threshold. Except from that, taking as ref-erence the TDHF spectrum (which is close to the exactspectrum), the LDA spectrum is a reasonable approx-imation to the photoionization spectrum and, again asnoted in Ref. 11, a reasonable continuous approximationto the photoexcitation spectrum. In comparison to LDA,TDLDA [75] gives smaller and less accurate oscillatorstrengths in the lower-energy part of the spectrum but,the TRK sum rule having to be preserved, larger oscilla-tor strengths in the higher-energy part of the spectrum,resulting in an accurate high-energy asymptotic behavioras seen in Fig. 7.Fig. 6b shows the spectra calculated with RSH, RSH-EXX, and TDRSH using for the range-separation param-eter the value µ = 1 .
115 bohr − which imposes the exactionization energy, as explained in Sec. IV B. The RSHspectrum is similar to the HF spectrum and does notrepresent a photoexcitation/photoionization spectrum.By contrast, the RSH-EXX spectra is qualitatively cor-rect for a photoexcitation/photoionization spectrum. Asshown in Table I, in comparison with TDHF, RSH-EXXgives more accurate Rydberg excitation energies, with alargest error of about 0.003 Ha (or 0.08 eV), but lessaccurate oscillator strengths which are significantly over-estimated. The TDRSH method also gives a correct pho-toexcitation/photoionization spectrum, with the advan-tage that it gives Rydberg excitation energies as accu-rate as the RSH-EXX ones and corresponding oscillatorstrengths as accurate as the TDHF ones. As shown inFig. 7, TDRSH also gives a slightly more accurate pho-toionization cross-section profile than TDHF. V. CONCLUSIONS
We have investigated the performance of the RSHscheme for calculating photoexcitation/photoionizationspectra of the H and He atoms, using a B-spline basisset in order to correctly describe the continuum part ofthe spectra. The study of these simple systems allowedus to quantify the influence on the spectra of the errorscoming from the short-range exchange-correlation LDAand from the missing long-range correlation in the RSHscheme. For the He atom, it is possible to choose a valuefor the range-separation parameter µ for which these er-rors compensate each other so as to obtain the exactionization energy.We have studied the differences between using thelong-range HF exchange nonlocal potential and the long-range EXX local potential. Contrary to the former,the latter supports a series of Rydberg states and thecorresponding RSH-EXX scheme, even without apply-ing linear-response theory, gives reasonable photoexcita-tion/photoionization spectra. Nevertheless, the most ac-curate spectra are obtained with linear-response TDRSH(or TDRSH-EXX since they are equivalent for one- andtwo-electron systems). In particular, for the He atom atthe optimal value of µ , TDRSH gives slightly more ac-curate photoexcitation and photoionization spectra thanstandard TDHF.The present work calls for further developments. First,the merits of TDRSH (and/or TDRSH-EXX) for calcu-lating photoexcitation/photoionization spectra of largeratoms and molecules, where screening effects are impor-tant, should now be investigated. Second, it would be in-teresting to test the effects of going beyond the LDA forthe short-range exchange-correlation functional [76, 77]and adding long-range wave-function correlation [18, 78,79]. Third, time-propagation TDRSH could be imple-mented to go beyond linear response and tackle strong-field phenomena, such as high-harmonic generation andabove-threshold ionization [80]. [1] E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997(1984). [2] E. K. U. Gross and W. Kohn, Phys. Rev. Lett. , 2850 (1985).[3] M. E. Casida, in Recent Advances in Density FunctionalMethods, Part I , edited by D. P. Chong (World Scientific,Singapore, 1995), p. 155.[4] M. Petersilka, U. J. Gossmann, and E. K. U. Gross, Phys.Rev. Lett. , 1212 (1996).[5] M. E. Casida, C. Jamorski, K. C. Casida, and D. R.Salahub, J. Chem. Phys. , 4439 (1998).[6] A. Dreuw, J. L. Weisman, and M. Head-Gordon, J.Chem. Phys. , 2943 (2003).[7] R. van Leeuwen and E. J. Baerends, Phys. Rev. A ,2421 (1994).[8] D. J. Tozer and N. C. Handy, J. Chem. Phys. , 10180(1998).[9] M. E. Casida and D. R. Salahub, J. Chem. Phys. ,8918 (2000).[10] P. R. T. Schipper, O. V. Gritsenko, S. J. A. van Gis-bergen, and E. J. Baerends, J. Chem. Phys. , 1344(2000).[11] A. Wasserman, N. T. Maitra, and K. Burke, Phys. Rev.Lett. , 263001 (2003).[12] A. Wasserman and K. Burke, Phys. Rev. Lett. , 163006(2005).[13] Y. Tawada, T. Tsuneda, S. Yanagisawa, T. Yanai, andK. Hirao, J. Chem. Phys. , 8425 (2004).[14] T. Yanai, D. P. Tew, and N. C. Handy, Chem. Phys. Lett. , 51 (2004).[15] M. J. G. Peach, T. Helgaker, P. Salek, T. W. Keal, O. B.Lutnæs, D. J. Tozer, and N. C. Handy, Phys. Chem.Chem. Phys. , 558 (2006).[16] E. Livshits and R. Baer, Phys. Chem. Chem. Phys. ,2932 (2007).[17] R. Baer, E. Livshits, and U. Salzner, Annu. Rev. Phys.Chem. , 85 (2010).[18] E. Fromager, S. Knecht, and H. J. A. Jensen, J. Chem.Phys. , 084101 (2013).[19] E. Rebolini, A. Savin, and J. Toulouse, Mol. Phys. ,1219 (2013).[20] I. C. Gerber and J. G. ´Angy´an, Chem. Phys. Lett. ,100 (2005).[21] T. Tsuneda, J.-W. Song, S. Suzuki, and K. Hirao, J.Chem. Phys. , 174101 (2010).[22] A. Zangwill and P. Soven, Phys. Rev. A , 1561 (1980).[23] Z. H. Levine and P. Soven, Phys. Rev. A , 625 (1984).[24] M. Stener, P. Decleva, and A. Lisini, J. Phys. B, , 4973(1995).[25] M. Stener, G. De Alti, G. Fronzoni, P. Decleva, Chem.Phys. , 197 (1997).[26] M. Stener and P. Decleva, J. Phys. B, , 4481 (1997).[27] M. Stener and P. Decleva, J. Chem. Phys. , 10871(2000).[28] M. Stener, P. Decleva, and A. G¨orling, J. Chem. Phys. , 7816 (2001).[29] M. Stener, G. Fronzoni, and P. Decleva, J. Chem. Phys. , 234301 (2005).[30] M. Stener, D. Toffoli, G. Fronzoni, and P. Decleva, J.Chem. Phys. , 114306 (2006).[31] D. Toffoli, M. Stener, and P. Decleva, Phys. Rev. A ,042704 (2006).[32] M. Stener, D. Toffoli, G. Fronzoni, and P. Decleva, Theor.Chem. Acc. , 943 (2007).[33] Z. Zhou and S.-I. Chu, Phys. Rev. A , 053412 (2009).[34] K. Lopata and N. Govind, J. Chem. Theory Comput. , 4939 (2013).[35] R. G. Fernando, M. C. Balhoff, and K. Lopata, J. Chem.Theory Comput. , 646 (2015).[36] A. Sissay, P. Abanador, F. Mauger, M. Gaarde, K. J.Schafer, and K. Lopata, J. Chem. Phys. , 094105(2016).[37] J. Toulouse, E. Rebolini, T. Gould, J. F. Dobson, P. Seal,and J. G. ´Angy´an, J. Chem. Phys. , 194106 (2013).[38] A. Savin, in Recent Developments of Modern DensityFunctional Theory , edited by J. M. Seminario (Elsevier,Amsterdam, 1996), pp. 327–357.[39] J. Toulouse, F. Colonna, and A. Savin, Phys. Rev. A ,062505 (2004).[40] J. G. ´Angy´an, I. C. Gerber, A. Savin, and J. Toulouse,Phys. Rev. A , 012510 (2005).[41] S. Paziani, S. Moroni, P. Gori-Giorgi, and G. B. Bachelet,Phys. Rev. B , 155111 (2006).[42] J. D. Talman and W. F. Shadwick, Phys. Rev. A , 36(1976).[43] A. G¨orling and M. Levy, Phys. Rev. A , 196 (1994).[44] A. G¨orling and M. Levy, Int. J. Quantum Chem. Symp. , 93 (1995).[45] J. Toulouse, P. Gori-Giorgi, and A. Savin, Int. J. Quan-tum Chem. , 2026 (2006).[46] J. Toulouse and A. Savin, J. Mol. Struct. (Theochem) , 147 (2006).[47] C. Filippi, C. J. Umrigar, and X. Gonze, Phys. Rev. A , 4810 (1996).[48] A. G¨orling, Phys. Rev. Lett. , 5459 (1999).[49] S. Ivanov, S. Hirata, and R. J. Bartlett, Phys. Rev. Lett. , 5455 (1999).[50] A. G¨orling, Phys. Rev. A , 3433 (1998).[51] A. G¨orling, Int. J. Quantum Chem. , 69 (1998).[52] C. de Boor, A practical Guide to Splines (Springer Books,1978).[53] H. Bachau, E. Cormier, P. Decleva, J. E. Hansen, andF. Mart´ın, Rep. Prog. Phys. , 1815 (2001).[54] R. D. Cowan, The Theory of Atomic Structure and Spec-tra , Los Alamos Series in Basic and Applied Sciences(University of California Press, Ltd., Berkeley, 1981).[55] O. ˇCert´ık, Ph.D. thesis, University of Nevada, Reno(2012).[56] Y. Qiu and C. Froese Fischer, J. Comput. Phys. , 257(1999).[57] J. G. ´Angy´an, I. Gerber, and M. Marsman, J. Phys. A , 8613 (2006).[58] S. L. Marshall, J. Phys.: Condens. Matter , 3175(2002).[59] R. Dick, Advanced Quantum Mechanics – Materials andPhotons (Springer, New York, 2012).[60] A. Mac´ıas, F. Mart´ın, A. Riera, and M. Y´anez, Int. J.Quantum Chem. , 279 (1988).[61] C. Umrigar, A. Savin, and X. Gonze, in Electronic Den-sity Functional Theory , edited by J. F. Dobson, G. Vig-nale, and M. P. Das (Springer US, 1998), pp. 167–176.[62] T. Stein, L. Kronik, and R. Baer, J. Am. Chem. Soc. , 2818 (2009).[63] T. Stein, L. Kronik, and R. Baer, J. Chem. Phys. ,244119 (2009).[64] H. Friedrich,
Theoretical Atomic Physics, 2nd ed. (Springer-Verlag, Berlin, Ltd., 1998).[65] Z.-h. Yang, M. van Faassen, and K. Burke, J. Chem.Phys. , 114308 (2009). [66] A. I. Al-Sharif, R. Resta, and C. J. Umrigar, Phys. Rev.A , 2466 (1998).[67] M. van Faassen and K. Burke, J. Chem. Phys. ,094102 (2006).[68] M. van Faassen and K. Burke, Phys. Chem. Chem. Phys. , 4437 (2009).[69] E. Wigner, Phys. Rev. , 1002 (1948).[70] H. R. Sadeghpour, J. L. Bohn, M. J. Cavagnero, B. D.Esryk, I. I. Fabrikant, J. H. Macek, and A. R. P. Rau, J.Phys. B, , R93 (2000).[71] H. A. Bethe and E. E. Salpeter, Quantum Mechanics ofone- and two-electron atoms (Springer, Berlin, 1957).[72] A. R. P. Rau, J. Astrophys. Astr. , 113 (1996).[73] M. Venuti, P. Decleva, and A. Lisini, J. Phys. B , 5315(1996).[74] A. Kono and S. Hattori, Phys. Rev. A , 2981 (1984).[75] Contrary to our Fig. 6a, the TDLDA spectrum of the Heatom shown in Fig. 6 of Ref. 11 has a larger maximum than the LDA spectrum. This discrepancy is due to thefact that the TDLDA spectrum shown in Ref. 11 comesin fact from Ref. 28 where it was calculated by replacingthe LDA 1s orbital energy by the opposite of the exactionization energy. We have checked that this results notonly in an energy shift of the spectrum but also to largeroscillator strengths. The true TDLDA spectrum of theHe atom is thus the one shown in the present Fig. 6a.[76] J. Toulouse, F. Colonna, and A. Savin, J. Chem. Phys. , 014110 (2005).[77] E. Goll, H.-J. Werner, H. Stoll, T. Leininger, P. Gori-Giorgi, and A. Savin, Chem. Phys. , 276 (2006).[78] E. D. Hedeg˚ard, F. Heiden, S. Knecht, E. Fromager, andH. J. A. Jensen, J. Chem. Phys. , 184308 (2013).[79] E. Rebolini and J. Toulouse, J. Chem. Phys. , 094107(2016).[80] M. Labeye, F. Zapata, E. Coccia, V. V´eniard,J. Toulouse, J. Caillat, R. Ta¨ıeb, and E. Luppi, J. Chem.Theory Comput.14