## Accurate optical spectra of solids from pure time-dependent density-functional theory

Sarah Cavo, 1 J. A. Berger, 2 and Pina Romaniello 1

1 Laboratoire de Physique Th´ eorique, Universit´ e de Toulouse, CNRS, UPS and and European Theoretical Spectroscopy Facility, France 2 Laboratoire de Chimie et Physique Quantiques, Universit´ e de Toulouse, CNRS, UPS and and European Theoretical Spectroscopy Facility, France (Dated: February 21, 2020)

We present accurate optical spectra of semiconductors and insulators within a pure Kohn-Sham time-dependent density-functional approach. In particular, we show that the onset of the absorption is well reproduced when comparing to experiment. No empirical information nor a theory beyond Kohn-Sham density-functional theory, such as GW , is invoked to correct the Kohn-Sham gap. Our approach relies on the link between the exchange-correlation kernel of time-dependent density functional theory and the derivative discontinuity of ground-state density-functional theory. We show explicitly how to relate these two quantities. We illustrate the accuracy and simplicity of our approach by applying it to various semiconductors (Si, GaP, GaAs) and wide-gap insulators (C, LiF, Ar).

Time-dependent-density functional theory (TDDFT) [1] has become, over the years, one of the few well-established first-principles' approaches to describe time-dependent phenomena for a large variety of systems, both in the linear-response regime and beyond (see, e.g., Refs [2, 3] and references therein). In the last two decades TDDFT has been increasingly applied to solids, and in particular to the calculation of the optical absorption spectra. Optical experiments in general are very useful tools to investigate and characterize condensed-matter systems; it is hence desirable to develop efficient and reliable theoretical approaches to complement experiment.

Sham valence and conduction bands are available. Indeed, optical spectra obtained with all currently available xc kernels show absorption at energies close to the Kohn-Sham band gap, thereby severely underestimating the absorption onset.

Within TDDFT the description of optical spectra depends crucially on the exchange-correlation (xc) kernel f xc which relates the response of the Kohn-Sham system to a small perturbation to the response of the true system. Thanks to the numerical efficiency of TDDFT, it is desirable that a simple but accurate xc kernel is available for the calculation of optical spectra. It is well-known that traditional xc kernels, such as the random-phase approximation (RPA), i.e., f xc = 0, and the adiabatic local-density approximation (ALDA) [4], fail to describe two important features of optical spectra: 1) excitonic effects and 2) the absorption onset. While excitonic effects can nowadays be described accurately for various systems with a several xc kernels [5-10], the correct description of the absorption onset within TDDFT remains an unsolved problem.

The starting point for a TDDFT calculation is the Kohn-Sham band structure. As is well known, the KohnSham band gap is, in general, estimated to be much smaller than the fundamental gap, i.e., the difference between the ionization potential and the electron affinity. [11-13]. Therefore, since TDDFT should give the exact absorption spectra, the TDDFT xc kernel has the difficult task to ensure that there is no absorption below the optical gap even though transitions between Kohn-

The standard approach to circumvent this problem is to add a scissors operator [14] to the Kohn-Sham Hamiltonian. The shift parameter is either obtained from experiment or from a method that goes beyond KS-DFT, such as GW [15-19], or generalised KS-DFT with a hybrid functional, such as those based on a screened Coulomb interaction [20-23]. The former approach is unsatisfactory because it is empirical while the latter approach is unsatisfactory both conceptually and numerically because i) one has to go to a theory outside of the KS-DFT framework, i.e., GW (or a hybrid functional), in order to obtain spectra that are comparable with those observed in experiment, and ii) a GW (or a hybrid functional) calculation is much more expensive than a pure KS-DFT calculation. Eliminating the intermediate GW step yields a fully coherent theory. Moreover, it leads to a large speed up of calculations, making the optical spectra of larger systems accessible. We note that one could use an ad hoc meta-GGA functional specifically constructed to obtain a KS gap similar to the fundamental gap. However, this is not in agreement with fundamental DFT theory because the KS gap is not equal to the fundamental gap, the difference between the two being the derivative discontinuity.

The TDDFT xc kernel f xc can be written exactly as [24]

