---
abstract: 'We present a framework for obtaining reliable solid-state charge and optical excitations and spectra from optimally-tuned range-separated hybrid density functional theory. , allows for accurate prediction of exciton binding energies. We demonstrate our approach through calculations of one- and two-particle excitations in pentacene, a molecular semiconducting crystal, where our work is in excellent agreement with experiments and prior computations. We further show that with one adjustable parameter, , this method accurately predicts band structures and optical spectra of silicon and lithium flouride, prototypical covalent and ionic solids. Our findings indicate that for a broad range of extended bulk systems, this method may provide a computationally inexpensive alternative to many-body perturbation theory, opening the door to studies of materials of increasing size and complexity.'
author:
- 'Sivan Refaely-Abramson'
- Manish Jain
- Sahar Sharifzadeh
- 'Jeffrey B. Neaton'
- Leeor Kronik
bibliography:
- 'tdrsh.bib'
title: |
    Solid-state optical absorption from optimally-tuned time-dependent\
    range-separated hybrid density functional theory
---

Many solid-state systems exhibit strong excitonic effects, notably an optical excitation spectrum that is affected substantially by interaction between excited electron and hole quasiparticle states. The nature of this electron-hole, or [*excitonic*]{}, interaction is of central importance for a variety of applications in, e.g., optoelectronics and photovoltaics [@Savoie2014]. Nevertheless, its accurate theoretical prediction remains a challenging task. It is common to account for such interactions within the framework of ab initio many-body perturbation theory, in which single-particle excitations are well-predicted from Dyson’s equation, typically solved within the GW approximation [@Hedin1965; @Hybertsen1986], and two-particle excitations are well-predicted using the Bethe-Salpeter equation (BSE) [@Rohlfing1998; @*Rohlfing2000; @Strinati1982; @*Strinati1984].

Current GW-BSE calculations are highly demanding and therefore presently impose significant practical limits on the calculated system size and complexity. Density functional theory (DFT), in both its time-independent [@DreizlerGross; @ParrYang] and time-dependent (TDDFT) [@Marques2012; @Casida1995; @Burke2005; @Ullrich_book] forms, is considerably more efficient computationally. However, common (semi-)local approximations to both DFT and TDDFT suffer from serious deficiencies which have precluded their use as a viable alternative to GW-BSE in the prediction of excitonic properties [@Onida2002]. First, quasi-hole and quasi-electron excitation energies are generally underestimated and overestimated, respectively, by the DFT Kohn-Sham eigenvalue spectrum [@Kummel2008; @Kronik2012]. While the same functionals often perform better in the prediction of optical excitation energies of isolated molecular systems, the Kohn-Sham gap is typically similar to the optical gap [@Salzner1997; @Chong2002; @Baerends2013; @Kronik2012; @Kronik2014]. In any case, they still fail in the solid-state limit [@Onida2002; @Ullrich_book; @Izmaylov2008; @Sharma2014; @Ullrich2014]. Therefore, neither one- nor two-particle excitations are well-predicted in the solid-state, and hence the nature of excitons or their binding energies are not obtained.

The failure of semi-local functionals in predicting solid-state absorption spectra has been traced back to an incorrect description of the long-range electron-electron and electron-hole interaction, manifested by the absence of a $1/q^2$ contribution [@Gonze1995; @*Ghosez1997; @Kim2002] to the interaction, where $q$ is a wavevector in the periodic system. Several ingenious schemes for overcoming this deficiency have been suggested, including the use of an exchange-correlation kernel of the form $f_{xc}(r,r')=-\alpha/(4\pi |r-r'|)$, where $\alpha$ is a system-dependent empirical parameter [@Reining2002; @Botti2004]; a static approximation to the exchange-correlation kernel based on a jellium-with-gap model [@Trevisanutto2013]; a “bootstrap” parameter-free kernel, achieved using self-consistent iterations of the random phase approximation (RPA) dielectric function [@Sharma2011; @Sharma2014]; a related “guided iteration” RPA-bootstrap kernel [@Rigamonti2015]; and the Nanoquanta kernel [@Onida2002; @Reining2002; @Sottile2003; @Marini2003; @*Adragna2003], derived by constructing the exchange-correlation kernel from an approximate solution to the BSE. Each correction provides a major step forward. However, none is a fully DFT-based solution, as single quasiparticle excitations are obtained from GW, RPA, a DFT+U approach, or a scissors-shift correction.

