Linear-response range-separated density-functional theory for atomic photoexcitation and photoionization spectra

TL;DR

This paper evaluates the effectiveness of range-separated hybrid density-functional theory, especially the TDRSH scheme, in accurately predicting atomic photoexcitation and photoionization spectra, highlighting its advantages over traditional methods.

Contribution

It demonstrates that the linear-response TDRSH scheme provides more accurate spectra for atomic systems than standard methods, emphasizing the potential of range-separated DFT for optical property calculations.

Findings

TDRSH yields more accurate spectra than linear-response HF.

Using local long-range exchange supports Rydberg states.

Range-separated DFT effectively models continuum states.

Abstract

We investigate the performance of the range-separated hybrid (RSH) scheme, which combines long-range Hartree-Fock (HF) and a short-range density-functional approximation (DFA), for calculating photoexcitation/photoionization spectra of the H and He atoms, using a B-spline basis set in order to correctly describe the continuum part of the spectra. The study of these simple systems allows us to quantify the influence on the spectra of the errors coming from the short-range exchange-correlation DFA and from the missing long-range correlation in the RSH scheme. We study the differences between using the long-range HF exchange (nonlocal) potential and the long-range exact exchange (local) potential. Contrary to the former, the latter supports a series of Rydberg states and gives reasonable photoexcitation/photoionization spectra, even without applying linear-response theory. The most accurate spectra are obtained with the linear-response time-dependent range-separated hybrid (TDRSH) scheme. In particular, for the He atom at the optimal value of the range-separation parameter, TDRSH gives slightly more accurate photoexcitation and photoioniza-tion spectra than standard linear-response time-dependent HF. More generally, the present work shows the potential of range-separated density-functional theory for calculating linear and nonlinear optical properties involving continuum states.

Figures

Tables

Table 1. Table 1: Excitation energies ( ω n subscript 𝜔 𝑛 \omega_{n} in Ha) and oscillator strengths ( f n subscript 𝑓 𝑛 f_{n} ) of the first discrete transitions calculated with different methods for the He atom. The ionization energy is also given.

Transition	xxxxxExactaxxxxx			xxxxxTDHFxxxxx			RSH-EXX ($\mu =1.115$)			TDRSH ($\mu =1.115$)
	${\omega }_n$	$f_n$		${\omega }_n$	$f_n$		${\omega }_n$	$f_n$		${\omega }_n$	$f_n$
11S $\to$ 21P	0.7799	0.2762		0.7970	0.2518		0.7766	0.3303		0.7827	0.2547
11S $\to$ 31P	0.8486	0.0734		0.8636	0.0704		0.8474	0.0857		0.8493	0.0708
11S $\to$ 41P	0.8727	0.0299		0.8872	0.0291		0.8721	0.0344		0.8729	0.0292
11S $\to$ 51P	0.8838	0.0150		0.8982	0.0148		0.8835	0.0172		0.8839	0.0148
11S $\to$ 61P	0.8899	0.0086		0.9042	0.0087		0.8897	0.0100		0.8899	0.0087
Ionization energy	0.9036			0.9180			0.9036			0.9036
aFrom Ref. Kono and Hattori, 1984.

Equations

$$w_{\text{ee}}(r)=w_{\text{ee}}^{\text{lr}}(r)+w_{\text{ee}}^{\text{sr}}(r),$$

$$w_{\text{ee}}^{\text{lr}}(r)=\frac{\mathrm{erf}(\mu r)}{r},$$

$$w_{\text{ee}}^{\text{sr}}(r)=\frac{\mathrm{erfc}(\mu r)}{r}.$$

$$\left(-\frac{1}{2}{\nabla }^2+v_{\text{ne}}(\mathbf{r})+v_{\text{H}}(\mathbf{r})+v_{\text{xc}}^{\text{sr}}(\mathbf{x})\right){\phi }_p(\mathbf{x})$$

$$+\int v_{\text{x}}^{\text{lr},\text{HF}}(\mathbf{x},{\mathbf{x}}^{\prime }){\phi }_p({\mathbf{x}}^{\prime })\text{d}{\mathbf{x}}^{\prime }={\epsilon }_p{\phi }_p(\mathbf{x}),$$

$$v_{\text{H}}(\mathbf{r})=\int n({\mathbf{x}}^{\prime })w_{\text{ee}}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)\text{d}{\mathbf{x}}^{\prime },$$

$$v_{\text{x}}^{\text{lr},\text{HF}}(\mathbf{x},{\mathbf{x}}^{\prime })=-\sum _i^{\text{occ}}{\phi }_i^{\ast }({\mathbf{x}}^{\prime }){\phi }_i(\mathbf{x})w_{\text{ee}}^{\text{lr}}(|\mathbf{r}-{\mathbf{r}}^{\prime }|),$$

$$v_{\text{xc}}^{\text{sr}}(\mathbf{x})=\frac{\delta {\bar{E}}_{\text{xc}}^{\text{sr}}}{\delta n(\mathbf{x})},$$

$$\left(-\frac{1}{2}{\nabla }^2+v_{\text{ne}}(\mathbf{r})+v_{\text{H}}(\mathbf{r})+v_{\text{xc}}^{\text{sr}}(\mathbf{x})+v_{\text{x}}^{\text{lr},\text{EXX}}(\mathbf{x})\right){\phi }_p(\mathbf{x})$$

$$={\epsilon }_p{\phi }_p(\mathbf{x}),$$

$$v_{\text{x}}^{\text{lr},\text{EXX}}(\mathbf{x})=\frac{\delta E_{\text{x}}^{\text{lr}}}{\delta n(\mathbf{x})},$$

$$v_{\text{x}}^{\text{lr},\text{EXX}}(\mathbf{x})=-v_{\text{H}}^{\text{lr}}(\mathbf{r}),$$

$$v_{\text{x}}^{\text{lr},\text{EXX}}(\mathbf{x})=-\frac{1}{2}{v}_{\text{H}}^{\text{lr}}(\mathbf{r}),$$

$$f_{ia}^0=\frac{2}{3}{\omega }_{ia}^0\sum _{\nu =x,y,z}|d_{\nu ,ia}{|}^2,$$