$$\mathfrak { h } { n } { - \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathf { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c r } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c r } \int _ { \mathfrak { c r } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c r } \int _ { \mathfrak { c c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c r } \int _ { \mathfrak {c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } { \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \ int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \chr { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { | \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \comfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfra \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int _ { \mathfrak { c } \int$$

where χ KS and χ 0 = -iGG are the Kohn-Sham and independent quasiparticle polarizability, respectively, and G (1 , 2) and Σ(1 2) , are the one-body Green function and the self-energy, respectively. The collective index (1) = ( x , t ) = ( r , s, t ) contains the space, spin and time coordinates. The xc kernel written in Eq. (1) clearly exhibits two distinct parts. The first part f (1) xc = χ -1 KS -χ -1 0 only involves independent (quasi-)particles, and, therefore, is responsible for the shift of the Kohn-Sham band gap to the fundamental gap, while the second part, f (2) xc , which includes the electron-hole interaction, accounts for the excitonic effects. In general, both terms are required to guarantee a correct onset of the absorption, unless the exciton binding energy is small, in which case f (1) xc is sufficient. Although Eq. (1) clearly distinguishes these two parts it is not useful in practical applications since it would require the calculation of G . From the above discussion one would expect a link between f (1) xc and the derivative discontinuity of ground-state DFT [25, 26], which is defined as the difference between the fundamental gap and the Kohn-Sham gap. One of the goals of this work is to make this link explicit.

In order to obtain this formal link we first generalize the two-point KS polarizability to four points: [27]

$$^ { 4 } \chi _ { \text{KS} } ( \text{x} _ { 1 }, \text{x} _ { 2 }, \text{x} _ { 3 }, \text{x} _ { 4 }, \omega ) & = \sum _ { i, j } ( f _ { j } - f _ { i } ) \\ & \times \frac { \phi _ { i } ( \text{x} _ { 1 } ) \phi _ { j } ( \text{x} _ { 2 } ) \phi _ { j } ^ { * } ( \text{x} _ { 3 } ) \phi _ { i } ^ { * } ( \text{x} _ { 4 } ) } { \omega - ( \epsilon _ { i } - \epsilon _ { j } ) + i \eta }, \quad ^ { 4 } \chi _ { 0 }$$

/negationslash where φ i is a KS spinorbital, /epsilon1 i its energy, f i its occupation (0 and 1 for unoccupied and occupied orbitals, respectively), and η is a positive infinitesimal that ensures causality. For the solids we study here i and a are multi-indices composed of a band index (comprising the spin) and a Bloch vector, k and k ′ , respectively. We note that, although the final goal is the description of optical absorption for which k ′ → k , the discussion below is completely general, i.e., k = k ′ . The usual two-point KS polarizability is retrieved from the diagonal part of χ (4) KS as

$$\chi _ { K S } ( x _ { 1 }, x _ { 2 }, \omega ) = { ^ { 4 } } \chi _ { K S } ( x _ { 1 }, x _ { 2 }, x _ { 1 }, x _ { 2 }, \omega ). \quad ( 3 ) \quad. \.$$

We can express χ (4) KS in the Kohn-Sham basis by using the following basis transformation

$$^ { 4 } \chi _ { K S [ n _ { 1 } n _ { 3 } ] } ^ { [ n _ { 4 } n _ { 2 } ] } ( \omega ) & = \int d { \mathbf x } _ { 1 } { \mathbf x } _ { 2 } { \mathbf x } _ { 3 } { \mathbf x } _ { 4 } ^ { 4 } \chi _ { K S } ( { \mathbf x } _ { 1 }, { \mathbf x } _ { 2 }, { \mathbf x } _ { 3 }, { \mathbf x } _ { 4 }, \omega ) \quad \text{$\quad$Th} \, \\ & \times \phi _ { n _ { 1 } } ^ { * } ( { \mathbf x } _ { 1 } ) \phi _ { n _ { 2 } } ^ { * } ( { \mathbf x } _ { 2 } ) \phi _ { n _ { 3 } } ( { \mathbf x } _ { 3 } ) \phi _ { n _ { 4 } } ( { \mathbf x } _ { 4 } ). \quad ( 4 ) \quad \text{matrix} \, \quad \, 4 f _ { x c, 0 } ^ { ( 1 ) }$$

This yields a 2 M × 2 M diagonal matrix with M the num- ber of KS excitations. It is schematically given by [28]

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00808v3-figure-b373ee4c-bcbc-471b-b9c2-7c5ca0d5aff4.jpg)

where ω i is a KS excitation energy, i.e., a pole of χ KS ( ω ), and the matrix elements are arranged in order of increasing excitation energy, i.e., ω 1 ≤ ω 2 ≤ ω 3 , etc. Note that the χ KS ( ω ) matrix representation has two blocks, one representing the resonant part (excitation energies) and the other one the antiresonant part (de-excitation energies). In particular, the lowest KS excitation energy is given by ω 1 = /epsilon1 CBM -/epsilon1 VBM , with /epsilon1 CBM and /epsilon1 VBM the KS energy of the conduction band minimum (CBM) and the KS energy of the valence band maximum (VBM), respectively.

We now assume that also χ 0 is diagonal in the KohnSham basis. This is an approximation, but it is in accordance with various practical calculations, in particular those based on the GW method, in which χ 0 is built with KS orbitals [17, 29]. We can thus write

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00808v3-figure-a509b24d-c450-4db3-907b-8f9943fe454e.jpg)

where Ω i are quasiparticle energy differences, i.e., differences of ionization potentials and electron affinities. In particular, the lowest excitation energy is given by Ω 1 = I -A , where I is the first ionization potential and A is the first electron affinity. Here we assume that the lowest excitation energy Ω 1 corresponds to the same matrix element as the lowest KS excitation energy. The fourpoint kernel 4 f (1) xc = 4 χ -1 KS -4 χ -1 0 has hence the following simple frequency-independent matrix representation

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00808v3-figure-929009f0-59dc-4333-b62f-91891d473131.jpg)

The absorption onset is determined by the head of the matrix 4 f (1) xc , which we will refer to in the following as 4 f (1) xc, 00 . It is given by

$$\ ^ { 4 } f _ { x c, 0 0 } ^ { ( 1 ) } = I - A - ( \epsilon _ { \text{CBM} } - \epsilon _ { \text{VBM} } ) \quad \ \ ( 8 )$$

It can be shown that the ionization potential is exactly equal to minus the KS energy at the VBM, i.e., I = -/epsilon1 VBM [30-32]. Since an equivalent relation holds for the N + 1 system, i.e., the system with one additional electron, and the fact that A should be equal to the first ionization potential of the N +1 system, one can deduce that A = -/epsilon1 N +1 VBM , where /epsilon1 N +1 VBM is the KS energy at the VBM of the N + 1 system. Therefore, we can rewrite Eq. (8) as

$$\L ^ { 4 } f _ { x c, 0 0 } ^ { ( 1 ) } = \epsilon _ { \text{VBM} } ^ { N + 1 } - \epsilon _ { \text{CBM} }. \quad \quad ( 9 ) \quad \begin{matrix} 3 2 \\ \text{sta} \end{matrix}$$

Although Eq. (9) seems a simple expression, it is not easy to calculate in practice. The problem arises from the fact that /epsilon1 N +1 VBM is difficult to evaluate in solids since they are usually described within the thermodynamic limit, which implies an infinite number of electrons from the outset. However, /epsilon1 N +1 VBM -/epsilon1 CBM is equal to ∆, the difference between the fundamental gap E g and the KS gap E KS [26, 32, 33]. This difference is also known as the derivative discontinuity [25, 26]. Therefore we arrive at the following relation,

$$\ ^ { 4 } f _ { x c, 0 0 } ^ { ( 1 ) } = \Delta. \quad \ \ ( 1 0 ) \quad \ \ a _ { F }$$

This is one of the main results of this work.

Assuming a rigid shift of the conduction bands, 4 f (1) xc can thus be approximated by 4 f (1) xc = ∆ I, where 4 4 I is the four-point identity matrix in transition space. In practice, however, we prefer to use two-point quantities for numerical efficiency. Therefore, it might be tempting to rotate 4 f (1) xc to a two-point kernel. However, this would lead to a very complicated frequency-dependent quantity and would hence be very difficult to apply in practice [34]. We will avoid this problem by including the effect of 4 f (1) xc on the four-point response function in transition space and only afterwards rotating the latter to its twopoint real-space representation.

We thus introduce a modified four-point Kohn-Sham polarizability 4 χ (1) KS ( ω ) defined by

$$[ ^ { 4 } \chi ^ { ( 1 ) } _ { K S } ] ^ { - 1 } ( \omega ) = [ ^ { 4 } \chi _ { K S } ] ^ { - 1 } ( \omega ) - \Delta ^ { 4 } I, \quad ( 1 1 ) \quad \text{$\ u c a $} \\ \text{$\ o r $} \, \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \text{$\ } \end{$$

where the superscript (1) indicates that χ (1) KS ( ω ) contains ∆. Using the inverse of the basis set transformation in Eq. (4), and taking the diagonal part of the resulting expression we obtain the following expression for the modified two-point KS polarizability,

$$\chi ^ { ( 1 ) } _ { K S } ( \mathbf x _ { 1 }, \mathbf x _ { 2 }, \omega ) = \sum _ { i, j } \frac { ( f _ { j } - f _ { i } ) \phi _ { i } ( \mathbf x _ { 1 } ) \phi _ { j } ( \mathbf x _ { 2 } ) \phi _ { j } ^ { * } ( \mathbf x _ { 1 } ) \phi _ { i } ^ { * } ( \mathbf x _ { 2 } ) } { \omega - ( \epsilon _ { i } - \epsilon _ { j } ) - s \mathbf g n ( \epsilon _ { i } - \epsilon _ { j } ) \Delta + i \eta }, \quad & O ( N ^ { 6 } ) \\ ( 1 2 ) \quad & \text{$rowin} \\ \text{with each can be odd, and in practice, to $two$} \quad \text{$1\mathcal{ }C\mathcal{ }T$}$$

which can be easily applied in practice. The true response function can then be written in terms of χ (1) KS as

$$\chi ( \omega ) = \chi _ { K S } ^ { ( 1 ) } ( \omega ) + \chi _ { K S } ^ { ( 1 ) } ( \omega ) \left [ v _ { c } + f _ { x c } ^ { ( 2 ) } ( \omega ) \right ] \chi ( \omega ), \quad ( 1 3 ) \quad \text{ample}$$

with v c the Coulomb potential. From χ ω ( ) one can readily obtain the inverse dielectric function ε -1 ( ω ) = 1 + v χ ω c ( ). The optical spectra are obtained from the imaginary part of /epsilon1 M ( ω ), the macroscopic part of /epsilon1 ( ω ):

$$\varepsilon ^ { M } ( \omega ) = \varepsilon _ { 1 } ( \omega ) + i \varepsilon _ { 2 } ( \omega ). \text{ \quad \ \ } ( 1 4 )$$

In order to apply our approach in practice we have to use an approximation for ∆. As first proposed by Kuisma et al. [35] and further discussed by Baerends [32], ∆ can be approximated in terms of simple groundstate KS-DFT quantities according to

$$\Delta & = K _ { x } \sum _ { i = 1 } ^ { N } \left [ \sqrt { \epsilon _ { \text{CBM} } - \epsilon _ { i } } - \sqrt { \epsilon _ { \text{VBM} } - \epsilon _ { i } } \right ] \\ & \quad \times \langle \phi _ { \text{CBM} } | \frac { | \phi _ { i } | ^ { 2 } } { \rho _ { 0 } } | \phi _ { \text{CBM} } \rangle \quad \quad ( 1 5 ) \\. & \quad \cdot \quad. \quad \cdot \quad. \quad. \quad. \quad. \quad.$$

where φ CBM is the KS spinorbital corresponding to the CBM, ρ 0 is the ground-state density and K x = 8 √ 2 3 π 2 ≈ 0 382. .

The expression in Eq. (15) can be obtained from the GLLB (Gritsenko-van Leeuwen-van Lenthe-Baerends) approximation to the ground-state xc potential derived in Ref. [36]. The GLLB functional is an approximation to the response part of the exact exchange optimized effective potential. A detailed derivation of Eq. (15) is given by Baerends [32]. The constant K x can be obtained from the uniform electron gas, where the GLLB exchange response potential becomes exact. [62] Fundamental gaps calculated using the derivative discontinuity in Eq. (15) have been reported for a large number of solids. [32, 35, 37-39]. In general, the results are excellent. We note, however, that for gapped materials for which the KS band gap is zero, the correction in Eq. (15) is zero as well. One finds the same problem when calculating a G W 0 0 band gap on top of a KS metallic band structure.

In practice, we use a slight generalization of TDDFT, namely TD-current-DFT (TDCDFT) [40-43]. The practical details of how we solve the KS equations within TDCDFT can be found elsewhere [44-49]. In particular our method works in real space with two-point quantities and we calculate the macroscopic polarization induced by a macroscopic electric field in terms of the macroscopic current. This allows us to have optical spectra with the correct optical gap and excitonic effects at the cost of an RPA calculation, i.e., O N ( 3 ) scaling instead of the O N ( 6 ) for methods working in transition space. We approximate f (2) xc with the polarization functional (PF) of Ref. [8] which accurately describes the excitonic effects in various systems. We will refer to the full kernel, i.e., f (1 ,GLLB ) xc + f (2 ,PF ) xc , as the Pure kernel to highlight the fact that it is based on pure Kohn-Sham theory. We implemented our approach in a modified version of the Amsterdam Density Functional (ADF) code [53-55]. We

TABLE I: GLLB-SC values for ∆ and fundamental gap E g

| system   |   ∆ [eV] |   E g [eV] |
|----------|----------|------------|
| Si       |     0.38 |      1.12  |
| GaP      |     0.76 |      2.48  |
| GaAs     |     0.27 |      1.03  |
| C        |     1.31 |      5.449 |
| LiF      |     4.11 |     14.91  |
| Ar       |     4.63 |     14.83  |

illustrate our approach by applying it to the calculation of the optical spectra of two very different classes of solids, standard semiconductors (Si, GaP, GaAs) and wide-gap insulators (C, Ar, LiF). We use the following lattice parameters: 5.43 ˚ for Si, 5.42 ˚ for GaP, 5.654 ˚ for GaAs, A A A 3.534 ˚ for diamond, 4.026 ˚ for LiF, and 5.26 ˚ for solid A A A argon. Moreover, we use the TZ2P (tripleζ + 2 polarization functions) and QZ4P (quadrupleζ + 4 polarization functions) basis sets provided by ADF for bulk silicon and LiF and for GaP, GaAs, diamond, and solid argon, respectively. The k -space integrals are done analytically using a Lehmann-Taut tetrahedron scheme [56]. The ground-state calculations are done within the GLLB-SC xc potential [32, 35, 36], which is based on the PBEsol [57] correlation potential and uses the GLLB approximation to the exchange optimized effective potential. For GaP and GaAs we also include scalar relativistic effects. The GLLB-SC values we obtained for ∆ and the fundamental gap are reported in Table I. We note that the approximation f (2 ,PF ) xc induces a gap renormalization [24, 58]; this is clear from the example of solid argon for which we have an optical gap of about 11.6 eV, which is about 3 eV smaller than the (direct) fundamental gap.

In Figs. 1, 2, and 3 we report the dielectric function of bulk Si, GaP, and GaAs, respectively, calculated with the Pure functional at 0 Kelvin and compare it to the RPA spectrum as well as to the experimental spectrum obtained at 30 Kelvin for Si, 15 Kelvin for GaP, and 22 Kelvin for GaAs.

The Pure kernel yields an absorption spectrum that is in good agreement with the experimental measurements for the three systems. In particular, the peak positions and the excitonic effects are well reproduced. Also, the real part of the dielectric function obtained with the Pure kernel compares well to the experiment. Instead, the RPA spectrum exhibits the well-known shortcomings mentioned before, i.e., the underestimation of the absorption onset and the absence of excitonic effects. In case of GaAs the RPA spectra are already in overall good agreement with experiment, and the Pure kernel only slightly improves the exciton position and intensity. Note that the splitting of the first peak in the experimental spectrum of GaAs is not reproduced in the calculated spectrum, because it is due to spin-orbit coupling, which is

FIG. 1: (Color online) Real ( /epsilon1 1 ( ω )) and imaginary ( /epsilon1 2 ( ω )) parts of the dielectric function of bulk silicon. Solid line (black): Pure functional (Pure); Dashed line (red): RPA; Dotted line (blue): experiment from Ref. [50].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00808v3-figure-d0993728-3b1e-4400-a908-8396fc04aa91.jpg)

neglected in our calculations. We note that the theoretical spectra have more structure than the experimental spectrum because it is calculated at 0 Kelvin and no broadening parameter is used to simulate temperature effects. In Figs. 4, 5, and 6 we show the dielectric function of diamond, LiF, and solid argon, respectively, calculated with the Pure kernel at 0 Kelvin and compare it to the RPA spectrum as well as to the experimental spectrum obtained at room temperature. We see that the onset of the absorption is well-reproduced as is the full spectrum, except for an overestimation of the intensity of the bound exciton in LiF and solid Ar. We note that this overestimation is common to similar kernels derived to describe excitonic effetcs [5, 7]. Unfortunately, to the best of our knowledge, there is no experimental data of /epsilon1 1 ( ω ) of solid Ar.

In conclusion, we have made explicit the link between the derivative discontinuity of ground-state DFT and the xc kernel of TDDFT. Using this link we proposed the Pure kernel, which combines the derivative discontinuity and the polarization functional, to describe optical spectra. We showed that it yields optical spectra in good agreement with experiment for typical examples of standard semiconductors and wide-gap insulators. The central issue here is that these results were obtained within a pure KS approach without resorting to empirical data or approaches that go beyond TDDFT. Finally we note that the kernel we propose in this work is an approximation to the TDDFT kernel f xc , it hence can in principle be used

FIG. 2: (Color online) Real ( /epsilon1 1 ( ω )) and imaginary ( /epsilon1 2 ( ω )) parts of the dielectric function of bulk GaP. Solid line (black): Pure functional (Pure); Dashed line (red): RPA; Dotted line (blue): experiment from Ref. [51].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00808v3-figure-41b9d263-0c6f-4e53-bece-423044a3972e.jpg)

to calculate all the properties TDDFT can access. In particular, it could be useful in total energy calculations.

This work has been supported through the EUR grant NanoX ANR-17-EURE-0009 in the framework of the 'Programme des Investissements d'Avenir' and by ANR (project no. ANR-18-CE30-0025-01).

- [1] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52 , 997 (1984).
- [2] C. Ullrich, Time-Dependent Density-Functional Theory: Concepts and Applications (Oxford University Press, 2012).
- [3] N. T. Maitra, J. Chem. Phys. 144 , 220901 (2016).
- [4] A. Zangwill and P. Soven, Phys. Rev. A 21 , 1561 (1980).
- [5] S. Sharma, J. K. Dewhurst, A. Sanna, and E. K. U. Gross, Phys. Rev. Lett. 107 , 186401 (2011).
- [6] P. E. Trevisanutto, A. Terentjevs, L. A. Constantin, V. Olevano, and F. Della Sala, Phys. Rev. B 87 , 205143 (2013).
- [7] S. Rigamonti, S. Botti, V. Veniard, C. Draxl, L. Reining, and F. Sottile, Phys. Rev. Lett. 114 , 146402 (2015).
- [8] J. A. Berger, Phys. Rev. Lett. 115 , 137402 (2015).
- [9] Z.-h. Yang, F. Sottile, and C. A. Ullrich, Phys. Rev. B 92 , 035202 (2015).
- [10] A. V. Terentjev, L. A. Constantin, and J. M. Pitarke, Phys. Rev. B 98 , 085123 (2018).
- [11] R. W. Godby, M. Schl¨ter, and L. J. Sham, Phys. Rev. u Lett. 56 , 2415 (1986).
- [12] R. W. Godby, M. Schl¨ter, and L. J. Sham, Phys. Rev. u B 36 , 6497 (1987).
- [13] M. Gr¨ uning, A. Marini, and A. Rubio, J. Chem. Phys. 124 , 154108 (2006).
- [14] Z. H. Levine and D. C. Allan, Phys. Rev. Lett. 63 , 1719 (1989).
- [15] L. Hedin, Phys. Rev. 139 , A796 (1965).
- [16] F. Aryasetiawan and O. Gunnarsson, Rep. Prog. Phys. 61 , 237 (1998).
- [17] W. G. Aulbur, L. J¨nsson, and J. W. Wilkins, in o Solid State Physics (Academic, New York, 2000), vol. 54, p. 1.
- [18] L. Reining, Wiley Interdisciplinary Reviews: Computational Molecular Science 8 , e1344 (2018).
- [19] D. Golze, M. Dvorak, and P. Rinke, Frontiers in Chemistry 7 , 377 (2019).
- [20] J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118 , 8207 (2003).
- [21] J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 124 , 219906 (2006).
- [22] Y.-i. Matsushita, K. Nakamura, and A. Oshiyama, Phys. Rev. B 84 , 075205 (2011).
- [23] J. H. Skone, M. Govoni, and G. Galli, Phys. Rev. B 93 , 235106 (2016).
- [24] F. Sottile, V. Olevano, and L. Reining, Phys. Rev. Lett. 91 , 056402 (2003).
- [25] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Phys. Rev. Lett. 49 , 1691 (1982).
- [26] J. P. Perdew and M. Levy, Phys. Rev. Lett. 51 , 1884 (1983).
- [27] G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. 74 , 601 (2002).
- [28] M. Rohlfing and S. G. Louie, Phys. Rev. B 62 , 4927 (2000).

FIG. 3: (Color online) Real ( /epsilon1 1 ( ω )) and imaginary ( /epsilon1 2 ( ω )) parts of the dielectric function of bulk GaAs. Solid line (black): Pure functional (Pure); Dashed line (red): RPA; Dotted line (blue): experiment from Ref. [52].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00808v3-figure-7fc68426-1abe-4f49-9fbd-bd6bf1ad16ec.jpg)

FIG. 4: (Color online) Real ( /epsilon1 1 ( ω )) and imaginary ( /epsilon1 2 ( ω )) parts of the dielectric function of diamond. Solid line (black): Pure functional (Pure); Dashed line (red): RPA; Dotted line (blue): experiment from Ref. [59].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00808v3-figure-6c5f6a8f-835f-408e-aa87-5822511634fb.jpg)