A different path for enabling TDDFT calculations in the solid state is the use of (global or range-separated) hybrid functionals. These are still well within density functional theory, using the generalized Kohn-Sham (GKS) framework [@Seidl1996; @Kummel2008; @Kronik2012], and their non-local Fock-like exchange component assists in the inclusion of long-range contributions. Although the time-dependent GKS equations have yet to be formally derived, hybrid functionals are already widely used for calculating optical properties. For gas-phase molecules, hybrid functionals can improve optical excitation energies, although standard hybrids still do not provide for accurate single-particle excitation energies [@Salzner1997; @Kummel2008; @Kronik2012; @Kronik2014]. TDDFT using the Heyd-Scuseria-Ernzerhof (HSE) short-range hybrid functional [@Heyd2003; @*Heyd2006], where non-local exchange is introduced only in the short-range, can improve the absorption spectra of semiconductors and insulators [@Paier2008], although some discrepancies remain. However, the HSE functional still does not provide the desired long-range non-local contribution. The B3LYP hybrid functional [@Becke1993; @*Stephens1994], in which a global 20% fraction of exact-exchange is used, was recently shown to yield TDDFT optical spectra for semiconductors in good agreement with experiment [@Bernasconi2011; @*Tomic2014], . Although in this case a non-local contribution to the kernel tail does exist, it is global and parameterized for a finite set of small organic molecules. However, global and short-range hybrid functionals were shown to be insufficient predictors of band-structures in solid-state systems [@Jain2011], notably for molecular crystals [@Refaely-Abramson2013] where excitonic effects are strong. Recently, Yang et al. [@Yang2015] suggested a screened exact-exchange (SXX) approach, in which the local part of the hybrid calculation is set to zero and the time-dependent Hartree-Fock exchange is scaled down non-empirically per system by using the inverse of the dielectric constant, based on a ground state obtained from a scissor-corrected local density approximation (LDA) calculation. Again, this led to improved performance for more strongly bound excitons.

![(Color online) (a) Band-structure (left) and density of states (right) of the pentacene solid, calculated using LDA (gray, dashed lines), G$_0$W$_0$@LDA (red, dashed lines), and OT-SRSH (black, solid lines). For all methods, the middle of the bandgap is shifted to zero. (b) The imaginary part of the dielectric function of the pantacene solid, with incident light polarization averaged over the $a$, $b$, and $c$ main unit-cell axes, calculated using TDLDA (gray, dashed lines), G$_0$W$_0$/BSE (red, dashed lines), and TD-OT-SRSH (black, solid lines). For visualization purposes, the leading absorption feature (between 0.5 to 2.5 eV) was multiplied by a factor of 10 with all computational methods used. The OT-SRSH and TD-OT-SRSH results were obtained using the parameters $\gamma=0.16$ Bohr$^{-1}$, $\alpha=0.2$, and $\varepsilon=3.6$. For computational details and convergence information, see the SI.[]{data-label="fig_pen"}](pen_DOS_abs_2h-averageLS-wan.pdf)

Ideally, we seek a DFT-based method where accurate one- and two-particle excitations can be read directly off of the eigenvalues of the time-independent (G)KS and the linear-response time-dependent (G)KS equations, respectively, using a single exchange-correlation functional, from which a consistent exchange-correlation potential and kernel are derived. This challenge is not met by any of the above-surveyed methods. Recently, it was met for gas-phase systems, using the optimally-tuned range separated hybrid functional (OT-RSH) approach [@Baer2010; @Stein2010; @Kronik2012], where the long- and short- range fraction of Fock-exchange is tuned non-empirically so as to obey rigorous physical constraints. This approach, elaborated further below, was shown to yield excellent fundamental and optical gaps for molecules [@Refaely-Abramson2011; @RSHpapers]. More recently, it has been generalized so as to provide accurate single-particle excitations for both molecules [@Refaely-Abramson2012; @Egger2014; @Nguyen2015] and molecular crystals [@Refaely-Abramson2013; @Lueftner2014], and it was shown to capture gap renormalization in molecular solids [@Refaely-Abramson2013]. Can this approach, then, resolve the long-standing challenge of providing an accurate one- and two-particle excitation spectrum in solid-state systems fully within the framework of (TD)DFT?

