Assessment of long-range-corrected exchange-correlation kernels for solids: accurate exciton binding energies via an empirically scaled Bootstrap kernel

TL;DR

This paper evaluates various long-range-corrected exchange-correlation kernels in time-dependent density-functional theory for solids, identifying limitations and proposing an empirically scaled Bootstrap kernel that accurately predicts exciton binding energies across different materials.

Contribution

The paper introduces a simple, universal, empirically scaled Bootstrap kernel that improves exciton binding energy predictions in solids with low computational cost.

Findings

No existing LRC kernel accurately predicts both spectra and binding energies.

The scaled Bootstrap kernel provides accurate exciton binding energies for all tested materials.

The study highlights the limitations of current kernels in simultaneously modeling optical spectra and exciton binding energies.

Abstract

In time-dependent density-functional theory, a family of exchange-correlation kernels, known as long-range-corrected (LRC) kernels, have shown promise in the calculation of excitonic effects in solids. We perform a systematic assessment of existing static LRC kernels (empirical LRC, Bootstrap, and jellium-with-a-gap model) for a range of semiconductors and insulators, focusing on optical spectra and exciton binding energies. We find that no LRC kernel is capable of simultaneously producing good optical spectra and quantitatively accurate exciton binding energies for both semiconductors and insulators. We propose a simple and universal, empirically scaled Bootstrap kernel which yields accurate exciton binding energies for all materials under consideration, with low computational cost.

Figures

Tables

Table 1. Table 1: LRC kernel parameters α 𝛼 \alpha and exciton binding energies E b subscript 𝐸 b E_{\mathrm{b}} (in meV) obtained from the RPA-Bootstrap kernel using different types of the scissors shift.

Scissors shift		GaAs		Ne
$E_{c𝐤}\to E_{c𝐤}+\Delta$	$\hat{p}\to {\hat{p}}_{\mathrm{renorm}}$	$\alpha$	$E_b$	$\alpha$	$E_b$
Yes	Yes	0.116	0.601	37.5	6000
Yes	No	0.284	0.246	66.7	7420
No	No	0.073	0.344	30.9	2400

Table 2. Table 2: Experimental and calculated exciton binding energies E b subscript 𝐸 b E_{\mathrm{b}} (in meV).

			GaAs	$\alpha$-GaN	$\beta$-GaN	AlN	MgO	LiF	Ar	Ne
Exp.			3.27	20.4	26.0	48.0	80.0	1600	1900	4080
scaled Boot	${ϵ}_{\mathrm{LDA}\left(\mathrm{DFPT}\right)}^{-1}$	Eq. (27)	3.30	23.1	21.4	97.4	90.1	1790	1230	5220
scaled Boot	${ϵ}_{\mathrm{LDA}\left(\mathrm{DFPT}\right)}^{-1}$	Eq. (28)	3.30	23.0	21.4	97.4	92.3	1790	1230	5190
scaled Boot	${ϵ}_{\mathrm{RPA}\left(\mathrm{SOS}\right)}^{-1}$	Eq. (27)	3.24	22.2	22.1	90.4	97.2	1710	1220	5410
scaled Boot	${ϵ}_{\mathrm{RPA}\left(\mathrm{SOS}\right)}^{-1}$	Eq. (28)	3.24	22.1	22.0	90.6	102	1720	1210	5350
RPA-Boot			0.344	1.06	1.01	0.00	2.12	94.7	96.0	2400
0-Boot			0.293	0.919	0.829	0.00	1.72	43.2	13.7	612
Boot			0.278	0.735	0.649	0.00	1.20	14.8	9.14	101
JGM			0.141	0.438	0.279	0.00	0.397	12.1	17.1	5.96
LRC			0.670	1.33	1.32	0.00	0.855	1.89	1.54	1.06

Table 3. Table 3: Fitting parameters for the scaling factor A 𝐴 A .

$x$	$a_1$	$a_2$	$a_3$	$a_4$	$b_1$	$b_2$	$b_3$	$b_4$
${ϵ}_{\mathrm{LDA}\left(\mathrm{DFPT}\right)}^{-1}$	25.9	-0.159	0.161	1.16	6.89	1.11	0.166	1.16
${ϵ}_{\mathrm{RPA}\left(\mathrm{SOS}\right)}^{-1}$	11.6	-0.00239	0.148	1.10	5.56	1.25	0.155	1.11

Equations

$$f_{\mathrm{xc}}^{\mathrm{LRC}}=-\frac{\alpha }{{\mathbf{q}}^2},$$

$${\chi }_{\mathbf{G}{\mathbf{G}}^{\prime }}(\mathbf{q},\omega )={\chi }_{\mathbf{G}{\mathbf{G}}^{\prime }}^{(0)}(\mathbf{q},\omega )+\sum _{{\mathbf{G}}_1{\mathbf{G}}_2}{\chi }_{\mathbf{G}{\mathbf{G}}_1}^{(0)}(\mathbf{q},\omega )$$

$$\times \left[V_{{\mathbf{G}}_1}(\mathbf{q}){\delta }_{{\mathbf{G}}_1{\mathbf{G}}_2}+f_{\mathrm{xc},{\mathbf{G}}_1{\mathbf{G}}_2}(\mathbf{q})\right]{\chi }_{{\mathbf{G}}_2{\mathbf{G}}^{\prime }}(\mathbf{q},\omega ),$$

$${\chi }_{\mathbf{G}{\mathbf{G}}^{\prime }}^{(0)}(\mathbf{q},\omega )=\frac{2}{\mathcal{V}}\sum _{nm\mathbf{k}}(f_{m\mathbf{k}+\mathbf{q}}-f_{n\mathbf{k}})$$

$$\times \frac{⟨m\mathbf{k}+\mathbf{q}|e^{i(\mathbf{k}+\mathbf{G})\cdot \mathbf{r}}|n\mathbf{k}⟩⟨n\mathbf{k}|e^{-i(\mathbf{k}+{\mathbf{G}}^{\prime })\cdot {\mathbf{r}}^{\prime }}|m\mathbf{k}+\mathbf{q}⟩}{E_{m\mathbf{k}+\mathbf{q}}-E_{n\mathbf{k}}-(\omega +i\eta )},$$

$${\chi }_{00}^{(0)}(\mathbf{q}\to 0,0)$$

$${ϵ}_{\mathrm{M}}(\omega )$$

$$=\underset{\mathbf{q}\to 0}{\lim }\frac{1}{1+V_0(\mathbf{q}){\chi }_{00}(\mathbf{q},\omega )},$$

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

$$A_{vc\mathbf{k},v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime }}$$

$$B_{vc\mathbf{k},v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime }}$$

$$F_{vc\mathbf{k},v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime }}^{\mathrm{H}}$$

$$F_{vc\mathbf{k},v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime }}^{\mathrm{xc}}$$

$$\times ⟨v^{\prime }{\mathbf{k}}^{\prime }|e^{-i(\mathbf{q}+{\mathbf{G}}^{\prime })\cdot \mathbf{r}}|c^{\prime }{\mathbf{k}}^{\prime }⟩.$$

