# Comparison of long-range corrected kernels and range-separated hybrids for
excitons in solids

Rita Maji  Dipartimento di Scienze e Metodi dell’Ingegneria, Università di
Modena e Reggio Emilia,Via Amendola 2 Padiglione Tamburini , I-42122 Reggio
Emilia, Italy  Elena Degoli  Dipartimento di Scienze e Metodi dell’Ingegneria,
Università di Modena e Reggio Emilia and Centro Interdipartimentale
En$\&$Tech, Via Amendola 2 Padiglione Morselli, I-42122 Reggio Emilia, Italy
Centro S3, Istituto Nanoscienze-Consiglio Nazionale delle Ricerche (CNR-
NANO),Via Campi 213/A, 41125 Modena, Italy
Centro Interdipartimentale di Ricerca e per i Servizi nel settore della
produzione, stoccaggio ed utilizzo dell’Idrogeno H$2$–MO.RE., Via Università
4, 41121 Modena, Italy  Monica Calatayud  Laboratoire de Chimie Théorique,
Sorbonne Université and CNRS F-75005 Paris, France  Valérie Véniard
<EMAIL_ADDRESS>Laboratoire des Solides Irradiés CNRS,
CEA/DRF/IRAMIS, École Polytechnique, Institut Polytechnique de Paris, F-91128
Palaiseau, France and European Theoretical Spectroscopy Facility (ETSF)
Eleonora Luppi<EMAIL_ADDRESS>Laboratoire de Chimie
Théorique, Sorbonne Université and CNRS F-75005 Paris, France

###### Abstract

The most accurate theoretical method to describe excitons is the solution of
the Bethe-Salpeter equation in the GW approximation (GW-BSE). However, because
of its computation cost time-dependent density functional theory (TDDFT) is
becoming the alternative approach to GW-BSE to describe excitons in solids.
Nowadays, the most efficient strategy to describe optical spectra of solids in
TDDFT is to use long-range corrected exchange-correlation kernels on top of GW
or scissor-corrected energies. In recent years, a different strategy based on
range-separated hybrid functionals started to be developed in the framework of
time-dependent generalised Kohn-Sham density functional theory (TDGKSDFT).
Here, we compare the performance of long-range corrected kernels with range-
separated hybrid functionals for the description of excitons in solids. This
comparison has the purpose to weight the pros and cons of using range-
separated hybrid functionals, giving new perspectives for theoretical
developments of these functionals. We illustrate the comparison for the case
of Si and LiF, representative of solid state excitons.

##  I Introduction

Excitons play a central role in the optical properties of materials for
optoelectronics, photovoltaics and photocatalysis applications.
PhysRevB.95.035125 ; PhysRevLett.116.066803 ; acs.jpclett.1c00543 ; Dong2020 ;
science.abm8511  Excitons are usually described as an electron-hole pair and
they are classified as Frenkel excitons (bound) localised at the atomic sites
and Mott-Wannier excitons (continuum) delocalised over the atomic unit cells.
onid+02rmp

An accurate description of the excitonic effect in the optical properties of
materials is still very challenging, and nowadays the most accurate
theoretical method to describe excitons is the solution of the Bethe-Salpeter
equation (BSE) in the GW approximation (GW-BSE). However, the computational
cost of GW-BSE can be very high. reboCRC2013 ; onid+02rmp

Time-dependent density functional theory (TDDFT) is an alternative approach to
BSE to describe excitons. TDDFT is mathematically simpler than BSE which makes
TDDFT computationally more efficient. The key quantity of TDDFT is the
exchange-correlation kernel $f_{\text{xc}}$ which needs to be approximated.
Nowadays, none of the proposed approximations for $f_{\text{xc}}$ reaches the
BSE accuracy with the only exception of the Nanoquanta exchange-correlation
kernel which, however, makes TDDFT as expensive as BSE. PhysRevLett.91.056402
; PhysRevB.68.165108 ; PhysRevLett.91.256402

The most efficient strategy to describe optical spectra of solids in TDDFT is
to use long-range corrected exchange-correlation kernels ($1/q^{2}$ in the
long wavelength limit) on top of GW or scissor-corrected energies.
SottIJQC2005 ; gaurJCP2019  The first long-range corrected (LRC) kernel was
derived by Reining et al. rein+02prl  through a comparison with BSE. LRC is an
empirical kernel which requires a material-dependent parameter. For a large
class of semiconductors, the parameter depends on the inverse dielectric
constant in a simple way. bott+04prb  This kernel demonstrated to correctly
describe only continuum excitons. Since then, a number of nonempirical
exchange-correlation kernels, corrected for the long-range interactions, have
been proposed in literature. PhysRevLett.114.146402 ; PhysRevLett.107.186401 ;
PhysRevB.87.205143 ; PhysRevB.95.205136 ; PhysRevLett.107.186401 ;
PhysRevB.95.205136 ; PhysRevLett.127.077401  Different efficient bootstrap
kernels have been developed, which describe both continuum and strong excitons
in insulators and semiconductors. PhysRevLett.114.146402 ;
PhysRevLett.107.186401  A kernel based on the jellium-with-gap model (JGM) was
proposed PhysRevB.87.205143  to describe both continuum and strong excitons in
different materials. However, despite the success of these kernels to describe
optical properties, it has been found that they cannot predict accurate
exciton binding energies. PhysRevB.95.205136  In order to recover both
properties simultaneously an empirically scaled bootstrap kernel has been
proposed. PhysRevB.95.205136

In recent years, a different strategy based on range-separated hybrid
functionals started to develop. PhysRevB.78.121201 ; PhysRevB.92.081204 ;
PhysRevB.92.035202 ; PhysRevResearch.2.013091 ; B812838C ; ZapJCP2019 ;
Rebo2013MolPhys  Range-separated hybrid functionals rely on the splitting of
the Coulomb electron-electron interaction $w_{\text{ee}}=1/r$ into a long-
range ($w_{lr}$) and a short-range ($w_{sr}$) contributions by a tunable
parameter $\mu$ which controls the range separation. Starting from this simple
idea different schemes exist which simulate the $w_{sr}$ by nonlocal Hartree-
Fock (HF) exchange and the $w_{lr}$ by (semi-)local density-functional theory
exchange functional or vice versa. HeyScu-JCP-04b ; AleScu16JPCL ; B812838C ;
PhysRevB.78.121201  This methodology has been developed in the framework of
the time-dependent generalised Kohn-Sham density functional theory (TDGKSDFT).
An important property of TDGKSDFT is that the exchange-correlation potential
and the kernel are fully consistent with the choice of the exchange-
correlation energy, being its first and second functional derivative with
respect to the density. This consistency is not provided by TDDFT with long-
range corrected kernels.

