# QS$G\hat{W}$: Quasiparticle Self consistent $GW$ with ladder diagrams in $W$

Brian Cunningham  School of Mathematics and Physics, Queen’s University
Belfast, Belfast BT7 1NN, Northern Ireland, United Kingdom  Myrta Grüning
School of Mathematics and Physics, Queen’s University Belfast, Belfast BT7
1NN, Northern Ireland, United Kingdom  European Theoretical Spectroscopy
Facilities (ETSF)  Dimitar Pashov  King’s College London, London WC2R 2LS,
United Kingdom  Mark van Schilfgaarde  King’s College London, London WC2R 2LS,
United Kingdom

###### Abstract

We present an approach to calculate the electronic structure for a range of
materials using the quasiparticle self-consistent $GW$ method with vertex
corrections included in the screened Coulomb interaction $W$. This is achieved
by solving the Bethe-Salpeter equation for the polarization matrix at all
$\bm{k}$-points in the Brillouin zone. We refer to this method as
QS$G\hat{W}$. We show that including ladder diagrams in $W$ can greatly reduce
the band gap overestimation of RPA-based QS$GW$. The resultant discrepency of
the calculated band gap in this method is then attributed mostly to the fact
that electron-phonon contributions to $W$ are neglected; which would allow one
to then obtain an estimate for the size of this effect. We present results for
a range of systems from simple $sp$ semiconductors to the strongly correlated
systems NiO and CoO. Results for systems where the RPA-based QS$GW$ band gap
is larger than expected are investigated, and an estimate for the Frölich
contribution to the gap is included in a few polar compounds where QS$GW$ can
overestimate the gap by as much as 2 eV. The improvement over QS$GW$ for the
dielectric constants is also presented.

###### pacs:

42.25.Bs,11.10.St,71.15.-m,78.20.-e

##  I Introduction

The electronic band structure of crystalline solids provides essential
information on material properties. For example, besides determining whether a
material is insulating or metallic, it is the starting point to understand and
compute the response to electric and magnetic fields and transport properties.
As a consequence, one of the challenges within electronic structure is to
develop methods that accurately and efficiently predict the band structure for
a broad range of systems. The main objective of this work is to present an
accurate first-principles approach based on many-body perturbation theoryHedin
(1965a) (MBPT) to compute the band structure of solids.

The on-hand first-principles approach to compute band structures is density-
functional theoryHohenberg and Kohn (1964) (DFT). In spite of being a ground-
state theory, in its standard Kohn-Sham (KS) formulationKohn and Sham (1965)
DFT outputs the band structure of the auxiliary KS system. In some cases the
latter provides a reasonable approximation to the band structure of the real
system. On the other hand, the calculated KS band structure is often
unsatisfactory: the best-known shortcoming is the large underestimation of the
fundamental band gap, but there are also instances in which DFT fails to
predict e.g. the insulating behavior of a material. At present, the state-of-
the-art approach to compute band structures is the so-called $GW$
approximationHybertsen and Louie (1986) (GWA) in which the quasiparticle (QP)
bandstructure is obtained perturbatively starting from a KS calculation. The
GWA is based on MBPT which in principle contains the physics to describe the
QP bandstructure. However, in its practical implementation the performance of
the GWA depends critically on the KS starting point.van Schilfgaarde  _et al._
(2006a) Then – even if in many cases the GWA satisfactorily corrects the KS
bandstructure, providing a good agreement with results from direct/inverse
photoemission experiments – when the behavior of the physical system deviates
significantly from that of an independent-particle system, the GWA breaks
down.111See, for example, the discussion and results for NiO in Ref. Kotani
_et al._ , 2007.

Advanced methods based on the GWA implement some form of self-consistency so
to ‘forget’ the KS starting point. For example, the quasiparticle self-
consistent $GW$ (QS$GW$) methodFaleev  _et al._ (2004); van Schilfgaarde  _et
al._ (2006b); Kotani  _et al._ (2007)—which is the method of choice in this
work—determines, self-consistently, the optimum starting point Hamiltonian.
While self-consistency schemes remedy to the most apparent failures of the
standard GWA, in general they do not necessarily improve over the GWA in spite
of the additional level of sophistication and computational cost.

In fact, introducing the self-consistency uncovers that a large part of the
GWA success relies on the cancellations of errors: the overestimation of the
dielectric screening—that stems from taking the KS electronic structure as the
starting point and contributes to underestimating the band gap—is compensated
by the neglect of the excitonic and electron-phonon contributions to the
dielectric response, which effectively overestimates the band gap. In this
work, the excitonic contributions are taken into account by including ladder
diagrams into the screened Coulomb interaction $W$ through use of the Bethe-
Salpeter equation (BSE)Sottile  _et al._ (2003); Marini  _et al._ (2003) for
the polarization. We also included a correction accounting for the Frölich
contribution from lattice polarization effects using the method outlined in
Ref. Lambrecht  _et al._ , 2017.

Other works have considered the effect of vertex corrections to the dielectric
screening on the band gap. For example, Refs. Franchini  _et al._ , 2010;
Shishkin  _et al._ , 2007; Bruneval  _et al._ , 2005 included ladder diagrams
through an effective nonlocal static kernel constructed within time-dependent
functional theory to mimic the BSE. More recently, Kutepov proposed several
schemes for the self-consistent solution of Hedin’s equations including vertex
corrections. Kutepov (2016) In particular, for selected semiconductors and
insulators, Kutepov (2017) he included the vertex correction for the
dielectric screening at the BSE level together with a so-called first-order
approximation for the vertex in the self-energy. In all the cases, an
improvement over the QS$GW$ results has been observed. With respect to these
previous works, we include vertex corrections to the dielectric screening only
at the BSE level, but introducing the usual static approximation for the BSE
kernel, which was lifted in Refs Kutepov, 2017, 2016. We also include the
lattice polarization effects, which can give corrections that are as large as
those from the vertex. These corrections were either neglected or estimated in
the previous works. Kutepov included a vertex for the polarization, and also a
first-order vertex for the self-energy $\Sigma$. We only include the former;
nevertheless we show it yields quite satisfactory results. Another difference
between our work and Kutepov is that our construction of the quasiparticlized
self-energy ${\Sigma_{ij}^{0}}$ coupling states $i$ and $j$ is evaluated from

$\displaystyle{\Sigma_{ij}^{0}}=\frac{1}{2}\sum_{ij}|\psi_{i}\rangle\left\\{{{\rm
Re}[\Sigma({\varepsilon_{i}})]_{ij}+{\rm
Re}[\Sigma({\varepsilon_{j}})]_{ij}}\right\\}\langle\psi_{j}|.$ (1)

Kutepov instead uses $\Sigma(\omega{=}0)$ and folds in an effective energy
dependence through $\Sigma^{\prime}(\omega{=}0)$. There are formal reasons why
Eq. 1 should be better: it satisfies a minimization condition, where a norm
defining the difference between the interacting and quasiparticled Green’s
function is minimized. It has also been shown to satisfy a variational
principle Ismail-Beigi (2017). On the other hand, a detailed comparison
between the two schemes has not been made.

Starting from the Hedin’s equations (Sec. II.1), we outline how the original
QS$GW$ approximation (Sec. II.2) is modified to include excitonic
contributions (Sec. II.3). The resulting approach, which is referred to as
QS$G\hat{W}$ (where the substitution $W{\rightarrow}\hat{W}$ implies that
vertex corrections are included in $W$), was numerically implemented within an
all-electron framework using a linear muffin-tin orbital basis set (Sec. II.4)
in the Questaal packageque  and applied to calculate the fundamental band gap
and static dielectric constants of a broad range of materials: Si, C, Ge, LiF,
LiCl, NaCl, MgO, GaAs, GaSb, InSb, InP, InAs, CuCl, CdTe, CdO, ZnO, CuAlO2,
hexagonal-BN, wurzite InN and GaN; the perovskite SrTiO3; the Bismuth-
Telleride halides BiTeBr and BiTeI; the transition metal dichalcogenide bulk
MoS2; and the strongly correlated transition metal monoxides NiO and CoO. For
the polar materials we then included the contribution from lattice
polarization effects. We obtain accurate band gaps and static dielectric
(optical) constants for the considered systems, and we quantify both the
contributions from excitonic effects and lattice polarization effects to the
band gap.

##  II Theory and numerical implementation

###  II.1 The GWA from the Hedin’s equations

The approach described in this work, the standard GWA and the QSGW, are all
derived from the many-body perturbative approach developed by HedinHedin
(1965a). In this method, the following set of closed coupled equationsHedin
(1965b); Schwinger (1951); Martin and Schwinger (1959) are to be solved
iteratively:

