Exact exchange-correlation kernels for optical spectra of model systems

TL;DR

This paper calculates the exact exchange-correlation kernels for two model systems in time-dependent density functional theory, revealing their frequency and spatial dependence and emphasizing the importance of the exchange-correlation potential's accuracy.

Contribution

It provides the first exact calculations of $f_{xc}$ for prototype systems, highlighting their simple form near the first excitation and complexity at higher energies.

Findings

Exact $f_{xc}$ has weak frequency dependence near first excitation

Spatial form of $f_{xc}$ is simple and non-local at low energies

Higher excitations show complex spatial and frequency behavior

Abstract

For two prototype systems, we calculate the exact exchange-correlation kernels $f_{\mathrm{xc}}(x,x',\omega)$ of time-dependent density functional theory. $f_{\mathrm{xc}}$, the key quantity for optical absorption spectra of electronic systems, is normally subject to uncontrolled approximation. We find that, up to the first excitation energy, the exact $f_{\mathrm{xc}}$ has weak frequency-dependence and a simple, though non-local, spatial form. For higher excitations, the spatial behavior and frequency-dependence become more complex. The accuracy of the underlying exchange-correlation potential is of crucial importance.

Figures

Equations

$$\sigma(\omega)=-\frac{4\pi\omega}{c}\int\int\text{Im}\big{(}\chi(x,x^{\prime},\omega)\big{)}\ x\ x^{\prime}\ dx\ dx^{\prime}.$$

$$\chi(x,x^{\prime},\omega)=\sum_{m}\bigg{[}\frac{\langle 0|\hat{n}(x)|m\rangle\langle m|\hat{n}(x^{\prime})|0\rangle}{\omega-(E_{m}-E_{0})+i\eta}+c.c.(-\omega)\bigg{]},$$

$${\chi }_0(x,x^{\prime },\omega )=\sum _{i,j}(f_i-f_j)\frac{{\varphi }_i^{\ast }(x){\varphi }_j(x){\varphi }_j^{\ast }(x^{\prime }){\varphi }_i(x^{\prime })}{\omega -({\epsilon }_j-{\epsilon }_i)+i\eta },$$

$$f_{\mathrm{xc}}={\chi }_0^{-1}-{\chi }^{-1}-u,$$

Full text

Time-dependent Kohn-Sham density functional theory Runge and Gross (1984); Gross et al. (1996) (TDDFT) is in principle an exact and efficient theory of the excited-state properties of many-electron systems, including a wide variety of important spectroscopies such as optical absorption spectra of molecules and solids. However its application is restricted by the limitations of the available approximate functionals for electron exchange and correlation – in particular the exchange-correlation kernel, $f_{\mathrm{xc}}$, the functional derivative of the exchange-correlation potential with respect to the electron density. To assist the construction of more powerful approximations for $f_{\mathrm{xc}}$, we calculate the exact $f_{\mathrm{xc}}$ for small prototype systems, and analyze its character, including key aspects in which it differs from the common approximations.

In the Runge-Gross formulation Runge and Gross (1984) of TDDFT the real system of interacting electrons is mapped onto an auxiliary system of noninteracting electrons moving in an effective local Kohn-Sham (KS) potential $v_{\mathrm{KS}}=v_{\mathrm{ext}}+v_{\mathrm{H}}+v_{\mathrm{xc}}$, with both systems having the same electron density $n$ at all points in space and time. Many TDDFT calculations are done within the framework of linear response theory, which describes how a system responds upon application of a weak, time-dependent external perturbation. The induced density is described by the interacting density-response function, the functional derivative $\chi=\delta n/\delta v_{\mathrm{ext}}$. $\chi$ is related to the non-interacting density-response function of the KS system, $\chi_{0}=\delta n/\delta v_{\mathrm{KS}}$, by the Dyson equation 111Matrix multiplication for the spatially non-local quantities $\chi$, $\chi_{0}$ and $f_{\mathrm{xc}}$, and $\omega$-dependence, are implied. Petersilka et al. (1996) $\chi=\chi_{0}+\chi_{0}(u+f_{\mathrm{xc}})\chi$, where $u$ is the bare Coulomb interaction. $\chi_{0}$ is to be obtained from a ground-state DFT calculation. $\chi$ can then be used to compute, for example, the optical absorption spectrum of the system,

[TABLE]

In practice, both $v_{\mathrm{xc}}$ and its functional derivative $f_{\mathrm{xc}}$ must be approximated. While there have been some successes, the commonly used adiabatic TDDFT functionals, such as the adiabatic local density approximation Zangwill and Soven (1980); Gross and Kohn (1985) (ALDA), often fail in extended systems. For example, the optical absorption spectra of many semiconductors and insulators are not even qualitatively described, with excitonic effects and many-electron excitations omitted Jamorski et al. (1996); Maitra et al. (2004), and the optical gap underestimated. Here, approximations for $f_{\mathrm{xc}}$ achieve little improvement over the random phase approximation (RPA), in which $f_{\mathrm{xc}}$ is neglected entirely Gavrilenko and Bechstedt (1997). Attempts to improve approximations for $f_{\mathrm{xc}}$ include exact-exchange methods Kim and Görling (2002a, b); Görling (1998); Hellgren and von Barth (2009); Ipatov et al. (2010); Görling (1999), diagrammatic expansions using perturbative methods Tokatly and Pankratov (2001); Tokatly et al. (2002), and adding long-range contributions Botti et al. (2004); Sharma et al. (2011); Trevisanutto et al. (2013); Rigamonti et al. (2015). Another approach involves calculations of the homogeneous electron gas (HEG) Nazarov et al. (2010); Qian and Vignale (2002); Richardson and Ashcroft (1994); Conti et al. (1997); Panholzer et al. (2018). Kernels derived from the Bethe-Salpeter equation Onida et al. (2002); Reining et al. (2002); Sottile et al. (2003); Adragna et al. (2003); Marini et al. (2003) have had some success, but require a relatively expensive many-body perturbation theory calculation as their input, and are outside the KS TDDFT framework.

There have been a limited number of studies conducted on analyzing the character of the exact $f_{\mathrm{xc}}$, all of which focus on its frequency-dependence. One approach has been to calculate the exact adiabatic $f_{\mathrm{xc}}$ for model systems Thiele and Kümmel (2009), in order to investigate its performance upon application and deduce when memory effects become important. This approach has been used in simple Hubbard systems Aryasetiawan and Gunnarsson (2002); Carrascal et al. (2018) and extended by analyzing additional properties, such as the frequency-dependence of the full $f_{\mathrm{xc}}$ around double excitations. Other research has explored how this frequency-dependence of $f_{\mathrm{xc}}$ turns the single-particle quantities of exact KS TDDFT into many-body excitations Thiele and Kümmel (2014) and its behaviour for long-range excitations has been analyzed in order to develop approximate kernels Maitra and Tempel (2006).

In this paper, we explore the properties of exact xc kernels, including full spatial and frequency-dependence, in order to inform the development of improved approximate functionals. We employ our iDEA code Hodgson et al. (2013) which solves the many-electron Schrödinger equation exactly for small, one-dimensional prototype systems 222We perform calculations for systems of two spinless electrons interacting via the appropriately softened Coulomb repulsion Gordon et al. (2005) $u(x,x^{\prime})=(|x-x^{\prime}|+1)^{-1}$, and work in Hartree atomic units: $m_{\mathrm{e}}=\hbar=e=4\pi\varepsilon_{0}=1$. 333See Supplemental Material at ‘Link’ for the parameters of the model systems, and details on our calculations to obtain converged results.. From the many-electron eigenstates of the system we calculate the exact $\chi$ using the Lehmann representation,

[TABLE]

where $|0\rangle$, $E_{0}$, $|m\rangle$ and $E_{m}$ are the ground state and its energy, and the $m-$th excited state and its energy, respectively, $\hat{n}$ is the density operator in the Heisenberg picture and $\eta$ is a positive infinitesimal. $\chi$ has poles at the excitation energies of the system, $E_{m}-E_{0}$. It is convenient to calculate Im$(\chi)$, with Re$(\chi)$ following from the Kramers-Kronig relations Pines (1964); Marder (2010). As is customary, we replace $\eta$ with a small positive number, to broaden the absorption peaks for ease of viewing.

We then determine the exact KS potential through our reverse-engineering algorithm Ramsden and Godby (2012). From the exact KS orbitals, we calculate the exact non-interacting density-response function,

[TABLE]

where the $\phi_{i}$, $\varepsilon_{i}$ are the exact solutions to the Kohn-Sham equations of ground-state DFT, and $f_{i}$ is the Fermi occupation (0 or 1) of $\phi_{i}$.

The Dyson equation may be manipulated into an expression for $f_{\mathrm{xc}}$,

[TABLE]

but the inverses of $\chi$ and $\chi_{0}$ are not well defined. For instance, a spatially uniform perturbation of any angular frequency induces no change in density, so both $\chi$ and $\chi_{0}$ have a zero eigenvalue and therefore a zero determinant. To overcome this, we find a pseudo-inverse of $\chi$ using truncated singular-value decomposition, discarding those eigenvectors with eigenvalues smallest in magnitude, which we term the eigenvalue cutoff. This procedure is repeated for $\chi_{0}$, discarding the same number of eigenvectors. From the modified response functions a kernel $f_{\mathrm{xc}}$ is now well defined. We confirm the validity of this procedure by verifying that the calculated $f_{\mathrm{xc}}$, together with the unmodified $\chi_{0}$, closely reproduces the unmodified $\chi$ via the Dyson equation. Additionally, we ensure that the zero-force sum rule is obeyed – a well known property of the exact $f_{\mathrm{xc}}$ 444See Supplemental Material for more details..

We begin by considering a system of two interacting electrons confined to a harmonic well potential ($\omega_{0}=0.25$ a.u.), where $\omega_{0}$ is the angular frequency of the well (inset of Fig. 1). We compute the exact optical absorption spectrum of the system 555For this harmonic well system, at the level of linear response theory, only one excitation appears in the absorption spectrum., with the first excitation at $\omega=\omega_{0}$ (Fig. 1). Additionally, we compute the absorption spectrum of the exact Kohn-Sham system, in which the absorption frequency is slightly too low ($\approx 0.01$ a.u.). We also calculate the RPA and ALDA absorption spectra, in which the RPA and ALDA Entwistle et al. (2018) kernels are combined with the exact $\chi_{0}$. This last point provides a strong reminder of the challenge of $f_{\mathrm{xc}}$: starting from the exact Kohn-Sham orbitals, a much better absorption peak is obtained by ignoring the induced changes in the Hartree and xc potentials ($\chi_{0}$) than by accounting for the first exactly and either neglecting (RPA) or approximating (ALDA) the second. This highlights the importance of obtaining a good approximation to the ground-state xc potential $v_{\mathrm{xc}}$, which leads to $\chi_{0}$.

We now turn to the spatial characteristics of $f_{\mathrm{xc}}$ (Fig. 2). Typically, several different choices of the eigenvalue cutoff yield kernels $f_{\mathrm{xc}}$ with varying degrees of spatial structure, all of which essentially yield the correct $\chi$ from the exact $\chi_{0}$ as set out above. Of these, we select the eigenvalue cutoff with the largest magnitude, resulting in the smoothest possible spatial structure without detriment to the exact absorption spectrum (Fig. 1). We observe that while $f_{\mathrm{xc}}$ has real and imaginary parts (see later), the real part alone is sufficient to reproduce the position and weight of the first excitation $(\omega=\omega_{0})$. Fig. 2(a) and Fig. 2(b) show Re$(f_{\mathrm{xc}})$ at $\omega=0$ and $\omega_{0}$, respectively. The behavior of $f_{\mathrm{xc}}$ away from the diagonal, $x\neq x^{\prime}$, represents the kernel’s non-locality, and it is evident that this is fairly simple in nature; analysis (Fig 2(c)) shows it to be similar to the negative of the Coulomb interaction, with which it therefore tends to cancel in the expression for $\chi$. The $\omega$-dependence of $f_{\mathrm{xc}}$ up to the first excitation is seen to be extremely weak, as observed in other model systems Aryasetiawan and Gunnarsson (2002); Carrascal et al. (2018); Thiele and Kümmel (2014). We analyze this more closely in the inset to Fig. 2(c).

To gain insight into these observations, we analyze the exact $\chi$ and $\chi_{0}$. Fig. 3 shows Re$(\chi)$ and Re$(\chi_{0})$; up to the first excitation, these exhibit strong – but closely similar – $\omega$-dependence. The similarity arises in part from the exact many-electron wavefunction being well approximated by the exact Kohn-Sham wavefunction 666We define this as a Slater determinant of the occupied KS orbitals., which reflects the dominance of exchange (including self-interaction correction) in the harmonic potential system 777E. Richardson, private communication. Therefore $\chi^{-1}$ and $\chi_{0}^{-1}$ largely cancel, so that $f_{\mathrm{xc}}$ is similar to $-u$, with weak $\omega$-dependence.

This can be demonstrated succinctly through a simple model, in which we take Eq. (2) and Eq. (3), and replace the spatially-dependent numerators (the oscillator strengths) with scalars. Specifically, we consider a system with a single excitation at $\omega=1$, set the numerator equal to 1, and let $\eta=0.05$ (top inset of Fig. 4). We do the same for the Kohn-Sham system, but choose the excitation to occur at $\omega=0.9$. By taking their inverses, we calculate Re$(\chi_{0}^{-1}-\chi^{-1})$, which is the $\omega$-dependent part of Re$(f_{\mathrm{xc}})$ in Eq. (4), and find this to be small in amplitude and have a fairly weak $\omega$-dependence up to the first excitations (bottom inset of Fig. 4). The inclusion of higher excitations, and taking the limit $\eta\rightarrow 0$, change little at these low frequencies.

Including higher excitations in the model $\chi$ causes Re$(\chi)$ to pass through zero between excitations. At these points Re$(\chi^{-1})$ also passes through zero, and Im$(\chi^{-1})$ peaks sharply (Fig. 4). As Im$(f_{\mathrm{xc}})$ = Im$(\chi_{0}^{-1}-\chi^{-1})$, we find that the $f_{\mathrm{xc}}$ in our simple model only has an imaginary component when $\chi$ or $\chi_{0}$ passes through zero between excitations, and hence is completely real up to the first excitations (as $\eta\rightarrow 0$). This supports our finding in the harmonic well system, in which Im$(f_{\mathrm{xc}})$ was very small up to the first excitations, and Re$(f_{\mathrm{xc}})$ was sufficient to reproduce the peak in the absorption spectrum.

We now consider a system whose absorption spectrum includes higher excitations – two interacting electrons in a softened atomic-like potential (top inset of Fig. 5). As in the harmonic well system, the absorption spectrum of the exact Kohn-Sham system is slightly too low for the first excitation (Fig. 5). Again, we find $f_{\mathrm{xc}}$ to be dominated by its real part and nearly $\omega$-independent, while exhibiting relatively simple spatial structure, up to and including the first excitation (Fig. 6(a)). The second excitation does not appear in the absorption spectrum, and so we move to the third excitation, which is much smaller in amplitude than the first, and once more observe the peak in the Kohn-Sham system to be slightly below but still very close to the exact (bottom inset of Fig. 5). Again, the closeness between the two peaks arises from the strong similarity between $\chi$ and $\chi_{0}$. In order to reproduce this excitation, a smaller eigenvalue cutoff is required, leading to higher spatial frequencies in $f_{\mathrm{xc}}$ 888As expected for higher energy excited states. (Fig. 6(b)).

For the atom-like system, we have investigated the extent to which local kernel approximations for $f_{\mathrm{xc}}$ may be meaningful. As we have observed $f_{\mathrm{xc}}$ to largely cancel with $u$ at low $\omega$, we choose to focus on the Hartree exchange-correlation kernel $f_{\mathrm{Hxc}}=f_{\mathrm{xc}}+u$ which is more local. We incorporate the nonlocal parts of $f_{\mathrm{Hxc}}$ by projecting them onto a local kernel 999We fold the nonlocal $f_{\mathrm{Hxc}}$ with an envelope function that suppresses the more distant nonlocal parts and projects the remainder onto the diagonal $x=x^{\prime}$.. We find this largely corrects the difference $\chi_{0}-\chi$, and hence the position of the peak in the absorption spectrum, for the first excitation, but fails to correct the height of the peak. Such a local kernel is completely inadequate to describe the third excitation.

In summary, we have calculated the exact $f_{\mathrm{xc}}(x,x^{\prime},\omega)$ for two prototype systems. At low $\omega$, we find the imaginary component of $f_{\mathrm{xc}}$ to be small, with the real part alone sufficient to reproduce the first excitation. Up to and including the first excitation, Re$(f_{\mathrm{xc}})$ exhibits strikingly weak $\omega$-dependence, stemming from strong, but closely similar $\omega$-dependence between the interacting and non-interacting density-response functions – boding well for the applicability of adiabatic kernels. Additionally, Re$(f_{\mathrm{xc}})$ here has a rather simple spatial form, which is similar to the negative of the Coulomb interaction $u$, indicating that approximations to $f_{\mathrm{Hxc}}$ may be more appropriate than those for $f_{\mathrm{xc}}$ alone. For higher excitations, $f_{\mathrm{xc}}$ exhibits both additional spatial structure and stronger $\omega$-dependence, indicating that more sophisticated approximations are needed. Throughout, the absorption spectrum of the exact Kohn-Sham system provides a very good starting point, signifying the crucial importance of an accurate approximation for the ground-state $v_{\mathrm{xc}}$.

We thank Phil Hasnip for helpful discussions. Data created during this research is available from the York Research Database 101010M. T. Entwistle and R. W. Godby, http://dx.doi.org/10.15124/56f576ee-b2de-4ca5-9251-831bfc3cae6f (2019)..