The Heyd-Scuseria-Ernzerhof (HSE) range-separated hybrid functional
PhysRevB.78.121201  was used to reproduce optical spectra of semiconductors
and insulators. The spectra improve with respect to semilocal functionals for
semiconductors showing a very good agreement with experiments, except for
large gap insulators. A Coulomb attenuating method (CAM) range-separation was
also proposed PhysRevB.92.081204  showing excellent agreement with both
semiconductors and insulators. CAM belongs to range-separated hybrid (RSH)
functionals. The main difficulty of these range-separated approaches is to
find a general criterium valid for different types of materials for the range-
separation parameter $\mu$ and for those parameters that control the weight of
nonlocal HF exchange and DFT exchange. To solve this problem, in
Ref.(PhysRevB.92.081204 ) the authors optimally tuned the parameters in order
to reproduce physical constraints.

Another promising approach is to screen a fraction of the nonlocal HF exchange
with the inverse dielectric constant and not to include semilocal exchange-
correlation functional. PhysRevB.92.035202  In this case, however the
calculation is not fully consistent as it is obtained from a scissor-corrected
local density approximation (LDA) calculation. The same approach was also
proposed combining the full range nonlocal HF exchange with a fraction of
local exchange functional and correlation. PhysRevResearch.2.013091

The use of range-separated hybrid functionals seems to be very promising and
open new perspectives for the calculation of optical spectra of solids.
However, the calculation of nonlocal HF exchange is computationally more
demanding than the standard TDDFT approach with long-range corrected kernels.
Moreover, a general rule valid for any materials concerning the choice of the
parameters needed in the calculations with hybrid functionals has not been
found yet.

In this paper, we compare the performance of long-range corrected kernels with
range-separated hybrid functionals for the description of excitons in solids.
This comparison has the purpose to weight the pros and cons of using range-
separated hybrid functionals, giving new perspectives for theoretical
developments of these functionals. Concerning long-range corrected kernels we
studied the LRC rein+02prl , scalar RPA bootstrap (RPA-BO)
PhysRevLett.114.146402  and JGM. PhysRevB.87.205143  Concerning hybrid
functionals, we investigated the short-range nonlocal HF exchange with and
without semilocal exchange-correlation PBE functional PhysRevLett.77.3865 . We
call this two schemes respectively TDHFsr,μ;α and
TDHFsr,μ;αXC${}^{\text{PBE}}$. Moreover, in the discussion we also include the
hybrid scheme presented in Refs. PhysRevB.92.081204 ;
PhysRevMaterials.3.064603  which has a long-range nonlocal HF exchange
component. We illustrate the comparison for the case of Si and LiF,
representative of solid state excitons.

In Section II we compare the kernels of BSE, TDDFT with long-range corrected
and TDGKSDFT with range-separated hybrid functionals. Section III is devoted
to computational details, while in Section IV we present and discuss the
results. Conclusions are in Section V.

##  II Theory

The macroscopic dielectric tensor $\epsilon_{\text{M}}(\omega)$ is

$\varepsilon_{\text{M}}(\omega)=\lim_{{\bf q}\to
0}\frac{1}{\varepsilon^{-1}_{{\bf G}_{1}=0,{\bf G}_{2}=0,}({\bf q},\omega)}$
(1)

where $\varepsilon^{-1}_{{\bf G}_{1},{\bf G}_{2}}({\bf q},\omega)$ is the
inverse microscopic dielectric matrix written in terms of the reciprocal-space
lattice vectors ${\bf G}_{1}$ and ${\bf G}_{2}$ for a given wave-vector ${\bf
q}$ and frequency $\omega$. The case ${\bf G}_{1}={\bf G}_{2}=0$ indicates the
head element of the inverse microscopic dielectric matrix. Through the
calculation of
$\varepsilon_{2}(\omega)=\text{Im}[\varepsilon_{\text{M}}(\omega)]$ the
absorption spectrum is obtained. onid+02rmp

There exist two different approaches which can be used in GW-BSE, TDDFT and
TDGKSDFT to obtain $\varepsilon_{\text{M}}(\omega)$. One approach is the
solution of the Dyson equation and the other approach is the solution of the
Casida’s equations. Converged spectra are identical within the two approaches.
PhysRevB.95.205136  Within these approaches different kernels are used
depending on the level of theory.

The kernel of the BSE in the GW approximation written in reciprocal space and
expressed in terms of the 4-space indices ($v{\mathbf{k}}$, $c{\mathbf{k}}$,
$v^{\prime}{\mathbf{k}}^{\prime}$ and $c^{\prime}{\mathbf{k}}^{\prime}$) of
the transition space is onid+02rmp ; PhysRevB.78.121201

$\Xi^{\text{GW-
BSE}}_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}=w_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}-W_{c{\mathbf{k}}c^{\prime}{\mathbf{k}}^{\prime},v{\mathbf{k}}v^{\prime}{\mathbf{k}}^{\prime}}(\omega).$
(2)

The first term is the Hartree contribution

$\displaystyle
w_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}=\lim_{{\bf q}\to
0}\sum_{{\bf G}_{1}{\bf G}_{2}}\frac{4\pi}{|{\bf q}+{\bf
G}_{1}|^{2}}\delta_{{\bf G}_{1}{\bf G}_{2}}\langle c{\mathbf{k}}|e^{i({\bf
q}+{\bf G}_{1})\cdot{\bf r}_{1}}|v{\mathbf{k}}\rangle\langle
v^{\prime}{\mathbf{k}}^{\prime}|e^{-i({\bf q}+{\bf G}_{2})\cdot{\bf
r}_{2}}|c^{\prime}{\mathbf{k}}^{\prime}\rangle,$ (3)

and the second term is the screened Coulomb interaction

$\displaystyle
W_{c{\mathbf{k}}c^{\prime}{\mathbf{k}}^{\prime},v{\mathbf{k}}v^{\prime}{\mathbf{k}}^{\prime}}(\omega)=4\pi\lim_{{\bf
q}\to 0}\sum_{{\bf G}_{1}{\bf G}_{2}}\frac{\varepsilon^{-1}_{{\bf G}_{1}{\bf
G}_{2}}({\bf q},\omega)}{|{\bf q}+{\bf G}_{1}|^{2}}\langle
c{\mathbf{k}}|e^{i({\bf q}+{\bf G}_{1})\cdot{\bf
r}_{1}}|c^{\prime}{\mathbf{k}}^{\prime}\rangle\langle
v^{\prime}{\mathbf{k}}^{\prime}|e^{-i({\bf q}+{\bf G}_{2})\cdot{\bf
r}_{2}}|v{\mathbf{k}}\rangle.$ (4)

where ${\bf q}$ is a vector in the first Brillouin zone, ${\bf G}_{1}$ and
${\bf G}_{2}$ are vectors of the reciprocal lattice.

Most of the GW-BSE calculations neglect the frequency dependence of the
inverse dielectric function $\varepsilon^{-1}_{{\bf G}_{2}{\bf G}_{1}}({\bf
q},\omega=0)$ which cause the GW-BSE kernel to be static.