$$\left(\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ -{\mathbf{B}}^{\ast }& -{\mathbf{A}}^{\ast }\end{array}\right)\left(\begin{array}{c}{\mathbf{X}}_n\\ {\mathbf{Y}}_n\end{array}\right)={\omega }_n\left(\begin{array}{c}{\mathbf{X}}_n\\ {\mathbf{Y}}_n\end{array}\right),$$

$$A_{ia,jb}=({\epsilon }_a-{\epsilon }_i){\delta }_{ij}{\delta }_{ab}+K_{ia,jb},$$

$$B_{ia,jb}=K_{ia,bj},$$

$$K_{ia,jb}$$

$$f_{\text{x}}^{\text{lr},\text{EXX}}({\mathbf{x}}_1,{\mathbf{x}}_2;\omega )=-f_{\text{H}}^{\text{lr}}({\mathbf{r}}_1,{\mathbf{r}}_2),$$

$$f_{\text{x}}^{\text{lr},\text{EXX}}({\mathbf{x}}_1,{\mathbf{x}}_2;\omega )=-\frac{1}{2}{f}_{\text{H}}^{\text{lr}}({\mathbf{r}}_1,{\mathbf{r}}_2),$$

$$f_n=\frac{2}{3}{\omega }_n\sum _{\nu =x,y,z}{\left|d_{\nu ,ia}(X_{n,ia}+Y_{n,ia})\right|}^2.$$

$${\phi }_p(\mathbf{r})=R_{n_p{l}_p}(r)Y_{l_p}^{m_p}(\Omega ),$$

$$R_{n_p{l}_p}(r)=\sum _{\alpha =1}^{N_{\text{s}}}{c}_{\alpha }^{n_p{l}_p}\frac{B_{\alpha }^{k_{\text{s}}}(r)}{r},$$

$$w_{\text{ee}}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)=\frac{1}{{\left(|\mathbf{r}{|}^2+|{\mathbf{r}}^{\prime }{|}^2-2|\mathbf{r}||{\mathbf{r}}^{\prime }|\cos \gamma \right)}^{1/2}},$$

$$w_{\text{ee}}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)=\sum _{k=0}^{\infty }\left[\frac{r_{<}^k}{r_{>}^{k+1}}\right]\sum _{m_k=-k}^k(-1{)}^{m_k}{C}_{-m_k}^k(\Omega )C_{m_k}^k({\Omega }^{\prime }),$$

$$⟨pq|w_{\text{ee}}|tu⟩$$

$$R^k(p,q;t,u)$$

$$R^k(\alpha ,\lambda ;\beta ,\nu )$$

$$w_{\text{ee}}^{\text{sr}}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)$$

$$S^k(r_{>},r_{<};\mu )=\mu {\Phi }^k(\Xi ,\xi ),$$

$${\Phi }^k(\Xi ,\xi )$$

Full text

I Introduction

Nowadays, time-dependent density-functional theory (TDDFT) Runge and Gross (1984), applied within the linear-response formalism Gross and Kohn (1985); Casida (1995); Petersilka et al. (1996), is a widely used approach for calculating photoexcitation spectra (transitions from bound to bound states) 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-density approximation (LDA) and generalized-gradient approximations (GGAs), for the exchange-correlation potential and its associated exchange-correlation kernel do not correctly describe long-range electronic transitions, such as those to Rydberg Casida et al. (1998) and charge-transfer Dreuw et al. (2003) states in atomic and molecular systems. A better description of Rydberg excitations can be obtained with exchange-correlation potential approximations having the correct $-1/r$ long-range asymptotic decay R. van Leeuwen and E. J. Baerends (1994); Tozer and Handy (1998); Casida and Salahub (2000); P. R. T. Schipper, O. V. Gritsenko, S. J. A. van Gisbergen, and E. J. Baerends (2000), even though it has been shown that accurate Rydberg excitation energies and oscillator strengths can in fact be extracted from LDA calculations in small atoms Wasserman et al. (2003); Wasserman and Burke (2005). A more general solution for correcting both Rydberg and charge-transfer excitations is given by range-separated TDDFT approaches Tawada et al. (2004); Yanai et al. (2004); Peach et al. (2006); Livshits and Baer (2007); Baer et al. (2010); Fromager et al. (2013); Rebolini et al. (2013) which express the long-range part of the exchange potential and kernel at the Hartree-Fock (HF) level. These range-separated approaches also give reasonably accurate values for the ionization energy threshold Yanai et al. (2004); Gerber and Ángyán (2005); Tsuneda et al. (2010).

Linear-response TDDFT has also been used for calculating photoionization spectra (transitions from bound to continuum states) of atoms and molecules Zangwill and Soven (1980); Levine and Soven (1984); Stener et al. (1995); M. Stener, G. De Alti, G. Fronzoni, P. Decleva (1997); Stener and Decleva (1997, 2000); Stener et al. (2001, 2005, 2006); Toffoli et al. (2006); Stener et al. (2007); Zhou and Chu (2009). These calculations are less standard in quantum chemistry since they involve spatial grid methods or B-spline basis sets for a proper description of the continuum states. In this case as well, usual (semi-)local DFAs provide a limited accuracy and asymptotically corrected exchange-correlation potential approximations give more satisfactory results. More accurate still, but less common, are photoionization spectra calculated with the exact-exchange (EXX) potential Stener et al. (2001) or the localized HF exchange potential and its associated kernel Zhou and Chu (2009). Recently, range-separated approximations have been successfully used for calculating photoexcitation and photoionization spectra of molecular systems using time-propagation TDDFT with Gaussian basis sets together with an effective lifetime model compensating for the missing continuum states Lopata and Govind (2013); Fernando et al. (2015); Sissay et al. (2016). However, to the best of our knowledge, range-separated approximations have not yet been used in frequency-domain linear-response TDDFT calculations of photoionization spectra.