$\displaystyle\Sigma(1,2)$ $\displaystyle={\rm i}\int{\rm
d}(34)~{}G(1,3^{+})W(1,4)\Lambda(3,2,4)$ (2) $\displaystyle G(1,2)$
$\displaystyle=G_{0}(1,2)+\int{\rm d}(34)~{}G_{0}(1,3)\Sigma(3,4)G(4,2)$ (3)
$\displaystyle W(1,2)$ $\displaystyle=v(1,2)+\int{\rm
d}(34)v(1,3)P(3,4)W(4,2)$ (4) $\displaystyle P(12)$ $\displaystyle=-{\rm
i}\int{\rm d}(34)G(1,3)G(4,1^{+})\Lambda(3,4,2)$ (5)
$\displaystyle\Lambda(1,2,3)$ $\displaystyle=\delta(1,2)\delta(1,3)+$
$\displaystyle~{}~{}~{}\int{\rm d}(4567)\frac{\delta\Sigma(1,2)}{\delta
G(4,5)}G(4,6)G(7,5)\Lambda(6,7,3)$ (6)

where $G$ is the Green’s function,
$v(\bm{r},\bm{r}^{\prime})=1/|\bm{r}-\bm{r}^{\prime}|$ is the bare Coulomb
interaction, $W$ is the screened Coulomb interaction, $\Lambda$ is the
irreducible vertex function, $P$ is the irreducible polarizability (the
functional derivative of the induced density with respect to the total
potential) and $\Sigma$ is the self-energy operator. In Eqs. (2) and (5), the
$+$ superscript implies $t^{\prime}=t+\eta$, with $\eta\rightarrow 0^{+}$.

In the standard GWA, also known as one-shot $GW$ or $G_{0}W_{0}$, Eq. (6) is
approximated as:

$\Lambda(1,2,3)\approx\delta(1,2)\delta(1,3)$ (7)

in both the expressions for the self-energy [Eq. (2)] and the irreducible
polarization [Eq. (5)]. In addition, $\Sigma$ and $P$ are both evaluated for
$G=G_{0}$, the independent-particles Green’s function,222This is equivalent to
perform one iteration, starting with $\Sigma=0$ in Eq. (2) that in the
frequency-domain takes the form

$G_{0}(\bm{r},\bm{r}^{\prime},\omega)=\sum_{n}\frac{\psi_{n}(\bm{r})\psi_{n}^{*}(\bm{r}^{\prime})}{\omega-\varepsilon_{n}\pm{\rm
i}\eta}.$ (8)