In this Letter, we present a solid-state OT-RSH approach that achieves just that. . We non-empirical calculations for pentacene, a prototypical molecular crystal provides excellent agreement with GW-BSE calculations. Furthermore, with one empirical parameter - set to reproduce the known fundamental gap - we again achieve results that are comparable with both GW-BSE and experiment for bulk silicon and LiF. The approach therefore emerges as promising for photoelectron and optical properties accurately and efficiently for a broad range of extended systems.

In the range-separated hybrid approach, the Coulomb interaction is range-partitioned [@Leininger1997], e.g., via: [@Yanai2004] $$\frac{1}{r}=\frac{\alpha+\beta\mathrm{erf}(\gamma r)}{r}+\frac{1-[\alpha+\beta\mathrm{erf}(\gamma r)]}{r}, \label{alphabetagamma}$$ where $\alpha$, $\beta$ and $\gamma$ are parameters, and $r$ is the inter-electron coordinate. The exchange expression corresponding to the first term on the right-hand side of Eq. (\[alphabetagamma\]) is then treated as in Hartree-Fock theory; the exchange expression corresponding to the second term is treated within the Kohn-Sham framework, typically using the LDA or the generalized gradient approximation (GGA). In Eq. (\[alphabetagamma\]), $\gamma$ is the range-separation parameter, i.e., it controls the range at which each of the terms dominates; $\alpha$ and $\beta$ dictate the limiting behavior of the Fock-like exchange, which tends to $\alpha/r$ for $r \rightarrow 0$ and to $(\alpha + \beta )/r$ for $r \rightarrow \infty$. The resulting exchange-correlation energy is of the form: $$\begin{split}   
E_{xc}^{RSH}=(1-\alpha)E_{\mathrm{KS}x}^{SR}+\alpha E_{xx}^{SR}+
(1-(\alpha+\beta))E_{\mathrm{KS}x}^{LR}+\\
(\alpha+\beta)E_{xx}^{LR}+E_{\mathrm{KS}c},
\label{SRSH}
\end{split}$$ where $\mathrm{KS}x$ and $\mathrm{KS}c$ denote (semi-)local KS exchange and correlation respectively, and $xx$ is a Fock-like exchange. SR and LR label short- and long-range terms, in which the Coulomb interaction is scaled using error functions. Within the GKS framework, the potential corresponding to the (semi-)local energy components is then obtained as a functional derivative, whereas the potential corresponding to the $xx$ energy components is obtained as a non-local Fock-like operator.