In this work, we explore the performance of the linear-response time-dependent range-separated hybrid (TDRSH) scheme Rebolini et al. (2013); Toulouse et al. (2013) for calculating photoexcitation and photoionization spectra of the H and He atoms using a B-spline basis set to accurately describe the continuum part of the spectra. The TDRSH scheme allows us to treat long-range exchange effects at the HF level 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 effect of replacing the long-range HF exchange nonlocal potential by the long-range EXX local potential (resulting in a scheme that we refer to as RSH-EXX). Second, oscillator strengths directly computed with the RSH and the RSH-EXX orbitals are compared with oscillator strengths obtained with the linear-response TDRSH scheme. The study of the H atom allows us to quantify the residual self-interaction error coming from the short-range exchange-correlation DFA, and the study of the He atom permits to quantify the effect of the missing long-range correlation in the RSH scheme. This work constitutes a first step for applying range-separated TDDFT to strong-field phenomena, such as high-harmonic generation or above-threshold ionization, where long-range effects 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 the RSH-EXX variant, and, secondly, we review the linear-response TDRSH method. In Sec. III, the basis set of B-spline functions is defined, and we indicate how the range-separated two-electron integrals are computed using an exact spherical harmonic expansion for the range-separated interaction. In Sec. IV results are presented and discussed. Firstly, we show the performance of the B-spline basis set for describing the density of continuum states of the H atom within the different methods. Secondly, the dependence of the orbital energies of the H and He atoms on the range-separation parameter is analyzed. Thirdly, different calculated photoexcitation/photoionization spectra for the H and He atoms are discussed and compared with exact results. In Sec. V, conclusions and perspectives are given. Unless otherwise indicated, Hartree atomic units are used throughout the paper.

II Range-separated density-functional theory

II.1 Range-separated hybrid scheme

Range-separated density-functional theory (see, e.g., Refs. Savin, 1996; Toulouse et al., 2004) is based on the splitting of the Coulomb electron-electron interaction $w_{\text{ee}}(r)=1/r$ into long-range (lr) and short-range (sr) contributions

[TABLE]

and the most common forms for the long-range and short-range interactions are

[TABLE]

and

[TABLE]

where $\operatorname{erf}$ and $\operatorname{erfc}$ are the error function and the complementary error function, respectively, and $\mu$ is a tunable range-separation parameter controlling the range of the separation. Using this decomposition, it is possible to rigorously combine a long-range wave-function approach with a complementary short-range DFA.

The simplest approach in range-separated density-functional theory consists in using a single-determinant wave function for the long-range interaction. This leads to the RSH scheme Ángyán et al. (2005) which spin orbitals $\{\varphi_{p}(\mathbf{x})\}$ (where $\mathbf{x}=(\mathbf{r},\sigma)$ are space-spin coordinates) and orbital energies $\varepsilon_{p}$ can be determined for a given system by the following eigenvalue problem,

[TABLE]

where $v_{\text{ne}}(\mathbf{r})$ is the nuclei-electron potential, $v_{\text{H}}(\mathbf{r})$ is the Hartree potential for the Coulomb electron-electron interaction,

[TABLE]

where $n(\mathbf{x})=\sum_{i}^{\text{occ}}|\varphi_{i}(\mathbf{x})|^{2}$ are the spin densities ($i$ refers to occupied spin orbitals), $v_{\text{x}}^{\text{lr},\text{HF}}(\mathbf{x},\mathbf{x}^{\prime})$ is the nonlocal HF exchange potential for the long-range electron-electron interaction,

[TABLE]

and $v_{\text{xc}}^{\text{sr}}(\mathbf{x})$ is the short-range exchange-correlation potential

[TABLE]

where $\bar{E}_{\text{xc}}^{\text{sr}}$ is the complement short-range exchange-correlation density functional. In this work, we use the short-range spin-dependent LDA exchange-correlation functional of Ref. Paziani et al., 2006 for $\bar{E}_{\text{xc}}^{\text{sr}}$. The long-range and short-range potentials, $v_{\text{x}}^{\text{lr},\text{HF}}(\mathbf{x},\mathbf{x}^{\prime})$ and $v_{\text{xc}}^{\text{sr}}(\mathbf{x})$, explicitly depend on the range-separation parameter $\mu$, and consequently the spin orbitals, the orbital energies, and the density also implicitly depend on it. For $\mu=0$, $v_{\text{x}}^{\text{lr},\text{HF}}(\mathbf{x},\mathbf{x}^{\prime})$ vanishes and $v_{\text{xc}}^{\text{sr}}(\mathbf{x})$ becomes the usual full-range LDA exchange-correlation potential, and thus the RSH scheme reduces to standard Kohn-Sham LDA. For $\mu\to\infty$, $v_{\text{x}}^{\text{lr},\text{HF}}(\mathbf{x},\mathbf{x}^{\prime})$ becomes the usual full-range HF exchange potential and $v_{\text{xc}}^{\text{sr}}(\mathbf{x})$ vanishes, and thus the RSH scheme reduces to standard HF.

In the present paper, we also consider the following variant of the RSH scheme,

[TABLE]

in which the long-range nonlocal HF exchange potential has been replaced by the long-range local EXX Talman and Shadwick (1976); Görling and Levy (1994, 1995) potential

[TABLE]

where $E_{\text{x}}^{\text{lr}}$ is the long-range exchange density functional Toulouse et al. (2006); Toulouse and Savin (2006). We will refer to this scheme as RSH-EXX. The calculation of the EXX potential is involved Filippi et al. (1996); Görling (1999); Ivanov et al. (1999), with the exception of one- and two-electron systems. Indeed, for one-electron systems, the long-range EXX potential is simply

[TABLE]

and for systems of two electrons in a single spatial orbital, it is

[TABLE]

where $v^{\text{lr}}_{\text{H}}(\mathbf{r})=\int n(\mathbf{x}^{\prime})w_{\text{ee}}^{\text{lr}}(|\mathbf{r}-\mathbf{r}^{\prime}|)\text{d}\mathbf{x}^{\prime}$ is the long-range Hartree potential. For these one- and two-electron cases, it can be shown that Eqs. (4) and (8) give identical occupied orbitals but different unoccupied orbitals. More generally, for systems with more than two electrons, the HF and EXX exchange potentials give similar occupied orbitals but very different unoccupied orbitals.

Once orbitals and orbital energies are obtained from Eqs. (4) and (8), the bare oscillator strengths can be calculated. They are defined as