In Eq. (8), $\psi_{n}$ and $\varepsilon_{n}$ are the single-particle
wavefunctions and energies (the index $n$ contains band, spin and wavevector
indices and the $+(-)$ is for unoccupied(occupied) bands.

The approximation for the irreducible polarization obtained by neglecting the
vertex is referred to as the Random-Phase-Approximation (RPA).Onida  _et al._
(2002) When using Eq. (8) in frequency space it takes the form:

$P^{\rm
RPA}(\bm{r},\bm{r}^{\prime};\omega)=\sum_{n_{1}n_{2}}(f_{n_{2}}-f_{n_{1}})\frac{\psi_{n_{2}}^{*}(\bm{r})\psi_{n_{1}}(\bm{r})\psi_{n_{1}}^{*}(\bm{r}^{\prime})\psi_{n_{2}}(\bm{r}^{\prime})}{\varepsilon_{n_{2}}-\varepsilon_{n_{1}}-\omega+{\rm
i}(f_{n_{2}}-f_{n_{1}})\eta},$ (9)

where $f_{n}$ are the single-particle occupations.

The Green’s function in Eq. (8) can be constructed from the Kohn-Sham
electronic structure, which is obtained from the self-consistent solution of
Schrödinger-like equations with the single-particle Hamiltonian

$H_{0}(\bm{r})=-\frac{1}{2}\bm{\nabla}^{2}+V_{\rm
ext}(\bm{r})[\rho]+V_{H}(\bm{r})[\rho]+V_{\rm XC}(\bm{r})[\rho].$ (10)

In Eq. (10), $V_{\rm ext}(\bm{r})$ is the external potential due to nuclei and
any external fields, $V_{H}(\bm{r})$ is the Hartree potential describing the
classical mean-field electron-electron interaction, and $V_{\rm XC}(\bm{r})$
is the exchange-correlation potential describing correlation effects missing
in $V_{H}(\bm{r})$. $\rho$ is the electronic density calculated as $\sum_{\rm
occ}|\psi_{n}(\bm{r})|^{2}$, that is constructed from the eigensolutions of
$H_{0}$ (from which the self-consistent solution is obtained). Because of
$V_{\rm XC}$, the corresponding single-particle $G_{0}$ effectively contains
many-body effects. Then in Eq. 3 for the many-body Green’s function, the self-
energy is replaced by $\Delta\Sigma=\Sigma-V_{\rm XC}$ to avoid double
counting. Eq. 3 is rewritten as a nonlinear equation for the quasiparticle
energies $E_{n\bm{k}}$

$E_{n\bm{k}}=\varepsilon_{n\bm{k}}+\langle\psi_{n\bm{k}}|\Sigma(E_{n\bm{k}})-V_{\rm
XC}|\psi_{n\bm{k}}\rangle,$ (11)

and solved after it is linearized:

$E_{n\bm{k}}=\varepsilon_{n\bm{k}}+Z_{n\bm{k}}\langle\psi_{n\bm{k}}|\Sigma(\varepsilon_{n\bm{k}})-V_{\rm
XC}|\psi_{n\bm{k}}\rangle,$ (12)

where the renormalization factor $Z_{n\bm{k}}$ is:

$Z_{n\bm{k}}=(1-\displaystyle\partial\Sigma(\omega)/\partial\omega|_{\omega=\varepsilon_{n\bm{k}}})^{-1}.$
(13)

###  II.2 Self-consistency: QS$GW$

In the DFT based $G_{0}W_{0}$ approximation, the $E_{n\bm{k}}$ are obtained as
a first-order correction of the Kohn-Sham single-particle energies. As
mentioned in the introduction, the GWA works when the Kohn-Sham system gives a
qualitatively correct description of the physical system, i.e., when the Kohn-
Sham single-particle energies are not ‘too far’ from the quasiparticle
energies. When this is not the case, the GWA does not give accurate results. A
solution to that is to replace the corrected energy $E_{n}$ [Eq. (12)] either
in the Green’s function [Eq. (8)], or in the RPA polarization [Eq. (9)]
entering the screened potential $W$, or in both. Another solution is to choose
a better starting point — as for example using a hybrid DFT approach,Schleife
_et al._ (2009) LDA+U Anisimov  _et al._ (1997) or the Coulomb-hole Screened
exchange approximation (COHSEX)Hedin (1965a).

Here we use the QS$GW$ approach in which the starting point is chosen to
effectively minimize $\Delta\Sigma$, which is the difference between the self-
energy and the exchange-correlation part of the single-particle potential. In
practice, once the self-energy has been calculated within the $GW$
approximation, a new effective single-particle static potential is determined
by Eq. 1. Then, by substituting $V_{\rm XC}$ with $\Sigma^{0}$ in Eq. (10), a
new set of single-particle energies and wavefunctions can be determined. In
turn, those can be used to re-calculate the $GW$ self-energy, and the whole
procedure can be repeated until self-consistency in the energies and
eigenvalues is achieved. In this procedure the resulting electronic structure
does not depend on the quality of the Kohn-Sham DFT electronic structure for
the system.

###  II.3 Ladder diagram corrections to $W$

In previous works employing QS$GW$ van Schilfgaarde  _et al._ (2006b); Kotani
_et al._ (2007); Faleev  _et al._ (2004); Cunningham  _et al._ (2018) the RPA
for $W$ was used. This—along with other assumptions made—leads to, e.g., a
significant band gap overestimation. Here we go beyond the RPA and include
ladder diagram corrections in $W$ through the BSE for the
polarization.Salpeter and Bethe (1951) To include the vertex in Hedin’s
equations we need to determine the interaction kernel $\delta\Sigma/\delta G$
in Eq. 6. Suppose we have adopted the $GW$ approximation for $\Sigma$333The
vertex in Eq. (2) can been shown to effectively cancel with the $Z$-factor,
see for example Appendix A in Ref. Kotani  _et al._ , 2007 and assume that
$\delta W/\delta G$ is negligible,Onida  _et al._ (2002) then
$\delta\Sigma(12)/\delta G(45)={\rm i}W(12)\delta(1,4)\delta(2,5)$, where $W$
is determined in the $GW$ approximation. This expression is then inserted in
Eq. 6.

Before we present the BSE for the polarization we will introduce the expansion
of the two-point polarization to its four-point counterpart:
$P(12)=P(1122)=P(1324)\delta(1,3)\delta(2,4)$. We are now able to present an
expression for the polarization that goes beyond the RPA using Eqs 5 and 6 and
adopting the expression for the interaction kernel from above;

$P(12)=P^{\rm RPA}(12)-\int{\rm d}(34)P^{\rm RPA}(1134)W(34)P(3422),$ (14)

where $P^{\rm RPA}(1234)=-{\rm i}G(13)G(42)$. The $W$ that appears in the
interaction kernel, $\delta\Sigma/\delta G$, is calculated at the level of the
RPA and this is usually assumed to be static,444In few works this
approximation has been relaxed, see e.g. Ref. Marini and Del Sole, 2003. i.e.,
$\delta\Sigma/\delta G={\rm i}W^{RPA}(\omega{=}0)$. To avoid confusion with
$W$ in Eq. 2 we will refer to the $W$ in Eq. 14 from here on as $K$.

The Dyson-like equation for the polarizability, Eq. 14, can be transformed to
an eigenproblem for an effective 2-particle Hamiltonian by introducing the
basis of single particle eigenfunctions that diagonalize the RPA polarization.
Using the completeness of the eigenfunctions, any 4-point quantity can be
expanded as

$S(\bm{r}_{1},\bm{r}_{2},\bm{r}_{3},\bm{r}_{4})=\sum_{\tiny
n_{1}n_{2}n_{3}n_{4}\normalsize}S_{n_{1}n_{2}n_{3}n_{4}}\psi_{n_{1}}(\bm{r}_{1})\psi_{n_{2}}^{*}(\bm{r}_{2})\psi_{n_{3}}^{*}(\bm{r}_{3})\psi_{n_{4}}(\bm{r}_{4}),$
(15)

where we have again combined band, spin and wavevector indices, and
$S_{n_{1}n_{2}n_{3}n_{4}}=\int{\rm
d}(\bm{r}_{1}\bm{r}_{2}\bm{r}_{3}\bm{r}_{4})S(\bm{r}_{1},\bm{r}_{2},\bm{r}_{3},\bm{r}_{4})\times$
$\psi_{n_{1}}^{*}(\bm{r}_{1})\psi_{n_{2}}(\bm{r}_{2})\psi_{n_{3}}(\bm{r}_{3})\psi_{n_{4}}^{*}(\bm{r}_{4})$.

Inserting the expression for the RPA polarization from Eq. (9) in Eq. (14),
one arrives at the following expression for the polarization

$P_{\scriptsize\begin{array}[]{l}n_{1}n_{2}\bm{k}\\\
n_{3}n_{4}\bm{k}^{\prime}\end{array}}(\bm{q},\omega)=\left[H(\bm{q})-\omega\right]_{\scriptsize\begin{array}[]{l}n_{1}n_{2}\bm{k}\\\
n_{3}n_{4}\bm{k}^{\prime}\end{array}}^{-1}(f_{n_{4}\bm{k}^{\prime}+\bm{q}}-f_{n_{3}\bm{k}^{\prime}}),$
(16)

whereby the conservation of momentum we have
$\bm{k}_{2(4)}=\bm{k}_{1(3)}+\bm{q}$; and555Note that if we are determining
the reducible polarizabililty $\chi$, such that $W=v+v\chi v$ then the Kernel
becomes $K=W-V$ Onida  _et al._ (2002)

$\begin{array}[]{rl}H_{\scriptsize\begin{array}[]{l}n_{1}n_{2}\bm{k}\\\
n_{3}n_{4}\bm{k}^{\prime}\end{array}}(\bm{q})=&(\varepsilon_{n_{2}\bm{k}^{\prime}+\bm{q}}-\varepsilon_{n_{1}\bm{k}^{\prime}})\delta_{n_{1}n_{3}}\delta_{n_{2}n_{4}}\delta_{\bm{k}\bm{k}^{\prime}}+(f_{n_{2}\bm{k}+\bm{q}}-f_{n_{1}\bm{k}})K_{\scriptsize\begin{array}[]{l}n_{1}n_{2}\bm{k}\\\
n_{3}n_{4}\bm{k}^{\prime}\end{array}}(\bm{q}),\\\ &\\\
K_{\scriptsize\begin{array}[]{l}n_{1}n_{2}\bm{k}\\\
n_{3}n_{4}\bm{k}^{\prime}\end{array}}(\bm{q})=&\int{\rm d}\bm{r}_{1}{\rm
d}\bm{r}_{2}\psi_{n_{1}}(\bm{r}_{1})\psi^{*}_{n_{3}}(\bm{r}_{1})W^{\rm
RPA}(\bm{r}_{1},\bm{r}_{2};\omega=0)\psi^{*}_{n_{2}}(\bm{r}_{2})\psi_{n_{4}}(\bm{r}_{2}).\end{array}$
(17)

The expression $(H-\omega)^{-1}$ can be expressed in the spectral
representation as:

$\left[H(\bm{q})-\omega\right]_{ss^{\prime}}^{-1}=\sum_{\lambda\lambda^{\prime}}\frac{A_{s}^{\lambda}(\bm{q})N_{\lambda,\lambda^{\prime}}^{-1}(\bm{q})A_{s^{\prime}}^{*\lambda^{\prime}}(\bm{q})}{E_{\lambda}(\bm{q})-\omega\pm{\rm
i}\eta},$ (18)

where $A_{s}^{\lambda}(\bm{q})$ is element $s=n_{1}n_{2}\bm{k}$ of the
eigenvector of $H(\bm{q})$ with corresponding eigenvalue $E_{\lambda}(\bm{q})$
and $N$ is the overlap matrix. When the Tamm-Dancoff approximation (TDA) is
adopted, Grüning  _et al._ (2009) $H$ is Hermitian and Eq. (18) reduces to
$\displaystyle\sum_{\lambda}\displaystyle\frac{A_{s}^{\lambda}(\bm{q})A_{s^{\prime}}^{*\lambda}(\bm{q})}{E_{\lambda}(\bm{q})-\omega\pm{\rm
i}\eta}$.

The polarization in Eq. 16 can then be expressed in real-space according to
Eq. 15 and contracted to its two-point form. This two-point polarization is
then used in Eq. 4 to obtain $W$ with ladder diagrams included. The updated
$W$ is then used in the expression for the self-energy with the vertex
$\Lambda$ neglected (this has been shown to effectively cancel with the
renormalization $Z$-factor, as discussed in Appendix A of Ref. Kotani  _et
al._ , 2007). The role of $\Lambda$ in the expression for the self-energy may
indeed play an important role in certain systems, for example, as we will
discuss in CuCl.

###  II.4 Numerical evaluation of the kernel matrix elements

Our numerical implementation of the BSE relies on a generalization of the
linear muffin-tin orbital basis.Methfessel  _et al._ (2000); Kotani  _et al._
(2007); Kotani and Van Schilfgaarde (2010) The eigenfunctions are expanded in
Bloch-summed muffin-tin orbitals in spheres around atom centers. The radial
part of the eigenfunctions in these spheres is expanded by numerical solutions
of the radial Schrödinger equation. In the region between the spheres, the
eigenfunctions are then expanded in either smoothed Hankel functionsKotani and
Van Schilfgaarde (2010) and/or plane waves. Expanding the interstitial in
plane waves, the eigenfunctions are

$\Psi_{n\bm{k}}(\bm{r})=\sum_{\bm{R}u}\alpha_{\bm{R}u}^{\bm{k}n}\varphi_{\bm{R}u}^{\bm{k}}(\bm{r})+\sum_{\bm{G}}\beta_{\bm{G}}^{\bm{k}n}P_{\bm{G}}^{\bm{k}}(\bm{r}),$
(19)

where $\bm{R}$ denotes the atomic site and $u$ is a composite index that
contains the angular momentum of the site along with an index that denotes
either: a numerical solution of the radial Schrödinger equation at some
representative energy; its energy derivative (since the energy dependence has
been linearized by expanding in a Taylor series about the representative
energyAndersen (1975)); or a local orbital which is a solution at an energy
well above or below the representative energy. In $GW$ and the BSE a basis is
required that expands the product of eigenfunctions. Expanding the
interstitial in plane waves, the product eigenfunctions will also be expanded
in plane waves, and within the spheres the basis is expanded by
$\varphi_{Ru}(\bm{r})\times\varphi_{Ru^{\prime}}(\bm{r})$. This mixed product
basis (MPB) is denoted $M_{I}^{\bm{k}}(\bm{r})$.

Using the notation in Ref. Kotani  _et al._ , 2007, the kernel $K$ in the MPB
is read as

$\displaystyle K_{{\scriptsize\begin{array}[]{l}n_{1}n_{2}\bm{k}\\\
n_{3}n_{4}\bm{k}^{\prime}\end{array}\normalsize}}(\bm{q})=\displaystyle\sum_{I,J}\langle\psi_{n_{3},\bm{k}^{\prime}}|\psi_{n_{1},\bm{k}}\widetilde{M}_{I}^{\bm{k}^{\prime}-\bm{k}}\rangle
W_{IJ}^{\rm
RPA}(\bm{k}^{\prime}-\bm{k};\omega=0)\times\langle\widetilde{M}_{J}^{\bm{k}^{\prime}-\bm{k}}\psi_{n_{2},\bm{k}+\bm{q}}|\psi_{n_{4},\bm{k}^{\prime}+\bm{q}}\rangle,$
(20)

where the matrix elements and $W_{IJ}^{\mathrm{R}PA}$ are calculated as in
Ref. Kotani  _et al._ , 2007.

Due to the huge computational demands of the BSE only a subset of transitions
that occur between bands within a selected energy range about the Fermi level
are considered. Contributions from transitions not included in the BSE are,
however, included at the level of the RPA. To include such contributions, we
effectively have a matrix $H$ that is diagonal except for a block
corresponding to the coupled transitions discussed above. To calculate the
polarization in this case, the RPA contribution from the subset of states that
are treated at the level of the BSE are not included in the full $P^{RPA}$
[calculated according to Eq. (32) in Ref. Kotani  _et al._ , 2007] and the
contribution from $P^{BSE}$ is added to $P^{RPA}$. The corrected polarization
is then transformed into the MPB and the dielectric matrix $\epsilon=1-vP$,
and hence $W=\epsilon^{-1}v$, are thus obtained. The so-obtained
$W_{IJ}(\bm{q},\omega)$ is used as in Eq. (34) of Ref. Kotani  _et al._ , 2007
to calculate the correlation part of the self-energy.

###  II.5 Divergence of the macroscopic dielectric function at $\bm{q}=0$

The macroscopic ($\bm{G}=\bm{G}^{\prime}=0$) dielectric function (head of the
dielectric matrix) is constructed from the divergent bare Coulomb interaction
($4\pi/|\bm{k}|^{2}$) and polarization function. Since the dielectric matrix
contains a three dimensional integral over $\bm{k}$, the dielectric matrix for
$\bm{k}=0$ itself remains finite but angular dependent; resulting in the
dielectric tensor. In this work we employ the offset $\Gamma$ method,Kotani
_et al._ (2007); Kotani (2014); Friedrich  _et al._ (2010) to treat the
divergent part of $W$, where an auxillary mesh is introduced that is shifted
from the original $\Gamma$ centered mesh. The averaged macroscopic dielectric
function calculated in a small cell around $\Gamma$ is then used to calculate
the macroscopic part of the screened Coulomb interaction for
$\bm{k}\rightarrow 0$, as in Ref. Kotani, 2014.

###  II.6 Including the Frölich contribution to the band gap

To correct the value for the band gap in this method due to the neglection of
electron-phonon interactions we can include an approximation for the
contribution from the Frölich contribution to the Fan term, which – in polar
insulators such as LiF—should be the dominant part. We include lattice
polarization corrections (LPC) using the method outlined in Ref. Lambrecht
_et al._ , 2017. The energy shift is determined from the following equation

$\Delta E_{n\bm{k}}=\frac{e^{2}}{4}\sqrt{\frac{2\omega_{\rm LO}m_{\rm
eff}}{\hbar}}\left(\epsilon_{\infty}^{-1}-\epsilon_{0}^{-1}\right),$ (21)

where $\epsilon_{\infty}$ is the ion-clamped static (optical) dielectric
constant and $\epsilon_{0}$ contains effects accounting for nuclear
relaxations. $\omega_{\rm LO}$ is the frequency for the longitudinal phonon
for which we consider the electron-phonon contribution. The values for
$\epsilon_{0}$, $\epsilon_{\infty}$ and $\omega_{\rm LO}$ used in this work
are taken from Refs Lambrecht  _et al._ , 2017; Bhandari  _et al._ , 2018. The
$m_{\rm eff}$ in Eq. 21 is the effective mass. Questaal’sque  built-in band-
edge program can be used for determining effective masses, as in Refs
Lambrecht  _et al._ , 2017; Bhandari  _et al._ , 2018.

##  III Results and Discussion

We apply the QS$G\hat{W}$ to calculate the electronic structure and dielectric
properties of a broad range of materials that were listed in the introduction.
We first discuss in detail a few prototypical materials: Si, SrTiO3 (Sec.
III.2), LiF (Sec. III.3), NiO and CoO (Sec. III.4). We then analyze the
performance of QS$G\hat{W}$ in predicting the band gap and the static
dielectric constant (Sec. III.5).

###  III.1 Computational Details

All results have been obtained using Questaal que . Table 4 contains the
relevant parameters used in the calculations. The $l$ cut-off for partial
waves inside muffin-tin spheres was set to $4$ and an $spdf$ $spd$ (double-
Kappa type) basis was used in all calculations. Local orbitals were also used
in some systems as indicated in table 4. Empty sites were used as placements
for additional site-centered Hankel functions without augmentation, to improve
the basis in systems with large interstital voids. When calculating the
polarization within the RPA, the tetrahedron methodKotani  _et al._ (2007) is
employed for integration over the Brillouin zone. In the BSE implementation, a
broadening was applied according to Eq. 18 and set to 0.01 Ha for vertex
calculations. The RPA and BSE approaches were used for calculating the optical
absorption spectra of LiF and NiO, the broadening was set to match the
experimental broadening. These further details are given in the corresponding
Figures’ caption.

|  $a$(Å) |  $c$(Å) |  $G_{\rm Max}$(Ry) | $|\bm{q}+\bm{G}|_{\psi}$ | $|\bm{q}+\bm{G}|_{V}$ |  $\Sigma_{\rm cut}$(Ry) | $N_{\bm{k}}$ | $N_{v}$ | $N_{c}$ | $\varphi_{\rm loc}$
---|---|---|---|---|---|---|---|---|---|---
Si | 5.43 | - | 5.0 | 3.4 | 2.8 | 3.50 | 8 | 4 | 8 | 4$p$, 4$d$
C | 3.57 | - | 7.8 | 5.4 | 4.2 | 7.50 | 8 | 3 | 5 | -
Ge | 5.66 | - | 4.8 | 3.3 | 2.7 | 3.30 | 8 | 4 | 5 | -
LiF | 4.03 | - | 9.8 | 4.6 | 3.7 | 5.98 | 8 | 4* | 4 | -
LiCl | 5.13 | - | 9.3 | 3.6 | 2.9 | 3.90 | 8 | 4* | 4 | -
NaCl | 5.62 | - | 9.0 | 3.3 | 2.7 | 3.30 | 8 | 4* | 4 | -
MgO | 4.22 | - | 8.7 | 3.5 | 2.8 | 3.70 | 8 | 4* | 4 | -
GaAs | 5.65 | - | 8.1 | 3.3 | 2.7 | 3.30 | 8 | 3 | 11 | 3$d$ Ga
GaSb | 6.09 | - | 8.5 | 3 | 2.5 | 2.90 | 6 | 9 | 11 |  | 3$d$ Ga
---
4$d$ Sb
InSb | 6.48 | - | 9.4 | 2.7 | 2.3 | 2.60 | 10 | 4 | 5 |  | 5$d$, 6$p$ 1/2 In
---
4$d$, 6$p$ 1/2 Sb
InP | 5.87 | - | 9.2 | 3.1 | 2.6 | 3.08 | 8 | 3 | 5 | 5$d$ In
InAs | 6.05 | - | 9.2 | 2.7 | 2.3 | 2.93 | 8 | 3 | 8 |  | 5$d$, 6$p$ 1/2 In
---
3$d$, 6$p$ 1/2 As
CuCl | 5.64 | - | 9.2 | 3.4 | 2.8 | 3.54 | 8 | 8 | 4 |  | 4$d$ Cu
---
4$s$, 4$p$ Cl
CdTe | 6.48 | - | 9.2 | 2.7 | 2.3 | 2.62 | 8 | 3 | 5 |  | 5$d$ Cd
---
4$d$ Te
CdO | 4.70 | - | 10.0 | 3.9 | 3.2 | 3.50 | 10 | 3 | 2 |  | 4$p$, 5$d$ Cd
---
3$s$, 3$p$ O
ZnO | 3.25 | 5.21 | 10.2 | 3.7 | 3.0 | 4.02 | 9,9,6 | 5 | 5 |  | 5$s$, 5$p$, 4$d$ Zn
---
3$s$, 3$p$, 4$d$ O
CuAlO2 | 3.47 | - | 9.2 | 4.3 | 3.5 | 4.00 | 6 | 10 | 4 | 4$d$ Cu
h-BN | 2.50 | 6.64 | 9.3 | 3.5 | 2.9 | 5.65 | 12,12,6 | 6 | 8 | -
w-InN | 3.53 | 5.69 | 9.2 | 3.4 | 2.7 | 3.46 | 9,9,6 | 6 | 5 |  | 5$d$ In
---
3$s$, 3$p$ N
w-GaN | 3.18 | 5.19 | 9.2 | 3.7 | 3 | 4.11 | 9,9,6 | 6 | 5 |  | 4$d$ Ga
---
3$s$, 3$p$ N
SrTiO3 | 3.9 | - | 9.2 | 3.5 | 2.8 | 3.20 | 6 | 9 | 11 |  | 4$s$, 4$p$ Sr
---
3$p$, 4$d$ Ti
NiO | 4.17 | - | 11 | 4.4 | 3.6 | 4.40 | 8 | 16 | 16 | 3$p$, 4$d$ Ni
CoO | 4.25 | - | 9.2 | 3.5 | 2.8 | 3.60 | 6 | 12 | 6 | 4$d$ Co
Table 1: Parameters used in the calculations: lattice constant $a$ (and $c$
where applicable); energy cut-off for the plane wave basis set $G_{\rm MAX}$;
cut-off for the eigenfunctions in units of $2\pi/a$; cut-off for the Coulomb
interaction in units of $2\pi/a$; the cut-off for the self-energy333The self-
energy is approximated with an average diagonal value for energies above
$\Sigma_{\rm cut}-0.5~{}{\rm Ry}$.; the number of $\bm{k}$-points along the
reciprocal lattice vectors $N_{\bm{k}}$ (where 3 numbers correspond to
different divisions along the reciprocal lattice vectors); the number of
valence444Note that ∗ indicates that all valence states are used in the BSE.
$N_{v}$ and conduction $N_{c}$ states used in the BSE. Note that states not
included in the BSE are included in the polarization at the level of the RPA.
$\varphi_{\rm loc}$ are the local orbitals used (see text). Waves marked as
$p\,1/2$ are local orbitals computed using the Dirac equation, as discussed in
Ref. Pashov  _et al._ (2020). This improves the matrix elements for spin-orbit
coupling.

Bulk Silicon is one of the most widely studied semiconductors and usually is
one of the prototypical systems when testing and developing novel methods and
codes for electronic structure calculations of bulk solids. QS$G\hat{W}$
predicts band gaps of 3.31 eV ($\Gamma$-$\Gamma$), 2.22 eV ($\Gamma$-L), 1.281
eV ($\Gamma$-$X$) and 1.47 eV (fundamental). These gaps are very close to
experimentOrtega and Himpsel (1993); Kittel (2005) and about 0.04 eV below the
values determined from QS$GW$ (see Fig. 3 in Ref. Kotani  _et al._ , 2007), as
one would anticipate. The value for $\epsilon_{\infty}$ determined within the
BSE and with the QS$G\hat{W}$ electronic structure is 10.2 about 10% less than
11.4 observed 0 KVuye  _et al._ (1993), and reported in Ref. Shishkin  _et
al._ , 2007. Kutepov reported a slighly larger value for the fundamental gap,
1.41 eV, close to what Shishkin et. al reportedShishkin  _et al._ (2007).
There could be many reasons for this: Kutepov for example implements a
somewhat different expression than Eq. 1, as previously noted. This
discrepancy seems to be in keeping with the general tendency for Kutepov’s
reported gaps to be slightly larger than our results. Further, there is a
linearization error in the augmentation spheres first noted by Friedrich and
coworkers Friedrich  _et al._ (2006) which in ZnO increases the gap by 0.25
eV, and in that case is not adequately accounted for by a single local
orbital. Finally we use an atom-centered basis, which has an essentially
infinite plane-wave cutoff, but unlike the LAPW method, the different
$\mathbf{G}$ components do not vary fully independently. We have found that
results change only slightly as we continue to enlarge the basis beyond the
sizes we have selected van Schilfgaarde  _et al._ (2006a), but it does not
entirely rule out the possibility of basis set incompleteness in the
interstitial region, in addition to the augmentation region. Finally there are
complexities dealing with the product basis and the coulomb singularity in the
polarization as $q{\rightarrow}0$, which are difficult to assess. In any case
the differences are on the order of 0.1 eV, which is approximately the
uncertainty of our implementation. Part of the overestimation is due to the
approximations adopted555 e.g. the use of the non-interacting Green’s
function; adopting the $GW$ approximation when determining the interaction
kernel [along with ignoring the contribution from $\delta W/\delta G$ in the
interaction kernel and also assuming this to be static]., and part is due to
the finite $\bm{k}$-grid.

In fact, the convergence of the dielectric response within the BSE with the
number of $\bm{k}$-points in the Brillouin zone is very slow for Si, mainly
because one needs to account for excitations between nearly parallel bands.
This is shown in Fig. 1 for the ion-clamped static (optical) dielectric
constant, $\epsilon_{\infty}$, determined within the BSE using, for
simplicity, the LDA electronic structure as input. As we can see from the
inset in Fig. 1 the value for $\epsilon_{\infty}$ converges for Si when using
just 3 valence and 4 conduction bands in the BSE (the remaining valence bands
and 60 extra conduction bands are included at the level of the RPA, see Sec.
II.3). On the other hand, from the main panel it is clear that a dense
$\bm{k}$-grid is needed to converge $\epsilon_{\infty}$.

For the QS$G\hat{W}$ calculation we chose a $12\times 12\times 12$ regular
$\bm{k}$-mesh, which is fine enough to produce a value that is relatively
close to the $1/N_{k}$ limit, whilst making the simulation feasible (thus
easily reproducible) without the need for huge computational resources. From
the figure, this grid overestimates the extrapolated value of
$\epsilon_{\infty}$ (i.e. for $1/N_{k}\rightarrow 0$) by about 0.4.

Figure 1: Main panel: ion-clamped static (optical) dielectric constant,
$\epsilon_{\infty}$, as a function of the total number of $\bm{k}$-points in
the Brillouin zone using a regular $\Gamma$-centered(blue) and shifted(violet)
$\bm{k}$-grid. Inset: $\epsilon_{\infty}$ as a function of: (black) the number
of conduction bands in the BSE with all 4 valence states used; and (red) the
number of valence bands used with 4 conduction bands. Note the BSE is solved
using the LDA electronic structure for simplicity.

We expect that the overestimation of $\epsilon_{\infty}$ due to the
underconverged $\bm{k}$-mesh leads to band gaps that are slightly
underestimated with respect to the $1/N_{k}$ limit. In fact, the vertex
correction (that is the difference between the band gap obtained with
QS$G\hat{W}$ and QS$GW$) decreases with increasing $N_{k}$, so that—although
the QS$GW$ gaps decrease with increasing $N_{k}$ (Fig. 3 in Ref. Kotani  _et
al._ , 2007)—the shift due to the inclusion of ladders also decreases,
resulting in band gaps that tend to converge from below as $N_{k}$ is
increased.

In most systems investigated the values for $\epsilon_{\infty}$ tend to
converge with fewer $\bm{k}$-points than for Si; Si is chosen here as an
extreme case. Considering the convergence behavior in Fig. 1, we determine a
numerically accurate value for $\epsilon_{\infty}$ by choosing a relatively
coarse $\bm{k}$-mesh and perform the QS$\bm{GW}_{\Lambda}$ calculation. Then,
using the converged $H_{0}$ from QS$\bm{GW}_{\Lambda}$ another calculation is
performed to obtain the dielectric function using a slightly finer mesh. We
then extrapolate to—in principle—the value on an infinitely dense
$\bm{k}$-grid.

###  III.2 SrTiO3

SrTiO3 is a perovskite material that may exist in the usual different
perovskite phases: cubic, tetragonal and orthorhombic. The experimental band
gap is 3.25 eV (indirect) and 3.75 eV (indirect),van Benthem  _et al._ (2001)
and according to Bhandari et al. in Ref. Bhandari  _et al._ , 2018 is almost
independent of the structure. We investigate SrTiO3 in the cubic phase with
lattice parameter $3.9~{}{\rm\AA}$. As in previous works we include some of
the semi-core states in the valence; the Sr 4s, 4p and Ti 3p and 4d states are
included as valence states here.666We choose a $6\times 6\times
6~{}\bm{k}$-grid, which is finer than the mesh used in Ref. Bhandari  _et al._
, 2018 since QS$G\hat{W}$ usually requires a finer mesh than QS$GW$ for the
same level of convergence. The basis set parameters and energy cut-offs used
in the calculation (Table 4) are similar to those discussed in Sec. IV of Ref.
Bhandari  _et al._ , 2018.

We found a LDA band gap of 1.77eV which is increased to 4.24 eV (4.67 eV at
$\Gamma$) with QS$GW$. Within QS$G\hat{W}$, the inclusion of ladder-diagrams
reduces the band gap to 3.77 eV (4.19 eV at $\Gamma$). Finally, the LPC 777The
values for $\epsilon_{\infty}$, $\epsilon_{0}$, $\omega_{\rm LO}$ and $m_{\rm
eff}$ used in Eq. 21 were taken from Ref. Bhandari  _et al._ , 2018, which
uses Questaal’s band-edge program to determine the effective masses. As is the
case with Si, the value for $\epsilon_{\infty}$ calculated in this work is
overestimated—comparing to the value quoted in Ref. Bhandari  _et al._ , 2018:
6.76 vs 5.52. Note that the effective masses determined using the QS$G\hat{W}$
band structure are the same, which can be seen since the ladders effectively
shift the bands down, and the band edges – in particular – have effectively
the same shape. gives a band gap renormalization of -0.404 eV, resulting in
band gaps of 3.37 eV (indirect) and 3.79 eV (at $\Gamma$), in very good
agreement with experiment. van Benthem  _et al._ (2001)

Figure 2: Band structure of cubic SrTiO3 for QS$GW$ (blue dashes) and
QS$G\hat{W}$ (solid lines). The bands are positioned so that the valence band
maxima are aligned in the two cases.

Figure 2 compares the band stuctures for QS$GW$ and QS$G\hat{W}$. The main
effect of the inclusion of the ladder diagrams are: an almost rigid shift
downwards of the conduction bands and a narrowing of the valence bandwidth.
The latter is also seen in band narrowing beyond QS$GW$ in highly correlated
systems, such as Sr2RuO4 Acharya  _et al._ (2019). One can observe as well
more subtle effects, such as the removal of the accidental degeneracy in the
lowest lying valence band at $M$, and the feature occuring in the lowest lying
conduction band at $X$, which, interestingly, is affected in a similar way
when lattice polarization corrections are included to the QS$GW$electronic
structure. Bhandari  _et al._ (2018) $W{\rightarrow}\hat{W}$ scarcely modifies
the shape of the highest lying valence bands in this system, which is
important for accurate determination of valence band offsets Chambers  _et
al._ (2004).

###  III.3 LiF

The macroscopic dielectric function of the polar insulator LiF was recently
calculated in Ref. Cunningham  _et al._ , 2018 within the BSE applied on top
of the QS$GW$ electronic structure. Since the vertex corrections are neglected
in QS$GW$ the screening of the exchange contribution was underestimated.
Together with the neglection of electron-phonon effects, this results in a
greatly overestimated band gap of $\sim$16.2 eV, i.e. about 2 eV larger than
the experimental value. The underscreening also caused an overestimation of
about 0.5–1 eV of the exciton binding energy. As a result of the partial
cancellation of these errors, producing the optical absorption spectrum using
the BSE with the QS$GW$ electronic structure results in a blue shift of $\sim$
0.9 eV with respect to experiment.

Figure 3: Imaginary part of the macroscopic dielectric function for LiF. The
experimental data (blue squares) Roessler and Walker (1967) is compared with
the results from the QS$G\hat{W}$ method. The spectrum is red-shifted by 0.483
eV to account for lattice polarization effects (see text). The spectral
broadening was increased linearly to match that of the experimental spectrum:
the broadening at the first excitonic peak is 0.15 eV and at the X exciton
peak ($\sim$22 eV) is 1.25 eV.

Here, we repeat the BSE calculation of the optical spectrum but on top of the
QS$G\hat{W}$ electronic structure and also consider the lattice polarization
effects. Including ladder-diagram vertex-corrections produces a fundamental
band gap of 14.71 eV; a reduction of over 1.4 eV, in agreement with the vertex
correction calculated in Ref. Shishkin  _et al._ , 2007. Including the 0.483
eV polaron shift correction from Ref. Lambrecht  _et al._ , 2017 gives then a
band gap in excellent agreement with the experimental value of 14.2$\pm$0.2
eV.Piacentini  _et al._ (1976) The exciton binding energy is around 2 eV, also
in agreement with the experimental valuePiacentini  _et al._ (1976) and as a
consequence Fig. 3 shows an excellent overall agreement between the
theoretical and experimental spectra.

| RPA | BSE |
---|---|---|---
| $G_{0}W_{0}$ | QS$GW$ | $G_{0}W_{0}$ | QS$GW$ | QS$G\hat{W}$ | Exp.
$\epsilon_{\infty}$ | 1.61 | 1.71 | 1.76 | 1.84 | 1.96 | 1.90-1.95
Table 2: Electronic part of the static dielectric constant,
$\epsilon_{\infty}$, for LiF. Values at different levels of theory are
compared with the experimental result. Lines (1990)

Table 2 reports the values for the ion-clamped static dielectric constant,
$\epsilon_{\infty}$, at the various levels of theory as in Ref. Cunningham
_et al._ , 2018. Calculating the macroscopic dielectric function at the BSE
level improved noticeably the agreement with the experimental value when
comparing to the RPA value. When the electronic structure is calculated at the
QS$GW$ rather than $G_{0}W_{0}$ level, the agreement with experiment is
improved. When the electronic structure is then calculated with the
QS$G\hat{W}$ electronic structure, the value is slightly overestimated ($1.96$
versus the experimental value of $1.90-1.95$).

###  III.4 NiO and CoO

The transition metal monoxides NiO and CoO are antiferromagnetic with
intermediate size band gaps and relatively large magnetic moments on the
transition metal site. The experimental values for the band gaps and magnetic
moments—along with values calculated for these quantities at the different
levels of theory—are given in table 3. Experimental values are taken from
Refs. Faleev  _et al._ , 2004; Sawatzky and Allen, 1984; Khler, 1993; Bredow
and Gerson, 2000 and references therein.

| Magnetic moment ($\mu_{\rm B}$) | Bandgap (eV)
---|---|---
| LDA | QS$GW$($\hat{W}$) | Exp | LDA | QS$GW$($\hat{W}$) | Exp
NiO | 1.23 | 1.71 (1.67) | 1.64–1.9 | 0.4 | 5.07 (4.02) | 3.8 - 4.3
CoO | $\sim 2.4$ | 2.72 (2.71) | 2.38–2.8 | 0 | 4.1 (3.1) |  $\sim$ 2.9
Table 3: Values calculated for the fundamental band gaps and magnetic moments
for the strong correlated antiferromagnetic transition metal monoxides NiO and
CoO along with experimental values for these quantities. The LDA, QS$GW$ and
vertex corrected (QS$G\hat{W}$) values are presented and compared with
experiment.Faleev  _et al._ (2004); Sawatzky and Allen (1984); Khler (1993);
Bredow and Gerson (2000)

Mean-field approaches, such as the LDA, fail to accurately predict band gaps
in these materials and poorly describe even qualitative features such as
whether the material is metallic or insulating. Conduction band dispersions
and the nature of the fundamental band gap (e.g., direct or indirect) are
poorly described. This was reinforced in Ref. Kotani  _et al._ , 2007 that
investigated the NiO electronic structure in depth within QS$GW$. In that
work, the calculated LDA fundamental band gap of NiO was about 0.4 eV, in
agreement with that found here; severely underestimated with respect to the
experimental value of 3.8–4.3 eV.Faleev  _et al._ (2004); Sawatzky and Allen
(1984); Bredow and Gerson (2000) Further, the analysis in that work of the LDA
and $GW$ (with eigenvalue-only self-consistency) bandstructures highlighted
the need for self-consistency in the eigenfunctions as well as the
eigenvalues. However, the selfconsistency at the QS$GW$ level gives an
indirect bandgap for NiO of 4.86 eV (direct gap of 5.56 eV at $\Gamma$),
overestimated by about 0.5–1.0 eV with respect to experiment (see table 3).
The main sources of errors at the QS$GW$ level are the missing vertex
corrections and temperature effects. The QS$G\hat{W}$ calculated bandgaps
(including vertex corrections) are 4.27 eV (fundamental) and 4.85 eV (at
$\Gamma$) in excellent agreement with the experimental value.

Similarly, for CoO, the LDA predicts a semi-metallic behaviour; i.e., a band
gap equal to zero in contrast with the experimental value of about 2.9 eV. At
the QS$GW$ level the calculated indirect gap is 4.08 eV (direct gap of 5.24 eV
at $\Gamma$), overestimated by about 1.5 eV. At the QS$G\hat{W}$ level, the
overestimation is reduced substantially: the calculated band gaps are 3.1 eV
(fundamental) and 4.08 eV (at $\Gamma$).

Along with the calculated and experimental band gaps, Table 3 also presents
calculated magnetic moments. As for the band gaps, the LDA underestimates the
experimental magnetic moments. The underestimation is corrected at the QS$GW$
level: $1.23~{}\mu_{\rm B}$ versus $1.71~{}\mu_{\rm B}$ at QS$GW$ for NiO and
$2.4~{}\mu_{\rm B}$ versus $2.72~{}\mu_{\rm B}$ for CoO. In both cases the
QS$GW$ calculated value fall into the range of the experimental results. In
contrast with the significant differences observed for the band gaps, the
differences between the magnetic moments calculated within QS$G\hat{W}$ and
QS$GW$ are negligible, of the order 0.01 (less than 1%).

Figure 4 shows the orbital character analysis based on a Mulliken
decomposition and compares the QS$GW$ and QS$G\hat{W}$ band structures for NiO
(top) and CoO (bottom). For the NiO, the Ni $d$-like bands (red) are the most
affected by the vetex corrections; the shift on the lowest $d$-like conduction
bands are almost 1 eV, whereas the $s$\- and $p$-like bands are shifted by
around 0.5 eV. The same trend is observed in CoO: 1.1 eV for $sp$-bands versus
1.9 eV for Co $d$-like bands.

Figure 5 shows the optical absorption spectrum of NiO calculated within the
BSE using either the QS$GW$ (from Ref. Cunningham  _et al._ , 2018) and the
QS$G\hat{W}$ electronic structures (this work). Calculations based on the
QS$G\hat{W}$ electronic structure show a marked improvement over the QS$GW$,
by correctly blue-shifting by $\sim$1 eV the position of the onset and first
peak.

Figure 4: Left panels: QS$G\hat{W}$ band structures with orbital contributions
(based on a Mulliken decomposition) for NiO (top) and CoO (bottom). Red is
used for bands with a predominantly Ni(Co) $d$-orbital character; green for
bands with a Ni(Co) and O $p$-like character. Right panels: comparison of the
band structures at the QS$GW$ (dotted blue lines) and QS$G\hat{W}$ (solid
lines) levels for NiO (top) and CoO (bottom). By comparing the left and right
panels for each material one can see the orbital character of the bands that
are most affected by the vertex corrections  Figure 5: Imaginary part of the
macroscopic dielectric function for NiO. The experimental data (blue
squares)Powell and Spicer (1970) is presented along with spectra calculated
within the BSE using a quasiparticlized self-energy at the level of QS$GW$ and
QS$G\hat{W}$. Also shown is the RPA dielectric function within QS$GW$. The
onset is blue-shifted, as is typical when vertex corrections are omitted: the
electron-hole attraction redshifts the peak by 1.7 eV, and is redshift by an
additional 1.3 eV in the QS$G\hat{W}$ case. We use the method outlined in Ref.
Cunningham  _et al._ , 2018 for calculating the $\bm{q}\rightarrow 0$
transition dipole matrix elements with a contribution from the non-local self-
energy. Figure 6: Calculated electronic band gaps versus experimental band
gaps for: (left) materials with band gaps in the range 0 - 2 eV; and (right)
materials with band gaps in the range 2.5 - 10 eV. Different levels of theory
are considered, QS$GW$ (green circles), QS$G\hat{W}$ (magenta squares) and
QS$G\hat{W}$ with a further correction accounting for lattice polarization
effects (black crosses). Both vertex corrections (included at the QS$G\hat{W}$
level) and the lattice polarization corrections tend to reduce the gap. In
general, considering these corrections reduces the overestimation observed for
the QS$GW$ results and bring the results closer to the experimental value.

###  III.5 Overall performance

Figure 6 compares the fundamental band gaps calculated at different levels of
theory with the experimental values. With the exception of CuCl, the QS$GW$
consistently overestimates the experimental band gap. This overestimation is
(partially) corrected at the QS$G\hat{W}$ bringing the results closer to the
experimental results. With few exceptions, the QS$G\hat{W}$ slightly
overestimates the experiemental band gaps (to within an acceptable tolerance).
Finally, for some polar materials the LPC was evaluated. This corrections is
negative and further reduces the calculated band gap. In general the
QS$G\hat{W}$+LPC values are in close agreement with the experimental value.
This agreement shows that, as previously argued, the overestimation in the
QS$GW$ is due mostly to the neglection of the vertex in $W$ and of phonon
effects.

CuCl is the only material considered for which QS$G\hat{W}$ significantly
underestimates the experimental value. We argue that this can be traced to
omitting the vertex corrections in the self-energy, Eq. 2. In fact, as a
consequence of this omission occupied states are systematically too shallow.
Further, this error increases with the localization of the state.Grüneis  _et
al._ (2014) Then in CuCl, the too shallow Cu $d$ state intermixes too strongly
with the Cl $p$ states. The Cl $p$ states are then pushed to a too high
energy, reducing the gap.

Note that for simplicity we did not include the error in the experimental
values which comes either from the measurement and from the sample
variability. As an example, the band gap of InN varies greatly from sample to
sample and thus it is not straightforward to determine experimentally. While
current estimates are between 0.6–0.7 eV, old estimates were as high as 1.9
eV!Tansley and Foley (1986) The value we found at the QS$G\hat{W}$+LPC level,
0.65 eV, agrees well with the band gap estimates from the some recent
experimentsYu  _et al._ (2005); Gallinat  _et al._ (2006); Higashiwaki and
Matsui (2004).

Figure 7 presents the static dielectric constants, $\epsilon_{\infty}$,
calculated at the BSE level using the QS$G\hat{W}$ electronic structure.
Results are plotted against the experimental values. In general, we found a a
very good agreement between the calculated and experimental values.

Figure 7: $\epsilon_{\infty}$ calculated from the the BSE starting from a
QS$G\hat{W}$ one-particle hamiltonian.

##  IV Conclusions

We proposed an approach based on the so-called QS$GW$ in which we include
vertex corrections to $W$, calculated at the level of the the BSE, and
lattice-polarization effects including the correction from Ref. Lambrecht  _et
al._ , 2017. This approach, QS$G\hat{W}$+LPC, sought to correct the systematic
overestimation of the bandgap in QS$GW$ which has been indeed attributed to
the omission of both vertex-corrections and lattice-polarization. We assessed
the performance of this approach on an extensive sample of materials. The
sample includes the common $sp$ semiconductors and insulators—for which state-
of-the-art $G_{0}W_{0}$ performed well—as well as materials, such as NiO and
CoO, for which $G_{0}W_{0}$ performed poorly.

The results confirm that the omission of vertex corrections and lattice
polarization is responsible for most of the error in the QS$GW$ band gaps and
that including these effects as proposed in this work produces an accurate
electronic structure from which the band gap, as well as dielectric and
optical properties can be reliably estimated. As the only exception we
presented CuCl for which the gap is significantly underestimated by
QS$G\hat{W}$+LPC. For this system we think it is important to include vertex
corrections to the self-energy as well. As discussed in Sec. III.5, the Cu $d$
state is too shallow and intermixes too strongly with the Cl $p$ states,
resulting in these being pushed up too high and reducing the band gap. The
method presnted in this work could be improved with including dynamical
effects in the BSE kernelMarini and Del Sole (2003), as well as extending the
method to go beyond the TDAGrüning  _et al._ (2009) (note that the
calculations performed in this work are already incredibly computationally
expensive and going beyond the TDA for some of the systems studied here would
not currently be practical). Finally, a more accurate inclusion of electron-
phonon effects in general systems (i.e., not just polar insulators) is
required to further enhance the method, and this is an area of active
research.

###### Acknowledgements.

The authors would like to thank all those involved in the CCP flagship
project: Quasiparticle Self-Consistent $GW$ for Next-Generation Electronic
Structure, especially Scott Mckechnie for his help. We are grateful for
support from the Engineering and Physical Sciences Research Council, under
grant EP/M011631/1. MvS was supported in part by the Simons Foundation. We are
grateful to the UK Materials and Molecular Modelling Hub for computational
resources, which is partially funded by EPSRC (EP/P020194/1)

## References

  * Hedin (1965a) L. Hedin, Phys. Rev. 139, A796 (1965a).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Hybertsen and Louie (1986) M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 (1986).
  * van Schilfgaarde  _et al._ (2006a) M. van Schilfgaarde, T. Kotani, and S. V. Faleev, Phys. Rev. B 74, 245125 (2006a).
  * Kotani  _et al._ (2007) T. Kotani, M. van Schilfgaarde, and S. V. Faleev, Phys. Rev. B 76, 165106 (2007).
  * Faleev  _et al._ (2004) S. V. Faleev, M. van Schilfgaarde, and T. Kotani, Phys. Rev. Lett. 93, 126406 (2004).
  * van Schilfgaarde  _et al._ (2006b) M. van Schilfgaarde, T. Kotani, and S. Faleev, Phys. Rev. Lett. 96, 226402 (2006b).
  * Sottile  _et al._ (2003) F. Sottile, V. Olevano, and L. Reining, Phys. Rev. Lett. 91, 056402 (2003).
  * Marini  _et al._ (2003) A. Marini, R. Del Sole, and A. Rubio, Phys. Rev. Lett. 91, 256402 (2003).
  * Lambrecht  _et al._ (2017) W. R. L. Lambrecht, C. Bhandari, and M. van Schilfgaarde, Phys. Rev. Materials 1, 043802 (2017).
  * Franchini  _et al._ (2010) C. Franchini, A. Sanna, M. Marsman, and G. Kresse, Phys. Rev. B 81, 085213 (2010).
  * Shishkin  _et al._ (2007) M. Shishkin, M. Marsman, and G. Kresse, Phys. Rev. Lett. 99, 246403 (2007).
  * Bruneval  _et al._ (2005) F. Bruneval, F. Sottile, V. Olevano, R. Del Sole, and L. Reining, Phys. Rev. Lett. 94, 186402 (2005).
  * Kutepov (2016) A. L. Kutepov, Phys. Rev. B 94, 155101 (2016).
  * Kutepov (2017) A. L. Kutepov, Phys. Rev. B 95, 195120 (2017).
  * Ismail-Beigi (2017) S. Ismail-Beigi, J. Phys.: Condens. Matter 29, 385501 (2017).
  * (18) “Questaal code website,” https://www.questaal.org, accessed: 2017-07-04. Our $GW$ implementation was adapted from the original ecalj package, now at https://github.com/tkotani/ecalj/.
  * Hedin (1965b) L. Hedin, Phys. Rev. 139, A796 (1965b).
  * Schwinger (1951) J. Schwinger, Proceedings of the National Academy of Sciences 37, 452 (1951), http://www.pnas.org/content/37/7/452.full.pdf .
  * Martin and Schwinger (1959) P. C. Martin and J. Schwinger, Phys. Rev. 115, 1342 (1959).
  * Onida  _et al._ (2002) G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. 74, 601 (2002).
  * Schleife  _et al._ (2009) A. Schleife, F. Fuchs, C. Rdl, J. Furthmller, and F. Bechstedt, physica status solidi (b) 246, 2150 (2009).
  * Anisimov  _et al._ (1997) V. I. Anisimov, F. Aryasetiawan, and A. Lichtenstein, Journal of Physics: Condensed Matter 9, 767 (1997).
  * Cunningham  _et al._ (2018) B. Cunningham, M. Grüning, P. Azarhoosh, D. Pashov, and M. van Schilfgaarde, Phys. Rev. Materials 2, 034603 (2018).
  * Salpeter and Bethe (1951) E. E. Salpeter and H. A. Bethe, Phys. Rev. 84, 1232 (1951).
  * Marini and Del Sole (2003) A. Marini and R. Del Sole, Phys. Rev. Lett. 91, 176402 (2003).
  * Grüning  _et al._ (2009) M. Grüning, A. Marini, and X. Gonze, Nano Letters 9, 2820 (2009), pMID: 19637906, http://dx.doi.org/10.1021/nl803717g .
  * Methfessel  _et al._ (2000) M. Methfessel, M. van Schilfgaarde, and R. A. Casali, “A full-potential lmto method based on smooth hankel functions,” in  _Electronic Structure and Physical Properies of Solids: The Uses of the LMTO Method Lectures of a Workshop Held at Mont Saint Odile, France, October 2–5,1998_, edited by H. Dreyssé (Springer Berlin Heidelberg, Berlin, Heidelberg, 2000) pp. 114–147.
  * Kotani and Van Schilfgaarde (2010) T. Kotani and M. Van Schilfgaarde, Physical Review B (Condensed Matter and Materials Physics) 81, 125117 (2010).
  * Andersen (1975) O. K. Andersen, Phys. Rev. B 12, 3060 (1975).
  * Kotani (2014) T. Kotani, Journal of the Physical Society of Japan 83, 094711 (2014), http://dx.doi.org/10.7566/JPSJ.83.094711 .
  * Friedrich  _et al._ (2010) C. Friedrich, S. Blügel, and A. Schindlmayr, Phys. Rev. B 81, 125102 (2010).
  * Bhandari  _et al._ (2018) C. Bhandari, M. van Schilgaarde, T. Kotani, and W. R. L. Lambrecht, Phys. Rev. Materials 2, 013807 (2018).
  * Pashov  _et al._ (2020) D. Pashov, S. Acharya, W. R. L. Lambrecht, J. Jackson, K. D. Belashchenko, A. Chantis, F. Jamet, and M. van Schilfgaarde, Comp. Phys. Comm. 249, 107065 (2020).
  * Ortega and Himpsel (1993) J. E. Ortega and F. J. Himpsel, Phys. Rev. B 47, 2130 (1993).
  * Kittel (2005) C. Kittel,  _Introduction to solid state physics_, 8th ed. (Hoboken, NJ Wiley, 2005) includes bibliographical references and index. ”Chapter 18, Nanostructures, was written by Professor Paul McEuen …”.
  * Vuye  _et al._ (1993) G. Vuye, S. Fisson, V. Nguyen Van, Y. Wang, J. Rivory, and F. Abelès, Thin Solid Films 233, 166 (1993).
  * Friedrich  _et al._ (2006) C. Friedrich, A. Schindlmayr, S. Blügel, and T. Kotani, Physical Review B (Condensed Matter and Materials Physics) 74, 045104 (2006).
  * van Benthem  _et al._ (2001) K. van Benthem, C. Elsässer, and R. H. French, Journal of Applied Physics 90, 6156 (2001), https://doi.org/10.1063/1.1415766 .
  * Acharya  _et al._ (2019) S. Acharya, D. Pashov, C. Weber, H. Park, L. Sponza, and M. van Schilfgaarde, Communications Physics 2, 163 (2019).
  * Chambers  _et al._ (2004) S. Chambers, T. Droubay, T. Kaspar, M. Gutowski, and M. van Schilfgaarde, Surface Science 554, 81 (2004).
  * Roessler and Walker (1967) D. M. Roessler and W. C. Walker, J. Opt. Soc. Am. 57, 835 (1967).
  * Piacentini  _et al._ (1976) M. Piacentini, D. W. Lynch, and C. G. Olson, Phys. Rev. B 13, 5530 (1976).
  * Lines (1990) M. E. Lines, Phys. Rev. B 41, 3372 (1990).
  * Sawatzky and Allen (1984) G. A. Sawatzky and J. W. Allen, Phys. Rev. Lett. 53, 2339 (1984).
  * Khler (1993) J. Khler, Angewandte Chemie 105, 477 (1993), https://onlinelibrary.wiley.com/doi/pdf/10.1002/ange.19931050352 .
  * Bredow and Gerson (2000) T. Bredow and A. R. Gerson, Phys. Rev. B 61, 5194 (2000).
  * Powell and Spicer (1970) R. J. Powell and W. E. Spicer, Phys. Rev. B 2, 2182 (1970).
  * Grüneis  _et al._ (2014) A. Grüneis, G. Kresse, Y. Hinuma, and F. Oba, Phys. Rev. Lett. 112, 096401 (2014).
  * Tansley and Foley (1986) T. L. Tansley and C. P. Foley, Journal of Applied Physics 59, 3241 (1986), https://doi.org/10.1063/1.336906 .
  * Yu  _et al._ (2005) K. M. Yu, Z. Liliental-Weber, W. Walukiewicz, W. Shan, J. W. Ager, S. X. Li, R. E. Jones, E. E. Haller, H. Lu, and W. J. Schaff, Applied Physics Letters 86, 071910 (2005), https://aip.scitation.org/doi/pdf/10.1063/1.1861513 .
  * Gallinat  _et al._ (2006) C. S. Gallinat, G. Koblmüller, J. S. Brown, S. Bernardis, J. S. Speck, G. D. Chern, E. D. Readinger, H. Shen, and M. Wraback, Applied Physics Letters 89, 032109 (2006), https://doi.org/10.1063/1.2234274 .
  * Higashiwaki and Matsui (2004) M. Higashiwaki and T. Matsui, Journal of Crystal Growth 269, 162 (2004), proceedings of the First ONR International Indium Nitride Workshop.