$$\underset{\mathbf{q}\to 0}{\lim }{f}_{\mathrm{xc},\mathbf{G}{\mathbf{G}}^{\prime }}(\mathbf{q})=\left(\begin{array}{cccc}\frac{{\kappa }_{00}}{{\mathbf{q}}^2}& \frac{{\kappa }_{01}}{\mathbf{q}}& \frac{{\kappa }_{02}}{\mathbf{q}}& \cdots \\ \frac{{\kappa }_{10}}{\mathbf{q}}& {\kappa }_{11}& {\kappa }_{12}& \cdots \\ \frac{{\kappa }_{20}}{\mathbf{q}}& {\kappa }_{21}& {\kappa }_{22}& \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right),$$

$$f_{\mathrm{xc}}^{\mathrm{LRC}(\mathrm{d})}=-\frac{\alpha }{4\pi }V=\left(\begin{array}{cccc}-\frac{\alpha }{{\mathbf{q}}^2}& 0& 0& \cdots \\ 0& -\frac{\alpha }{{\mathbf{G}}_1^2}& 0& \cdots \\ 0& 0& -\frac{\alpha }{{\mathbf{G}}_2^2}& \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right)$$

$$f_{\mathrm{xc}}^{\mathrm{LRC}(\mathrm{h})}=-\frac{\alpha }{4\pi }{V}_0=\left(\begin{array}{cccc}-\frac{\alpha }{{\mathbf{q}}^2}& 0& 0& \cdots \\ 0& 0& 0& \cdots \\ 0& 0& 0& \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right).$$

$${\alpha }_{\mathrm{LRC}}=\frac{C_1}{{ϵ}_{\infty }}-C_2,$$

$$\displaystyle F^{\mathrm{Hxc}}_{vc\mathbf{k},v^{\prime}c^{\prime}\mathbf{k^{\prime}}}=\frac{2}{\cal V}\Big{(}\sum_{\mathbf{G}\neq 0}\frac{4\pi-\bar{\alpha}}{|\mathbf{G}|^{2}}\langle c\mathbf{k}|e^{i\mathbf{G}\cdot\mathbf{r}}|v\mathbf{k}\rangle\langle v^{\prime}\mathbf{k^{\prime}}|e^{-i\mathbf{G}\cdot\mathbf{r}}|c^{\prime}\mathbf{k^{\prime}}\rangle$$

$$\displaystyle-\alpha_{0}\frac{\langle c\mathbf{k}|\hat{p}+i[V_{\mathrm{NL}},\hat{r}]|v\mathbf{k}\rangle}{E_{c\mathbf{k}}-E_{v\mathbf{k}}}\frac{\langle c^{\prime}\mathbf{k^{\prime}}|\hat{p}+i[V_{\mathrm{NL}},\hat{r}]|v^{\prime}\mathbf{k^{\prime}}\rangle^{*}}{E_{c^{\prime}\mathbf{k^{\prime}}}-E_{v^{\prime}\mathbf{k^{\prime}}}}\Big{)},$$

$$f_{\mathrm{xc},\mathbf{G}{\mathbf{G}}^{\prime }}^{\mathrm{B}}(\mathbf{q},\omega )=\frac{V_{\mathbf{G}}^{1/2}(\mathbf{q}){ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-1}(\mathbf{q},0)V_{{\mathbf{G}}^{\prime }}^{1/2}(\mathbf{q})}{1-{ϵ}_{\mathrm{RPA},00}(\mathbf{q},0)},$$

$$f_{\mathrm{xc}}^{\mathrm{B}}=\left(\begin{array}{cccc}\frac{{\beta }_{00}}{{\mathbf{q}}^2}& \frac{{\beta }_{01}}{|\mathbf{q}||{\mathbf{G}}_1|}& \frac{{\beta }_{02}}{|\mathbf{q}||{\mathbf{G}}_2|}& \cdots \\ \frac{{\beta }_{10}}{|{\mathbf{G}}_1||\mathbf{q}|}& \frac{{\beta }_{11}}{{\mathbf{G}}_1^2}& \frac{{\beta }_{12}}{|{\mathbf{G}}_1||{\mathbf{G}}_2|}& \cdots \\ \frac{{\beta }_{21}}{|{\mathbf{G}}_2||\mathbf{q}|}& \frac{{\beta }_{21}}{|{\mathbf{G}}_2||{\mathbf{G}}_1|}& \frac{{\beta }_{22}}{{\mathbf{G}}_2^2}& \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right),$$

$${\beta }_{\mathbf{G}{\mathbf{G}}^{\prime }}=\frac{4\pi {ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-1}(\mathbf{q},0)}{1-{ϵ}_{\mathrm{RPA},00}(\mathbf{q},0)}.$$

$$f_{\mathrm{xc}}^{\mathrm{B}(\mathrm{h})}=\left(\begin{array}{cccc}\frac{{\beta }_{00}}{{\mathbf{q}}^2}& 0& 0& \cdots \\ 0& 0& 0& \cdots \\ 0& 0& 0& \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right).$$

$${\alpha }_{\mathrm{B}}=\frac{4\pi {ϵ}_{00}^{-1}(0,0)}{{ϵ}_{\mathrm{RPA},00}(0,0)-1}.$$

$${\alpha }_{0B}=\frac{4\pi {ϵ}_{\mathrm{RPA},00}^{-1}(0,0)}{{ϵ}_{\mathrm{RPA},00}(0,0)-1}.$$

$${\alpha }_{\mathrm{RPAB}}=\frac{4\pi {ϵ}_{\mathrm{RPA},00}^{-1}(0,0)}{1/{ϵ}_{\mathrm{RPA},00}^{-1}(0,0)-1}$$

$$f_{\mathrm{xc}}^{\mathrm{JGM}}(\mathbf{q};n,E_{\mathrm{g}})$$

$$k_{n,E_{\mathrm{g}}}^{\prime }=k_n+\frac{E_{\mathrm{g}}^2}{4\pi n{\mathbf{q}}^2}\left(\frac{1+E_{\mathrm{g}}}{B(n)+E_{\mathrm{g}}}\right).$$

$$f_{\mathrm{xc}}^{\mathrm{sB}}=A(x)f_{\mathrm{xc}}^{\mathrm{RPAB}}=-A(x)\frac{4\pi x}{(1/x-1){\mathbf{q}}^2},$$

$$A(x)$$

$$=b_1{e}^{-x^{b_2}/b_3}+b_4.$$

Full text

I Introduction

The optical properties of insulators and semiconductors in the energy range close to the gap are strongly influenced by excitons. The accurate and efficient calculation of excitonic properties is an important task of computational materials science, since it is a key requirement in the design of novel photovoltaic materials of desired properties. For example, low exciton binding energies in perovskite solar cells promote the electron-hole separation and thereby enhance power conversion efficiencies.Miyata et al. (2015)

Many-body perturbation theory is a standard theoretical method for excitonic effects in solids: accurate exciton binding energies $E_{\mathrm{b}}$ and optical absorption spectra of semiconductors and insulators are obtained by solving the Bethe-Salpeter equation (BSE).Rohlfing and Louie (1998); Onida et al. (2002); Martin et al. (2016) However, the BSE is computationally too expensive to be applied to large systems. Time-dependent density-functional theory (TDDFT)Ullrich (2012); Ullrich and Yang (2015) provides alternatives to the BSE which are computationally much cheaper.