[TABLE]

where $i$ and $a$ refer to occupied and unoccupied spin orbitals, respectively, $\omega_{ia}^{0}=\varepsilon_{a}-\varepsilon_{i}$ are the bare excitation energies and $d_{\nu,ia}=\int\varphi_{i}^{*}(\mathbf{x})r_{\nu}\varphi_{a}(\mathbf{x})\text{d}\mathbf{x}$ are the dipole-moment transition integrals. We will consider these bare excitation energies $\omega_{ia}^{0}$ and oscillator strengths $f_{ia}^{0}$ for a first approximation to photoexcitation/photoionization spectra.

II.2 Linear-response time-dependent range-separated hybrid scheme

In the time-dependent extension of the RSH scheme within linear response (referred to as TDRSH) Rebolini et al. (2013); Toulouse et al. (2013); Fromager et al. (2013), one has to solve the following pseudo-Hermitian eigenvalue equation

[TABLE]

whose solutions come in pairs: excitation energies $\omega_{n}>0$ with eigenvectors $(\mathbf{X}_{n},\mathbf{Y}_{n})$, and de-excitation energies $\omega_{n}<0$ with eigenvectors $(\mathbf{Y}_{n}^{*},\mathbf{X}_{n}^{*})$. The elements of the matrices $\mathbf{A}$ and $\mathbf{B}$ are

[TABLE]

[TABLE]

where $i,j$ and $a,b$ refer to occupied and unoccupied RSH spin orbitals, respectively, and the coupling matrix $\mathbf{K}$ contains the contributions from the Hartree kernel $f_{\text{H}}(\mathbf{r}_{1},\mathbf{r}_{2})=w_{\text{ee}}(|\mathbf{r}_{1}-\mathbf{r}_{2}|)$, the long-range HF exchange kernel $f_{\text{x}}^{\text{lr},\text{HF}}(\mathbf{x}_{1},\mathbf{x}_{2};\mathbf{x}_{1}^{\prime},\mathbf{x}_{2}^{\prime})=-w_{\text{ee}}^{\text{lr}}(|\mathbf{r}_{1}-\mathbf{r}_{2}|)\delta(\mathbf{x}_{1}-\mathbf{x}_{2}^{\prime})\delta(\mathbf{x}_{1}^{\prime}-\mathbf{x}_{2})$, and the adiabatic short-range exchange-correlation kernel $f_{\text{xc}}^{\text{sr}}(\mathbf{x}_{1},\mathbf{x}_{2})=\delta v_{\text{xc}}^{\text{sr}}(\mathbf{x}_{1})/\delta n(\mathbf{x}_{2})$

[TABLE]

where $\langle aj|w_{\text{ee}}|ib\rangle$ and $\langle aj|w_{\text{ee}}^{\text{lr}}|bi\rangle$ are the two-electron integrals associated with the Coulomb and long-range interactions, respectively, and $\langle aj|f_{\text{xc}}^{\text{sr}}|ib\rangle=\iint\varphi_{a}^{*}(\mathbf{x}_{1})\varphi_{j}^{*}(\mathbf{x}_{2})f_{\text{xc}}^{\text{sr}}(\mathbf{x}_{1},\mathbf{x}_{2})\varphi_{i}(\mathbf{x}_{1})\varphi_{b}(\mathbf{x}_{2})\text{d}\mathbf{x}_{1}\text{d}\mathbf{x}_{2}$. Since we use the short-range LDA exchange-correlation density functional, for $\mu=0$ the TDRSH scheme reduces to the usual linear-response time-dependent local-density approximation (TDLDA). For $\mu\to\infty$, the TDRSH scheme reduces to standard linear-response time-dependent Hartree-Fock (TDHF).

The time-dependent extension of the RSH-EXX variant within linear response (referred to as TDRSH-EXX) leads to identical equations with the exception that the long-range HF exchange kernel $f_{\text{x}}^{\text{lr},\text{HF}}(\mathbf{x}_{1},\mathbf{x}_{2};\mathbf{x}_{1}^{\prime},\mathbf{x}_{2}^{\prime})$ is replaced by the long-range frequency-dependent EXX kernel Görling (1998, 1998) $f_{\text{x}}^{\text{lr},\text{EXX}}(\mathbf{x}_{1},\mathbf{x}_{2};\omega)=\delta v_{\text{x}}^{\text{lr},\text{EXX}}(\mathbf{x}_{1},\omega)/\delta n(\mathbf{x}_{2},\omega)$. For one-electron systems, the long-range EXX kernel is simply

[TABLE]

and, for systems with two electrons in a single spatial orbital, it is

[TABLE]

where $f_{\text{H}}^{\text{lr}}(\mathbf{r}_{1},\mathbf{r}_{2})=w_{\text{ee}}^{\text{lr}}(|\mathbf{r}_{1}-\mathbf{r}_{2}|)$ is the long-range Hartree kernel. For these one- and two-electron cases, TDRSH and TDRSH-EXX give rise to identical excitation energies and oscillator strengths.

Finally, we can calculate the corresponding TDRSH (or TDRSH-EXX) oscillator strengths as

[TABLE]

In the limit of a complete basis set, the linear-response oscillator strengths in Eq. (19) always fulfill the Thomas-Reiche-Kuhn (TRK) sum rule, $\sum_{n}f_{n}=N$ where $N$ is the electron number. The bare oscillator strengths of Eq. (12) fulfill the TRK sum rule only in the case where the orbitals have been obtained from an effective local potential, i.e. for LDA and RSH-EXX but not for HF and RSH (see Ref. Toulouse et al., 2013).

III Implementation in

a B-spline basis set

In practice, each spin orbital is decomposed into a product of a spatial orbital and a spin function, $\varphi_{p}(\mathbf{x})=\varphi_{p}(\mathbf{r})\delta_{\sigma_{p},\sigma}$ where $\sigma_{p}$ is the spin of the spin orbital $p$, and we use spin-adapted equations. As we investigate atomic systems, the spatial orbitals are written in spherical coordinates,

[TABLE]