The screened interaction $W$ has a long-range behaviour $1/q^{2}$ and is
attractive, opposite to the Hartree $w^{\text{ee}}$ interaction which is
repulsive. We observe that by neglecting $W$ in the kernel we recover the
random-phase approximation (RPA).

The kernel of the TDDFT is

$\Xi^{\text{TDDFT}}_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}=w_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}+f^{\text{xc}}_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}},$
(5)

where the first term is still the Hartree contribution and the second term is
the exchange-correlation kernel in the reciprocal and transition space defined
as

$\displaystyle
f_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}^{\text{xc}}=\lim_{{\bf
q}\to 0}\sum_{{\bf G}_{1}{\bf G}_{2}}f_{\text{xc},{\bf G}_{1}{\bf G}_{2}}({\bf
q})\langle c{\mathbf{k}}|e^{i({\bf q}+{\bf G}_{1})\cdot{\bf
r}_{1}}|v{\mathbf{k}}\rangle\langle v^{\prime}{\mathbf{k}}^{\prime}|e^{-i({\bf
q}+{\bf G}_{2})\cdot{\bf r}_{2}}|c^{\prime}{\mathbf{k}}^{\prime}\rangle.$ (6)

This quantity is expected to describe the electron correlations that in the
BSE are described by $W$. Therefore a good mathematical approximation for the
$f^{\text{xc}}$ should include the $1/q^{2}$ long-range behaviour and the
screening. Note also that, by comparing the matrix elements of $W$ in Eq.(4)
and $f^{\text{xc}}$ in Eq.(6), there is an exchange between $vk$ and
$c^{\prime}k^{\prime}$ between $W$ and $f^{\text{xc}}$. The structure of
$f^{\text{xc}}$ is similar to the Hartree contribution. In the case of
nonlocal HF exchange we obtain exactly the same structure of the matrix
elements of $W$ as it will be shown later.

The correct long-range behaviour and screening are described by the long-range
corrected kernels. The first of this type of kernels presented in literature
is the LRC kernel rein+02prl ; bott+04prb  defined as

$\displaystyle f^{\text{LRC}}_{\text{xc},{\bf G}{\bf G}^{\prime}}({\bf
q})=-\frac{\alpha^{\text{LRC}}}{|{\bf q}+{\bf G}^{\prime}|^{2}}\delta({\bf
G}^{\prime},{\bf G}).$ (7)

For materials with a small inverse dielectric constant $\varepsilon^{-1}_{0}$,
the parameter $\alpha^{\text{LRC}}$ can be approximated by
$\alpha^{\text{LRC}}=4.651\varepsilon^{-1}_{0}-0.213$. bott+04prb  This kernel
has demonstrated to be able to simulate continuum excitons but not strong
excitons rein+02prl ; bott+04prb ; PhysRevB.95.205136 .

Another long-range corrected kernel is the RPA-BO PhysRevLett.117.159702
which is defined as

$\displaystyle f^{\text{RPA-BO}}_{\text{xc}}({\bf
q})=\frac{\varepsilon^{-1}_{\text{RPA},00}w({\bf
q})}{1-1/\varepsilon^{-1}_{\text{RPA,00}}({\bf q},0)}.$ (8)

This kernel is scalar $({\bf G}=0$ and ${\bf G}^{\prime}=0$) and the screening
is given by the RPA inverse dielectric constant
$\varepsilon^{-1}_{\text{RPA},00}$. The RPA-BO kernel gives good results for
both continuum and strong excitons in semiconductors and insulators.
PhysRevLett.117.159702

The JGM kernel PhysRevB.87.205143 ; gaurJCP2019  based on the jellium-with-gap
model is another kernel with the correct $1/q^{2}$ behaviour and is defined as

$\displaystyle f^{\text{JGM}}_{\text{xc},{\bf G}{\bf G}^{\prime}}({\bf
q})=-4\pi\frac{B^{\prime}({\bf G}-{\bf G}^{\prime})}{|{\bf q}+{\bf
G}^{\prime}|^{2}}+4\pi\frac{H({\bf G}-{\bf G}^{\prime},{\bf
G}^{\prime})}{|{\bf q}+{\bf G}^{\prime}|^{2}}-\frac{D^{\prime}({\bf G}-{\bf
G}^{\prime})}{1+1/|{\bf q}+{\bf G}^{\prime}|^{2}}$ (9)

where $B^{\prime}$, $H$ and $D^{\prime}$ depend on the density and on the
electronic gap. The precise definition of the quantity is given in
Ref.PhysRevB.87.205143 . The JGM kernel gives also good results for both
continuum and strong excitons in semiconductors and insulators.
PhysRevB.87.205143 ; gaurJCP2019

The TDHFsr,μ;α kernel is

$\Xi^{\text{TDHF${}^{sr,\mu;\alpha}$}}_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}=w_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}-\alpha\,w^{\text{HF${}^{sr,\mu}$}}_{c{\mathbf{k}}c^{\prime}{\mathbf{k}}^{\prime},v{\mathbf{k}}v^{\prime}{\mathbf{k}}^{\prime}}$
(10)

and contains the Hartree term and the short-range nonlocal HF exchange

$\displaystyle
w^{\text{HF${}^{sr,\mu}$}}_{c{\mathbf{k}}c^{\prime}{\mathbf{k}}^{\prime},v{\mathbf{k}}v^{\prime}{\mathbf{k}}^{\prime}}=\lim_{{\bf
q}\to 0}\sum_{{\bf G}_{1},{\bf G}_{2}}\frac{4\pi}{|{\bf q}+{\bf
G}_{1}|^{2}}\delta_{{\bf G}_{1},{\bf G}_{2}}\langle
c{\mathbf{k}}v^{\prime}{\mathbf{k}}^{\prime}|e^{i({\bf q}+{\bf
G}_{1})\cdot({\bf r}_{1}-{\bf r}_{2})}\text{erfc}(({\bf r}_{1}-{\bf
r}_{2})\mu)|v{\mathbf{k}}c^{\prime}{\mathbf{k}}^{\prime}\rangle.$ (11)

screened by a parameter $\alpha$. In the case $\mu=0$ we obtain that
$w^{\text{HF${}^{sr,0}$}}$ is equal to nonlocal HF exchange $w^{\text{HF}}$
which is defined as

$\displaystyle
w^{\text{HF}}_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}=\lim_{{\bf
q}\to 0}\sum_{{\bf G}_{1}{\bf G}_{2}}\frac{4\pi}{|{\bf q}+{\bf
G}_{1}|^{2}}\delta_{{\bf G}_{1},{\bf G}_{2}}\langle c{\mathbf{k}}|e^{i({\bf
q}+{\bf G}_{1})\cdot{\bf r}_{1}}|c^{\prime}{\mathbf{k}}^{\prime}\rangle\langle
v^{\prime}{\mathbf{k}}^{\prime}|e^{-i({\bf q}+{\bf G}_{2})\cdot{\bf
r}_{2}}|v{\mathbf{k}}\rangle,$ (12)