The main challenge for TDDFT lies in finding approximations to the exchange-correlation (xc) kernel $f_{\mathrm{xc}}$ which yield accurate excitonic properties. The random-phase approximation (RPA) ($f_{\mathrm{xc}}=0$), the local-density approximation (LDA), and generalized gradient approximations (GGAs) fail to capture excitonic effects in solids due to their inadequate long-range behavior. The so-called “nanoquanta kernel”,Reining et al. (2002); Sottile et al. (2003); Adragna et al. (2003); Marini et al. (2003) constructed by reverse-engineering the BSE, yields very good optical spectra of solids and thus provides an important proof of principle; however, it is computationally as expensive as the BSE.

Hybrid xc functionals (mixtures of semilocal xc functionals with a fraction of nonlocal Fock exchange) are very widely used in TDDFT. The B3LYP hybrid functional Stephens et al. (1994) gives reasonably good optical spectra for systems whose gap is not too large.Bernasconi et al. (2011); Tomić et al. (2014) For organic molecular crystals, the so-called optimally tuned range-separated hybrids produce excellent results.Refaely-Abramson et al. (2015) A scaled exact exchange approach was recently shown to yield good excitonic binding energies for a wide variety of materials.Yang et al. (2015) However, the nonlocal exchange contribution adds to the computational cost of the hybrid methods; it is therefore desirable to work with purely local xc functionals.

A simple nonlocal model kernel, which is known as the long-range-corrected (LRC) kernel,Reining et al. (2002); Botti et al. (2004, 2007)

[TABLE]

where ${\bf q}$ is the momentum transfer in the first Brillouin zone (BZ), can account for bound excitons in solids, but it requires a material-dependent parameter $\alpha$, a positive scalar. A number of xc kernels proposed in the literature, such as the empirical LRC, Bootstrap, RPA-Bootstrap, and jellium-with-gap-model (JGM) kernels, Botti et al. (2004, 2007); Sharma et al. (2011); Rigamonti et al. (2015); Trevisanutto et al. (2013) report that the long-range part in them gives the most important contribution to their results, and we hence refer to them as the family of LRC-type kernels. These kernels have been applied to simple bulk semiconductors and insulators, with some degree of success. However, there also were reports of conflicting results, giving rise to some recent controversies in the literature.Sharma et al. (2016); Rigamonti et al. (2016)

Testing the performance of the various LRC-type kernels is a complex task which depends on many choices. For instance, the xc kernel, which is formally a matrix in reciprocal space, can be implemented as head-only, diagonal, or a full matrix. Local-field effects can be fully or partially included, or completely ignored. The calculated optical spectra depend on the input band structure (LDA with or without scissors correction, GGA, LDA+U, hybrids, or GW) and on the method (such as all-electron versus pseudopotential-based). And, last but not least, the selection of the materials is important. Given the large number of choices that have to be made, an unbiased assessment and a comparison between different LRC methods is challenging, and conflicting results can arise.

In this paper, we will perform a systematic assessment of the various existing static LRC-type kernels (i.e. we do not assess dynamical LRC-type kernels such as those proposed in Refs. Botti et al., 2005; Berger, 2015), for a variety of materials ranging from small-gap semiconductors to large-gap insulators, comparing calculated optical spectra and exciton binding energies to experimental data. The main finding is that the existing LRC-type kernels, while often producing good-looking optical spectra for semiconductors, all fail to yield consistently good exciton binding energies. We propose an empirical scaling approach, to be used in conjunction with the RPA-Bootstrap method, which gives accurate $E_{\rm b}$ for all materials under study, but the resulting optical spectra may have unsatisfactory distributions of oscillator strength.

This paper is organized as follows. In Section II, we give an overview of the formal framework of linear-response TDDFT, comparing two approaches to describe optical properties of solids: the Dyson-equation approach and the Casida equation. We then review the existing static LRC-type xc kernels, and the different choices for their implementation. We also discuss some computational details. Section III then presents our results. We demonstrate the sensitivity of the optical spectra to the choice of the $\alpha$-parameter, and then propose a scaled RPA-Bootstrap kernel which gives accurate exciton binding energies. Section IV contains our conclusions.

II Background and methodology

II.1 Linear-response TDDFT for solids: Dyson equation vs Casida equation approach

There are several ways to calculate optical absorption spectra of periodic systems using linear-response TDDFT.Ullrich and Yang (2015) The most common approach is based on the interacting density-density response function $\chi_{{\bf G}{\bf G}^{\prime}}(\mathbf{q},\omega)$, where ${\bf G}$ and ${\bf G}^{\prime}$ are reciprocal lattice vectors, and $\omega$ is the frequency. The response function is obtained from the following Dyson-type equation:

[TABLE]

where $\chi^{(0)}$ is the noninteracting response function and $V_{\bf G}({\bf q})=4\pi/|\mathbf{q}+\mathbf{G}|^{2}$ is the Coulomb interaction. It is convenient to write $V=V_{0}+\bar{V}$, where $V_{0}$ is the long-range ($\mathbf{G}=0$) part of the Coulomb interaction, and $\bar{V}$ is the Coulomb interaction without the long-range part. $f_{{\rm xc},{\bf G}{\bf G}^{\prime}}(\mathbf{q})$ is the xc kernel in the adiabatic approximation, i.e., independent of $\omega$. $\chi^{(0)}$ is explicitly given by Gajdoš et al. (2006)

[TABLE]

where $\mathbf{k}$ lies within the first BZ, $n$ and $m$ are band indices, $E_{n\mathbf{k}}$ and $E_{m\mathbf{k+q}}$ are the associated Kohn-Sham single-particle energies, $f=1(0)$ for occupied (unoccupied) states, the factor of 2 accounts for the spin (we here only consider non-spin-polarized systems), $\mathcal{V}$ is the crystal volume, and $\eta$ is an infinitesimal. In the optical limit ($\mathbf{q}\to 0$), the head ($\mathbf{G}=\mathbf{G^{\prime}}=0$) of $\chi^{(0)}$ at $\omega=0$ becomesBaroni and Resta (1986)

[TABLE]

where $v$ and $c$ are valence and conduction band indices, respectively, $\hat{p}$ is the momentum operator, $\hat{r}$ is the position operator, and $V_{\mathrm{NL}}$ is the non-local part of the pseudopotential. The ${\bf q}^{2}$-dependence will be important for the construction of the Bootstrap kernels, see below. It is also important that $\chi^{(0)}_{00}(\mathbf{q}\to 0,0)$ is always negative.

The optical spectrum is obtained from the imaginary part of the macroscopic dielectric function $\epsilon_{\mathrm{M}}$:

[TABLE]

where $\epsilon^{-1}$ is the inverse dielectric function.Onida et al. (2002) We shall refer to this method as the Dyson approach; it has a moderate computational cost, and is therefore the method of choice for calculating optical spectra. However, the drawback of the Dyson-equation approach is that fine details of the spectra, in particular the binding energies of weakly bound excitons, cannot be obtained, because the spectral broadening washes out any subtle features of the order of a few meV (see also Sec. II.3).

As an alternative which is strictly equivalent to the Dyson equation, optical spectra and exciton binding energies can be obtained from the Casida equation:Casida (1995)

[TABLE]

where $\mathbf{A}$ and $\mathbf{B}$ are excitation and de-excitation matrices, respectively, $X_{n}$ and $Y_{n}$ are $n$th eigenvectors, and $\omega_{n}$ is the $n$th excitation energy. The matrix elements of $\mathbf{A}$ and $\mathbf{B}$ are given by

[TABLE]