where $Y_{l_{p}}^{m_{p}}(\Omega)$ are the spherical harmonics ($\Omega$ stands for the angles $\theta,\phi$) and the radial functions $R_{n_{p}l_{p}}(r)$ are expressed as linear combinations of B-spline functions of order $k_{\text{s}}$,

[TABLE]

where $N_{\text{s}}$ is the dimension of the basis. To completely define a basis of B-spline functions, a non-decreasing sequence of $N_{\text{s}}+k_{\text{s}}$ knot points (some knot points are possibly coincident) must be given de Boor (1978). The B-spline function $B^{k_{\text{s}}}_{\alpha}(r)$ is non zero only on the supporting interval $[r_{\alpha},r_{\alpha+k_{\text{s}}}]$ (containing $k_{\text{s}}+1$ consecutive knot points) and is a piecewise function composed of polynomials of degree $k_{\text{s}}-1$ with continuous first $k_{\text{s}}-m$ derivatives across each knot of multiplicity $m$. We have chosen the first and the last knots to be $k_{\text{s}}$-fold degenerate, i.e. $r_{1}=r_{2}=\cdots=r_{k_{\text{s}}}=R_{\text{min}}$ and $r_{{N_{\text{s}}+1}}=r_{{N_{\text{s}}+2}}=\cdots=r_{{N_{\text{s}}+k_{\text{s}}}}=R_{\text{max}}$, while the multiplicity of the other knots is unity. The spatial grid spacing was chosen to be constant in the whole radial space between two consecutive non-coincident points and is given by $\Delta r=R_{\text{max}}/(N_{\text{s}}-k_{\text{s}}+1)$. In the present work, the first and the last B-spline functions were removed from the calculation to ensure zero boundary conditions at $r=R_{\text{min}}$ and $r=R_{\text{max}}$. The results presented in this paper have been obtained using the following parameters: $k_{\text{s}}=8$, $N_{\text{s}}=200$, $R_{\text{min}}=0$, and $R_{\text{max}}=100$ bohr. Moreover, we need to use only s and pz spherical harmonics.

Working with such a B-spline representation, one must compute matrix elements involving integrals over B-spline functions. The principle of the calculation of one-electron and two-electron integrals over B-spline functions are well described by Bachau et al. in Ref. Bachau et al., 2001. We will now briefly review the computation of the standard Coulomb two-electron integrals over B-spline functions, and then we will present the calculation of the long-range or short-range two-electron integrals over B-spline functions, the latter being original to the present work.

III.1 Coulomb two-electron integrals

The Coulomb electron-electron interaction is given by

[TABLE]

where $\mathbf{r}$ and $\mathbf{r}^{\prime}$ are electron vector positions and $\gamma$ is the angle between them. The multipolar expansion for this interaction is

[TABLE]

where $r_{<}=\mathrm{min}(|\mathbf{r}|,|\mathbf{r}^{\prime}|)$ and $r_{>}=\mathrm{max}(|\mathbf{r}|,|\mathbf{r}^{\prime}|)$ and $C_{m_{k}}^{k}(\Omega)=\left(4\pi/(2k+1)\right)^{1/2}Y_{k}^{m_{k}}(\Omega)$ are the renormalized spherical harmonics. The Coulomb two-electron integrals, in the spatial orbital basis, can then be expressed as the sum of products of radial integrals and angular factors

[TABLE]

where $R^{k}(p,q;t,u)$ are the two-dimensional radial Slater integrals 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 coefficients Cowan (1981); Čertík (2012). The coefficient $c^{k}(l,m,l^{\prime},m^{\prime})$ is non zero only if $|l-l^{\prime}|\leq k\leq l+l^{\prime}$ and if $l+l^{\prime}+k$ is an even integer, which makes the sum over $k$ in Eq. (24) exactly terminate. The Slater integrals are defined as

[TABLE]

where $R^{k}(\alpha,\lambda;\beta,\nu)$ are the Slater matrix elements given by

[TABLE]

In order to compute the Slater matrix elements $R^{k}(\alpha,\lambda;\beta,\nu)$, we have implemented the integration-cell algorithm developed by Qiu and Froese Fischer Y. Qiu and C. Froese Fischer (1999). This algorithm exploits all possible symmetries and B-spline properties to evaluate efficiently the integrals in each two-dimensional radial region on which the integrals are defined. Gaussian quadrature is used to compute the integrals in each cell.

III.2 Long-range and short-range two-electron integrals

A closed form of the multipolar expansion of the short-range electron-electron interaction defined in Eq. (3) was determined by Ángyán et al. Ángyán et al. (2006), following a previous work of Marshall Marshall (2002) who applied the Gegenbauer addition theorem to the Laplace transform of Eq. (3). This exact expansion is

[TABLE]

where the $\mu$-dependent radial function is written in terms of the scaled radial coordinates $\Xi=\mu\;r_{>}$ and $\xi=\mu\;r_{<}$ as

[TABLE]

with

[TABLE]

and the introduced auxiliary functions

[TABLE]

and

[TABLE]

In order to arrive at a separable expression in $\Xi$ and $\xi$, Ángyán et al. Ángyán et al. (2006) also introduced a power series expansion of the radial function $\Phi^{k}(\Xi,\xi)$ in the smaller reduced variable $\xi$. However, the range of validity of this expansion truncated to the first few terms is limited to small values of $\xi$, i.e. $\xi\lesssim 1.5$, and higher-order expansions show spurious oscillations. After some tests, we decided to use the exact short-range radial function $\Phi^{k}(\Xi,\xi)$ without expansion in our work.

The expression of the short-range two-electron integrals $\langle pq|w^{\text{sr}}_{\text{ee}}|tu\rangle$ is then identical to the one in Eq. (24) with the simple difference that the radial term is not given by the standard Slater matrix elements. Now, the radial kernel in Eq. (26) is changed to that of Eq. (28). Due to the fact that the radial kernel is not multiplicatively separable in the variables $r_{>}$ and $r_{<}$, the integration-cell algorithm is modified in order to calculate all integrals as non-separable two-dimensional integrals. In a second step, the long-range two-electron integrals can be simply obtained by difference

[TABLE]

IV Results and discussion