- [29] R. M. Martin, L. Reining, and D. M. Ceperley, Interacting Electrons: Theory and Computational Approaches (Cambridge University Press, 2016).
- [30] M. Levy, J. P. Perdew, and V. Sahni, Phys. Rev. A 30 , 2745 (1984).
- [31] C.-O. Almbladh and U. von Barth, Phys. Rev. B 31 , 3231 (1985).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00808v3-figure-c67884bc-c152-4de8-b32c-5ed1507a9962.jpg)

ε (ω)

FIG. 5: (Color online) Real ( /epsilon1 1 ( ω )) and imaginary ( /epsilon1 2 ( ω )) parts of the dielectric function of bulk LiF. Solid line (black): Pure functional (Pure); Dashed line (red): RPA; Dotted line (blue): experiment from Ref. [60].

- [44] F. Kootstra, P. L. de Boeij, and J. G. Snijders, J. Chem. Phys. 112 , 6517 (2000).
- [45] J. A. Berger, P. L. de Boeij, and R. van Leeuwen, Phys. Rev. B 71 , 155104 (2005).
- [46] P. Romaniello and P. L. de Boeij, Phys. Rev. B 71 , 155108 (2005).
- [32] E. J. Baerends, Phys. Chem. Chem. Phys. 19 , 15639 (2017).
- [47] J. A. Berger, P. L. de Boeij, and R. van Leeuwen, Phys. Rev. B 75 , 035116 (2007).
- [48] J. A. Berger, P. Romaniello, R. van Leeuwen, and P. L. de Boeij, Phys. Rev. B 74 , 245117 (2006).
- [33] L. J. Sham and M. Schl¨ter, u Phys. Rev. Lett. 51 , 1888 (1983), URL https://link.aps.org/doi/10.1103/PhysRevLett.51.1888 . [49]
- [34]
- I. V. Tokatly and O. Pankratov, Phys. Rev. Lett. 86 , 2078 (2001), URL
- P. Romaniello and P. L. de Boeij, The Journal of Chemical Physics 127 , 174111 (2007), https://doi.org/10.1063/1.2780146, URL https://doi.org/10.1063/1.2780146 .
- https://link.aps.org/doi/10.1103/PhysRevLett.86.2078 .50] [35] M. Kuisma, J. Ojanen, J. Enkovaara, and T. T. Rantala, Phys. Rev. B 82 , 115106 (2010). [ [51]
- [36] O. Gritsenko, R. van Leeuwen, E. van Lenthe, and E. J. Baerends, Phys. Rev. A 51 , 1944 (1995).
- [37] I. E. Castelli, T. Olsen, S. Datta, D. D. Landis, S. Dahl, K. S. Thygesen, and K. W. Jacobsen, Energy Environ. Sci. 5 , 5814 (2012).
- [38] F. H¨ser, T. Olsen, and K. S. Thygesen, Phys. Rev. B u 87 , 235132 (2013).
- [39] F. A. Rasmussen and K. S. Thygesen, J. Phys. Chem. C 119 , 13169 (2015).
- [40] A. K. Dhara and S. K. Ghosh, Phys. Rev. A 35 , 442 (1987).
- [41] S. K. Ghosh and A. K. Dhara, Phys. Rev. A 38 , 1149 (1988).
- [42] G. Vignale, Phys. Rev. B 70 , 201102(R) (2004).
- [43] D. Sangalli, J. A. Berger, C. Attaccalite, M. Gr¨ uning, and P. Romaniello, Phys. Rev. B 95 , 155203 (2017).
- P. Lautenschlager, M. Garriga, L. Vina, and M. Cardona, Phys. Rev. B 36 , 4821 (1987).
- S. Zollner, M. Garriga, J. Kircher, J. Huml´ ıˇek, c M. Cardona, and G. Neuhold, Phys. Rev. B 48 , 7915 (1993), URL https://link.aps.org/doi/10.1103/PhysRevB.48.7915
- [52] P. Lautenschlager, M. Garriga, S. Logothetidis, and M. Cardona, Phys. Rev. B 35 , 9174 (1987), URL https://link.aps.org/doi/10.1103/PhysRevB.35.9174
- [53] G. te Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. A. van Gisbergen, J. G. Snijders, and T. Ziegler, J. Comp. Chem. 22 , 931 (2001).
- [54] C. Fonseca Guerra, J. G. Snijders, G. te Velde, and E. J. Baerends, Theor. Chem. Acc. 99 , 391 (1998).
- [55] ADF2018 (modified version), SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands, www.scm.com.
- [56] G. Lehmann and M. Taut, Phys. Status Solidi B 54 , 469 (1972).