and which has the same matrix form of the unscreened $W$ of the GW-BSE in
Eq.(4). In this case the kernel is
$\Xi^{\text{TDHF${}^{sr,0;\alpha}$}}=w-\alpha\,w^{\text{HF}}$. This kernel has
been proposed in Refs. PhysRevB.92.035202  under the name of screened-exact
exchange (SXX). In the case of $\mu\to\infty$ we obtain that
$w^{sr,\mu}_{\text{HF}}$ is equal to zero and
$\Xi^{\text{TDHF${}^{sr,\mu\to\infty;\alpha}$}}=w$ reduces to RPA.

The kernel TDHFsr,μ;αXC${}^{\text{PBE}}$ is

$\Xi^{\text{TDHF${}^{sr,\mu;\alpha}$XC${}^{\text{PBE}}$}}_{cv{\bf
k},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}=w_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}-\alpha\,w^{\text{HF${}^{sr,\mu}$}}_{c{\mathbf{k}}c^{\prime}{\mathbf{k}}^{\prime},v{\mathbf{k}}v^{\prime}{\mathbf{k}}^{\prime}}+(1-\alpha)f^{\text{x,PBE}}_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}+f^{\text{c,PBE}}_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}.$
(13)

The same kernel is proposed in Ref.PhysRevResearch.2.013091 . In the case
$\mu=0$ we obtain
$\Xi^{\text{TDHF${}^{sr,0;\alpha}$XC${}^{\text{PBE}}$}}=w-\alpha\,w^{\text{HF}}+(1-\alpha)f^{\text{x,PBE}}+f^{\text{c,PBE}}$.
In the case of $\mu\to\infty$ we obtain
$\Xi^{\text{TDHF${}^{sr,\mu\to\infty;\alpha}$XC${}^{\text{PBE}}$}}=w+(1-\alpha)f^{\text{x,PBE}}+f^{\text{c,PBE}}$.

In the discussion, we also show the comparison with the range-separated CAM
proposed in Refs. PhysRevMaterials.3.064603 ; PhysRevB.92.081204  and which
also includes a fraction of nonlocal long-range HF exchange. In this case the
CAM kernel is

$\displaystyle\Xi^{\text{TDCAM}^{sr,\mu;\alpha,\beta}}_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}=w_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}-\alpha\,w^{\text{HF${}^{sr,\mu}$}}_{c{\mathbf{k}}c^{\prime}{\mathbf{k}}^{\prime},v{\mathbf{k}}v^{\prime}{\mathbf{k}}^{\prime}}-(\alpha+\beta)w^{\text{HF${}^{lr,\mu}$}}_{c{\mathbf{k}}c^{\prime}{\mathbf{k}}^{\prime},v{\mathbf{k}}v^{\prime}{\mathbf{k}}^{\prime}}$
$\displaystyle+(1-\alpha)f^{\text{x,PBE${}^{sr,\mu}$}}_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}+(1-\alpha-\beta)f^{\text{x,PBE${}^{lr,\mu}$}}_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}+f^{\text{c,PBE}}_{cv{\mathbf{k}},c^{\prime}v^{\prime}{\mathbf{k}}^{\prime}}.$
(14)

This approach requires an additional parameter $\beta$ calculated as
$\alpha+\beta=1/\varepsilon_{0}$ where $\varepsilon_{0}$ is the material’s
dielectric constant. PhysRevMaterials.3.064603 ; PhysRevB.92.081204  In the
case $\mu=0$ we obtain
$\Xi^{\text{TDCAM${}^{sr,0;\alpha,\beta}$}}=w-\alpha\,w^{\text{HF}}+(1-\alpha)f^{\text{x,PBE}}+f^{\text{c,PBE}}$.
In the case of $\mu\to\infty$ we obtain
$\Xi^{\text{TDCAM${}^{sr,\mu\to\infty;\alpha,\beta}$}}=w-(\alpha+\beta)w^{\text{HF}}+(1-\alpha-\beta)f^{\text{x,PBE}}+f^{\text{c,PBE}}$.

##  III Computational Details

The TDDFT optical spectra with long-range corrected exchange-correlation
kernels have been calculated with DP DP  and 2light lupp+10jcp ; luppi_ab_2010
codes interfaced with the norm-conserving (NC) pseudopotentials and plane-wave
basis set ABINIT code. abinit1 ; abinit2

The optical spectra in GW-BSE and TDGKSDFT with range-separated hybrid
functionals have been calculated with the plane-wave based Vienna Ab initio
Simulation Package (VASP) with projector augmented-wave (PAW)
pseudopotentials. Hafner ; Kresse

In the case of NC pseudopotentials we used an energy cutoff of 10 Ha for Si
and 40 Ha for LiF, while for PAW pseudopotentials we used an energy cutoff of
9 Ha for Si and 16 Ha for LiF.

All the calculations have been performed using the experimental lattice
parameter 5.430Å for Si and 4.026Å for LiF. We used the experimental lattice
parameter in order to be consistent between the different theoretical methods
for the spectra comparison.

The convergence parameters for the optical spectra are reported in Table 1.
Note that for TDDFT calculations we used shifted k-points grids, while for GW-
BSE and TDGKSDFT we averaged the dielectric function over multiple k-points
shifted grids.PhysRevB.78.121201 ; Hafner ; Kresse  111In the case of GW-
BSE/TDGKSDFT, the construction of the Hamiltonian scales as
$N_{k}N^{2}_{v}N^{2}_{c}N_{G}$ where $N_{k}$ is the number of k-points in the
Brillouin zone, $N_{v}$ is the number of valence bands, $N_{c}$ is the number
of conduction bands and $N_{G}$ is the number of G-vectors. However, except
for very large systems, the main limiting factor usually comes from the
diagonalization of the Hamiltonian which scales cubically with the matrix rank
$N_{rank}^{3}$ and $N_{rank}=N_{k}N_{v}N_{c}$. TDDFT requires the
diagonalization for each k-point of a matrix of rank $N_{rank}=(N_{v}+N_{c})$
and the evaluation of the response function scales as
$N_{k}N_{v}N_{c}N_{G}^{2}$.

A broadening of 0.05 eV for all optical spectra have been used.

Table 1: Si and LiF convergence parameters.  | Material | k-points | empty bands | G-vectors
---|---|---|---|---
TDDFT (NC) | Si | 30$\times$30$\times$30 (shifted) | 4 | 89
| LiF | 32$\times$32$\times$32 (shifted) | 26 | 89
TDGKSDFT (PAW) | Si | 8$\times$8$\times$8 (29 shifted) | 16 | 163
| LiF | 8$\times$8$\times$8 (29 shifted) | 32 | 294
GW-BSE (PAW) | Si | 8$\times$8$\times$8 (29 shifted) | 12/128 | 150
| LiF | 8$\times$8$\times$8 (29 shifted) | 12/160 | 270