In this section, photoexcitation and photoionization spectra for the H and He atoms are presented. Photoexcitation and photoionization processes imply transitions from bound to bound and from bound to continuum states, respectively. For this reason, we first check the density of continuum states obtained with our B-spline basis set. After that, we show how orbital energies for the H and He atoms are influenced by the range-separation parameter $\mu$. Finally, having in mind these aspects, we discuss the different calculated spectra. All the studied transitions correspond to dipole-allowed spin-singlet transitions from the Lyman series, i.e. $1\text{s}\rightarrow n\text{p}$.

IV.1 Density of continuum states

In Fig. 1, the radial density of states (DOS) of a free particle in a spherical box is compared with the radial DOS of the continuum p orbitals of the H atom computed with the exact Hamiltonian or with the HF or LDA effective Hamiltonian using the B-spline basis set. The radial DOS of a free particle is given by Bachau et al. (2001) $\rho(\varepsilon)=R_{\text{max}}/\pi\sqrt{2\varepsilon}$ where $R_{\text{max}}$ is the radial size of the box, while for the different Hamiltonians using the B-spline basis set (with the same $R_{\text{max}}$) the radial DOS is calculated by finite differences as $\rho(\varepsilon_{p})=2/(\varepsilon_{p+1}-\varepsilon_{p-1})$ where $\varepsilon_{p}$ are positive orbital energies.

As one can observe, the radial DOS computed with the LDA or the HF Hamiltonian is essentially identical to the DOS of the free particle. This can be explained by the fact that since the unoccupied LDA and HF orbitals do not see a $-1/r$ attractive potential they are all unbound and they all contribute to the continuum, similarly to the free-particle case. By contrast, for the exact Hamiltonian with the same B-spline basis set, one obtains a slightly smaller DOS in the low-energy region. This is due to the presence of the $-1/r$ attractive Coulomb potential which supports a series of bound Rydberg states, necessarily implying less unoccupied orbitals in the continuum for a given basis.

We have checked that, by increasing the size of the simulation box, together with the number of B-spline functions in the basis so as to keep constant the density of B-spline functions, the DOS of the exact Hamiltonian converges, albeit slowly, to the free-particle DOS. This must be the case since, for potentials vanishing at infinity, the global density of unbound states is independent of the potential for an infinite simulation box (only the local DOS depends on the potential, see e.g. Ref. Dick, 2012). From a numerical point of view, the computation of the DOS can be seen as a convergence test. With the present basis set, a huge energy range of the continuum spectrum is described correctly, and the difference between the DOS of the exact Hamiltonian and the free-particle DOS at low energies ($0.0-0.2$ Ha) is only about $10^{-4}$ Ha*-1*. This difference is small enough to fairly compare the different methods considered in this paper.

The calculation of the DOS is also important in order to compute proper oscillator strengths involving continuum states. Because of the use of a finite simulation box, the calculated positive-energy orbitals form, of course, a discrete set and not strictly a continuum. These positive-energy orbitals are thus not energy normalized as the exact continuum states should be. To better approximate pointwise the exact continuum wave functions, the obtained positive-energy orbitals should be renormalized. Following Macías et al. Macías et al. (1988), we renormalize the positive-energy orbitals by the square root of the DOS as $\tilde{\varphi}_{p}(\mathbf{r})=\sqrt{\rho(\varepsilon_{p})}\varphi_{p}(\mathbf{r})$.

IV.2 Range-separated orbital energies

In Fig. 2 we show the 1s and the low-lying p orbital energies for the H atom calculated with both the RSH and RSH-EXX methods as a function of the range-separation parameter $\mu$.

As one observes in Fig. 2a, with the RSH method only the 1s ground state is bound, and the energy of this state is strongly dependent on $\mu$. At $\mu=0$, the self-interaction error introduced by the LDA exchange-correlation potential is maximal. But, when $\mu$ increases, the long-range HF exchange potential progressively replaces the long-range part of the LDA exchange-correlation potential and the self-interaction error is gradually eliminated until reaching the HF limit for $\mu\to\infty$, where one obtains the exact 1s orbital energy. The p orbitals (and all the other unoccupied orbitals) are always unbound and their (positive) energies are insensible to the value of $\mu$. One also observes that the approximate continuum of p orbitals has a DOS correctly decreasing as the energy increases, as previously seen in Fig. 1.

In Fig. 2b, one sees that the 1s orbital energy computed with the RSH-EXX method is identical to the 1s orbital energy obtained by the RSH scheme, as expected. However, a very different behavior is observed for the unoccupied p orbitals. Starting from the LDA limit at $\mu=0$ where all unoccupied orbitals are unbound, when the value of $\mu$ increases one sees the emergence of a series of bound Rydberg states coming down from the continuum. This is due to the introduction of an attractive $-1/r$ term in the long-range EXX potential, which supports a Rydberg series. For $\mu\to\infty$, we obtain the spectrum of the exact hydrogen Hamiltonian calculated with the B-spline basis set. Necessarily, with the finite basis used, the appearance of the discrete bound states is accompanied by a small reduction of the density of continuum states, as we already observed in Fig. 1 with the exact Hamiltonian.

Another interesting aspect that can be observed in Fig. 2b is the fact that the different bound-state energies reach their exact $\mu\to\infty$ values at different values of $\mu$. Thus, for a fixed small value of $\mu$, each bound-state energy is affected differently by the self-interaction error. For the compact 1s orbital, the self-interaction error is eliminated for $\mu\gtrsim 1$ bohr*-1*. For the more diffuse 2p Rydberg state, the self-interaction error is essentially eliminated with $\mu\gtrsim 0.5$ bohr*-1*. When we continue to climb in the Rydberg series, the orbitals become more and more diffuse and the self-interaction error is eliminated from smaller and smaller values of $\mu$.