where $F^{\mathrm{Hxc}}=F^{\mathrm{H}}+F^{\mathrm{xc}}$ is the Hartree-exchange-correlation (Hxc) matrix.Ullrich (2012) In the optical limit, $F^{\mathrm{H}}$ and $F^{\mathrm{xc}}$ are given by

[TABLE]

For the elements of $F^{\mathrm{xc}}$ in Eq. (11) to remain finite (i.e., neither vanishing nor diverging) in the $\mathbf{q}\to 0$ limit, the head (${\bf G}={\bf G}^{\prime}=0$) of $f_{\mathrm{xc}}$ should be proportional to $\mathbf{q}^{-2}$, the wings (${\bf G}=0,{\bf G}^{\prime}\neq 0$ or vice versa) should be proportional to $\mathbf{q}^{-1}$, and the body (${\bf G},{\bf G}^{\prime}\neq 0$) should be independent of $\mathbf{q}$. In other words, the most general form is

[TABLE]

where the $\kappa_{\mathbf{GG^{\prime}}}$ are constants (in general, they are functionals of the density). We will discuss various approximations of the xc kernel in the following subsection.

The excitation energy spectrum $\omega_{n}$ of the Casida equation (7) for periodic solids with a gap has discrete levels, which correspond to bound excitons, and a continuous part, which corresponds to the unbound particle-hole excitations. For the adiabatic xc kernels considered here, only one excitonic level is found, which can be identified as the lowest bound exciton (to obtain an excitonic Rydberg series with a scalar xc kernel requires the kernel to be frequency-dependent).Yang et al. (2012); Yang and Ullrich (2013) We calculate the exciton binding energy as that energy which separates this discrete level from the onset of the continuum. Since no artificial spectral broadening is involved, exciton binding energies can be calculated in principle with arbitrary precision. However, the Casida-equation approach is computationally expensive because it requires building and diagonalizing a large matrix.

Note that by using a very small broadening width and a very fine frequency grid, one may be able to obtain $E_{\mathrm{b}}$ of semiconductors from the Dyson-equation optical spectrum, but the broadening width and the frequency grid spacing always cause an error that may be greater than $E_{\mathrm{b}}$ of interest. Note also that Ref. Rigamonti et al., 2015 proposed a method to “read” $E_{\mathrm{b}}$ from the real part of $\epsilon_{\mathrm{M}}^{\mathrm{RPA}}(\omega)$, but this approach works only for head-only kernels and only for wide-gap insulators (i.e. one cannot obtain small $E_{\mathrm{b}}$ on the order of a few meV), and it has a moderate precision ($\sim$0.1 eV). By contrast, the Casida equation works for all forms of the xc kernel and for all materials, and it has a high precision ($\sim$0.01 meV).

A widely used approach to simplify the Casida equation is the so-called Tamm-Dancoff approximation (TDA), which decouples excitations and de-excitations by setting $\mathbf{B}$ to zero in Eq. (7). However, we have found Byun and Ullrich (2017) that the TDA underestimates LRC $E_{\mathrm{b}}$ of insulators significantly (i.e. by more than 100%) (e.g. TDA and full Casida equations using the RPA-Bootstrap kernel without the scissors shift yield $E_{\mathrm{b}}=666$ and 2400 meV, respectively, for solid Ne), so we will only use the full Casida equation in this work.

The local-field effect (LFE) is determined by the number of $\mathbf{G}$ vectors included, and has different forms in the Dyson and Casida equations. In the Dyson approach, the LFE means including not only the head, but also the wings and body of the matrix in $\mathbf{G,G^{\prime}}$, which leads to $\epsilon_{00}\neq 1/\epsilon^{-1}_{00}$. The Dyson equation is used to calculate optical spectra and Bootstrap-type kernel parameters (more about this later). In the Dyson equation for optical spectra, the LFE is not a matter of choice and should be included. However, in the Dyson equation for Bootstrap-type kernel parameters, the LFE is a matter of choice because of the freedom of defining Bootstrap-type kernels. In this work, we will include the LFE when calculating Bootstrap-type kernel parameters, following the convention adopted in the literature.Sharma et al. (2011)

In the Casida equation, the LFE means including not only the head, but also other terms in the summation of $F^{\mathrm{Hxc}}$ matrix elements in Eqs. (11) and (16). Mathematically, the LFE in the Dyson equation is exactly transformed into the summation in the Casida equation. Therefore, if the LFE is included in the Dyson equation, it should be included in the Casida equation, too.

II.2 LRC-type xc kernels

In reciprocal space, the xc kernels $f_{{\rm xc},{\bf G}{\bf G}^{\prime}}({\bf q})$ are matrices in ${\bf G}$ and ${\bf G}^{\prime}$, see Eq. (12). In the following, we list the xc kernels we have tested, paying particular attention to distinguish between head-only, diagonal or full matrix forms. In all expressions, the optical limit $({\bf q}\to 0)$ is understood.

II.2.1 Empirical LRC kernel

The empirical LRC kernel was originally designed as a kernel for optical spectra of semiconductors. Reining et al. (2002) The diagonal and the head-only versions of the empirical LRC kernel are defined, respectively, as

[TABLE]

and

[TABLE]

Here, $\alpha$ is an empirical parameter, given byBotti et al. (2004)

[TABLE]

where $C_{1}=4.615$, $C_{2}=0.213$, and $\epsilon_{\infty}$ is the high-frequency dielectric constant. Note that here we use $\epsilon_{\mathrm{RPA}}^{-1}$ instead of the experimental $1/\epsilon_{\infty}$, where $\epsilon_{\mathrm{RPA}}^{-1}$ is greater than $1/\epsilon_{\infty}$ by $\sim$10%. Also note that the empirical LRC kernel used calculated lattice parameters, while we take experimental ones. Due to these differences, empirical parameters $C_{1}$ and $C_{2}$ should be re-fitted to our choices, but it turns out that such differences make little effect on LRC results for semiconductors (see below).

In general, when a head-only or diagonal LRC kernel is used, $F^{\mathrm{Hxc}}$ simplifies to

[TABLE]

where $\alpha=\alpha_{0}$ and $\bar{\alpha}=0$ for the head-only kernel $f_{\mathrm{xc}}^{\rm LRC(h)}$, and $\alpha=\alpha_{0}=\bar{\alpha}$ for the diagonal kernel $f_{\mathrm{xc}}^{\rm LRC(d)}$. Note that head-only or diagonal LRC kernels reduce the exciton Hamiltonian building time drastically because this removes the double loop over $\mathbf{G,G^{\prime}}$ in Eq. (11).

It turns out that the body of $f_{\mathrm{xc}}^{\rm LRC(d)}$ has a negligible effect on optical spectra of semiconductors such as Si:Botti et al. (2004) this is because $\bar{\alpha}\approx 0.2\ll 4\pi$ in Eq. (16). However, $f_{\mathrm{xc}}^{\rm LRC(h)}$ and $f_{\mathrm{xc}}^{\rm LRC(d)}$ can produce very different results for insulators, and one needs to state clearly which version, (h) or (d), of the xc kernel is used.

II.2.2 Bootstrap kernels

The original Bootstrap kernel is a parameter-free xc kernel for optical spectra of semiconductors and insulators. Sharma et al. (2011) The original Bootstrap kernel is defined as

[TABLE]

where $\epsilon^{-1}$ is the self-consistent (i.e. bootstrapped) inverse dielectric function. In matrix form, the bootstrap kernel is given by