##  IV Results and Discussion

The goal of this work is to compare TDDFT and TDGKSDFT to describe optical
spectra of solids. As in our calculations we used both NC and PAW
pseudopotentials, we have first analysed the electronic structures of Si and
LiF.

For TDGKSDFT, the electronic structure was calculated with the HSE exchange-
correlation functional
$E^{\text{HSE}^{\mu;\alpha}}_{\text{xc}}=\alpha\,E^{\text{HF}^{\text{sr},\mu}}_{\text{x}}+(1-\alpha)E^{\text{PBE}^{\text{sr},\mu}}_{\text{x}}+E^{\text{PBE}^{\text{lr},\mu}}_{\text{x}}+E_{\text{c}}^{\text{PBE}}$,
where the long-range PBE exchange functional is
E${}^{\text{PBE}^{\text{lr},\mu}}_{\text{x}}$ = E${}_{\text{x}}^{\text{PBE}}$
\- E${}^{\text{PBE}^{\text{sr},\mu}}_{\text{x}}$. HeyScu-JCP-04b  Considering
$\mu\to\infty$ we have $E^{\text{HSE}^{\mu\to\infty;\alpha}}_{\text{xc}}$ =
$E^{\text{PBE}}_{\text{xc}}$ and, instead, considering $\mu=0$ we have
$E^{\text{HSE}^{0;\alpha}}_{\text{xc}}$ =
$\alpha\,E^{\text{HF}}_{\text{x}}+(1-\alpha)E^{\text{PBE}}_{\text{x}}+E^{\text{PBE}}_{\text{c}}$.
Following the optimally tuned strategy, we chose the parameters $\mu$ and
$\alpha$ in order to have a good agreement with the GW gaps. For Si we used as
($\mu$;$\alpha$) : (0.2;0.25), (0.3;0.25), (0.3;0.3) and (0.0;0.125), while in
the case of LiF we used (0.0;0.4), (0.0;0.45) and (0.0;0.5).

In Table (2) we report the Si gaps calculated in PBE with NC and PAW, together
with GW, HSE0.2;0.25, HSE0.3;0.25, HSE0.3;0.3 and HSE0.0;0.125 gaps calculated
with PAW pseudopotentials. The HSE0.2;0.25 gives the closest agreement with
the GW gaps. Increasing the value of $\mu$ keeping the value of $\alpha$
constant, as in HSE0.3;0.25, has the effect to lower the values of the gaps.
This is due to a smaller percentage of nonlocal HF exchange included in the
calculation. Instead, increasing the value of $\alpha$ keeping constant the
value of $\mu$, as in HSE0.3;0.3, has the effect to increase the gap values.
In this case a larger amount of nonlocal HF exchange is considered. The
HSE0.0;0.125 includes a full-range nonlocal HF exchange and the parameter
$\alpha$ acts as a screening. A value of $\alpha=0.125$ gives a good agreement
with GW gaps.

In Table (3) we report LiF gaps calculated in PBE with NC and PAW, together
with GW and HSE0.0;0.4 gaps calculated with PAW pseudopotentials. LiF is a
large gap insulator and it requires the correct long-range behaviour of the
nonlocal HF exchange. For this reason the HSE performs well only for $\mu=0$.
The value of $\alpha=0.4$ was found by imposing the constraints to recover the
GW gaps.

Table 2: Si gaps (eV). The use of norm-conserving pseudopotentials is indicated with the label NC, otherwise PAW pseudopotentials have been used.  Si | PBE${}^{\text{NC}}$ | PBE | GW | HSE0.2;0.25 | HSE0.3;0.25 | HSE0.3;0.3 | HSE0.0;0.125 | Exp
---|---|---|---|---|---|---|---|---
$\Gamma_{c}$ \- $\Gamma_{v}$ | 2.58 | 2.57 | 3.34 | 3.33 | 3.15 | 3.27 | 3.28 | 3.35111Reference PhysRevB.5.497 .
$X_{c}$ \- $\Gamma_{v}$ | 0.58 | 0.66 | 1.28 | 1.29 | 1.13 | 1.22 | 1.33 | 1.17222Reference kittel .
$L_{c}$ \- $\Gamma_{v}$ | 1.61 | 1.57 | 2.18 | 2.24 | 2.06 | 2.17 | 2.20 | 2.40333Reference PhysRevLett.54.142 ., 2.06444Reference HULTHEN19761341 .
Table 3: LiF gaps (eV). The use of norm-conserving pseudopotentials is indicated with the label NC, otherwise PAW pseudopotentials have been used.  LiF | PBE${}^{\text{NC}}$ | PBE | GW | HSE0.0;0.4 | Exp
---|---|---|---|---|---
$X_{c}$ \- $\Gamma_{v}$ | 9.21 | 9.12 | 14.21 | 13.99 | 14.20111Reference PhysRevB.13.5530 .
$X_{c}$ \- $\Gamma_{v}$ | 11.19 | 11.25 | 16.36 | 16.12 |
$L_{c}$ \- $\Gamma_{v}$ | 13.46 | 13.36 | 19.07 | 18.53 |

On top of the electronic structure we calculated the optical spectra starting
from the lowest level of theory, i.e. the independent-particle approximation
(IPA).

In Fig. (1) we show IPA-PBE for Si and LiF. We observe that the agreement is
excellent between NC and PAW pseudopotentials. PhysRevB.63.125108  This
implies that the differences we can observe in the spectra calculated with an
higher level of theory are only due to the relevance of the TDDFT and TDGKSDFT
kernels for the description of the excitons.

Figure 1: $\varepsilon_{2}$ for Si (top panel) and LiF (bottom panel)
calculated in IPA using NC and PAW pseudopotentials.

In Fig. (2) and in Fig. (3) we compare IPA-GW and IPA-HSE which are IPA
spectra calculated respecitvely on top of GW electronic structure and
HSEsr,μ;α electronic structure where $\mu$ and $\alpha$ values are those that
reproduce the GW gaps (see Table (2) and Table (3)).

For Si the calculations are consistent as shown in Fig. (2). The trend is the
same we observed for the electronic gaps of Table (2). In fact, IPA-
HSEsr,0.3;0.25 is slightly lower than IPA-HSEsr,0.2;0.25 due to a larger value
of $\mu$. Instead, using the same value of $\mu=0.3$ but an higher value of
$\alpha=0.3$ as in IPA-HSEsr,0.3;0.3, the spectrum shifts at higher energy due
to a larger percentage of nonlocal HF exchange. Also in the case of LiF we
found a good agreement as shown in Fig. (3), where we compared the IPA-GW
spectrum with IPA-HSEsr,0.0;0.4.

Figure 2: $\varepsilon_{2}$ for Si calculated in IPA using GW, IPA-
HSEsr,0.2;0.25, HSEsr,0.3;0.25, HSEsr,0.3;0.3 and HSEsr,0.0;0.125 and PAW
pseudopotentials.