In Fig. 3, the 1s and low-lying p orbital energies for the He atom are shown. Again, for the RSH method, one sees in Fig. 3a that only the occupied 1s orbital is bound and all the unoccupied p orbitals are in the continuum. Similarly to the case of the H atom, at $\mu=0$ the 1s orbital energy is too high, which can essentially be attributed to the self-interaction error in the LDA exchange-correlation potential. This error decreases when $\mu$ increases and the 1s orbital energy converges to its HF value for $\mu\to\infty$. However, contrary to the case of the H atom, for this two-electron system, the 1s HF orbital energy is not equal to the opposite of the exact ionization energy but is slightly too low due to missing correlation effects. In the spirit of the optimally tuned range-separated hybrids Livshits and Baer (2007); Stein et al. (2009a, b), the range-separation parameter $\mu$ can be chosen so that the HOMO orbital energy is equal to the opposite of the exact ionization energy, which gives $\mu=1.115$ bohr*-1* for the He atom.

As regards the RSH-EXX method, one sees again in Fig. 3b that, for this two-electron system, the 1s RSH-EXX orbital energy is identical to the 1s RSH orbital energy. As in the case of the H atom, the introduction of the long-range EXX potential generates a series of bound Rydberg states, whose energies converge to the Kohn-Sham EXX orbital energies for $\mu\to\infty$. For the Rydberg states of the He atom, it turns out that the Kohn-Sham EXX orbital energies are practically identical to the exact Kohn-Sham orbital energies Umrigar et al. (1998), implying that the Kohn-Sham correlation potential has essentially no effect on these Rydberg states. As we will see, contrary to the RSH case, the set of unoccupied RSH-EXX orbitals can be considered as a reasonably good first approximation for the computation of photoexcitation and photoionization spectra, even before applying linear-response theory.

IV.3 Photoexcitation and photoionization

spectra for the hydrogen atom

In Fig. 4, photoexcitation/photoionization spectra for the H atom calculated with different methods are shown. For the calculation using the exact Hamiltonian, the spectrum is correctly divided into a discrete and a continuum part, corresponding to the photoexcitation and photoionization processes, respectively. As already discussed in Sec. IV.1, for all calculations, the continuum states have been renormalized, or equivalently the oscillator strengths of the continuum part of the spectrum have been renormalized as $\tilde{f}_{1\text{s}\to n\text{p}}=\rho(\varepsilon_{n\text{p}})f_{1\text{s}\to n\text{p}}$ where $\rho(\varepsilon_{n\text{p}})$ is the DOS at the corresponding positive orbital energy $\varepsilon_{n\text{p}}$. Moreover, for better readability of the spectra, following Refs. Friedrich, 1998; Wasserman et al., 2003; Z.-h. Yang, M. van Faassen, and K. Burke, 2009, we have also renormalized the oscillator strengths of the discrete part of the spectrum as $\tilde{f}_{1\text{s}\to n\text{p}}=n^{3}f_{1\text{s}\to n\text{p}}$ where $n$ is the principal quantum number of the excited p orbital. This makes the transition between the discrete and the continuum part of the spectrum smooth. Another thing is, since we are working with a finite B-spline basis set principally targeting a good continuum, we obtain only a limited number of Rydberg states and the last Rydberg states near the ionization threshold are not accurately described. In particular, the corresponding oscillator strengths are overestimated (not shown). To fix this problem, we could for example use quantum defect theory in order to accurately extract the series of Rydberg states Al-Sharif et al. (1998); Friedrich (1998); van Faassen and Burke (2006, 2009). However, for the propose of the present work, we did not find necessary to do that, and instead we have simply corrected the oscillator strengths of the last Rydberg states by interpolating between the oscillator strengths of the first five Rydberg states and the oscillator strength of the first continuum state using a second-order polynomial function of the type $\tilde{f}_{n}=c_{0}+c_{1}\;\omega_{n}+c_{2}\;\omega_{n}^{2}$. This procedure was applied for all spectra having a discrete part.

Let us first discuss the spectra in Fig. 4a. The LDA spectrum, calculated using the bare oscillator strengths of Eq. (12), does not possess a discrete photoexcitation part, which was of course expected since the LDA potential does not support bound Rydberg states, as seen in the $\mu=0$ limit of Fig. 2. The ionization threshold energy, giving the onset of the continuum spectrum, is much lower than the exact value (0.5 Ha) due to the self-interaction error in the ground-state orbital energy. At the ionization threshold, the LDA oscillator strengths are zero, in agreement with the Wigner-threshold law Wigner (1948); Sadeghpour et al. (2000) for potentials lacking a long-range attractive $-1/r$ Coulomb tail. Close above the ionization threshold, the LDA spectrum has an unphysical large peak, which corresponds to continuum states with an important local character. However, as noted in Ref. Wasserman et al., 2003, at the exact Rydberg transition energies, the LDA continuum oscillator strengths are actually reasonably good approximations to the exact discrete oscillator strengths, which was explained by the fact that the LDA potential is approximately the exact Kohn-Sham potential shifted by a constant. Moreover, above the exact ionization energy, LDA reproduces relatively well the exact photoionization spectrum and becomes essentially asymptotically exact in the high-energy limit. This is consistent with the fact that, at a sufficiently high transition energy, the LDA continuum orbitals are very similar to the exact ones, at least in the spatial region relevant for the calculation of the oscillation strengths, as shown in Fig. 5.

The TDLDA spectrum differs notably from the LDA spectrum only in that the unphysical peak at around $0.3$ Ha, close above its ionization threshold, has an even larger intensity. This increased intensity comes from the contribution of the LDA exchange-correlation kernel (not shown). The LDA exchange-correlation kernel being local, its larger impact is for the low-lying LDA continuum orbitals having a local character. As the TRK sum rule must be satisfied, the higher peak in the TDLDA spectrum is followed by a decrease of the oscillator strengths faster than in the LDA spectrum, until they reach the same asymptotic behavior.

The HF spectrum in Fig. 4a not only has no discrete photoexcitation part, as expected since the unoccupied HF orbitals are unbound (see the $\mu\to\infty$ limit of Fig. 2a), but does not even look as a photoionization spectrum. The HF unoccupied orbitals actually represent approximations to the continuum states of the H- anion, and are thus much more diffuse than the exact continuum states of the H atom, as shown in Fig. 5. Consequently, the HF spectrum has in fact the characteristic shape of the photodetachment spectrum of the H- anion Bethe and Salpeter (1957); Rau (1996) (with the caveat that the initial state is the 1s orbital of the H atom instead of the 1s orbital of the H- anion). Finally, note that, for the H atom, linear-response TDHF gives of course the exact photoexcitation/photoionization spectrum.