[TABLE]

where

[TABLE]

Neglecting the wings and body of $f_{\mathrm{xc}}^{\mathrm{B}}$, which can be viewed as neglecting the LFE, yields a head-only Bootstrap kernel:

[TABLE]

Comparing $f_{\mathrm{xc}}^{\rm B(h)}$ with $f_{\mathrm{xc}}^{\rm LRC(h)}$, we define the LRC $\alpha$-parameter for the Bootstrap kernel as

[TABLE]

Whereas $f_{\mathrm{xc}}^{\rm LRC(d)}$ and $f_{\mathrm{xc}}^{\rm LRC(h)}$ give quite different results for insulators, we have found that $f_{\mathrm{xc}}^{\rm B}$ and $f_{\mathrm{xc}}^{\rm B(h)}$ make a relatively small difference for both semiconductors and insulators (see Table I in Supplemental MaterialSee Supplemental Material at http://… ), which is consistent with the findings of Refs. Sharma et al., 2011; Rigamonti et al., 2015. In Ref. Trevisanutto et al., 2013 the same trend was found in the JGM kernel (see below). Therefore, in view of the reduced computational effort, we use the head-only form for all kernels in the following unless stated otherwise. The only exception is when we verify the results of Dyson-equation optical spectra obtained from $f_{\mathrm{xc}}^{\mathrm{B}}$ using Casida-equation calculations (see Table I in Supplemental MaterialSee Supplemental Material at http://… ). We emphasize again that we only consider the ${\bf q}\to 0$ limit here; at finite ${\bf q}$, the matrix character of the Bootstrap kernel appears to play a more significant role.Sharma et al. (2012)

We also consider two simpler variations of the Bootstrap kernel. The first one, referred to as the 0-Bootstrap kernel,Sharma is the head-only Bootstrap kernel (20) without the built-in self-consistency (i.e., “0” means no iteration, similar to the G0W0 version of the GW approach) for optical spectra of semiconductors and insulators. The LRC $\alpha$-parameter for the 0-Bootstrap kernel is thus given by

[TABLE]

Note that $\alpha_{\mathrm{0B}}>\alpha_{\rm B}$ by about 10% because $\epsilon_{\mathrm{RPA},00}^{-1}(0,0)$ is greater than $\epsilon^{-1}_{00}(0,0)$ by about 10%.

The second simplified Bootstrap kernel is the RPA-Bootstrap kernel,Rigamonti et al. (2015) which is a head-only kernel with

[TABLE]

for exciton binding energies of insulators. Note that $\alpha_{\mathrm{RPAB}}>\alpha_{\mathrm{0B}}$ by about 10% because $1/\epsilon^{-1}_{\mathrm{RPA},00}(0,0)-1<\epsilon_{\mathrm{RPA},00}(0,0)-1$ by about 10%. Note also that without the LFE (i.e. when $\epsilon_{00}=1/\epsilon^{-1}_{00}$), the 0-Bootstrap and RPA-Bootstrap kernels become identical.

II.2.3 Jellium with a gap model

The JGM kernel is a parameter-free kernel for optical spectra of semiconductors and insulators.Trevisanutto et al. (2013) The JGM kernel is defined as

[TABLE]

with

[TABLE]

Here, $E_{\mathrm{g}}$ is the band gap, $n$ is the electron density, and $k_{\mathrm{F}}$ is the Fermi wave vector; $k_{n}$, $B(n)$, and $C(n)$ are defined in Ref. Constantin and Pitarke, 2007. $f^{\mathrm{JGM}}_{\mathrm{xc},\mathbf{GG^{\prime}}}(\mathbf{q};E_{\mathrm{g}})$ is defined as a full matrix, obtained from the Fourier transform in reciprocal space and the symmetrization in $\mathbf{G,G^{\prime}}$; however, we here use it in the head-only form. Whereas other LRC-type kernels depend on dielectric constants, the JGM kernel depends on band gaps.

II.3 Computational aspects

We used the Abinit code,Gonze et al. (2009) which is based on norm-conserving pseudopotentials, for calculating the Kohn-Sham band structures including scissors corrections, as well as GW band gaps within the LDA. Experimental lattice parameters were used for all materials. We used the dp codeOlevano et al. for calculating optical spectra from the imaginary part of the dielectric function. We calculated exciton binding energies from the Casida equation with our own homemade code.

Optical spectra were obtained with a Lorentzian broadening of 0.15 eV for GaAs, and 0.2 eV for all other materials. Note that a Lorentzian broadening smaller than 0.15 eV, which is an optimal value that makes calculated and experimental $E_{2}$ peaks have similar heights, generates an artificial $E_{1}$ peak in the excitonic region of the RPA and LRC spectra of GaAs, and makes RPA and LRC $E_{2}$ peaks higher than the experimental one, and thus should not be used (see the top panel of Fig. 2 for $E_{1}$ and $E_{2}$ peaks). The Lorentzian broadening has a physical meaning (i.e. it simulates the lifetime broadening and can be calculated from GW), Cazzaniga et al. (2011) so it should not be used as an adjustable parameter to improve the appearance of calculated optical spectra.

We used experimental band gaps, $E_{\mathrm{g}}^{\mathrm{exp}}$, as onsets of optical spectra instead of GW band gaps, $E_{\mathrm{g}}^{\mathrm{GW}}$, because there are differences on the order of 1 eV between GW and experimental band gaps in wide-gap insulators. As shown below, these differences are comparable to the exciton binding energies in the materials under consideration, and can therefore cause an artificial cancellation of the two errors in $E_{\mathrm{g}}$ and $E_{\mathrm{b}}$ when one compares the excitonic peak position in the calculated optical spectrum using $E_{\mathrm{g}}^{\mathrm{GW}}$ with the experimental optical spectrum.

In the Dyson equation for optical spectra, we used a 16$\times$16$\times$16 Monkhorst-Pack $\mathbf{k}$-point mesh, 4 valence bands, and 20 conduction bands. We found that TDDFT-LRC shows a slower convergence with respect to the number of conduction bands ($N_{\mathrm{c}}$) than the BSE (e.g. for LiF, $N_{\mathrm{c}}=6$ is enough for the BSE, Rohlfing and Louie (1998) while $N_{\mathrm{c}}\geq 12$ is needed for both Dyson and Casida equations). An insufficient number of conduction bands causes blueshifts of the excitonic peak (i.e. decreases the exciton binding energy) and reduces its oscillator strength in the LRC spectrum of wide-gap insulators significantly. This slow convergence also occurs for the real part of $\epsilon_{\mathrm{M}}^{\mathrm{LRC}}(\omega=0)$ (i.e., the LRC dielectric constant).

In the Dyson equation for Bootstrap-type kernel parameters, we used a 20$\times$20$\times$20 (20$\times$20$\times$10) $\Gamma$-centered $\mathbf{k}$-point mesh, 4 (8) valence bands, 20 (20) conduction bands, and 59 (73) $\mathbf{G}$ vectors for GaAs, $\beta$-GaN, MgO, LiF, solid Ar, and solid Ne ($\alpha$-GaN and AlN).

In the Casida equation, we used a 28$\times$28$\times$28 $\Gamma$-centered $\mathbf{k}$-point mesh, 3 valence bands, 2 conduction bands, and 59 $\mathbf{G}$ vectors for GaAs. The corresponding parameters for the other materials are: 16$\times$16$\times$16, 3, 6, 59 for $\beta$-GaN and MgO, 16$\times$16$\times$8, 6, 9, 73 for $\alpha$-GaN and AlN, and 8$\times$8$\times$8, 3, 24, 59 for LiF, solid Ar, and solid Ne.

To calculate $\alpha_{\mathrm{JGM}}$, we used a 8$\times$8$\times$8 (8$\times$8$\times$4) $\Gamma$-centered $\mathbf{k}$-point mesh and 59 (73) $\mathbf{G}$ vectors for GaAs, $\beta$-GaN, MgO, LiF, solid Ar, and solid Ne ($\alpha$-GaN and AlN).

All computational parameters listed here were chosen after performing systematic convergence tests.

II.4 Effect of the scissors shift on LRC results

The exact xc kernel can be written in the form $f_{\rm xc}=f_{\rm xc}^{\rm qp}+f_{\rm xc}^{\rm ex}$,Stubner et al. (2004); Ullrich (2012) where the “quasiparticle” part, $f_{\rm xc}^{\rm qp}$, is responsible for correcting the Kohn-Sham gap, and $f_{\rm xc}^{\rm ex}$ is the excitonic part. In the standard TDDFT approach for semiconductors and insulators,Onida et al. (2002); Ullrich and Yang (2015); Sharma et al. (2014) the quasiparticle part of $f_{\rm xc}$ is ignored, and any corrections to the Kohn-Sham gap are made directly in the input band structure, usually by means of GW or hybrid functionals; the remaining part of the xc kernel, $f_{\rm xc}^{\rm ex}$, is then approximated.

A standard and inexpensive method for correcting LDA band structures is by applying the so-called scissors shift.Levine and Allan (1989); Gonze and Lee (1997) There are several ways of applying the scissors shift to Dyson and Casida equations in Eqs. (3), (4), (11), and (16) and LRC-type kernels. The scissors shift can be applied to only the conduction bands (i.e. replacing $E_{c\mathbf{k}}$ by $E_{c\mathbf{k}}+\Delta$) or to the momentum operator as well (i.e. replacing $\hat{p}$ by $\hat{p}_{\mathrm{renorm}}=\{(E_{c\mathbf{k}}+\Delta-E_{v\mathbf{k}})/(E_{c\mathbf{k}}-E_{v\mathbf{k}})\}\hat{p}$, where $\hat{p}_{\mathrm{renorm}}$ is the renormalized momentum operator), Levine and Allan (1989); Del Sole and Girlanda (1993) where $\Delta$ is the difference between the experimental (or GW) and DFT band gaps.

Excitonic effects within the LRC approximation are quite sensitive to the particular implementation of the scissors shift. For example, Table 1 shows $\alpha$ and $E_{\mathrm{b}}$ of GaAs and solid Ne obtained from the RPA-Bootstrap kernel using different types of the scissors shift. We find that the scissors shift affects the LRC results significantly.

In the following, we choose not to apply the scissors shift to $E_{c\mathbf{k}}$ and $\hat{p}$, i.e., we calculate exciton binding energies from the Casida equation using the uncorrected LDA band structure as input. Optical spectra, obtained from the Dyson-equation approach, are also calculated based on the uncorrected LDA band structure, and then rigidly shifted to align them with the experimental band gap. We have chosen this approach for several reasons.

Firstly, the scissors shift is not related to excitons. The scissors shift is a matter of choice for the study of excitonic effects in solids. Our aim is to make the simplest choices (e.g. the LDA and the head-only kernel) and to focus on identifying the origin of conflicting results in existing kernels and designing a new kernel. Note that we applied the scissors shift to $E_{c\mathbf{k}}$ and $\hat{p}$ to reproduce the results of existing kernels, which are provided in the Supplemental Material.See Supplemental Material at http://…

Secondly, we were concerned about the unphysically strong influence of the scissors shift on the LRC results. The scissors shift increases $\epsilon_{00}^{-1}$ by $\sim$10%; this small increase in $\epsilon_{00}^{-1}$ affects the LRC results for wide-gap insulators significantly when the LRC-type kernel depends on the dielectric constant (see below). In other words, $f_{\mathrm{xc}}^{\mathrm{ex}}=f_{\mathrm{xc}}^{\mathrm{ex}}(f_{\mathrm{xc}}^{\mathrm{qp}})$, but this is not what $f_{\mathrm{xc}}^{\mathrm{qp}}$ and $f_{\mathrm{xc}}^{\mathrm{ex}}$ are meant to be. The big change in the LRC results due to the scissors shift is not associated with excitons.

Thirdly, it allows us to eliminate one source of conflicting results. Some kernels use $E_{\mathrm{g}}^{\mathrm{GW}}$ in the scissors shift, compare their optical spectra with experiment, and reproduce or predict the excitonic peak position for wide-gap insulators by interpreting the $\sim$1 eV error in $E_{\mathrm{g}}^{\mathrm{GW}}$ as $E_{\mathrm{b}}$. Sharma et al. (2011); Berger (2015) In addition, the small difference between $E_{\mathrm{g}}^{\mathrm{exp}}$ and $E_{\mathrm{g}}^{\mathrm{GW}}$ (or $E_{\mathrm{g}}^{\mathrm{GW}}$ obtained from different potential methods) makes a small difference in the scissors-shifted dielectric constant, which can cause a big difference in the LRC results for wide-gap insulators. By not using the scissors shift, we can avoid these unnecessary complications.

Lastly, by not using the scissors shift we can avoid expensive $E_{\mathrm{g}}^{\mathrm{GW}}$ calculations for unknown materials when we need only $E_{\mathrm{b}}$. When $E_{\mathrm{g}}^{\mathrm{GW}}$ is not calculated and the dielectric constant is calculated from density-functional perturbation theory (DFPT), Baroni et al. (2001); Gajdoš et al. (2006) which is computationally much cheaper than the sum-over-states (SOS) method (6) because conduction bands are not needed, large-scale or high-throughput screening exciton calculations become possible.

III Results and Discussion

III.1 Comparison of LRC $\alpha$-parameters

In the following, we will discuss our results for the excitonic properties of the bulk semiconductors GaAs, $\alpha$-GaN, and $\beta$-GaN, the narrow-gap insulators AlN and MgO, and the wide-gap insulators LiF, solid Ar, and solid Ne. The experimental exciton binding energies are obtained from Refs. Parenteau et al., 1992; As et al., 1997; Muth et al., 1997; Haensel et al., 1969; Roessler and Walker, 1967a; Saile and Koch, 1979; Leute et al., 2009; Roessler and Walker, 1967b. We point out again that all results shown below were obtained without using the scissors shift.

Let us begin with an assessment of the LRC $\alpha$-parameters for various materials. Figure 1 compares $\alpha_{\rm LRC}$, $\alpha_{\rm JGM}$, $\alpha_{\rm B}$, $\alpha_{\rm 0B}$, and $\alpha_{\rm RPAB}$ with the $\alpha$-parameter $\alpha_{\mathrm{exp}}$ which, when used in the head-only LRC kernel (14), reproduces the experimental exciton binding energy for each material under consideration. We see that $\alpha$ varies from $\sim$0.1 ($\alpha_{\mathrm{RPAB}}$ for GaAs) to $\sim$30 ($\alpha_{\mathrm{RPAB}}$ for solid Ne). All calculated $\alpha$-parameters are smaller than the experimentally fitted ones.

III.2 Sensitivity of optical spectra to changes in $\alpha$

Next, we examine the effects of the head-only LRC kernel on optical absorption spectra. Figure 2 shows calculated optical spectra of GaAs and solid Ne obtained from the Dyson equation using the LRC kernel with $\alpha=A\alpha_{\mathrm{RPAB}}$, where $A$ is a scaling factor, and compares them with experiment.Parenteau et al. (1992); Saile and Koch (1979) We chose GaAs and solid Ne because they are extreme examples of semiconductors with weakly bound Wannier-Mott excitons and insulators with strongly bound Frenkel excitons. In the case of GaAs, the optical spectrum shows two prominent peaks above the band gap; $E_{1}$ can be interpreted as a continuum exciton. The bound exciton below the gap is very weak, and not visible on the scale of this plot because $E_{\rm b}$ is much smaller than the line broadening. To see the bound Wannier-Mott exciton of GaAs, high-resolution spectroscopy at low temperatures is needed.Ulbrich (1985) On the other hand, for solid Ne the excitonic peak is very prominent and far from the gap, and it is easy to obtain $E_{\rm b}$ from the spectrum.

In the top panel of Fig. 2, we show calculated optical spectra of GaAs for a range of $A$ between 0.8 and 4.0. We find that the spectra are rather insensitive to the scaling: a 10% change in $\alpha$ has only a very small effect: in other words, $\alpha$ has a big margin for semiconductors. The RPA spectrum of GaAs is already quite similar to experiment, apart from the height of the $E_{1}$-peak. To obtain the experimental height of the $E_{1}$-peak, a scaling factor of $A\approx 4$ (i.e. $\alpha\approx 0.3$) is needed; however, this also increases the peak width, and the valley between the $E_{1}$ and the $E_{2}$ peak becomes too high.

The bottom panel of Fig. 2 shows the calculated spectra of solid Ne for a much smaller range of $A$, between 0.8 and 1.1. Here, the spectra are very sensitive to the change in $\alpha$: a 10% change shifts excitonic peaks by about 1 eV: in other words, $\alpha$ has a small margin for insulators. Clearly, the RPA spectrum of Ne is completely different from experiment, and the LRC kernel reshapes it significantly. Using $A\approx 1.1$ puts the excitonic peak at the right position; however, the peak height and width is now drastically overestimated.

The low sensitivity of LRC results for semiconductors to changes in $\alpha$ explains why there are so many LRC-type kernels. LRC-type kernels only slightly modify RPA spectra of semiconductors, which are already very close to experiment, and $\alpha$ has a big margin for semiconductors. Thus, all LRC-type kernels produce similar and seemingly good optical spectra of semiconductors even when they use different choices and yield very different $\alpha$ values (e.g. $\alpha_{\mathrm{LRC}}\approx 0.2$ and $\alpha_{\mathrm{B}}\approx 0.1$ for Si and GaAs).

The high sensitivity of LRC results for insulators to changes in $\alpha$ is consistent with the finding of Ref. Rigamonti et al., 2015. The idea of the RPA-Bootstrap kernel is to increase $E_{\mathrm{b}}$ for wide-gap insulators from $\sim$0.1 eV to $\sim$1 eV by increasing $\alpha_{\mathrm{B}}$ by $\sim$20% for all materials. However, the $\sim$20% increase in $\alpha_{\mathrm{B}}$ does not fix the problem of the Bootstrap kernel of not reproducing an excitonic peak in the optical spectrum of semiconductors such as Si, because of the low $\alpha$-sensitivity of LRC results for semiconductors.

We also point out that the LRC results show a similar sensitivity trend to diagonal versus head-only LRC kernels and full versus TDA Casida equations (i.e. insensitive for semiconductors, but sensitive for insulators). Byun and Ullrich (2017)

These two examples already indicate a general limitation of the LRC kernel that applies to all materials: it is impossible to obtain the correct position and the correct height and width of an excitonic peak in the LRC spectrum, for both semiconductors and insulators. We will provide more evidence for this conclusion and give more examples below. To reproduce a given excitonic feature for both semiconductors and insulators (e.g., the peak height or the peak position), it is clear that a nonuniform scaling factor for Bootstrap-type kernels will be needed: the scaling factor should be close to 1 for insulators, but much greater than 1 for semiconductors. Any method which nearly uniformly scales Bootstrap-type kernels for all materials [such as using different dielectric constants (e.g. bootstrapped vs not, scissors-shifted vs not, macroscopic vs microscopic, or RPA vs LDA, all of which are different from each other by $\sim$10%) in the numerator and/or the denominator of Eqs. (17), (21), (22), and (23) or using different band structures] is likely to fail to produce satisfactory results across the board.

III.3 Nonuniformly scaled Bootstrap kernel

In Fig. 1 we compared the $\alpha$ values from head-only LRC-type kernels for various materials, and found that for wide-gap insulators, $\alpha_{\mathrm{RPAB}}$ shows the most similar trend to $\alpha_{\mathrm{exp}}$ (e.g. Bootstrap and 0-Bootstrap kernels yield $E_{\mathrm{b}}$ of solid Ar that is smaller than that of LiF). We therefore choose it as the basis for constructing a new, scaled Bootstrap xc kernel.

Let us first define $f_{\rm xc}^{\rm B(h)}=A_{\rm B}f_{\rm xc}^{\rm RPAB}$ and $f_{\rm xc}^{\rm 0B}=A_{\rm 0B}f_{\rm xc}^{\rm RPAB}$. The values of $A_{\rm B}$ and $A_{\rm 0B}$ are plotted in Fig. 3 as a function of $\epsilon_{\rm RPA}^{-1}$ for various materials; we find that $A_{\mathrm{B}}\approx 0.8$ and $A_{\mathrm{0B}}\approx 0.9$ for all materials (i.e. Bootstrap-type kernels are nearly uniformly scaled to each other). On the other hand, if we define $f_{\rm xc}^{\rm exp}=A_{\rm exp}f_{\rm xc}^{\rm RPAB}$ (i.e. $\alpha_{\mathrm{exp}}=A_{\mathrm{exp}}\alpha_{\mathrm{RPAB}}$) as the head-only LRC xc kernel which reproduces the experimental exciton binding energy, we can see that $A_{\rm exp}$ varies strongly as a function of material, from $\sim$1.1 (solid Ne) to $\sim$5 (GaAs). This non-uniform variation is consistent with our observations from the optical spectra of GaAs and solid Ne, see Fig. 2.

The values of $A_{\rm exp}$ show a rather smooth behavior as a function of $\epsilon_{\rm RPA}^{-1}$, which suggests that a scaled Bootstrap kernel can be defined via a fit to the experimental data:

[TABLE]

where $x=\epsilon_{\mathrm{RPA}}^{-1}$ (alternatively, choosing $x=\epsilon_{\mathrm{LDA}}^{-1}$ would have been possible as well). Note that both SOS and DFPT methods yield the same $x$ value. Gajdoš et al. (2006) Among many $\epsilon^{-1}$, we used $\epsilon_{\mathrm{RPA(SOS)}}^{-1}$ and $\epsilon_{\mathrm{LDA(DFPT)}}^{-1}$ in this work because they can be easily obtained from the Abinit code. Among the two $\epsilon^{-1}$, we used $\epsilon_{\mathrm{RPA(SOS)}}^{-1}$ to obtain $E_{\mathrm{b}}$ and optical spectra in this work unless stated otherwise.

We found two fitting functions, which describe well the non-uniformity of $A_{\rm exp}$,

[TABLE]

The fitting parameters $a_{i}$ and $b_{i}$, where $i=1,2,3,4$, are summarized in Table 3. Note that these fitting parameters are appropriate for the specific choices made here: experimental lattice constant, pseudopotential method, LDA band structure, head-only LRC kernel, LFE, and no scissors shift. If other choices are made, such as an all-electron method or using the scissors shift, one needs to re-fit the parameters $a_{i}$ and $b_{i}$. This calibration is inevitable due to the high sensitivity of the LRC results for wide-gap insulators. We found that the two fitting functions yield almost the same result for $A$ (and thus $\alpha$ and $E_{\mathrm{b}}$) except for $\epsilon_{\infty}\gg 10$ (see Table 2). Among the two fitting functions, we chose to use Eq. (28) to obtain $E_{\mathrm{b}}$ and optical spectra.

Figure 4 and Table 2 show experimental and calculated $E_{\mathrm{b}}$ of various materials. Whereas other kernels underestimate $E_{\mathrm{b}}$ for all materials by $\sim$10 times, the scaled Bootstrap kernel yields accurate and consistent $E_{\mathrm{b}}$. The most significant deviations are for AlN (where all other kernels give zero exciton binding energy) and for solid Ar (where even the BSE underestimates $E_{\rm b}$ by $\sim$0.3 eV Sottile et al. (2007)).

Figure 5 shows experimental and calculated optical spectra of GaAs, MgO, LiF, and solid Ne. We included LiF because it is one of two extreme examples of wide-gap insulators. We also included MgO because the LRC spectrum of MgO is very different from experiment at all $\alpha$ values, so it is impossible to determine an optimal $\alpha$ value for MgO by varing $\alpha$ (i.e., no $\alpha$ exists that reproduces the correct excitonic peak height or position). Botti et al. (2004) Here, we report the LRC spectrum of MgO when $\alpha\approx\alpha_{\mathrm{exp}}$. Bootstrap-type kernels with similar $\alpha$ values produce very similar optical spectra of GaAs (a semiconductor) and MgO (a narrow-gap insulator), but very different ones of LiF and solid Ne (wide-gap insulators). As discussed earlier, this is due to the different sensitivity of LRC spectra to semiconductors and insulators.

Our scaled Bootstrap kernel, which is designed to reproduce $E_{\mathrm{b}}^{\mathrm{exp}}$, yields excitonic peaks with overestimated (i.e. higher and wider) oscillator strengths in optical spectra of GaAs and MgO, while other Bootstrap-type kernels, which underestimate $E_{\mathrm{b}}$ by $\sim$10 times, barely produce excitonic peaks. This indicates that the LRC kernel cannot produce correct exciton binding energies and optical spectra at the same time for all materials (i.e. for semiconductors and insulators). Our finding is consistent with the LRC spectrum of ZnO, in which the calculated excitonic peak is much higher and wider than the experimental one.Gori et al. (2010)

We emphasize that our kernel is empirical, but universal in that it works for all materials and all choices. In contrast, the RPA-Bootstrap kernel, one of parameter-free kernels, works only for wide-gap insulators under special conditions such as experimental lattice parameters, the head-only kernel, and the scissors shift. In principle, a parameter-free LRC-type kernel cannot be universal for all choices due to the high sensitivity of LRC results for wide-gap insulators; thus, a trade-off between parameter-free and universal is unavoidable.

Our scaling approach is not just another Bootstrap-type kernel or a correction to the RPA-Bootstrap kernel: it is a method to predict $E_{\mathrm{b}}$ of unknown materials using the experimental $E_{\mathrm{b}}$ of a few known materials as input. The RPA-Bootstrap kernel is merely used as a fitting function, which was chosen to demonstrate the problems of popular Bootstrap-type kernels Sharma et al. (2011); Rigamonti et al. (2015); Berger (2015) and to suggest a simple way to fix them. One has the full freedom to use any other LRC-type fitting functions for our method.

IV Conclusions

In this paper, we have carried out a systematic numerical assessment of the family of static long-range-corrected (LRC) xc kernels for solids. The main challenge faced by TDDFT for the optical spectral properties of semiconductors and insulators is to reproduce the excitonic peaks at the right position and with the correct strength. We have used two methods: the Dyson-equation approach, which yields optical spectra, and the Casida-equation approach, which allows a precise determination of exciton binding energies. The two methods are equivalent, i.e., they give, in principle, the same excitonic peak positions, but in their practical implementations they are very different: from the Dyson equation approach, and the resulting macroscopic dielectric function, one cannot extract the binding energies of weakly bound excitons. Hence, the Casida approach is a very useful method, complementing the standard Dyson approach.

We have studied a group of materials, ranging from small-gap semiconductors to large-gap insulators, with exciton binding energies between a few meV and several eV. For these materials, we have tested the empirical LRC kernel, several flavors of the Bootstrap kernel, and the jellium-with-a-gap model. Most of these methods produce decently-looking optical spectra for semiconductors, but the exciton binding energies are consistently underestimated. We proposed a new xc kernel, obtained via a material-dependent scaling of the RPA-Bootstrap kernel. The scaled Bootstrap kernel is designed to produce accurate exciton binding energies for all materials under study, at very low computational cost. However, there is a price to pay: it turns out that it is impossible to obtain accurate exciton binding energies and good optical spectra at the same time for all materials using any LRC method—if the exciton peak is at the right place, the oscillator strength (i.e., the peak height and width) tends to be exaggerated for both semiconductors and insulators.

In general, assessing the performance of xc kernels for excitonic properties is a delicate task, because there are many choices involved. Here, we chose to use LDA band structures obtained with a pseudopotential code, we included local-field effects, and we implemented the xc kernels in their head-only forms. These choices will affect the numerical results: whereas the spectra of semiconductors are relatively insensitive to the strength $\alpha$ of the head of the LRC kernel, the spectra of insulators are very sensitive. Hence, it is crucial that all choices made are clearly identified, in order to facilitate comparison between results obtained by different research groups.

The main outcome of our work is that we have developed a method which can produce accurate exciton binding energies at a low computational cost. In practice, the parameters for the scaling function should be re-fitted for each particular implementation, using a small test set of small- and large-gap materials. It should then be possible to obtain accurate exciton binding energies for other, more complicated materials. Such calculations are currently in progress.

The ultimate goal is to develop TDDFT approaches that yield both accurate exciton binding energies and spectral shapes. As we have seen, the LRC method is too restricted to achieve both. TDDFT is in principle exact; however, going beyond the LRC approach is very challenging: we may need to better understand the role of the wings and body of the xc kernel, and the frequency dependence of the xc kernel may have to be taken into account. Alternatives beyond pure TDDFT, such as hybrid functionals, are therefore very promising. Such methods are currently under development.

Acknowledgements.

This work was supported by NSF grant No. DMR-1408904. The computation for this work was performed on the high performance computing infrastructure provided by Research Computing Support Services at the University of Missouri-Columbia. We thank Lucia Reining and Sangeeta Sharma for discussions and for providing valuable comments on the manuscript.