Figure 3: $\varepsilon_{2}$ for LiF calculated in IPA using GW, IPA-
HSEsr,0.0;0.4 with PAW pseudopotentials.

The spectra of Fig. (2) and Fig. (3) do not include excitonic effects as IPA
is the lowest level of approximation for the calculation of optical spectra.

As already pointed out in rein+02prl , we show in Fig. (4) the excellent
agreement of GW-BSE with the experimental spectrum of Si. GW-TDPBE, as
expected, is not able to reproduce excitonic effects and it only slightly
improves the spectrum with respect to GW-RPA (see Refs. bott+04prb ;
PhysRevLett.91.056402 ; PhysRevB.87.205143 ). In fact, in TDPBE the exchange-
correlation kernel is PBE which has not the proper spatial nonlocality.
onid+02rmp  TDHF${}^{\text{sr},\mu;\alpha}$ optical spectra have a reasonable
shape but the intensity of the first peak around 3.5 eV is too low. In order
to increase the intensity of this peak, we need to increase the percentage of
the nonlocal HF exchange, as can be seen by comparing TDHFsr,0.3;0.25 with
TDHFsr,0.2;0.25. Otherwise, another strategy would be to increase the value of
the mixing parameter $\alpha$ as observed by comparing TDHFsr,0.3;0.25 with
TDHFsr,0.3;0.3. However, we believe that the use of short-range HF exchange
has not the necessary flexibility to improve further the spectrum. In fact,
increasing $\alpha$ or $\mu$ would change also the energy position of the
peaks.

TDHF${}^{\text{sr},0.0;0.125}$ contains the full range nonlocal HF exchange.
The $\alpha=0.125$ we have chosen, permits to be consistent with the previous
step, i.e. a correct electronic structure. However, this value of $\alpha$ is
still too small to correctly reproduce the experimental spectrum. Similar
results were also obtained by Yang et al. PhysRevB.92.035202  using for
$\alpha$ the value of the inverse RPA dielectric constant ($\sim$0.08). A
better description of the first peak could be done by increasing the $\alpha$
value, but also in this case this would cause a change in the position of the
energy peaks.

In the case of LiF, the GW-BSE reproduces an excitonic peak of 12.2 eV, which
is slightly lower PhysRevResearch.2.013091  than the experimental peak of
12.75 eV, as shown in Fig. (5). Instead, as expected, GW-TDPBE can not
reproduce the excitonic peak. TDHFsr,0.0;0.4 gives an excellent agreement with
the energy position of the experimental exciton. By increasing the value of
$\alpha$ we include more nonlocal HF exchange and therefore the exciton is
more strongly bound as we have shown for TDHFsr,0.0;0.45 and TDHFsr,0.0;0.5,
see Fig. (5).

The comparison between TDHFsr,μ;α and TDHFsr,μ;αXC${}^{\text{PBE}}$ is in Fig.
(6) for Si and in Fig. (7) for LiF. From Eq. (10) and Eq. (13) the difference
between these kernels is the addition to the nonlocal HF exchange of a
fraction of the semilocal exchange PBE $(1-\alpha)f^{\text{x,PBE}}$ and the
PBE correlation $f^{\text{c,PBE}}$. In the case of Si, adding a fraction of
semilocal exchange increases the intensity of the first peak around 3.5 eV,
therefore, improving the agreement with the experiment. The energy position of
the peak is not changed. Instead, in the case of LiF the energy position of
the peak is slightly shifted to lower energy and the intensity of the peak
changes. However, concerning the peak intensity we did not find a clear trend.

In Fig. (8) and Fig. (9) we finally present the TDHFsr,μ;αXC${}^{\text{PBE}}$
spectra which give the best agreement with experiment and we compare them to
the GW-BSE spectra.

Figure 4: Comparison of experimental $\varepsilon_{2}$ for Si with GW-BSE, GW-
TDPBE, TDHFsr,0.2;0.25, TDHFsr,0.3;0.25, TDHFsr,0.3;0.3 and TDHFsr,0.0;0.125
and PAW pseudopotential. Experiment is from Ref. PhysRevB.36.4821 . Figure 5:
Comparison of experimental $\varepsilon_{2}$ for LiF with GW-BSE, GW-TDPBE,
TDHFsr,0.0;0.4, TDHFsr,0.0;0.45 and TDHFsr,0.0;0.5 and PAW pseudopotentials.
Experiment is from Ref. Roessler:67 . Figure 6: $\varepsilon_{2}$ for Si :
effect of inclusion of a fraction of PBE exchange-correlation. PAW
pseudopotential has been used. Experiment is from Ref. PhysRevB.36.4821 .
Figure 7: $\varepsilon_{2}$ for LiF : effect of inclusion of a fraction of PBE
exchange-correlation. PAW pseudopotentials have been used. Experiment is from
Ref. Roessler:67 . Figure 8: Comparison of experimental $\varepsilon_{2}$ for
Si with GW-BSE and TDHFsr,0.3;0.25XC${}^{\text{PBE}}$. PAW pseudopotential has
been used. Experiment is from Ref. PhysRevB.36.4821 . Figure 9: Comparison of
experimental $\varepsilon_{2}$ for LiF with GW-BSE and
TDHFsr,0.0;0.4XC${}^{\text{PBE}}$. PAW pseudopotentials have been used.
Experiment is from Ref. Roessler:67 . Figure 10: Comparison of experimental
$\varepsilon_{2}$ for Si with GW-BSE, TDHFsr,0.3;0.3XC${}^{\text{PBE}}$ (PAW
pseudopotential) and with the long-range corrected kernels : LRC, RPA-BO and
JGM (NC pseudopotential). TDCAM is from Ref. PhysRevB.92.081204 . Experiment
is from Ref. PhysRevB.36.4821 . Figure 11: Comparison of experimental
$\varepsilon_{2}$ for LiF with GW-BSE, TDHFsr,0.3;0.3XC${}^{\text{PBE}}$ (PAW
pseudopotentials) and with the long-range corrected kernels (NC
pseudopotentials): LRC, RPA-BO and JGM. TDCAM is from Ref. PhysRevB.92.081204
. Experiment is from Ref.Roessler:67 .

Finally, for Si, we compare in Fig. (10) the selected
TDHFsr,0.3;0.25XC${}^{\text{PBE}}$ spectrum to TDDFT with long-range corrected
kernels with a scissor shift of 0.7 eV (see Table (2)).

We show the results for LRC ($\alpha^{\text{LRC}}=0.20$), RPA-BO
($\alpha^{\text{RPA-BO}}=0.13$) and JGM ($\alpha^{\text{JGM}}=0.12$) kernels.
gaurJCP2019  For the long-range corrected kernels a higher value of the
$\alpha$ parameter can be interpreted as if a higher nonlocal HF long-range
contribution was included. This contribution is higher in LRC than in RPA-BO
and JGM kernels. This is the reason for which LRC better describes the
spectrum around 3.5 eV.