Let us now discuss the spectra obtained with the range-separated methods in Fig. 4b. The common value of the range-separation parameter $\mu=0.5$ bohr*-1* has been used Gerber and Ángyán (2005). The RSH spectrum looks like the photodetachment spectrum of the H- anion. This is not surprising since the RSH effective Hamiltonian contains a long-range HF exchange potential. The RSH continuum orbitals are similarly diffuse as the HF continuum orbitals, as shown in Fig. 5. The RSH ionization threshold energy is slightly smaller than the exact value (0.5 Ha) due to the remaining self-interaction error in the 1s orbital energy stemming from the short-range LDA exchange-correlation potential at this value of $\mu$. The RSH-EXX ionization threshold is identical to the RSH one, but, contrary to the RSH spectrum, the RSH-EXX spectrum correctly shows a discrete photoexcitation part and a continuum photoionization part. Beside the small redshift of the spectrum, the self-interaction error at this value of $\mu$ manifests itself in slightly too small RSH-EXX oscillator strengths. The RSH-EXX continuum orbitals are very similar to the exact continuum orbitals, as shown in Fig. 5. Finally, at this value of $\mu$, TDRSH gives a photoexcitation/photoionization spectrum essentially identical to the RSH-EXX spectrum.

IV.4 Photoexcitation and photoionization

spectra for the helium atom

In Fig. 6, different photoexcitation/photoionization spectra for the He atom are shown. As in the H atom case, the oscillator strengths of the discrete part of the TDHF, RSH-EXX, and TDRSH spectra have been interpolated (using again the oscillator strengths of first five Rydberg states and of the first continuum state) to correct the overestimation of the oscillator strengths for the last Rydberg transitions. The excitation energies and the (non-interpolated) oscillator strengths of the first five discrete transitions are reported in Table 1 and compared with exact results. The photoionization part of some of the calculated spectra are compared with full configuration-interaction (FCI) calculations and experimental results in Fig. 7.

In Fig. 6a, one sees that the HF spectrum looks again like a photodetachment spectrum, corresponding in this case to the He- anion. By contrast, TDHF gives a reasonable photoexcitation/photoionization spectrum. In particular, for the first discrete transitions listed in Table 1, TDHF gives slightly too large excitation energies by at most about 0.02 Ha (or 0.5 eV) and slightly too small oscillator strengths by at most about 0.025. The ionization energy is also slightly too large by about 0.015 Ha, as already seen from the HF 1s orbital energy in the $\mu\to\infty$ limit of Fig. 3. As regards the photoionization part of the spectrum, one sees in Fig. 7 that TDHF gives slightly too large photoionization cross sections.

The LDA spectrum in Fig. 6a is also similar to the LDA spectrum for the H atom. The ionization threshold energy is much too low, and the spectrum lacks a discrete part and has an unphysical maximum close above the ionization threshold. Except from that, taking as reference the TDHF spectrum (which is close to the exact spectrum), the LDA spectrum is a reasonable approximation to the photoionization spectrum and, again as noted in Ref. Wasserman et al., 2003, a reasonable continuous approximation to the photoexcitation spectrum. In comparison to LDA, TDLDA Zap gives smaller and less accurate oscillator strengths in the lower-energy part of the spectrum but, the TRK sum rule having to be preserved, larger oscillator strengths in the higher-energy part of the spectrum, resulting in an accurate high-energy asymptotic behavior as seen in Fig. 7.

Fig. 6b shows the spectra calculated with RSH, RSH-EXX, and TDRSH using for the range-separation parameter the value $\mu=1.115$ bohr*-1* which imposes the exact ionization energy, as explained in Sec. IV.2. The RSH spectrum is similar to the HF spectrum and does not represent a photoexcitation/photoionization spectrum. By contrast, the RSH-EXX spectra is qualitatively correct for a photoexcitation/photoionization spectrum. As shown in Table 1, in comparison with TDHF, RSH-EXX gives more accurate Rydberg excitation energies, with a largest error of about 0.003 Ha (or 0.08 eV), but less accurate oscillator strengths which are significantly overestimated. The TDRSH method also gives a correct photoexcitation/photoionization spectrum, with the advantage that it gives Rydberg excitation energies as accurate as the RSH-EXX ones and corresponding oscillator strengths as accurate as the TDHF ones. As shown in Fig. 7, TDRSH also gives a slightly more accurate photoionization cross-section profile than TDHF.

V Conclusions

We have investigated the performance of the RSH scheme for calculating photoexcitation/photoionization spectra of the H and He atoms, using a B-spline basis set in order to correctly describe the continuum part of the spectra. The study of these simple systems allowed us to quantify the influence on the spectra of the errors coming from the short-range exchange-correlation LDA and from the missing long-range correlation in the RSH scheme. For the He atom, it is possible to choose a value for the range-separation parameter $\mu$ for which these errors compensate each other so as to obtain the exact ionization energy.

We have studied the differences between using the long-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 the corresponding RSH-EXX scheme, even without applying linear-response theory, gives reasonable photoexcitation/photoionization spectra. Nevertheless, the most accurate spectra are obtained with linear-response TDRSH (or TDRSH-EXX since they are equivalent for one- and two-electron systems). In particular, for the He atom at the optimal value of $\mu$, TDRSH gives slightly more accurate photoexcitation and photoionization spectra than standard TDHF.

The present work calls for further developments. First, the merits of TDRSH (and/or TDRSH-EXX) for calculating photoexcitation/photoionization spectra of larger atoms and molecules, where screening effects are important, should now be investigated. Second, it would be interesting to test the effects of going beyond the LDA for the short-range exchange-correlation functional Toulouse et al. (2005); Goll et al. (2006) and adding long-range wave-function correlation Fromager et al. (2013); Hedegård et al. (2013); Rebolini and Toulouse (2016). Third, time-propagation TDRSH could be implemented to go beyond linear response and tackle strong-field phenomena, such as high-harmonic generation and above-threshold ionization Labeye et al. (2018).