.

.

ε (ω)

FIG. 6: Color online) Real ( /epsilon1 1 ( ω )) and imaginary ( /epsilon1 2 ( ω )) parts of the dielectric function of solid argon. Solid line (black): pure functional (Pure); Dashed line (red): RPA; Dotted line (blue): experiment from Ref. [61].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00808v3-figure-58e7fb9f-f06c-4ac1-a758-f8e3cefd811b.jpg)

- [57] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 100 , 136406 (2008).
- [58] A. Marini, R. Del Sole, and A. Rubio, Phys. Rev. Lett. 91 , 256402 (2003), URL https://link.aps.org/doi/10.1103/PhysRevLett.91.256402
- [59] H. R. Phillip and E. A. Taft, Phys. Rev. 136 , A1445 (1964), URL https://link.aps.org/doi/10.1103/PhysRev.136.A1445

.

- [60] D. M. Roessler and W. C. Walker, J. Opt. Soc. Am. 57 , 835 (1967), URL
- http://www.osapublishing.org/abstract.cfm?URI=josa-57-6 [61] V. Saile, M. Skibowski, W. Steinmann, P. G¨rtler, E. E. u Koch, and A. Kozevnikov, Phys. Rev. Lett. 37 , 305 (1976).
- [62] Note that one could consider variation of K as a way to incorporate correlation effects (using a K xc ).

.
Sarah Cavo J. A. Berger Pina Romaniello
Accurate optical spectra of solids from pure time-dependent density-functional theory
We present accurate optical spectra of semiconductors and insulators within a
pure Kohn-Sham time-dependent density-functional approach. In particular, we
show that the onset of the absorption is well reproduced when comparing to
experiment. No empirical information nor a theory beyond Kohn-Sham
density-functional theory, such as $GW$, is invoked to correct the Kohn-Sham
gap. Our approach relies on the link between the exchange-correlation kernel of
time-dependent density functional theory and the derivative discontinuity of
ground-state density-functional theory. We show explicitly how to relate these
two quantities. We illustrate the accuracy and simplicity of our approach by
applying it to various semiconductors (Si, GaP, GaAs) and wide-gap insulators
(C, LiF, Ar).