However, despite the better behaviour of LRC kernel, we want to point out that
RPA-BO and JGM kernels do not require any adjustable parameter, which is an
enormous advantage as they can be applied to any kind of materials.

The TDHFsr,0.3;0.25XC${}^{\text{PBE}}$ gives a reasonable Si spectrum, similar
to the spectra from the RPA-BO and the JGM kernels.

We add to the comparison the result from the range-separated CAM proposed by
Rafaely-Abramson et al. PhysRevB.92.081204  which contains a fraction of long-
range nonlocal HF exchange. The approach of Rafaely-Abramson et al.
PhysRevB.92.081204  is in excellent agreement with experiment and seems also
to improve with respect to GW-BSE. However, this approach contains 3
parameters $\alpha$, $\beta$ and $\mu$. The parameter $\alpha$ is the amount
of short-range exact exchange, $\beta$ is calculated as
$\alpha+\beta=1/\varepsilon_{0}$ and $\mu$ is the range-separation parameter.
To obtain this excellent agreement, $\mu$ was optimally tuned in order to
reproduce the electronic gap, and $\alpha$ and $\beta$ are obtained from the
material’s dielectric constant $\varepsilon_{0}$. They used $\mu=0.11$ Bohr-1,
$\alpha=0.2$ and $\varepsilon_{0}=12$.

Comparing the method with a fraction of short-range nonlocal HF exchange
($\mu=0.3$ Bohr-1 and $\alpha=0.25$) to the one with long-range nonlocal HF
exchange ($\mu=0.11$ Bohr-1, $\alpha=0.2$), and considering that the mixing
parameter $\alpha$ is of the same order of magnitude, we observe that the
range-separation parameter $\mu$ is larger when a fraction of short-range is
used. This is reasonable as the role of $\mu$ is opposite between short and
long-range.

However, using TDHFsr,0.3;0.25XC${}^{\text{PBE}}$ (short-range), it is not
possible to obtain the same agreement with experiment that is reproduced when
a fraction of long-range nonlocal HF exchange is included. To obtain the same
performance of the long-range scheme, the value of $\mu$ should be increased.
However, this would cause a shift of the excitonic peaks to lower energy and
the spectrum will be wrong.

In Fig. (11), for LiF, we compare TDHFsr,0.0;0.4XC${}^{\text{PBE}}$ with TDDFT
with long-range corrected kernels with a scissor shift of 5.0 eV (see Table
(3)).

We used $\mu=0.0$ Bohr-1 and $\alpha=0.4$ as we need the full range nonlocal
HF exchange to reproduce the experimental spectrum. In fact, we did not find
any finite values of $\mu$ different from zero for which using only a fraction
of short-range nonlocal HF exchange it would be possible to reproduce the
experimental spectrum.

We show LRC ($\alpha^{\text{LRC}}=8.0$), RPA-BO ($\alpha^{\text{RPA-BO}}=8.8$)
and JGM ($\alpha^{\text{JGM}}=7.93$) kernels. gaurJCP2019  The
TDHFsr,0.0;0.4XC${}^{\text{PBE}}$ and RPA-BO are in excellent agreement with
the energy position of the excitonic peak. However, RPA-BO, as well as LRC and
JGM overestimate the peak intensity, which in Fig. (11) has been multiplied by
0.1 in order to compare the theoretical approaches. Furthermore, we observe
that the energy of the JGM peak is around 1 eV higher than the result
presented in the original work of Trevisanutto et al. PhysRevB.87.205143 .
This is due to the different scissor value taken to correct the energies.

We add to this comparison also the result from the range-separated CAM
proposed by Rafaely-Abramson et al. PhysRevB.92.081204  which is in excellent
agreement with the experiment and also improve with respect to the GW-BSE.
They used $\mu=0.58$ Bohr-1, $\alpha=0.2$ and $\epsilon_{0}=1.9$.

##  V Conclusion

We compared the performance of TDGKSDT range-separated hybrid functionals and
TDDFT long-range corrected kernels for the description of excitons in solids.
The comparison was illustrated for the case of Si and LiF, representative of
continuum and strong excitons.

We studied hybrid functionals with a fraction of short-range nonlocal HF
exchange. For Si, by optimally tuning $\mu$ and $\alpha$, it is possible to
reproduce the satisfactory experimental spectrum. Instead, for LiF we did not
find any finite values of $\mu$, different from zero, for which is possible to
reproduce the experimental spectrum. In the case of LiF we need to use the
(full range) nonlocal HF exchange ($\mu=0.0$) in order to satisfactory
reproduce the experiment. Therefore, exchange is much more important for
strong excitons than for weak ones.

We also studied the long-range corrected kernels: LRC rein+02prl , RPA-BO
PhysRevLett.114.146402  and JGM PhysRevB.87.205143 . These kernels perform
comparably to hybrid functionals with short-range nonlocal HF exchange. Except
that for LiF the intensity of the excitonic peak is strongly overestimated.

We included in our discussion the hybrid scheme of Refs. PhysRevB.92.081204 ;
PhysRevMaterials.3.064603  which has a long-range nonlocal HF exchange
component. This approach has an excellent agreement with experiment for both
Si and LiF and it also seems to improve with respect to GW-BSE. This approach
is the most flexible.

From this comparison it appears that the hybrid scheme with long-range
nonlocal HF exchange performs better than the hybrid scheme with short-range
nonlocal HF exchange. We believe that for Si, and therefore for weak excitons,
it is the lack of long-range component of the nonlocal HF exchange which
causes a not yet excellent description of the exciton around 3.4 eV. The
situation is even worse for LiF, and therefore for strong excitons, where the
short-range separation demonstrated not to work.

The main difficulty of using range-separated hybrid functionals is their
dependence on parameters that have to be chosen and which strongly depends on
the material. A general strategy to find these parameters is needed.

Moreover, despite the promising behaviour of range-separated schemes, long-
range kernels continue to be attractive. In fact, the computational cost is
lower and a kernel such as RPA-BO does not require any adjustable parameters.

###### Acknowledgements.

This work was performed using HPC resources from GENCI-IDRIS Grant
2021-x2021082131 and Grant 2021-x202109544.