Here we generalize this approach to the time-dependent case, within the usual linear-response formalism of Casida [@Casida1995; @Tretiak2003] and the Tamm-Dancoff approximation. This is achieved by coupling GKS electron-hole pairs via an exchange-correlation kernel. This leads to an eigenvalue equation in which the eigenvalues are related to optical excitation energies and the eigenvectors can be used to compute oscillator strengths, so that the complete optical absorption spectrum can be computed. The part of the TDDFT kernel originating from the Hartree and (semi-)local Kohn-Sham potential in the ground-state DFT can be expressed as: $${\langle ai|}\left[\frac{1}{|r-r'|}+(1-\alpha)f_{xc}^{SR}+(1-\alpha-\beta)f_{xc}^{LR}\right]{|bj\rangle} \label{ker-l}$$ where $f_{xc}^{SR}=\frac{\delta V_{xc}^{SR}}{\delta n(r)}\delta(r-r')$, $f_{xc}^{LR}=\frac{\delta V_{xc}^{LR}}{\delta n(r)}\delta(r-r')$, and $V_{xc}^{SR},V_{xc}^{LR}$ are the short- and long-range contributions of the (semi-)local KS exchange-correlation potential; and where $a,b$ and $i,j$ denote occupied and unoccupied states, respectively. The non-local exchange potential in the ground-state GKS leads to an additional term in the kernel, of the form $$-{\langle ab|}\left[ \alpha\frac{\mathrm{erfc(\gamma(|r-r'|)}}{|r-r'|}+(\alpha+\beta)\frac{\mathrm{erf(\gamma(|r-r'|)}}{|r-r'|} \right] {|ij\rangle} \label{ker-nl}$$ (see supplementary information (SI) for additional formal details).

As discussed above, optimal-tuning of the RSH parameters was shown to be crucial for achieving accurate description of molecular single-particle and optical excitations. This tuning procedure is challenging in the solid-state, as it involves a calculation of the system’s ionization potential and electron affinity from total energy differences, a problematic procedure for periodic systems (see [@Vlcek2015], and references therein). For the special case of molecular crystals, however, it was shown [@Refaely-Abramson2013] that predictive bandstructures can be achieved if $\gamma$ and $\alpha$ are optimally-tuned so as to obey the ionization potential theorem for an isolated molecule, with $\beta$ chosen such that $\alpha + \beta = 1/\varepsilon_0$, where $\varepsilon_0$ is the scalar dielectric constant, itself computed from first principles. We refer to this procedure as the optimally-tuned [*screened*]{}-RSH (OT-SRSH) approach.

To test the efficacy of this approach for optical properties, we examine pentacene, a molecular semiconducting crystal of extreme interest in organic electronics and photovoltaics. For pentacene, the optimal tuning parameters are found to be $\alpha=0.2$, $\beta=0.08$ (corresponding to $\varepsilon_0=3.6$), and $\gamma=0.16$ Bohr$^{-1}$ [@Refaely-Abramson2013]. With these parameters we construct the appropriate time-independent (Eq. (\[SRSH\])) and time-dependent (Eqs. (\[ker-l\]),(\[ker-nl\])) equations. These are then solved with the PARATEC [@Ihm1979] and BerkeleyGW [@Deslippe2012] codes, which we modified to handle range-separated hybrids (see SI for a complete discussion).

![image](Si_LiF_multiplot-midgap-v2.pdf)

Figure \[fig\_pen\] (a) shows the resulting band structure and density of states (DOS) calculated from LDA, OT-SRSH, and GW eigenvalues. The band gap predicted by OT-SRSH is 1.9 eV, in good agreement with the GW gap of 2.1 eV. As expected, both improve drastically on the LDA gap, which is 0.7 eV. Save for the slight difference in gap values, the OT-SRSH and GW band-structure and DOS are remarkably similar, demonstrating for this system that the eigenvalues of the OT-SRSH method are quantitatively useful approximations to single-quasiparticle excitation energies, as they are close to GW quasiparticle energies for several eV away from the band edges. Figure \[fig\_pen\] (b) presents the resulting imaginary part of the dielectric function, Im($\varepsilon$), calculated using TDLDA, TD-OT-SRSH and BSE, with the response averaged over incident light polarization along the $a$, $b$ and $c$ directions of the pentacene lattice (see SI for additional details). Importantly, our TD-OT-SRSH results do as well as GW-BSE in providing the expected picture of excitonic binding, i.e., comparison of the results with and without electron-hole interactions suggests an exciton binding energy of 0.45 eV from both methods. Furthermore, the BSE result agrees well with previously reported ones, e.g. those of Refs. [@Tiago2003; @Cudazzo2012; @Sharifzadeh2013; @*Sharifzadeh2012-2]. The first singlet excitation is predicted to be at 1.46 eV in TD-OT-SRSH and 1.64 eV in BSE. This small quantitative difference is essentially within the expected accuracy of either calculation. Moreover, this difference is consistent with that computed at the single-particle excitation level, and likely is inherited from it. In a similar manner, all presented optical excitation energies resulting from TD-OT-SRSH and BSE are very similar, with remaining differences within the desired accuracy. However, there are differences in the oscillator strength in part of the spectrum, notably at the first peak between 1.5 eV and 2.5 eV. These differences are primarily associated with the $a$-axis direction of incident light polarization. However, the scale of these differences is small when compared to the entire spectrum. (Note that this spectral range is enhanced by a factor of 10 in Fig.\[fig\_pen\](b) to be visible.) Overall, then, the TD-OT-SRSH spectrum, while not identical to that of BSE, is in very good agreement with it.

We now turn to two prototypical covalent and ionic bulk solids: Si and LiF, respectively. For these systems, we cannot use the original tuning procedure, because it was designed for systems where inter-molecular hybridization is small. However, it still possible to set $\alpha=0.2$ as a universally useful amount of short-range exact exchange [@Rohrdanz2009; @Egger2014] and demand $\alpha+\beta=1/\varepsilon_0$, as before. We are then left with only one parameter, the range-separation parameter, $\gamma$. Here, we simply choose $\gamma$ so as to obtain the fundamental gap. .

As shown in the top panel of Fig \[fig\_si\](a), the SRSH Si bandstructure is fully comparable to the canonical empirical pseudopotential work of Chelikowsky and Cohen [@CC1976], as well as to GW results, whereas with LDA the gap, as expected, is too small. From the bottom panel we see that, in agreement with previous work [@Reining2002; @Botti2004], the TDLDA spectrum does not reproduce experiment satisfactorily. However, once again the TDSRSH result does almost as well as GW-BSE in predicting the experimental results, without any scissor shift or other correction operators at the single-quasiparticle level. Qualitatively, the TDSRSH lineshape is indistinguishable from that of GW-BSE. Quantitatively, the peak positions are similar to within 0.05 eV. A similar picture emerges for LiF, as both single-particle and two-particle excitations are very close to those of GW/BSE. Specifically, both TDSRSH and BSE predict the first and largest excitation at 12.6 eV, very close to the experimental value at 12.75 eV, and the overall spectral shape is satisfying. Here too, then, SRSH and TDSRSH are shown to capture excitonic effects, even though the Si excitons are much more delocalized and weakly bound compared to the case of pentacene, and the LiF excitons are strongly-bound and known to be highly challenging for many of the TDDFT methods discussed in the introduction [@Rigamonti2015; @Sharma2014; @Yang2015; @Trevisanutto2013].

What are the physical origins of this success of the TDSRSH approach in the solid state? First, the ground-state calculation, being a [*generalized*]{} Kohn-Sham one, is capable in principle of describing quasi-particle excitations owing to the non-local potential operator. Range-separation combined with long-range dielectric screening allows us to fulfill asymptotic potential constraints while retaining the crucial balance of short-range exchange and correlation components, thereby making such prediction sufficiently accurate also in practice. Second, the asymptotic form of the TDSRSH kernel generates the desired non-local $1/q^2$ contribution by construction. Furthermore, it is already scaled correctly via the non-empirical $1/\varepsilon_0$ parameter. Third, the non-locality of the TDSRSH kernel also alleviates the need for frequency-dependence, which would be necessary for bound exciton prediction with (semi-)local exchange-correlation kernels [@Sharma2014]. For these reasons, highly accurate single-quasiparticle and two-quasiparticle properties of extended systems can be obtained in a predictive manner within the unifying framework of DFT at computationally-modest cost. Last but not least, Onida [*et al.*]{} [@Onida2002] have already noted some time ago that “both the Green‘s functions and the TDDFT approaches profit from mutual insight.’’ Here, we believe we achieve an important milestone towards that vision: on the one hand, the above work can be justified entirely from (TD)DFT reasoning. On the other hand, it is clear that the work has been motivated by the need to achieve the elegant quasi-particle picture obtained so naturally within many-body perturbation theory and that by achieving this goal we have created an effective simplified framework mimicking this picture.

In conclusion, we have presented a new approach for quantitative determination of single- and two-particle excitations in solids, based on range-separated hybrid density functional theory. The approach that are in quantitative agreement with those of many-body perturbation theory in the GW-BSE approximation, for three prototypical systems: Si, a covalent solid, LiF, an ionic solid, and pentacene, a molecular solid. In particular, it fully captures excitonic interactions for both strongly and weakly bound excitons. We envision that it could emerge as a useful low-cost substitute to many-body perturbation theory, as well as provoke further developments within density functional theory.

This work was supported by the European Research Council, the Israel Science Foundation, the United States-Israel Binational Science Foundation, the Helmsley Foundation, and the Wolfson Foundation. S.R.A. is supported by an Adams fellowship of the Israel Academy of Sciences and Humanities. S.S. was partially supported by the Scientific Discovery through Advanced Computing (SciDAC) Partnership program funded by US Department of Energy, Office of Science, Advanced Scientific Computing Research and Basic Energy Sciences. J.B.N. was supported by the US Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering (Theory FWP) under Contract No. DE-AC02- 05CH11231. The work performed at the Molecular Foundry was also supported by the Office of Science, Office of Basic Energy Sciences, of the US Department of Energy. We thank the National Energy Research Scientific Computing center for computational resources.