## References

  * (1) J. Koskelo et al.,  Phys. Rev. B 95, 035125 (2017).
  * (2) P. Cudazzo et al.,  Phys. Rev. Lett. 116, 066803 (2016).
  * (3) A. R. Kshirsagar, X. Blase, C. Attaccalite, and R. Poloni,  The Journal of Physical Chemistry Letters 12, 4045 (2021),  PMID: 33881873.
  * (4) Y. Dong et al.,  Nature Communications 11, 4617 (2020).
  * (5) E. Barré et al.,  Science 376, 406 (2022).
  * (6) G. Onida, L. Reining, and A. Rubio,  Rev. Mod. Phys. 74, 601 (2002).
  * (7) E. Rebolini, J. Toulouse, A. Savin, in ”Concepts and Methods in Modern Theoretical Chemistry, Vol. 1: Electronic Structure and Reactivity”, CRC Press, 367-390, (2013).
  * (8) F. Sottile, V. Olevano, and L. Reining,  Phys. Rev. Lett. 91, 056402 (2003).
  * (9) G. Adragna, R. Del Sole, and A. Marini,  Phys. Rev. B 68, 165108 (2003).
  * (10) A. Marini, R. Del Sole, and A. Rubio,  Phys. Rev. Lett. 91, 256402 (2003).
  * (11) F. Sottile et al.,  International Journal of Quantum Chemistry 102, 684 (2005).
  * (12) N. Gauriot, V. Véniard, and E. Luppi,  The Journal of Chemical Physics 151, 234111 (2019).
  * (13) L. Reining, V. Olevano, A. Rubio, and G. Onida,  Phys. Rev. Lett. 88, 066404 (2002).
  * (14) S. Botti et al.,  Phys. Rev. B 69, 155112 (2004).
  * (15) S. Rigamonti et al.,  Phys. Rev. Lett. 114, 146402 (2015).
  * (16) S. Sharma, J. K. Dewhurst, A. Sanna, and E. K. U. Gross,  Phys. Rev. Lett. 107, 186401 (2011).
  * (17) P. E. Trevisanutto, A. Terentjevs, L. A. Constantin, V. Olevano, and F. D. Sala,  Phys. Rev. B 87, 205143 (2013).
  * (18) Y.-M. Byun and C. A. Ullrich,  Phys. Rev. B 95, 205136 (2017).
  * (19) J. Sun, C.-W. Lee, A. Kononov, A. Schleife, and C. A. Ullrich,  Phys. Rev. Lett. 127, 077401 (2021).
  * (20) J. Paier, M. Marsman, and G. Kresse,  Phys. Rev. B 78, 121201 (2008).
  * (21) S. Refaely-Abramson, M. Jain, S. Sharifzadeh, J. B. Neaton, and L. Kronik,  Phys. Rev. B 92, 081204 (2015).
  * (22) Z.-h. Yang, F. Sottile, and C. A. Ullrich,  Phys. Rev. B 92, 035202 (2015).
  * (23) J. Sun, J. Yang, and C. A. Ullrich,  Phys. Rev. Research 2, 013091 (2020).
  * (24) B. G. Janesko, T. M. Henderson, and G. E. Scuseria,  Phys. Chem. Chem. Phys. 11, 443 (2009).
  * (25) F. Zapata, E. Luppi, and J. Toulouse,  The Journal of Chemical Physics 150, 234104 (2019).
  * (26) E. Rebolini, A. Savin, and J. Toulouse,  Molecular Physics 111, 1219 (2013).
  * (27) J. Heyd and G. E. Scuseria,  J. Chem. Phys. 121, 1187 (2004).
  * (28) A. J. Garza and G. E. Scuseria,  The Journal of Physical Chemistry Letters 7, 4165 (2016).
  * (29) J. P. Perdew, K. Burke, and M. Ernzerhof,  Phys. Rev. Lett. 77, 3865 (1996).
  * (30) D. Wing et al.,  Phys. Rev. Materials 3, 064603 (2019).
  * (31) S. Rigamonti et al.,  Phys. Rev. Lett. 117, 159702 (2016).
  * (32) Dp code, http://www.dp-code.org/.
  * (33) E. Luppi and H. Hübener and V. Véniard, J. Chem. Phys., accepted.
  * (34) E. Luppi, H. Hübener, and V. Véniard,  Phys. Rev. B 82, 235201 (2010).
  * (35) X. Gonze and B. Amadon and P.-M. Anglade and J.-M. Beuken and F. Bottin and P. Boulanger and F. Bruneval and D. Caliste and R. Caracas and M. Cote and T. Deutsch and L. Genovese and Ph. Ghosez and M. Giantomassi and S. Goedecker and D.R. Hamann and P. Hermet and F. Jollet and G. Jomard and S. Leroux and M. Mancini and S. Mazevet and M. J. T. Oliveira and G. Onida and Y. Pouillon and T. Rangel and G.-M. Rignanese and D. Sangalli and R. Shaltaf and M. Torrent and M. J. Verstraete and G. Zerah and J.W. Zwanziger, Computer Phys. Commun. 180, 2582-2615 (2009).
  * (36) X. Gonze and G.-M. Rignanese and M. Verstraete and J.-M. Beuken and Y. Pouillon and R. Caracas and F. Jollet and M. Torrent and G. Zerah and M. Mikami and Ph. Ghosez and M. Veithen and J.-Y. Raty and V. Olevano and F. Bruneval and L. Reining and R. Godby and G. Onida and D. R. Hamann and D.C. Allan, Zeit. Kristallogr. 220, 558-562 (2005).
  * (37) G. Kresse and J. Hafner,  Phys. Rev. B 47, 558 (1993).
  * (38) G. Kresse and J. Furthmüller,  Phys. Rev. B 54, 11169 (1996).
  * (39) In the case of GW-BSE/TDGKSDFT, the construction of the Hamiltonian scales as $N_{k}N^{2}_{v}N^{2}_{c}N_{G}$ where $N_{k}$ is the number of k-points in the Brillouin zone, $N_{v}$ is the number of valence bands, $N_{c}$ is the number of conduction bands and $N_{G}$ is the number of G-vectors. However, except for very large systems, the main limiting factor usually comes from the diagonalization of the Hamiltonian which scales cubically with the matrix rank $N_{rank}^{3}$ and $N_{rank}=N_{k}N_{v}N_{c}$. TDDFT requires the diagonalization for each k-point of a matrix of rank $N_{rank}=(N_{v}+N_{c})$ and the evaluation of the response function scales as $N_{k}N_{v}N_{c}N_{G}^{2}$.
  * (40) M. Welkowsky and R. Braunstein,  Phys. Rev. B 5, 497 (1972).
  * (41) C. Kittel, Introduction to Solid State Physics, 6th ed. (Wiley, New York, 1986) p. 185.
  * (42) D. Straub, L. Ley, and F. J. Himpsel,  Phys. Rev. Lett. 54, 142 (1985).
  * (43) R. Hulthén and N. Nilsson,  Solid State Communications 18, 1341 (1976).
  * (44) M. Piacentini, D. W. Lynch, and C. G. Olson,  Phys. Rev. B 13, 5530 (1976).
  * (45) B. Adolph, J. Furthmüller, and F. Bechstedt,  Phys. Rev. B 63, 125108 (2001).
  * (46) P. Lautenschlager, M. Garriga, L. Vina, and M. Cardona,  Phys. Rev. B 36, 4821 (1987).
  * (47) D. M. Roessler and W. C. Walker,  J. Opt. Soc. Am. 57, 835 (1967).