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

We present an extension of the quasiparticle self-consistent $GW$ approximation (QS$GW$) [Phys. Rev. B, 76 165106 (2007)] to include vertex corrections in the screened Coulomb interaction $W$. This is achieved by solving the Bethe-Salpeter equation for the polarization matrix at all $k$-points in the Brillouin zone. We refer to this method as QS$G\hat{W}$. QS$GW$ yields a reasonable and consistent description of the electronic structure and optical response, but systematic errors in several properties appear, notably a tendency to overestimate insulating bandgaps, blue-shift plasmon peaks in the imaginary part of the dielectric function, and underestimate the dielectric constant $\epsilon_{\infty}$. A primary objective of this paper is to assess to what extent including ladder diagrams in $W$ ameliorates systematic errors for insulators in the QS$GW$ approximation. For benchmarking we consider about 40 well understood semiconductors, and also examine a variety of less well characterized nonmagnetic systems, six antiferromagnetic oxides, and the ferrimagnet Fe$_3$O$_4$. We find ladders ameliorate shortcomings in QS$GW$ to a remarkable degree in both the one-body Green's function and the dielectric function for a wide range of insulators. New discrepancies with experiment appear, and a key aim of this paper is to establish to what extent the errors are systematic and can be traced to diagrams missing from the theory. One key finding of this work is to establish a relation between the bandgap and the dielectric constant $\epsilon_{\infty}$. Good description of both properties together provides a much more robust benchmark than either alone. We show how this information can be used to improve our understanding of the one-particle spectral properties in materials systems such as SrTiO$_3$ and FeO.


I. INTRODUCTION
The one-particle Green's function G(r, r ′ , ω) provides essential information about material properties.Besides having value in its own right, determining both ground state properties (total energy, charge and magnetic densities), and excitation energies, it is the starting point for transport and other two-particle properties, e.g.spin and charge response functions, and superconductivity.
As a consequence, knowledge of G is of the first importance, and a vast amount of effort has been dedicated to finding prescriptions to yield G both efficiently and with high-fidelity ab initio (without recourse to models or adjustable parameters).Density-functional theory [1] (DFT), where the electron density n replaces G as the fundamental variable, is an alternative, and indeed it is far more popular because of its efficiency and good scaling with system size.DFT is a ground-state theory, but it generates an auxilliary one-body H 0 , with fictitious eigenvalues and eigenfunctions.H 0 often provides a reasonable approximation to excitations of the real system, but it is often unsatisfactory, e.g. its notorious tendency to underestimate splitting between occupied and unoccupied levels.Wave function methods, widely used in quantum chemistry, but less so in materials physics, use the single-particle orbitals ψ i as the fundamental variable.They can provide high-fidelity solutions to the many-body Schrodinger equation.As Walter Kohn noted in his Nobel prize lecture [2], wave function methods contain more information than is needed or useful, but nevertheless require concomitant effort needed to compute observables.For that reason they are expensive and scale poorly with system size.Also, spectral properties are not readily computed.
Green's function (GF) methods lie between the two: G has more information than n but less than the wave functions.As with DFT, G-based methods create an effective one-body potential Σ(r, r ′ , ω), but differ in that Σ is nonlocal and energydependent.They are computationally more intensive than DFT; however they can be made to scale reasonably well with system size, and because of their better fidelity it is likely they will ultimately outphase DFT methods for many functional materials, particularly when excitations are involved.Thus, GF theories might be called the "Goldilocks" approach.GF methods possess a key advantage in another respect: dynamical screening becomes the predominant many-body effect for systems involving many atoms.Hedin's equations [3] can be expanded diagramatically in powers of the screened coulomb interaction W, and encapsulate this phenomenon in a natural way, even in the lowest order (GW).The traditional target applications are also different: quantum chemical methods focus mostly on ground state properties while GF methods focus on spectral properties, especially two-particle spectra.GF methods do not yet possess the fidelity of wave function methods, and to what extent their fidelity can eventually approach them remains a key open question. 1ne key aim of this paper is to provide a partial answer to this question.The quasiparticle self-consistent GW approximation (QSGW) provides an effective way to implement GW theory without relying on a lower level approximation as a starting point.This makes discrepancies with experimental data much more uniform, and it is essential to distinguish errors intrinsic to the theory itself, from accidents as a result of the starting point.We assess in some detail the extent to which discrepancies QSGW displays with experiment can be mitigated by adding the low-order diagram (ladder diagram) to the random-phase approximation (RPA) 2 for the bare polarizability.As noted, the GW approximation is the lowest order diagram in the manybody perturbation theory (MBPT) of Hedin [3], and while GW shows significant improvement over DFT (including functionals designed to surmount the well-gap underestimate [7]), it has well known problems.First, it is a perturbation theory, which typically starts from some reference noninteracting G 0 and generates a correction to it.The most common choice of G 0 is one based in DFT, but many kinds of choices have been made to improve on the final result.This situation is unsatisfactory in two respects: • G 0 can be (and often is) tuned to improve agreement with experiments.G 0 plays the role of a free parameter, and in this sense the theory isn't really ab initio any more.
• The errors inherent in low-order MBPT, e.g.GW, can be masked by the arbitrariness in G 0 .Sometimes qualitiatively wrong conclusions can be drawn, or good agreement with experiment found, but for the wrong reason.This is a quite common, albeit not well appreciated, difficulty with the theory; see for example Ref. [8].
By employing the GW approximation in the QSGW form, we can circumvent these difficulties.QSGW is a procedure where G 0 is determined self-consistently.Self-consistency is used not to minimize the total energy, but instead some measure (norm) of the difference between G −1 0 and G −1 [9].With a definition for optimal construction for G 0 , it surmounts the ambiguities from arbitrariness in the starting-point [10].It provides a good and systematic G 0 so that discrepancies with experiment that appear tend to be similarly systematic, making it possible to associate these discrepancies with diagrams missing from the theory.
There is no unique definition of the norm, but one intuitively appealing definition leads to a static (quasiparticlized) self-energy Σ 0 generated from the dynamical one as i and j are eigenstates of the one-particle hamiltonian.Ismail-Beigi showed this construction satisfies a variational principle, not for the total energy but its gradient [11].One other important consequence of Eq. ( 1) is that at self-consistency the poles of G and the poles of G 0 coincide: thus in contrast to DFT, the energy bands of H 0 generated by QSGW correspond to true excitations of the system.Results generated by QSGW may sometimes worsen agreement with experiment over other forms of GW.For example, ε ∞ generated from a Kohn-Sham band structure is often better than the QSGW one.We will argue that this stems from a fortutious cancellation of errors; see §III B 3. Fortutitous error cancellation in QSGW is much less pronounced, and as a result, discrepancies with experiment are better exposed, and moreover they are much more uniform.Several of the most salient discrepancies are connected to the inadequate description of the dielectric polarizability.This forms the primary motivation for the present work: to make a detailed assessment of how the simplest extension to the RPA polarizability improves both G and the dielectric response.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) [12,13] for the polarization.A high-fidelity G is essential for a good description of any response function, including the magnetic one, as shown for NiO [14] and for yttrium iron garnet [15]; and the particle-particle correlation function that governs superconductivity (see e.g.Ref. [16]).
Other works have considered the effect of vertex corrections to the dielectric screening on the band gap.For example, Refs.17-19 included ladder diagrams through an effective nonlocal static kernel constructed within time-dependent density functional theory to mimic the BSE.More recently, Kutepov proposed several schemes for the self-consistent solution of Hedin's equations including vertex corrections [20].In particular, for selected semiconductors and insulators [21], 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, Σ = iGW Γ.In all the cases, an improvement over LQSGW was observed (see §II D 2 for a comparison between LQSGW and QSGW ).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 20 and 21. 3 Here we omit the first-order vertex for Σ, in keeping with our present objective -to find the best single-Slater determinant construction.The QSGW philosophy incorporates this vertex in an approximate way, via a Ward identity in Γ that goes as 1/Z in the q→0, ω→0 limit, cancelling the Z factor that is the predominant difference between the quasiparticlized G 0 and the interacting G [22].
The approach described in this work, the standard GWA and the QSGW, are all derived from the many-body perturbative approach developed by Hedin [3].In this method, the following set of closed coupled equations [25][26][27] are to be solved iteratively: where G is the Green's function, v(r, r ′ ) = 1/|r − r ′ | is the bare Coulomb interaction, W is the screened Coulomb interaction, Γ is the irreducible vertex function, P is the irreducible polarizability (the functional derivative of the induced density with respect to the total potential) and Σ is the self-energy operator.In Eqs. ( 2) and ( 5), the indices subsume position and time and the + superscript implies t ′ = t + η, with η → 0 + .In the standard GWA, also known as one-shot GW or G 0 W 0 , Eq. ( 6) is approximated as: in both the expressions for the self-energy [Eq.(2)] and the irreducible polarization [Eq.(5)].In addition, Σ and P are both evaluated for G = G 0 , the independent-particles Green's function, 4 that in the frequency-domain takes the form In Eq. ( 8), ψ n and ε 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 variously as the independent particle approximation, the time-dependent Hartree approximation, and the Random-Phase-Approximation (RPA) [28].When using Eq. ( 8) in frequency space it takes the form: 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 In Eq. ( 10), V ext (r) is the external potential due to nuclei and external fields, V H (r) is the Hartree potential describing the classical mean-field electron-electron interaction, and V XC (r) is the exchange-correlation potential describing correlation effects missing in V H (r). ρ is the electronic density calculated as occ |ψ n (r)| 2 , that is constructed from the eigensolutions of H 0 (from which the self-consistent solution is obtained).Because of V 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 ∆Σ = Σ − V XC to avoid double counting.Eq. 3 is rewritten as a nonlinear equation for the quasiparticle energies E nk and solved after it is linearized: where the renormalization factor Z nk is: B. Ladder diagrams in W In previous works employing QSGW [9,22,29,30] the RPA was used to make W .This leads to errors noted in the introduction, 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 [31].To include the vertex in Hedin's equations we need to determine the interaction kernel δΣ/δG in Eq. 6. Suppose we have adopted the GW approximation for Σ 5 and assume that δW/δG is negligible [28], then δΣ(12)/δG (45) = iW (12)δ (1,4)δ (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)δ(1, 3)δ (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; where P RPA (1234) = −iG(13)G (42).The W that appears in the interaction kernel, δΣ/δG, is calculated at the level of the RPA and this is usually assumed to be static, 6 i.e., δΣ/δG = iW RPA (ω=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 where we have again combined band, spin and wavevector indices, and ). Inserting the expression for the RPA polarization from Eq. ( 9) in Eq. ( 14), one arrives at the following expression for the polarization whereby the conservation of momentum we have k 2(4) = k 1(3) + q; and7 The expression (H − ω) −1 can be expressed in the spectral representation as: where A λ s (q) is element s = n 1 n 2 k of the eigenvector of H(q) with corresponding eigenvalue E λ (q) and N (q) is the overlap matrix.When the Tamm-Dancoff approximation (TDA) is adopted, [33] H is Hermitian and Eq. ( 18) reduces to 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.In what follows we will denote W RPA with the symbol W , and refer to W with ladders included as Ŵ .The updated W or Ŵ is then used in the expression for the self-energy with the vertex Γ in the exact self-energy (iGW Γ) omitted.The justification for omission, and the consequences of it, is taken up in §II D.
In works that report optical absorption α(ω), we construct it using the relation where ε is calculated from the macroscopic part of the dielectric matrix ε = 1 − vP .
C. Self-consistency: QSGW In the DFT based G 0 W 0 approximation, the E nk 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.It is often improved in practice by choosing some other G 0 constructed, e.g. from a hybrid functional.Another route is to replace the corrected energy E n [Eq.(12)] either in the Green's function [34] [Eq.( 8)], or in the RPA polarization [Eq.( 9)] entering the screened potential W , or in both.
Here we use the QSGW approach in which the starting point is chosen to effectively minimize ∆Σ: the difference between the dynamical self-energy and (static) quasiparticlized one.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 XC in Eq. ( 10) with Σ 0 , Eq. 1, 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, and equally important, it removes the arbitrariness in starting point [10].Fig. 1 shows a flow chart of the process.These are used to construct the non-interacting Green's function G 0 , Coulomb interaction v, RPA polarization P 0 = − iG 0 G 0 , and W =(1 − vP 0 ) −1 v. W is used to make a vertex and better P via Eq.14, which gives the improved Ŵ .One cycle makes the static self-energy Σ 0 , that is passed to H 0 (green) and the cycle repeated to self-consistency.

Why self-consistency is important
Self-consistency is not typically performed in weakly correlated materials.LDA-based GW can do very well (see description of Bi 2 Te 3 , §III E 2) but self-consistency improves the theory and makes the discrepancies with experiment systematic.Recent work shows this to be the case even for simple sp metals such as Li, Na, and Mg [35] where errors in RPA for W are likely to be less important than in insulating systems.Similarly excellent agreement is found for the Fermi liquid regime of Fe, in a detailed study examining several properties [36].
Perhaps the first study applying GW to a correlated material was the work of Aryasetiawan and O. Gunnarsson [34], in which case a starting point better than the LDA becomes essential.This issue arises for many kinds of narrow-band systems, and even in weakly or moderately correlated ones the starting point can be important.Narrow-gap semiconductors in which the LDA has a negative gap offer one notable illustration of this.Using GW in the usual manner (correcting the reference eigenvalues via Eq.12) cannot correct the wrong topology of the starting point [37].InN is a classic example (Fig. 2).Even while the states at the k-point Γ have the correct ordering, the improper initial ordering leads to unphysical dispersions in the band structure in the vicinity of Γ.Other systems which fall into this class are Ge, PbTe, InAs, and InSb.(In PbTe, a gap appears at L within the LDA, but with L + 6 and L − 6 wrongly ordered; see Fig. 13 of Ref. [38]).Another effect of self-consistency is to modify the one-body part of H or G.This is because not only the eigenvalues but the density is significantly renormalized relative to the LDA.GW induces a corresponding change in the effective potential through the inverse of the susceptibility, χ −1 (x 1 , x 2 ) = δV (x 1 )/δn(x 2 ).Starting from the perturbation LDA ], and if we assume that χ −1 (x 1 , x 2 ) is adequately approximated by the LDA, we can estimate the change in n from δn = χ 0 δV 0 , and from this obtain the attendant screening potential as Here V Hxc is the combined Hartree + (LDA) V xc .In practice the Questaal codes execute an operation similar to this in the Low-energy band structure of wurzite InN at four different levels of approximation.Colors depict orbital character of the bands: red for N pxy character, green for N pz character, blue for In s character.LDA bands are shown in top left panel: the state of In s character at Γ lies below the three states of N p character, reflecting an inverted gap.Upper right panel shows the effect of GW treated perturbatively from the LDA, i.e.Eq. 12 with Z=1 for reasons explained in the text.GW rectifies the inverted gap at Γ, but without off-diagonal parts of Σ it cannot undo the wrong topology given by the starting point, and thus the bands cross near Γ.Bottom left panel is the classic QSGW result, Ref. 9. It provides a good description of the InN energy bands; however the gap is overestimated (1.01 eV) relative to experiment (0.68 eV).The QSGW dielectric constant ε∞ is calculated to be 6.1, about 3/4 of the measured value (8.4).Bottom right panel is the QSGW result with W augmented by ladder diagrams.The gap (0.74 eV) is slightly larger than experimental one, and differs by approximately the electron phonon interaction (estimated to be 0.07 eV [39]).
natural course of self-consistency: an internal loop is performed in the one-body code adding the fixed Σ QSGW as an external potential and making the density self-consistent.This accelerates convergence to self-consistency, but for the present we use that process to estimate the effect of δV scr on the bandgap.In Table I we compare the bandgaps for GW generated from the LDA, in various forms.It compares the usual Eq. ( 12) (with Z=1), GW including the full matrix structure of Σ without updates to the one-body hamiltonian [40]; an estimate for the change in one-body potential as just described; and finally QSGW.A key take-away is that the off-diagonal parts of Σ are unimportant only in the simplest nearly homogeneous systems, such as Ge.Even in SrTiO 3 , a simple d 0 transition metal compound, they are significant, modifying the eigenvalues both directly, and indirectly through changes in the density.Another important finding is that both direct and indirect contributions vary widely in both magnitude and sign, and indeed the change is often larger than the well-recognized need to account for the electron-phonon contribution [41].
One solution is to perform partial eigenvalue-only self-consistency: i.e. use Eq. ( 12) in a self-consistent manner by updating the eigenvalues without changing eigenfunctions.There is a simpler way to approximate eigenvalue-only self-consistency by simply omitting the Z factor in Eq. (12).This was shown to be rigorously true for a two-level system in the Appendix of Ref. 37. Eigenvalue-only self-consistency can significantly reduce the discrepancies with experiments, but it cannot resolve the topology problems or the modifications to the density noted above.Further, the off-diagonal parts of Σ can have nontrivial effects on the quasiparticle spectrum, as noted for example in the discussion around Fig. 2. Another solution is to choose a better starting point, e.g. based on an extension of LDA -as for example using a hybrid DFT approach [42] or LDA+U [43] -or the Coulomb-hole Screened exchange approximation (COHSEX) [3].Since the starting-point dependence can be chosen freely, the theory loses its ab initio flavor.This freedom is lost with QSGW, and errors that appear better reflect the nature of the approximations made.TABLE I. Dependence of the bandgap on various kinds of treatment in the off-diagonal parts of the self-energy.In all cases the starting hamiltonian is the LDA.EG(Σ diag ) is the outcome from a treatment similar to the usual way GW is employed (Eq.12 but with Z= 1).EG(n0) adds the full Σ − V LDA xc to the LDA hamiltonian, including the off-diagonal elements, but without updating the density.EG(n0+δn) is similar to EG(n0) but the density is updated in a "small" loop keeping Σ fixed, as described in the text around Eq. (20).QSGW is the quasiparticle self-consistent result (QSG Ŵ result in parenthesis).Values reported for Ge, GaSb, and TiSe2 are for the direct gap at Γ, with TiSe2 in the high-temperature P 3m1 phase.In each of these three cases, the valence and conduction band edge states are inverted in the LDA, similar to Fig. 2. QSGW is already a good approximation in many systems, but it is well known that discrepancies with experiment appear.They tend to be very systematic, and mostly related to the RPA approximation to W. Bandgaps being systematically overestimated, the high-frequency dielectric constant ε ∞ underestimated, and blue shifts in peaks in Im ε(ω) -all fairly universal with QSGW -are connected to the RPA approximation to W. It has long been known, starting from independent work in the groups of Louie [44] and of Reining [45], that if the RPA is extended to include ladder diagrams, optical response is significantly improved in simple semiconductors.
Our primary focus here is to determine to what extent ladder diagrams in W ameliorate these discrepancies.As we will show here, when W is extended to Ŵ and the cycle carried through to self-consistency, many of the systematic errors in the QSGW RPA self-energy are ameliorated to a remarkable degree for a wide range of weakly and strongly correlated insulators.(We restrict outselves to insulators since that is where ladders are most important [46].)While this is encouraging, some discrepancies remain, and these form a major focus of this paper.A very important feature of QSGW RPA has been that when discrepancies with experiment appear, the origin can often be clearly associated with a particular missing diagram, enabling the possibility for a systematic, hierarchical extension of the theory.We will show that this remains mostly true with Σ = iG Ŵ : QSG Ŵ improves on QSGW but systematic errors remain.The following omissions account for many of the shortcomings in results presented in this paper.

Shortcomings in QSG Ŵ
1.The electron-phonon interaction is a well-identified contribution to the self-energy and, if lattice vibrations are in the harmonic approximation, consists of two contributions (Fan and Debye Waller terms) [3].The diagram usually reduces insulating bandgaps; it also is needed to capture optical transitions between states of different wave numbers, e.g. in indirect gap semiconductors.For its effect on the index of refraction, see §III D.
2. Omission of Γ in the exact self-energy iGW Γ.As noted in the last paragraph of the introduction, there is a partial renormalization Z-factor connecting G to G 0 (see Appendix A of Ref. 9), which we rely on in the QSG Ŵ approximation.Typically this vertex pushes down all the states in an approximately uniform manner with a minimal effect on the bandgap [47].The effect is more pronounced for a nearly dispersionless d or f state, and when such a state comprises the valence band maximum, the gap is underestimated.Semicore d states in semiconductors such as CdTe and GaSb lie about 0.7 eV above photoemission experiments ( Fig 11).Also, bandgaps in materials systems whose valence band consist of a 3d state, or a strong admixture of it, tend to be too small (Table IV).An extreme manifestation of this is EuO: the valence band maximum consists of a nearly dispersionless, atomic-like f state, and as a result the QSGW gap is underestimated [48].From this calculation it was inferred that the Eu 4f state should be pushed down by ∼0.7 eV -a somewhat larger shift than for a flat 3d state (presumably it is even more atomic like).A shift in core-like d levels of order 0.5 eV was first explicitly demonstrated by Grüneis et al., who introduced a simple first-order vertex into G LDA W LDA [47].Very recently Kutepov added a first-order vertex in a somewhat more rigorous manner [49].
As regards the present work the most important error seems to occur with systems with shallow core-like levels, particularly when they occur near the valence band maximum.See §III D 2 and also the discussion around Table IV for instances where this neglect is important.
3. Higher order diagrams in the polarizability.The interaction kernel W (Eq. 14) is taken from the RPA and moreover it is assumed to be static.Other diagrams, have been considered in a few works, e.g. the second order screened exchange [50].This diagram when augmenting the RPA, was quite successfully used to predict total energies in chemical systems [6].We consider only one additional diagram, namely to use W BSE as the kernel in generating Eq. 14, and note its effects on a few systems ( §III C 1).
4. Inadequate treatment of spin fluctuations.In the theory presented here, the only spin contribution to the self-energy comes from the Fock exchange.We present some spectral properties of correlated antiferromagnetic insulators ( §III F), and show that even in such correlated cases, the response in the charge channel seems to be reasonably described.This is likely because, in contrast to spin fluctuations, charge fluctuations sense the long range coulomb interaction.The situation may be different when the gap closes or becomes small on the scale of spin excitations ( 0.1 eV).In such cases there may be cross coupling between spin and charge channels.Our solution to date has been to augment QSGW with Dynamical Mean Field Theory (DMFT).DMFT is a nonperturbative method and exact solutions are possible with e.g.Continuous Time Quantum Monte Carlo [51] that include all diagrams.However the vertex is assumed to be local, which is reasonable for spin fluctuations as the vertex is thought to predominantly reside on-site among the correlated orbitals where the fluctuations occur [52].Indeed the QSGW ++ framework ( ++ referring generally to extensions of QSGW), augmented either by ladder diagrams or by DMFT, does seem to have unprecedented predictive power in a number of strongly correlated materials [16,36,[53][54][55][56][57].Yet there are places where a nonlocal vertex may be important, e.g. to explain the nematic phase of FeSe.Some approaches have been formulated to improve on DMFT, e.g. the "DΓA" approximation -a nonperturbative, semilocal approach [58], but it is extremely demanding in practice.
In addition to the T matrix [52], somewhat more sophisticated low-order diagrams that treat spin fluctuations on the same footing as charge fluctuations have been proposed [59], but this has not been attempted yet in an ab initio context.As spin fluctuations tend to be low energy, many channels are possible, so whether a low-order theory is sufficient or not remains an open question.A low-order perturbation theory that could replace DMFT would be very advantageous, since DMFT has its own unique set of challenges.We do not consider such cases in the present work, but it should be noted that the claim that charge fluctuations are well described already in low order is not universal [60] and whether a low-order perturbative theory can be sufficient remains an open question in low-density and strongly correlated metals.
Perhaps surprisingly, this obvious deficiency does not seem to play a significant role in the systems we study here.The present work considers only systems with bandgaps, and the likely explanation is that spin wave frequencies are typically small energy compared to the optical gaps, which suppresses spin fluctuations.
Other discrepancies with experiment in zincblende semiconductors will be presented that do not appear to have a simple interpretation.Perhaps the most notable unexplained error are the errors in the k-dependent dispersion of the conduction band minimum in zincblende semiconductors ( §III C).Such systems are weakly correlated and the origin of the error is not readily explained.One significant possibility is that Questaal's present implementation does not include a scheme to make the basis set truly complete [61].This would not be a limitation of the theory itself, but in its implementation.
Finally, several outliers are noted often because the distinction between optical gaps and fundamental gaps is ignored, e.g. in ScN ( §III E 3), SrTiO 3 ( §III E 6), TiO 2 ( §III E 5), and CuAlO 2 ( §III E 7), or are likely artifacts of inaccurate experiments, e.g.hBN (see discussion around Table VII), and in correlated systems where the experiments are less reliable.FeO seems to be an extreme example of this ( §III F 4).The connection between discrepancies in one-particle properties and those in two-particle properties is discussed in §III D 2.

Relation between LQSGW and QSGW
We noted earlier that Kutepov constructed both a quasiparticlized scheme and a fully self-consistent one.His quasiparticlized scheme (LQSGW ) is somewhat different from QSGW .They are similar, but with the extensions to RPA presented here and in his work, the fidelity becomes high enough that the difference can be significant.To show this we present here a brief analysis of the relation between LQSGW and QSGW .
Kutepov quasiparticlizes the self-energy with Σ(ω=0), but folds in an effective energy dependence through Σ ′ (ω=0) (LQSGW ) while preserving the ability to construct a static hamiltonian.The Appendix derives a rough estimate for the expected difference in QP levels between LQSGW and QSGW , obtaining leading contribution from the omitted quadratic term; (Eq.31).In the QSGW scheme, the diagonal element of static (quasiparticlized) Σ 0 nn is by construction equal to dynamical self-energy at the quasiparticle level ω qp , so Σ nn (ω) = Σ QSGW nn at ω=ω qp .For LQSGW this is no longer true; thus the quasiparticle levels of the static G 0 do not coincide with the poles of G.We can estimate the difference between the LQSGW and QSGW QP levels from the difference between Σ nn (ω qp ) and the linear approximation to it, Σ nn (ω=0) + ω qp • Σ ′ nn (ω=0).This is depicted in Fig. 3 for the first band in Na, and the highest valence band and lowest conduction band in NiO.In Na, it corresponds to the QSGW -LQSGW change in bandwidth; in NiO, the change is the QSGW -LQSGW difference in the direct gap at X.For the Na case, according to the simple perturbative expression, Eq.31, LQSGW and QSGW should differ by 0.11 eV.A better estimate is the difference noted in the previous paragraph.The graphs of Fig. 3 indicate that the LQSGW bandwidth in Na should be slightly larger than QSGW , and the NiO bandgap also slightly larger.Numerically the difference in self-energies in the Na case, at the QSGW QP energy, is 0.17 eV.According to first order perturbation theory, this is the expected difference between LQSGW and QSGW QP energies.Indeed the discrepancy between LQSGW and QSGW appears to be of this order: one LQSGW and two QSGW calculations have been reported for Na8 (see Table in Fig. 3).0.17 eV is similar to the spread between QSGW and a recent photoemission measurement [63]; see the Table in Fig. 3.As GW is known to break down at sufficiently low densities, an accurate determination of the bandwidth in Na is important since it is one of the best realizations of a nearly homogeneous low-density metal.

Hybrid QSGW self-energies
As Questaal has no implementation of the electron-phonon vertex as yet, or the vertex modifying GW , we cannot evaluate its effect ab initio.However, by perturbing slightly the QSG Ŵ self-energy with an admixture of the QSGW Σ when the gap is underestimated or LDA V xc when it is overestimated, we can modify to Σ h to reach a target bandgap without affecting the eigenfunctions too severely.That permits us to assess the effect of the error in E G on ε ∞ .Table VII presents cases where both E G and ε ∞ are well known, and it establishes that discrepancies in the two are intimately connected for several systems.That provides an independent confirmation that the one-body hamiltonian would be of high fidelity if this perturbation were properly included.
To make a reasonable proxy to the QSG Ŵ self-energy, e.g. for the electron-phonon self-energy and missing vertex noted in §II D 1 (points 1 and 2), we will construct a hybrid one-body self-energy Σ h defined as This equation has often been used with α=0, β=0.8, γ=0.2, because ε ∞ computed from QSGW has been found empirically to be very nearly 80% of the true value for a wide range of semiconductors (see Fig. 4b).This formula has been empirically found to yield very good bandgaps in many kinds of materials sytems [64,65].In §III D 2 we use it to show how the errors in Σ (whatever their origin) are closely connected to discrepancies in the dielectric function.A caveat should be noted here: even while the bandgap can be rendered accurate with such a hybrid self-energy, ǫ ∞ computed from the RPA is not similarly improved, so the theory cannot capture both quantities in a consistent manner.§III B 3 discusses this at greater length.
We will also assume ths connection to hold in cases where the fundamental gap is uncertain while ε ∞ is better known.By aligning ε ∞ , or the frequency-dependent ε(ω) with measured data, we can make a reasonable estimate for the fundamental gap.This is done for several systems, e.g.TiO 2 ( §III E 5) and FeO ( §III F 4).

E. Numerical evaluation of the kernel matrix elements
Our numerical implementation of the BSE relies on a generalization of the linear muffin-tin orbital basis [9,66,67].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 functions [67] and/or plane waves.Expanding the interstitial in plane waves, the eigenfunctions are where 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 energy [68]); 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 ϕ Ru (r) × ϕ Ru ′ (r).This mixed product basis (MPB) is denoted M k I (r).Using the notation in Ref. 9, the kernel K in the MPB is read as where the matrix elements and W RPA IJ are calculated as in Ref. 9. Owing 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. 9] and the contribution from P BSE is added to P RPA .The corrected polarization is then transformed into the MPB and the dielectric matrix ǫ = 1 − vP , and hence W = ǫ −1 v, are thus obtained.The so-obtained W IJ (q, ω) is used as in Eq. ( 34) of Ref. 9 to calculate the correlation part of the self-energy.The macroscopic (G = G ′ = 0) dielectric function (head of the dielectric matrix) is constructed from the divergent bare Coulomb interaction (4π/|k| 2 ) and polarization function.Since the dielectric matrix contains a three dimensional integral over k, the dielectric matrix for k = 0 itself remains finite but angular dependent; resulting in the dielectric tensor.In this work we employ the offset Γ method [9,69,70], to treat the divergent part of W , where an auxillary mesh is introduced that is shifted from the original Γ centered mesh.The averaged macroscopic dielectric function calculated in a small cell around Γ is then used to calculate the macroscopic part of the screened Coulomb interaction for k → 0, as in Ref. 69.
The G=G'=0 component of the irreducible polarizability should vanish at q=0.Owing to numerical errors this is not exactly the case, so to correct for this its value is subtracted from the irreducible polarizability for all q.This adjustment stabilizes the calculations and also improves on the k convergence of the polarizability and self-energy.We performed careful checks for the k-convergence in ε ∞ in the RPA, and found for example in zincblende semiconductors an 8×8×8 mesh was reasonably good, and a 12×12×12 mesh converged ε ∞ to ∼1% in all cases but the smallest gap semiconductors.

G. 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. 71.The energy shift is determined from a P is the polaron length scale, which, in the effective mass approximation is computed from the optical mode phonon frequency ω LO and the effective mass m * .a P is different for electrons and holes, and we take an average of the electron and hole contributions, following Ref.71. ǫ ∞ is the ion-clamped static (optical) dielectric constant and ǫ 0 contains effects accounting for nuclear relaxations.The values for ǫ 0 , ǫ ∞ and ω LO used in this work are taken from Refs 71 and 72 for materials discussed there.For many of the systems studied here a more rigorous calculations of the gap shift has been published (Ref.[41]).Where available we use those results.
For a given shift, we use a proxy Σ h (Eq.21) to estimate the effect on the band structure and ǫ ∞ .

H. Effective Oscillator model for Index of Refraction
In Ref. [73] it was established that the frequency-dependent index of refraction of many compounds can be fit reasonably well by a single oscillators model.The model has the form where E 0 is the oscillator energy, and E d is a measure of the strength of interband optical transitions.Empirically, (n 2 −1) −1 has been found to be a mostly linear function of ( ω) 2 , for a wide range of ionic materials, which lends credence to the model.
In some experiments where n(ω) is tabulated, we use Eq. 26 to extrapolate to ω=0.

III. RESULTS AND DISCUSSION
A. Computational Details All results have been obtained using Questaal [24].Table XII contains the relevant parameters used in the calculations.The ℓ cut-off for partial waves inside muffin-tin spheres was set to 4 and an spdfg-spdf basis was used in all calculations, except in some lighter systems where the g orbitals were omitted.Local orbitals were also used in some systems as indicated in table XII.Empty sites were used as placements for additional site-centered Hankel functions (to ℓ max = 2 or 3) without augmentation, to improve the basis in systems with large interstital voids.When calculating the polarization within the RPA, the tetrahedron method is employed for integration over the Brillouin zone [9].In the BSE implementation, a broadening was applied according to Eq. 18 and set to 0.01 Ha for vertex calculations.
The TDA was also adopted due the huge increase in compution required to store, calculate and diagonalize the non-Hermitian matrix that has twice as many rows and columns as the Hermitian TDA one.We would, however, not expect going beyond the TDA too have to much of an impact on the systems investigated in this work [33].We did remove the TDA in a few cases, e.g.InSb, ScN, and MgO, and found the effect to be minor, as anticipated ( §III C 1).

Treatment of the Screened Coulomb interaction
For numerical reasons the Questaal codes compute Fock matrix elements not of the bare coulomb interaction 1/q 2 but a slightly screened one, 1/(q 2 + κ 2 ).A small value for κ is taken, between 10 −5 and 10 −4 .The QSGW self-energy is not usually sensitive to the value of κ; however, the dielectric constant ε ∞ can vary by a few percent for κ ranging between 10 −5 and 10 −4 .For that reason, we compute ε ∞ for three values of κ between 10 −5 and 3×10 −5 and extrapolate to κ=0.
Also, to avoid evaluating matrix elements at q = 0, we use the offset-Γ method [9], which requires generating the polarizability for small values of q near zero.For ε ∞ we evaluate ε(ω, q) for three small finite values of q and extrapolate to q=0.The direction of approach to q=0 gives us the orientation dependence of ε(q = 0).
Both kinds of extrapolations are done in one process.The difference between extrapolated and finite-(q, κ) values can differ by a few percent.Where available, this was taken from Ref. [41]; otherwise it was estimated from the Frölich expression, Eq. ( 24).(bottom) ε RPA ∞ calculated from G0 generated QSGW (blue circles) and ε BSE ∞ calculated from G0 generated QSG Ŵ (green diamonds).The dark dashed line corresponds to perfect agreement with experimental data; The light dashed line corresponds to ε th ∞ /ε exp ∞ = 0.8.For hBN, we used ε in the basal plane; see Table VI.Bencharking ε∞ in the antiferromagnetic oxides CuO, MnO, FeO, CoO, NiO and Cu is more complex.They're omitted here but discussed in §III F.

B. Survey of Results
We begin with a birds-eye view of some key results.Fig. 4 shows the fundamental bandgaps (E G ) and high-frequency dielectric constant (ε ∞ ) for a wide variety of materials, comparing classical QSGW results to QSG Ŵ .This figure elucidates the general trends: QSGW tends to overestimate bandgaps, and underestimate ε ∞ by an almost universal constant factor 0.8.As anticipated, addition of ladder diagrams ameliorates both of these discrepancies.Apart from a few exceptions (see discussion in §III D 2).QSG Ŵ greatly improves on QSGW .On the wide scale of the figure the ability of QSG Ŵ to predict optical properties (Fig. 4) looks stellar, but discrepancies appear on closer inspection.A main theme of this paper is to seek out these deviations, and associate them, where possible, with the missing diagrams noted in Sec II D 1. Fortunately, most of the discrepancies with measured data are fairly systematic, which opens the possibility that the shortcomings can be rectified with relatively simple low-order diagrams.

Index of refraction
Typically, ε ∞ is obtained by extrapolating the frequency-dependent index of refraction ε(ω) to zero using, e.g., Eq. 26.Its value is known only to a resolution of a few percent even in the best cases, and the uncertainty is often larger.An extreme case is AlN, where several values have been reported ranging between 3.8 [74] and 4.8 [75], and hBN is another instance ( §III D).When several experimental values are available for the compounds in Fig. 4, we use an average value.Reported values for ε ∞ for antiferromagnetic transition metal oxides (not shown in Fig. 4) also show variations, and calculations show larger deviations from the average value.They are discussed in §III F.
There is a small indeterminacy on the theory side also.Besides the extrapolation noted in §III A, care must also be taken to converge the uniform k mesh entering into numerical Brillouin zone integration: for narrow-gap semiconductors pushing this mesh beyond 12 × 12 × 12 divisions was not feasible, leading to a slight tendency to underestimate ε ∞ .These approximations lead to an uncertainty of a few percent.
The bulk of the remaining paper focuses on discrepancies where either E G or ε ∞ are outside the experimental uncertainty, which appear in some systems.One primary aim of this work is to draw a connection between E G and ε ∞ .Generally speaking, in well-characterized systems the discrepancies in E G and ε ∞ occur largely at the same time: when the gap is accurately described, ε ∞ is also.We believe this to be a significant finding, and it is taken up in §III B 3. §III D and §III D 2 present cases where E G deviates the most strongly from experiment.We use hybrid self-energies (Eq.21) to correct the gap, to see how the change in ε ∞ tracks it.[77], compared to results calculated at three levels of approximation: RPA@QSGW , BSE@QSGW , and BSE@QSG Ŵ .As is typical, the shoulder of RPA@QSGW is blue-shifted, by roughly 2 eV in this case.

NiO as archetype system
Adding ladders (BSE@QSGW ) shifts the shoulder towards the experiment, but it is still ∼1 eV too high, as a consequence of overestimate of the QSGW fundamental gap (Table IX).BSE@QSG Ŵ describes Im ǫ(ω) rather well, including peaks around 6 eV, 13.5 eV and 17 eV.However the shoulder around 3.5 eV is slightly red-shifted compared to experiment, indicating that the fundamental gap is underestimated.ε∞ is also overestimated.See §III F 1 for more details.
NiO is an archetype system that exhibits many of the phenomena that are the subject of this work.Fig. 5 shows in greater detail how ladder diagrams renormalize the QSGW self-energy in NiO.This manifests as shifts in QSG Ŵ energy bands and peaks in the density-of-states (DOS).DOS are compared to bremsstrahlung-isochromat-spectroscopy (BIS) and x-ray photo-emission (XPS) measurements in the left panel [78].
• BIS data exhibits three peaks between 0 and 9 eV, which the QSG Ŵ DOS captures quite well, except for a small underestimate of the fundamental gap seen in both BIS and optics (see §III F 1).This shows that ladders do an excellent job of capturing the frequency dependence of the local (k-integrated) spectral function.
• The corresponding QSGW peaks are blue shifted relative to experiment, but in varying amounts.Peak 1, which is composed almost entirely of flat Ni d states, is shifted about 1.5 eV while Peak 2, derived essentially of dispersive Ni sp states, is shifted by about half of that.This reflects a universal tendency: flat bands are affected by ladders more than dispersive ones.Fe 3 O 4 offers a particularly striking example ( §III F 7).
• QSG Ŵ significantly narrows the occupied Ni d bands relative to QSGW .Red bands (depicting O character) are almost unaffected, while there is a significant narrowing of the green bands relative to QSGW .This is a potentially important finding.
It is well known that the LDA severely overestimates d band widths in narrow-band transition metal compounds.Further, it has been shown in several works, e.g.Ref. [79], that QSGW narrows d bands relative to the LDA, but nevertheless continue to overestimate these bandwidths, especially in systems with strong spin fluctuations such as BaFe 2 As 2 and FeSe.In cases we have studied where experimental information is also available, e.g. in Sr 2 RuO4 4 [57], this overestimate is remedied very well by augmenting QSGW with DMFT, which includes vertices in both charge and spin channels.Whether the bandwidth can be captured entirely by a combination of low-order diagrams in both spin and charge channels remains an intriguing possibility.
To the extent it is true, this greatly simplifies the complexity of the electronic structure problem in correlated systems.This will be explored in a future work.
Some more details for NiO are presented in §III F 1. Also, there are some strong parallels with La 2 CuO 4 ; see §III F 6.

Consistency between one and two particle properties
The consistency between benchmarks for one-and two-particle quantities (E G and ε ∞ in Fig. 4) is striking.Apart from some outliers to be discussed in §III D, the calculated values ε ∞ agree with measured ones to within the available resolution.When this is not the case, usually there is a corresponding discrepancy in the fundamental gap: discrepancies in E G and ε ∞ occur largely at the same time: overestimate of E G yields underestimate of ε ∞ , and underestimate of E G yields overestimate of ε ∞ .
The internal consistency between one-and two-particle properties is a signature of consistency of the theory, since the same quantities (G and W) construct both ε(ω) and the potential Σ(ω) that makes G.
If we assume the fidelity of the theory is sufficient for this principle to be universally applicable, the extra information provides an ansatz to predict optical properties in materials with stronger correlations, where benchmarking is less simple.In such cases there is often a large uncertainty in the benchmark itself, not only owing to a wide variation in reported experimental data, but also the extraction from one-particle properties (e.g.fundamental gap) from two-particle response functions.This is reasonable for tetrahedral semiconductors where excitonic effects are small (see Fig. 9), but has less validity in general.For these more correlated cases our approach will be to compare optical experiments directly with calculated response functions.Combining such a comparison with the observed relation between calculated one-particle and two-particle properties, we can benchmark the theory, and sometimes provide values of quasiparticle levels where not well known, or new interpretation of accepted values.RPA@QSGW BSE RPA@LDA FIG. 6. ε∞ generated by two inconsistent approximations: ε RPA ∞ @QSG Ŵ (blue circles) and ε BSE ∞ @QSGW (green diamonds), to be compared against Fig. 4. Also shown are LDA results: ε RPA ∞ @LDA (red squares).Agreement is best for indirect gap semiconductors, where gaps for vertical transitions are relatively large.Fig. 4 was generated from two consistent approximations: ε RPA ∞ @QSGW , and ε BSE ∞ @QSG Ŵ .Consider by contrast two inconsistent approximations, ε BSE ∞ @QSGW and ε RPA ∞ @QSG Ŵ (Fig. 6).Both of these approximations show more randomness than either ε RPA ∞ @QSGW , or ε BSE ∞ @QSG Ŵ of Fig. 4. Yet there is a strong similarity between the green diamonds and the blue circles in the two figures.The green diamonds in Fig. 6 fall slightly below the ideal line, showing a modest but nonnegligible effect of improving the reference G 0 .The blue circles rise slightly above the 80% showing that the RPA continues to underestimate ε ∞ , even with a nearly ideal reference G 0 .This affirms that most-but not all-of the underestimate of ε ∞ originates from the RPA itself.This sheds light on the commonly observed fact that ε RPA ∞ , when computed from the G LDA , often provides a rather good estimate for ε expt ∞ , e.g. in sp semiconductors.The obvious, naive reason for this is a fortuitous error cancellation: LDA underestimates bandgaps, which tends to overestimate ε ∞ , while the RPA's neglect of electron-hole attraction tends to underestimate the screening, and thus tends to underestimate ε ∞ .However there has been some speculation that the good agreement is not accidental, but a consequence of characteristics inherent in the RPA and the LDA.In particular a recent work [80] asserts that ε RPA @LDA should be a good approximation for insulators, based on two arguments.First, ladders involve tunneling processes, and are effective at short range but not long range; thus the long-range screening that predominately controls ε(q=0) is well described by the RPA [81,82].The first argument is rather appealing, and consistent with prior work establishing that the largest corrections to the RPA occur at short distances [81][82][83].This argument can be rigorously checked by comparing blue circles in Fig. 6 to green diamonds in Fig. 4: both share the same eigenfunctions generating ε; the only difference being the presence or absence of ladder diagrams.Agreement is fairly good, differing by 15-20%, which explains in part why QSGW is a good and consistent theory.The argument of van Loon et al. [80] is only partially true: even if the vertex part of P 0 is short range, v is long range so ε=1 − vP 0 can have a long range contribution from the short ranged part of P 0 .Notably, the difference does not get smaller as the gap becomes large, as Ref. [80] asserted based on the tunneling argument.This is apparently because as the gap closes the screening becomes large, so the long range contribution from a short range vertex becomes relatively less important.It is nevertheless striking that the BSE correction is so insensitive to the bandgap.
The second argument of Ref. [80] is that the local vertex in insulators is approximately accounted for by using LDA eigenfunctions.The argument is based on a connection between the LDA derivative discontinuity and the missing vertex, which emerges in a model.To examine this proposition, ε RPA ∞ @LDA is also presented in Fig. 6 (red hexagons).The second argument is more difficult to assess quantitatively because screening modifies both G 0 and ε ∞ , but roughly speaking the difference between blue circles and green diamonds should be similar to the difference between red hexagons and blue circles.One might attribute the poor agreement to inadequacy of the LDA functional (distinct from the derivative discontinuity in the exact functional), but at least in a few systems where it has been tested, the primary gap error has been shown to originate largely from the derivative discontinuity, and not inadequacy of the functional [84].There is a distinct tendency for LDA to better predict ε ∞ for indirect gap tetrahedral semiconductors than for direct-gap ones: compare diamond, AlAs, AlSb, GaSb in Fig. 6 to ZnO, ZnS, CdTe, InP and GaAs.Since the derivative discontinuity does not vary wildly between direct-and indirect-gap materials, this is a hint that some other parameter controls the errors in ε RPA ∞ @LDA.Note also that in the one-oscillator model [73], §II H, the effective oscillator energy E 0 tends to better align with the smallest direct gap than the fundamental one.To disentangle the various effects, Fig. 7 plots the relative error in ε RPA ∞ @LDA against the relative error in the gap, which is a proxy for the derivative discontinuity.Excepting the d 0 and f 0 systems, the relationship between E LDA G /E exp G and ε RPA ∞ @LDA/ε exp ∞ is roughly linear.The sensitivity of ε RPA ∞ @LDA to the derivative discontinuity, together with the tendency of RPA to underestimate bandgaps established earlier, provides strong indication that ε RPA ∞ @LDA yields reasonable ε expt ∞ only sometimes, and since it generally produces values larger than experiment while the RPA underestimates the screening as we have shown, there is an additional hidden benefit from fortuitous error cancellation.

C. Benchmarks in weakly correlated semiconductors
The tetrahedrally coordinated sp 3 compounds form a good benchmark for weakly correlated systems in part because they are the best characterized of any family of materials, but also because weak correlations make it possible to well identify transitions between single-particle levels, especially associating peaks in ellipsometry measurements with them.The valence band maximum falls at or very near Γ for all tetrahedrally coordinated semiconductors, which simplifies the analysis.Besides the lowest Γ−Γ transition E 0 , the next Γ−Γ transition E ′ 0 has been measured for some materials.Ellipsometry also measures E 1 and E 2 shown in Fig. 8. E 0 and E 1 are easier to measure accurately, because there is larger volume in k space where the valence and conduction bands are parallel.E 2 has been measured for most semiconductors, but its determination is less certain (excepting compounds such as Si, C, and SiC where the global conduction band minimum lies near X).Some data for E ′ 0 are available, but their values are also less well known.
Energy bands in GaAs, depicting vertical transitions E0, E ′ 0 , E1 and E2 that can be measured by ellipsometry.Circles depict measurements of states at high-symmetry points.They have been determined by ARPES [85] to a resolution of about 0.1 eV.E0, E ′ 0 , E1 and E2 are reported by Lautenschlager et al, in Ref. [86].Combining this data provides one way to determine levels at X and L in the conduction band.Colored bands are taken from QSG Ŵ calculations, with red and green showing projections onto Ga and As, respectively.Gray lines show results of QSGW calculations.In the valence, QSG Ŵ and QSGW are nearly indistinguishable.QSG Ŵ and QSGW dispersions in the conduction band are very similar, with QSGW slightly higher in energy.
The wider conclusions we draw from the detailed analysis to be described below are as follows.
1. Bandgaps in light (and especially polar) materials are overestimated (MgO, LiF, LiCl, NaCl, TiO 2 , SrTiO 3 , C).The primary cause is the electron-phonon interaction ( §III D).A diagrammatic electron-phonon contribution to Σ has long been known [3] though historically speaking, reliably determining its magnitude has posed a challenge.A fairly high fidelity calculation of it has recently appeared (Ref.[41]), and we use their results to estimate this term where available.In other cases we make a simple estimate using the Frölich approach of Ref 72 ( §II G).The fundamental gaps with these adjustments are shown as black crosses in Fig. 4. See also §III D 2.
2. The gap in compounds with shallow, nearly dispersionless d levels are too small (Table IV), and semicore d levels are too shallow (Fig 11).This is a consequence of the imperfect Z factor cancellation noted in §II D 1, Point 2. To correct it would require the missing vertex Γ in the exact self-energy, GW Γ.Several instances of this are presented in §III D 2.
3. k-dispersions in the conduction bands of zincblende semiconductors show systematic errors of the order ±0.

eV (See discussion around Tables II, III
).There is no obvious diagram that explains this discrepancy.
Fig. 9 benchmarks E 0 , E ′ 0 , E 1 and E 2 transitions in zincblende semiconductors where ellipsometry data is available.E 1 shows close agreement, but E 0 and E 2 exhibit discrepancies with distinct patterns: Tables II and III establish that there is a systematic, k-dependent error in the conduction band in zincblende semiconductors on the order of 100 meV.The consequences can be significant: note for example, that QSG Ŵ predicts GaSb to have a global minimum at L, with E Γ −E L =0.05 eV, while experimentally at 0K it is a direct gap, with E Γ −E L = − 0.09 eV [90].Also, where gaps are overestimated, effective masses are too large (Fig. 10).
The k-dependent gap error is further discussed in §III C 1, but we can find no obvious explanation for it.One possibility is that Questaal's implementation of QSGW contains an error not inherent in QSGW itelf.In particular the incomplete basis noted by Betzinger et   Data for III-V semiconductors are taken from Ref. 87; other data taken from Adachi's compilation [88].CdTe, ZnTe, and GaAs are taken from a two-photon magnetoabsorption experiment [89] which is thought to be reliable.AlSb data is for the conduction band at Γ. so many systems is a fortuitous artifact of the implementation, or fortuitous cancellation of higher order diagrams.At all events there is no simple explanation that reconciles these inconsistencies.Fig. 10 benchmarks for effective masses and bandgaps in tetrahedral semiconductors.For the direct gap systems (circles and hexagons), the discrepancy in m * compared to the experimental value scales approximately in proportion to the discrepancy in E G (compare to the light dashed gray line).k•p theory predicts a m * to be proportional to E G , assuming a fixed matrix element coupling valence and conduction band, showing that errors in m * have the same origin as whatever causes the gap to be too large.
Core levels Fig. 11 presents d core levels at different levels of theory and compares to photoemission results.As is well known, their position is too shallow in the LDA.(GW ) LDA improves agreement with experiment, but levels remain too shallow.Selfconsistency (QSGW ) shows further improvements.QSG Ŵ fares slightly worse than QSGW on average, but the consistency improves with the level of theory.
Table IV shows some materials system where the valence band maximum is nearly flat and dispersionless.The bandgap is consistently underestimated in QSG Ŵ .(QSGW fares better, but this is an artifact of fortuitous error cancellation).
Shortcomings shown in Fig. 11 and Table IV have a common origin, the missing vertex ( §II D 1, point 2).Note that the discrepancy with experiment increases with distance from the Fermi level: from ∼0.5 eV for the valence states near E F , ∼0.9 eV TABLE II.Bandgaps at X or L for some zincblende semiconductors.All the semiconductors listed above have a global conduction band minimum near X, except for Ge and GaSb.When the electron-phonon interaction is taken into account, gaps at X are systematically underestimated by ∼0. 15  for cation levels around −10 eV, and ∼1.3 eV for the deeper anion levels.

Valence band parameters
The structure of the valence band around Γ provide less reliable benchmarks because of experimental uncertainty in the parameters.Key parameters are the effective masses, and in the wurtzite structure, the crystal field splitting arising from inequivalence of the z and xy directions.As regards the masses, the matter is considerably complicated by intermixing three states at the valence band maximum near the Γ point and the nontrivial role of spin-orbit coupling that splits the threefold degeneracy at Γ and pushes the band maximum slightly off it.
To encapsulate the many different masses, a Luttinger model is typically used, which has only three independent parameters.The Luttinger parameters can be generated from the following effective masses: m lh and m hh denote light-hole and heavy-hole masses.Table V shows Luttinger parameters for a few systems where they are best known.The range of values shown in the experimental columns correspond to the range collated from different measurements.In the two cases where the bandgap is close to experiment (Si and InP) the calculated Luttinger parameters fall within the range of experimental data.In the other two cases (Ge and GaAs) the parameters are underestimated for the same reason the conduction band effective masses are overestimated (see Fig. 10): the direct gap is somewhat overestimated; (see e.g.Table III) In Ge, for example, the conduction band mass at Γ was measured to be 0.037 [99], while our QSG Ŵ mass is 0.047.
Table V also show shows crystal-field splitting ∆ cr in the III-N semiconductors (splitting between states of p z and p xy character at Γ) in the absence of spin-orbit coupling.The QSG Ŵ result is within ∼0.01 eV of the measured values, which is quite satisfactory.This quantity is rather sensitive to find details of the potential.To obtain ∆ cr reliably, a fine k mesh of 9×9×6 divisions was needed: its value increased by 0.005 eV compared to the standard 6×6×4 mesh.Note that OEP-based GW reported in Ref. [100] yields quite different values for ∆ cr .[101] b compilation in Ref. [87] c Experimental data and G 0 W 0 results taken from Ref. [100].

Two extensions to the theory and their effect on zincblende semiconductors
As two possible sources of error on the QP levels of zincblende semiconductors, we first considered eliminating the Tamm-Dancoff approximation (TDA).Here we focus on InSb as it has the largest relative gap error.Removing the TDA reduced the QSG Ŵ E 0 gap by 0.03 eV -considerably less than the discrepancy with experiment.We also considered whether eliminating the TDA improves the k-dispersion, in particular the wrong prediction of the global minimum in GaSb noted above.Removing the TDA reduces the gap in GaSb by 0.03 eV (similar to InSb), but the shift was essentially independent of k and did not rectify this shortcoming.
We also considered the effect of using a better kernel in the BSE.In all the calculations presented here, we used W RPA for the kernel (Eq.23).It is possible that the dispersion errors in these compounds is a consequence of W RPA being too removed from the exact vertex.We can assess this effect by using a better kernel, namely BSE W as the kernel for the BSE.If we assume naively that the main effect of BSE is to reduce W (q = 0, ω = 0) (i.e. the change in ǫ ∞ ) then the substitution W RPA →W BSE in Eq. 23 would reduce the strength of the electron-hole attraction and shift the electronic structure to (e.g. the bandgap) something intermediate between QSGW and QSG Ŵ .This is roughly what happens in some cases, e.g.CrX 3 [102].In sp semiconductors, however, using a better W in the vertex causes the gap to decrease still further by a small amount, e.g. by 0.03 eV in InSb.This is another manifestation of vertex corrections being short ranged, as noted earlier.
To conclude, the combined effect of eliminating the TDA and better W in the vertex, are not sufficient to explain the tendency to overestimate the direct gap in small-gap zincblende semiconductors, or errors in the band dispersion.

Birefringence
Birefringence occurs when the refractive index depends on the polarization and propagation direction of light.It is normally measured as a difference in the principal axes of the ellipsoid's index of refraction, sometimes called the "ordinary" and "extraordinary" indices when there are two inequivalent ones.We consider a few materials where n in the basal plane differs from n normal to it: Birefringence is measured as the difference ∆n.

Relation between gap and dielectric function
Fig. 4(b) appears to predict ε ∞ very well, but there are discrepancies.Here we focus on systems for which ε BSE ∞ falls outside the uncertainty of experimental values (estimated by the variation in reported values), and show that these errors directly correlate with errors in the fundamental gap.
Several known potential sources of error in Σ were enumerated in §II D. Among them, the electron phonon interaction is significant for wide-gap, light-element compounds, especially polar ones where the narrow valence band enhances the Frölich interaction, Eq. 24.The electron-phonon interaction usually reduces gaps, by as much as 0.5 eV in an extreme case such as MgO.Table VII selects some materials where this reduction exceeds 0.3 eV.In such cases ε ∞ is slightly underestimated.As we noted previously, at present Questaal does not have the capability to incorporate the electron-phonon self-energy into the QSGW cycle; however, we can make a proxy by making a hybrid of the LDA and QSG Ŵ potentials to reduce the gap (Eq.21).We choose β=0 and pick the mixing parameter α to approximate the gap change from electron-phonon interaction calculated in Ref. [41].This should be a reasonable proxy for Σ e−ph since for these systems the LDA and QSG Ŵ bands differ mostly in a simple rigid shift of the conduction band.Materials in Table VII above the dividing line show systems for which the electron-phonon interaction exceeds 0.3 eV, and where both E G and ε ∞ are thought to be reliably known.Renormalization causes a modest increase in ε ∞ , and the systematic tendency to underestimate it is reduced to approximately the experimental uncertainty. 9or compounds in the bottom half of the table, benchmarking becomes murkier, because the electronic structure is not known or is poorly understood.hBN might have been put in the top half of the table, if so it would present a severe anomaly.Data for hBN in Tables VI and VII were computed from an average of Refs.[107] and [108].To suggest the possible source of the anomaly, Table VII also shows an entry where experimental data is taken only from Ref. [107], and by such a comparison the agreement is in line with other materials.Further experiments are needed to determine the true values (both ordinary and extraordinary) for ε ∞ in hBN.
For less well characterized systems, if we make the ansatz that the calculated ε ∞ should coincide with the experimental one when E G also coincides, we can assess the effect of the error in the fundamental gap if ε ∞ is better known (this is a common situation).We can estimate what E G should be by matching ε ∞ (more generally ε(ω)) to experiment.In later sections we apply this technique to several materials systems, e.g.CuAlO 2 , §III E 7 and FeO, §III F 4.

E. Band Structure and dielectric function in selected nonmagnetic materials
In this section we present a variety of selected materials.Where sufficient experimental information is available (e.g.LiF) those results are used to benchmark the theory.For most of the systems presented here, the available experimental information is partial, confused, or contradictory.For these systems we use a mix of theory and what experimental information seems sufficiently reliable, to arrive at a consistent picture where it seems reasonable to do so.In a few cases it is not fully possible (see CuAlO 2 , §III E 7).
The analysis relies on the ansatz stated in §III D 2, namely that if G 0 is good enough to well characterize one-particle properties, it also well characterizes two-particle properties provided an adequate theory for the vertex is used; and moreover, that ladder diagrams are sufficient for the vertex.This hypothesis was affirmed in nearly every case in the present study where reliable information is available.

LiF
The macroscopic dielectric function of the polar insulator LiF was recently calculated in Ref. 30 within the BSE using QSGW as the starting G 0 .Since the vertex corrections are omitted in QSGW the screening of the exchange was underscreened and the gap too large.Combined with the neglect of the electron-phonon self-energy, this results in a greatly overestimated band gap of ∼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 QSGW electronic structure results in a blue shift of ∼ 0.9 eV with respect to experiment.
Here, we repeat the BSE calculation of the optical spectrum but on top of the QSG Ŵ electronic structure and also consider the lattice polarization effects.Including ladder-diagram vertex-corrections produces a fundamental band gap of 14.7 eV; a reduction of over 1.4 eV, in agreement with the vertex correction calculated in Ref. 18. Including the 0.48 eV polaron shift correction from Ref. 71 gives then a band gap in excellent agreement with the experimental value of 14.2±0.2eV [132].The exciton binding energy is around 2 eV, also in agreement with the experimental value [132] and as a consequence Fig. 12 shows an excellent overall agreement between the theoretical and experimental spectra.The BSE ǫ ∞ (1.95) is close to the experimental one (1.96,Ref. [133] and 1.92, Ref. [134]).
TABLE VII.Estimated change to ε∞ induced by adjusting the QSG Ŵ self-energy according to Eq. ( 21), using α given in the table.∆EG<0 shows the change in fundamental bandgap, in eV: ∆EG<0 indicates the gap is reduced by taking β=0 and γ=1−α.∆EG>0 indicates the gap is increased by taking γ=0 and β=1−α.ε and ∆ε are defined in Eq. 28.Column "org" indicates the probable predominant physical origin of ∆EG: one of eph (electron-phonon); Γ (missing vertex in Σ), or * (unknown).Top box displays systems where both EG and ǫ∞ are fairly reliably known, or reliably known.Bottom box cotains entries where EG, and to some extent ǫ∞, are not well known.For CoO and MnO no adjustment was made owing to uncertainty in ǫ∞, and lack of information about the effect of the electron-phonon interaction.

Bi2Te3
Bi 2 Te 3 is widely studied because it has topological surface states protected by time-reversal symmetry [135].It is a narrow-gap system with reported energy gaps between 130 and 170 meV [136][137][138][139].The QSGW and QSG Ŵ bands are shown in Fig. 13, and are seen to be nearly identical.This system was studied previously [40] within the LDA and G LDA W LDA approximations (albeit including the off-diagonal parts of Σ).Fig. 13 is in close agreement with Fig. 1 of Ref 40.That the three many-body calculations (G LDA W LDA of Ref 40, QSGW and QSG Ŵ ) are so similar suggests that W is already well described by the LDA.Evidently ladder diagrams have almost no effect.This is perhaps not suprising, since the LDA and GW bands are also similar, with the LDA gap slightly smaller at 50 meV [40].QSGW and QSG Ŵ energy gaps are both 145 meV, slightly larger than 120 meV reported in Ref. [40] (presumably because of self-consistency), and within the range of reported experiments [136][137][138][139].

ScN
ScN is a material of considerable interest in opto-electronics applications, especially as a buffer layer.It is an indirect gap material with the conduction band minimum at X. Its bandgap has been controversial with many reported values ranging from 2.03-3.2eV for the direct gap and 0.9-1.5 eV for the indirect one.Theoretical predictions similarly vary, with predictions ranging between 1.82-2.59eV (direct) and 0.79-1.70(indirect) (see Ref. [140] for a summary and detailed discussion).
The most recent and detailed experimental study taking into account prior work (Ref.[140]) yields an optical indirect gap of 0.92±0.05eV (Table in Fig. 14).QSG Ŵ predicts a larger fundamental gap, 1.27 eV.The latter should be reduced by the electron-phonon interaction; Unfortunately no information is available in the literature, but it is likely to be similar to InN.The longitudinal and transverse mode phonon frequencies are similar, while the electron effective masses in ScN are heavier (0.34 eV for ScN, 0.07 for InN).Thus according to the Frölich formula, Eq. 24, the electron-phonon renormalization should be larger, by a factor between 1 and 2. A reasonable estimate is 0.1 eV, which is used in the Table in Fig. 14.A study of ǫ BSE (ω) yields an exciton at 2.19 eV; thus there is a spread of 0.14 eV between fundamental and optical gaps.This is apparent in the dielectric function (right panel of Fig. 14)).The shoulder between 2 and 2.3 eV are the subgap excitonic transitions.The dotted line is a guide to the eye, extrapolating the onset of the second shoulder to zero.
Combining this shift with an (admittedly crude) estimate of 0.1 eV for the electron-phonon interaction, brings the QSG Ŵ and ellipsometry direct optical gaps to within 0.1 eV (approximately the uncertainty in both theory and experiment).ε ∞ also agrees well with Ref. [140] (Table in Fig 14).As ε BSE (ω) does not include indirect couplings via the electron-phonon interaction, and we did not consider contribution from q =0 transitions, and cannot determine the exciton binding for the indirect transition, but it is reasonable to assume it is similar to the direct one.To summarize, a consistent picture emerges in close agreement with the recent work of Ref. [140], with the proviso that the one-particle and two-particle gaps must be distinguished.

CeO2
Electric conduction in CeO 2 takes place both by ionic and electronic conduction, and can be controlled by changing the O 2 pressure.Its unusual electrical properties make it a promising candidate anode in solid oxide fuel cells, or for intermediate temperature electrolytes [141].Its energy bandgap has been measured optically by a number of groups, by absorption [142,143] or reflectance [141,144,145].They vary in details, but all report optical bandgaps ranging between 3.1 and 3.3 eV.We compare against absorption data in Ref. [142] because the reported energy range was wide enough to show two peaks (Fig. 15).[142].Red curve is the BSE absorption α generated from the QSG Ŵ hamiltonian (α=4πn/λ, n 2 =dielectric function).The energy axis for the calculated function is redshifted by 0.5 eV as described in the text.
The QSG Ŵ bandgap is computed to be 4.24 eV (4.93 eV in QSGW ).The valence band is almost pure O-p character, the lowest conduction bands are nearly dispersionless Ce-4f bands (see left panel, Fig. 15).Above the narrow Ce-4f bands are states of mixed Ce-5d and Ce-6s character.LDA bands are not shown, but there is an orbital-selective shift: Ce-4f bands shift about 2.4 eV relative to the LDA, while the Ce-5d bands shift only 0.8 eV.(There is a smaller but orbital-selective shift of the opposite sign as ladders are added to QSGW , as happens for NiO, Fig. 5.) The LDA gap (1.8 eV) is not so far removed from the optical gap (∼3.2 eV), but this is largely fortuitous, for several reasons.

• The LDA badly understimates the shift in empty Ce-4f states
• There is a large renormalization of the Hartree part of the hamiltonian, which reduces the gap substantially (Table I), and this partially cancels the first error.
• The optical gap and fundamental gap apparently differ by ∼0.6 eV.Note the absorption spectra in the right panel of Fig. 15.In that figure, the BSE optical spectra were red-shifted by 0.5 eV to align them with the absorption data.Thus the QSG Ŵ optical gap is ∼3.6 eV, about 0.6 eV less than the fundamental gap.QSG Ŵ still overestimates the optical gap by ∼0.5 eV (Fig. 15); however, some portion of this difference can be attributed to the electron-phonon interaction, as explained below.
The global valence band maximum falls on the Σ line, about 2/3 between Γ and K (see Fig. 15), though other local maxima are nearly degenerate with it.
Comparing ǫ ∞ to experiment is not straightforward because of the wide dispersion in reported experimental data, as well as preparation conditions [146] and the crystallinity of the material.Reported values vary from 4.7 [147] to a range between 5.8 and 6.6 [148] on single crystals.A measurement of highly oriented crystalline films yields ǫ ∞ =6.1 [146].ǫ ∞ from QSG Ŵ +BSE is found to be 5.8, which fits comfortably within the range of reported experimental values.
An estimate of the zero-point motion can be made using Eq. ( 24).For this equation, a P is needed separately for the conduction band and valence band; however, the effective mass approximation is not meaningful for the almost flat conduction band, and we consider only the valence band here.Various experiments put ω LO =30±5 meV [149]; we compute the QSG Ŵ hole masses to be (0.86, 1.3, 1.8)m 0 , for an average mass of 1.27.Using these values, we obtain a P =18.9a 0 .The static dielectric constant is roughly ∼25 [150].Using ǫ ∞ =6, Eq. ( 24) predicts the valence band contribution to gap reduction to be 46 meV.Effective mass theory cannot be applied to the nearly dispersionless conduction band, but it is reasonable to expect its contribution to the total gap reduction to be several times larger.A factor of three larger conduction band contribution would give a total gap correction of ∼0.2 eV.This accounts in part, but it would seem not entirely, for the apparent ∼0.3-0.5 eV overestimate of the optical gap predicted by QSG Ŵ .It would seem a discrepancy of order 0.3 eV remains, but a better determination of the electron-phonon interaction is needed to know the discrepancy reliably.Assuming that for whatever reason, the QSG Ŵ gap is too large by 0.5 eV, the fundamental gap should be about 3.75 eV.
The RPA dielectric function calculated from the LDA [149] also shows a peak around 3 eV, but this is an artifact of error cancellation: RPA omits strong excitonic effects and the LDA gap is too small.TiO 2 , also known as titania, is a widely used commercial compound, as food coloring (E number E171) and as a pigment in paint for example.We consider here only the rutile phase.TiO 2 is a d 0 system with valence band essentially O 2p character and conduction Ti d character (Fig. 16).The optical absorption edge was measured to be ∼3.3 eV by Cardona and Harbeke from reflection measurements [105], and as 3.03 eV from optical transmission [151].This is well below the calculated QSG Ŵ fundamental gap of 3.88 eV.Part of the difference may be attributed to the electron-phonon interaction, which was calculated in Ref. 41 to reduce the gap by 0.34 eV.Reducing the gap by this amount puts it in line with a PES/BIS study, which reported a fundamental gap of 3.3±0.5 eV [152].
The first peak in Im ǫ(ω) is blue shifted about 0.2 eV compared to experiment [105], and correspondingly ǫ ∞ is about 5% smaller than experiment (Table VII).The corresponding QSG Ŵ birefringence is also slightly underestimated (Table VI): ∆n BSE = 0.22 compared to ∆n expt = 0.26.The corresponding RPA values from QSGW are about 80% of experiment, as is typical, with ∆n RPA = 0.16.Excitonic effects are strong in TiO 2 : compare Im RPA ǫ to Im BSE ǫ in Fig. 16 (both were generated from the same hamiltonian).The BSE redshifts the peak in Imǫ and significantly changes the shape.
Both Im BSE || ǫ(ω) and Im BSE ⊥ ǫ(ω) show reasonable resemblance to the experiment below 6 eV: peak at 4 eV, shoulder at 5 eV in both Im ε || and Im ε ⊥ .Three bound, weakly active excitons are found in the region (E c −0.45, E c −0.17 eV) below the fundamental gap, and several bright ones between E c −0.13 and E c .Here E c is the conduction band minimum.Thus Im BSE ǫ shows a strong peak close to the fundamental gap, with a tail extending below.
The electron-phonon interaction missing in QSG Ŵ well accounts for the blue shift in leading shoulders in Im BSE || ǫ(ω) and Im BSE ⊥ ǫ(ω) relative to experiment.As a proxy to account for it, we repeat the calculation with a hybrid G 0 , consisting of 90% QSG Ŵ and 10% LDA.This reduces the gap by the amount calculated in Ref. 41.With this shift, ǫ ∞ and the birefringence both align closely to available experiments (Table VII).
Thus, discrepancies in QSG Ŵ fundamental gap and experimental optical gap are fully explained in terms of a combination of excitonic effects, and the electron-phonon interaction.Adding the latter to the QSG Ŵ fundamental gap, we conclude that its true value is 3.5±0.1 eV, significantly larger than the widely accepted value of ∼3 eV.SrTiO 3 is a perovskite material that may exist in the usual different perovskite phases: cubic, tetragonal and orthorhombic.Like TiO 2 , it is a d 0 compound with valence band essentially O 2p character and conduction Ti d character (Fig. 17).In this respect it is very similar to TiO 2 , and care must be taken in interpreting the experiments to determine the fundamental bandgap.The experimental band gap is reported to be 3.25 eV (indirect) and 3.75 eV (direct) [153] and according to Bhandari et al. [72] it is almost independent of the structure.According to QSG Ŵ , the system has a indirect gap, of 4.06 eV, with the valence band maximum at R and conduction band minimum at Γ.The direct gap at Γ is 4.51 eV, larger than the indirect one by 0.45 eV; however, the QSG Ŵ fundamental gap and optical direct gap differ by about 0.75 eV (Fig. 17).Note the BSE code has no electron-phonon coupling and cannot detect indirect transitions.
The shoulder in QSG Ŵ dielectric function in Fig. 17 lies about 0.1 eV above the ellipsomentry measurement of Ref. [153].This suggests that the band structure is close to the true one.ǫ ∞ is slightly lower than the average experimental value (Table VIII), which suggests that the uncorrected QSG Ŵ gap is slightly too large.Using a hybrid functional to reduce the gap, ǫ ∞ moves close to the average experimental value for ǫ ∞ (Table VIII), but the shoulder in ε(ω) is slightly redshifted relative to the Benthem data.Thus there is a slight inconsistency.This excludes a precise determination of the fundamental gap, but we conclude it is 3.75±0.1 eV, which is about 0.5 eV larger than the reported optical gap.

CuAlO2
CuAlO 2 has received a great deal of attention because it is a p-type transparent conducting oxide (TCO); indeed, it seems to be the only known TCO that can be doped p-type.
,1/2,0), and (0,0,-1/2), respectively, as multiples of the reciprocal lattice vectors.(Right) dielectric function (average of x and z directions) as function of frequency ω (eV).Solid line shows Im ǫ BSE (ω), generated from QSG Ŵ .The green dashed line shows Im ǫ RPA also generated from the QSG Ŵ hamiltonian.Note Im ǫ BSE and Im ǫ RPA approach 0 near the fundamental direct gap at 4 eV; however, Im ǫ BSE has additional subgap peaks from strongly bound excitons.The blue dotted line shows Im ǫ RPA (ω) generated from QSGW .The RPA functions look similar except for a 1 eV shift, because of the difference in bandgaps.
CuAlO 2 was not included in Fig. 4 because reports of its magnitude vary widely.Reports of the lowest (indirect) gap ranges between 1.65 [155,156] and 2.99 eV [157,158], and a direct gap ranging between 3.3 and 4.2 eV [157,159,160].QSG Ŵ finds that the valence band is almost exclusively Cu 3d character, with its maximum at a low-symmetry point near S in Fig. 18.The conduction band minumum at Γ is mixed Cu-sd and O-p character.The fundamental (indirect) gap is found to be 3.5 eV, and the lowest direct gap at L is 4.0 eV.There is a large region of k space where the highest valence band is nearly dispersionless.A prior calculation, using a hybrid functional, found similar gaps but not the nearly dispersionless highest valence band [161].
The principal axes for hole mass are along low-symmetry directions, with a moderate mass (1.3) on one principal axis and large masses (2.7, ∼10) on the other two.The conduction band at Γ, by contrast, has much smaller masses and they are along the Cartesian axes (1.3 in xy, 0.41 along z).It has been argued that there should be a large gap renormalization from the electronphonon interaction [162]; however the measured difference ǫ ∞ (5.1) and ǫ 0 (7.7) [122] is fairly small, and the eigenstates at the band edges are mostly Cu-like instead of O (as it is in MgO), both of which reduce Σ e−ph in the Frölich model ( §II G). 10The fundamental gap is in line with the photoemission study of Ref. [155], which measures the DOS, a one-particle property.However, when two-particle properties are considered, QSG Ŵ predicts the situation to be more complicated.The right panel of Fig. 18 compares the BSE and RPA dielectric functions.Both approach 0 at the fundamental direct gap.However, the BSE shows strong peaks below the fundamental gap, around 3.2 eV; there are also several excitons between 3.7 eV and the fundamental gap at 4 eV.Such deep excitons are not typical in sp semiconductors, but it can be understood as an artifact of the nearly dispersionless Cu-like valence band, as well as a relatively small dielectric constant of 5.1.QSG Ŵ predicts ǫ ∞ rather well.If the strong correlation between the reliability of ǫ ∞ and the bandgap §III B 3 applies equally to CuAlO 2 , the gap should be close to the QSG Ŵ prediction of 3.5 eV.ǫ(ω) was computed without an electron-phonon contribution, so it can only measure direct transitions.Presumably there will be other excitons for bound electron-hole pairs coupling Γ and states in the valence band as well; thus the optical response will show some intensity in a spread below the peak at 3.2 eV, larger than what is shown in Fig. 18, possibly as much as the difference between the direct and indirect gap.Since most determinations of the gap are performed with optical measurements, much of the confusion in the literature likely originates from these deep excitons.These excitons cannot explain a gap as low as 1.8 eV, however; such a gap likely originates from a defect band, which explains why it is not always seen.Indeed, recent work [163] shows that the optical absorption edge is strongly dependent on preparation and post-annealing conditions.Defects apparently play an important role, which adds to the confusion about experimental reports on this materials system.

VO2
In the low-temperature monoclinic (M 1 ) phase, VO 2 has a gap approximately 0.7 eV [164,165].M 1 is a deformation of the high-symmetry rutile phase.The unit cell, consisting of four V atoms all equal in the rutile phase, dimerize into two pairs with short bond lengths.It is generally agreed that the V dimerization is what is responsible for the gap, splitting the V d manifold into a single occupied d bond per dimer, and a corresponding antibond (Peierls transition).Wentzcovitch et al. [166] calculated the energy band structure and suggested that despite the LDA yielding no gap, that the origin of the gap was more Peierls like than Mott-like.This picture is further supported by the observation that LDA augmented by single-site DMFT is also metallic [167], which would not happen in a simple Mott description.A cluster form of DMFT added to LDA does yield a gap [168].This indicates that the nonlocality of the self-energy is essential, and explains why it is too small in the LDA.Gatti et al. [169] employed GW to study this system, which captures the nonlocality quite well.While they found G LDA W LDA failed to open a gap, a self-consistent GW scheme within the COHSEX approximation did so.Counterbalancing this view, a DMFT work [170] argued the M 1 phase should be characterized as the Mott transition in the presence of strong intersite exchange.In our view Gatti's work is the most definitive, as it does not rely on the LDA, partitioning, or adjustable parameters.It also confirms the original Wentzcovitch conjecture: VO 2 is a simple band insulator.In further support of this conjecture, magnetism appears to play no role in this phase, as we show next.Here we computed the electronic structure of VO 2 in the QSGW and QSG Ŵ approximations.VO 2 has a nearly dispersionless core-like valence band maximum (Fig. 19).This makes it a prime candidate for the gap to be underestimated, owing to the missing vertex (see discussion around Table IV).Indeed QSGW predicts a rather good gap (E G =0.76 eV) owing to a fortuitous error cancellation, while QSG Ŵ underestimates the gap (0.43 eV), reminiscent of CuCl.
If Mott physics were involved, magnetism should play a role.We find within QSGW , that magnetism is totally suppressed in the M 1 phase: attempts to find a magnetic solution always reverted to a nonmagnetic one with self-consistency.The situation is very different in the metastable M 2 phase, where half of the V pairs dimerize and the other half do not.Nonmagnetic QSGW predicts a metallic phase.An insulating phase forms, however, if the system is allowed to be magnetic.To determine the magnetic structure, each of the four V atoms was assigned an arbitrary moment and the system driven self-consistent.We find that the magnetic moment on the dimerized pair vanishes, while spins on the undimerized pair becomes antiferromagnetically aligned with a local moment of 0.8µ B , which opens a gap of 0.7 eV 11 .The band structure looks remarkably similar to the M 1 phase, even though the physical basis for the gap is very different.Strikingly, one of the two states forming the upper valence band consists almost purely of dimerized V, while the other is almost purely undimerized V.There is very little hybridization between them, or between V and O.
That the physical basis for gap formation differ in the M 1 and M 2 phases was already pointed out in a comment to the Wentzcovitch paper [171].Their argument was based on NMR and EPR evidence for low-lying spin excitations in the M 2 phase, which is consistent with the present work.
The picture from QSG Ŵ is similar to that of the DMFT calculation of Ref. [170] for the M 2 phase, but differs for the M 1 phase.It finds a simple Peierls distortion accounts for the known properties, and magnetism plays essentially no role for the latter.Ref. [170] argues that the temperature-dependence of the bandgap is electronic in origin, and uses this as support for the Mott picture; however, QSGW calculations point to phonons playing an important role in controlling the bandgap at high temperature, with strong support from experimental data [172].Fig. 3 of that work also presented the conductivity derived from the BSE dielectric function, with QSGW as a reference hamiltonian.Agreement with ellipsometry data [173] is quite satisfactory.

F. Antiferromagnetic insulating oxides
The monoxide crystal structures MnO, FeO, CoO and NiO are all of rocksalt form.The magnetic structure consists of sheets of spins antiferromagnetically ordered, which doubles the size of the unit cell.According to the classic paper by Roth [174], the alternating sheets lie in the (111) plane, but the spin orientation depends on the monoxide.MnO and NiO are predicted to be band insulators even within the LDA.In these cases the spin orientation scarcely affects the electronic structure, and we assume the simpler [001] orientation.CoO and FeO are different: LDA predicts both to be metallic.Spins point in the [ 11 7] direction in CoO, and perpendicular to the (111) plane in FeO [174].For these sytems we orient the spin quantization axis along these directions and also do not assume time-reversal symmetry.
All of the rocksalt structure oxides have sizeable magnetic moments.By contrast, CuO is monoclinic with 8 formula units in the unit cell ( §III F 5) and a small local moment (Table IX).
Rödl and Bechstedt modeled the QP band structure of the rocksalt oxides with GW starting from a GGA+U functional [175], and later these authors used the BSE framework to examine the optical response, using a reference potential generated by a G HSE03 W HSE03 functional [176].
In each of these systems, (and probably Fe 3 O 4 §III F 7), QSGW significantly overestimates the bandgap, but not the local moment (Table IX), the difference being more pronounced than in nonmagnetic counterparts.QSG Ŵ greatly amelioriates this overestimate, sometimes slightly overcorrecting QSGW because of the missing vertex.Precise benchmarking is difficult owing to the large uncertainty in experimental data, especially in the strongly correlated cases.One measure of correlation is the Z factor, Eq. 13.Table IX presents a band-and kaveraged Z factor, namely the ratio of the interacting to non-interacting spectral functions A(ω)/A 0 (ω) at an energy just below the Fermi level.The degree of correlation differs in each case so each system is dealt with individually.
1. NiO Fig. 5 shows the (noninteracting) energy bands of NiO, and compares the DOS to BIS data.On the scale of the figure, agreement is excellent.However, QSG Ŵ apparently slightly underestimates the bandgap, which is apparent in both the BIS data and the optics data of Fig. 5. Also, the BSE value for ε ∞ , at 6.15, is outside the range of reported values (Table IX).Replacing Σ(QSG Ŵ ) with a hybrid 0.9 Σ(QSG Ŵ ) + 0.1 Σ(QSGW ) (Eq. 21), the gap increases by 0.17 eV, and ε ∞ decreases to 5.97.This is perhaps the best characterized correlated materials system, though even for NiO there is some spread in reported values for both the dielectric function and the fundamental gap.We can conclude that to within this experimental uncertainty, the close connection between gap and ε ∞ ( §III D 2) is affirmed.[76,117,120] e Optical absorption edge, Ref. [77]; BIS Refs [177].See §III F 2 f Ref. [120] g Optical absorption edge [147]; PES/BIS and XPS [178,179] h Reflectivity, Ref. [180], optical conductivity, Ref. [54] i Reststrahlen spectrum, Ref. [130].See §III F 3 j Refs [77,118,119] k Refs [77,128,129].See §III F 2 l Refs.[126,127] m Refs.[124,125] n reported in Ref. [180].Refs.[181,182] report anomalously large index of refraction (so that ǫ∞∼25-50), which is likely connected to excess holes in nominal La 2 CuO 4 .o values cited in Ref. [183], taken from Ref. [184] and Ref. [185] p spin moment from Ref. [186].Orbital moment estimated to be ∼1µ B .q Ref. [174] r Refs.[187,188] s A consensus value of 0.64µ B ±10% from several sources, Ref. [189] Some ARPES experiments on this correlated antiferromagnet have been published [190].Accordingly we generated the fully dynamical self-energy to compute the k resolved spectral function and compare to it (Fig. 20).The extent to which a particular band is broadened is strongly band-dependent.Agreement with ARPES is satisfactory, in light of the fact that that the matrix element and final state effects would have to be included for a direct comparison.Quite remarkably, the QSG Ŵ spectral function looks nearly identical to one generated by LDA+DMFT, Ref. [191].The close similarity between these two completely different approaches lends support to the thesis that both are characterizing the actual spectral function of NiO.

CoO
Our QSG Ŵ gap is 3.28 eV, and larger than a gap of 2.5 eV measured by a combination of XPS and BIS [177].Optical gap of similar size (∼2.6 eV) has been observed [77].These two experimental findings are consistent only if there are no excitonic effects to reduce the gap.Our QSG Ŵ calculations show, however, that there are a multiplicity of excitons throughout the gap at q=0.The deeper ones (ranging between 0.4 and 1.7 eV) are dark, but strongly active ones at 2.8, 2.96, 3.0, and 3.1 eV also appear.These are likely broadened somewhat, e.g.via some phonon-mediated transitions linking different q, which our calculation does not take into account.Finally, the QSG Ŵ dielectric constant (5.15), aligns well with the mean value of various experiments (5.43 [128], 5.29 [129], 4.75 [77]).If the consistency between gap and ε ∞ argued in §III B 3 can be relied on, it provides another indication that the QSG Ŵ fundamental gap is close to the true one.

MnO
J. van Elp et al. measured the fundamental gap of MnO by x-ray photoelectron and BIS spectroscopies, and obtained a gap of 3.9 eV [193].Three kinds of subgap transitions have been recorded by several groups, labelled as A, B, C transitions, and  [192] compared to BSE@QSG Ŵ .identified with the following symmetries 6 A 1g → 4 T 1g (A band), 6 A 1g → 4 T 2g (B band), 6 A 1g → 4 A 1g + 4 E 1g (C band) [194,195].These transitions are forbidden due to spin and parity selection rules, though significant oscillator strengths have been observed.Huffman et.al reported two additional peaks [196], the highest at ∼3.5 eV.[192] (circles connected by red dots) compared to BSE@QSG Ŵ .Also shown at low energy (circles connected by green dots) is a scaled absorption, α 2 /ω 2 , taken from Ref. [196].
The band structure depicted in the left panel of Fig. 22, is roughly similar to the one depicted by Rödl et al. [176].Unlike NiO, the conduction band is essentially pure Mn s character, the Mn d(t 2g ) appearing at 6-8 eV.The direct gap is at Γ, and is calculated to be 3.6 eV, slightly smaller than the XPS/BIS value (3.9 eV) reported in Ref. [193].
The BSE value for ε BSE ∞ is slightly smaller than observed in a Reststrahlen experiment, Ref. [130] (Table IX).This suggests an inconsistency with the gap being underestimated; however, the maximum value of Im ǫ derived from n and k presented in that experiment is about an order of magnitude too large, so it is not clear how reliable the measurement is.
We find a dark exciton at 3.07 eV and several bright ones at ∼3.5 eV, which can be seen from the shoulder in Im ǫ(ω) below the fundamental gap.These possibly correspond to the highest peaks observed by Huffman et al.. [196].We do not find the weak A and B excitons [194,195]; possibly these are associated with a phonon-assisted transition and an electronic part at finite q.The main shoulder in Im ǫ(ω) starts to rise about 0.3 eV earlier than the reflectance data of Ref. [192] (Fig. 22).This is consistent with the discrepancy in the XPS/BIS measurement of Ref. [193].Thus we tentatively conclude that the QSG Ŵ gap is ∼0.3 eV too low, though reliable experimental evidence is too limited to draw strong conclusions.

FeO
FeO poses one of the most challenging benchmarks in this study.Its highest valence state consists of a single, almost dispersionless d orbital whose m character changes with wave number (Fig. 23).The small Z factor (Tab. IX) provides a clear indication that FeO is strongly correlated.
Experimental information about FeO is sparse and somewhat inconsistent.Two values for ǫ ∞ have been reported: 9.24 [127] and 11.9 [126].The former may be more reliable, since the latter experiment was performed on Fe x O, with x deviating several percent from unity.Bowen et al. investigated the infrared absorption [197], which we use here to benchmark against the BSE.
Regarding the fundamental gap, there is a general expectation in both experimental and theoretical literature that it is of order 2.5 eV [175,183,191,198].Rödl et al. associated the sharp rise in α(ω) observed around 2.4 eV in Ref. [197], with the  [197] compared to BSE@QSG Ŵ (red) with absorption computed from Eq. 19.Also shown are BSE generated from a hybrid of QSG Ŵ and QSGW Σ, Eq.21, with β=0.3 (BSE1) and β=0.6 (BSE2); see Tab.X. fundamental gap.Hiraoka et al. [199] also assumed the fundamental gap was of this order but observed a peak in Im ǫ at ∼1 eV and tentatively assigned it to a defect band.Absorption data shows peaks at both 1.2 eV and 2.4 eV [197]; see Fig. 23.
Turning to theory, at the QSGW level the fundamental gap is found to be 1.9 eV, with the smallest direct gap 2.4 eV.In all the other antiferromagnetic oxides of this study, QSGW overestimates the gap by ∼1 eV (Tab.IX), so these gaps are likely too high.At the QSG Ŵ level, the fundamental gap is much smaller, 0.64 eV.This leads to a puzzle: why does QSG Ŵ yield such a gap so different from the accepted values in the literature?
Counterbalancing the experiments just mentioned, Zimmermann observed the one-particle spectral function XPS/BIS [120].He did not attempt to extract a bandgap, but based on his Fig. 15, it would be of order 1 eV.The XPS/BIS data and the optical measurements both point to the lowest excitation of order 1 eV, even though no deep excitons were found by QSG Ŵ to explain the absorption peak there.Thus our QSG Ŵ analysis suggests a different interpretation, namely that the observed peak in the absorption α(ω) around 1.2 eV (Fig. 23) corresponds to the true fundamental gap.
The QSG Ŵ prediction for E G is likely too small: FeO's practically dispersionless valence band strongly resembles that of VO 2 ( §III E 8), and we can expect the gap to be similarly underestimated in QSG Ŵ .Using materials in Table IV as a guide, the gap can be expected to be underestimated by ∼0.5 eV.Two other pieces of evidence point to the gap being underestimated: ε BSE ∞ is much larger than the two experiments noted earlier (Table X) and the peak in α(ω) falls at ∼0.8 eV, well below the 1.2 eV peak reported in Ref. [197] (Fig. 23).
To adjust for the probable QSG Ŵ gap underestimate, we consider a hybrid of QSG Ŵ and QSGW , Eq. ( 21), and benchmark both ǫ ∞ and the two peaks in absorption, for different admixtures β of QSGW into QSG Ŵ (Fig. 23).Table X shows the variation of E G and ǫ ∞ with β.Perfect alignment with the best available experimental value for ǫ ∞ corresponds to E G =1.05 eV.Thus if FeO has a fundamental gap 1.05-1.10eV, a consistent picture emerges.First the two peaks in α(ω) for BSE2 and BSE1 (Fig. 23 and Table X)) bracket the two experimental peaks from above and below.Second, the one-particle DOS is consistent with XPS/BIS [120].Finally, ǫ ∞ is consistent with the best available experimental data.

CuO
CuO has a monoclinic lattice structure of 4 formula units [200], while the magnetic structure is antiferromagnetic, and is a √ 2×1× √ 2 supercell of lattice with 8 formula units [187,201].The nominal configuration Cu 2+ O 2− would imply a single unpaired d electron; however the magnetic moment is substantially smaller than 1µ B /atom (Table IX).
The QSG Ŵ energy band structure is depicted in the left panel of DFT calculation of Filippetti and Fiorentini [201], who assigned the highest valence to O p (the LDA puts O p too high, so the O character will be overestimated).A band gap has been measured optically from the absorption edge [147], and also by PES/BIS [178] and XPS [179].All three measurements report bandgaps in 1.3-1.4eV range, slightly smaller than the 1.52 eV fundamental gap from QSG Ŵ (Table IX).However, QSG Ŵ overestimates ǫ ∞ (Table IX) by about 20%, which if the consistency between the gap and ε ∞ can be relied on ( §III B 3), the QSG Ŵ fundamental gap must be too small.Moreover, if we compare ǫ(ω) against ellipsometry measurements of Ito et al. [124], the peaks of Im ǫ(ω) are seen to fall ∼0.3 eV below the experimental data.From this we conclude it is likely that the fundamental gap is closer to 1.6 eV, assuming the dielectric data of Ito et al. [124] is reliable.This would mean the PES/BIS is underestimated.It is perhaps not unexpected since the BIS should be larger than the optical gap.
That experimentally Im ǫ(ω) is smoother than the QSG Ŵ one, can be attributed (at least in part) to ǫ(ω) being generated from a noninteracting G 0 (QSG Ŵ ).The frequency-dependence of Σ reduces quasiparticle weights (compare the DOS of the interacting G to that of G 0 , left panel of Fig. 25), and the imaginary part smears out the quasiparticle (right panel of Fig. 25).CuO is very strongly correlated: note the sharp reduction in the DOS around 2 eV.These dynamical effects do not shift the average position of the bands (owing to the QSGW construction) but will smooth out transitions between occupied and unoccupied states, and correspondingly, the imaginary part of the longitudinal dielectric function in the basal plane, Im ǫ xx .

La2CuO4
La 2 CuO 4 is the parent compound for one of the most widely studied superconductors.A gap forms because Cu d NiO.Such mass renormalization plays a critical role in the correlations of this orbital, which drives superconductivity.Spin fluctuations will narrow this band still further, but whether low-order perturbation theory will be sufficient to yield the true bandwidth remains an open question.
• The O p band is pushed down relative to QSGW , thus reducing the hybridization of O into the Cu d x 2 −y 2 state.The LDA often misaligns orbitals of different character, but it is notable that ladder diagrams not only reduce the gap, but induce a shift to the QSGW valence bands.
• For both QSG Ŵ and QSGW the La 4f states are pushed well above the Fermi level, something the LDA fails to do.
• Self-consistency plays a very important role in this system (compare G LDA W LDA to QSGW energy bands, and see Table I).
As with NiO, the G LDA W LDA bandgap is severely underestimated [34].The one-shot gap can be improved by using LDA+U or a hybrid functional instead of the LDA, but the resulting energy bands depend on the choice, as will other parts of the spectrum (e.g. the position of O 2p states).
The left panel of Fig. 27 shows two measurements of the dielectric function, Im ǫ xx , one inferred from reflectivity at 122 K [180], and the other from low-temperature optical conductivity [54].Ref [180] also shows results from a photoconductivity measurement, which looks similar to the blue squares in the figure but slightly blue-shifted.QSG Ŵ results are also shown: the peak in Im ǫ xx appears at slightly higher energy (0.1-0.2 eV) than the experimental data.The QSG Ŵ result has sharper structure, in particular there appears a pronounced sub-gap peak centered at ∼1.5 eV.A corresponding peak (albeit much weaker) is seen in the 122 K reflectivity data, though this peak is washed out as the temperature increases [180].QSG Ŵ predicts a spectrum of 30 or so subgap excitons, ranging between 1.2 eV and the fundamental gap with widely varying oscillator strengths.A particularly bright exciton appears at 1.45 eV; it is is partly responsible for the peak in Im ǫ xx there.As for the fundamental gap, QSG Ŵ predicts an indirect gap of 1.66 eV, but the lowest direct gap is ∼2.1 eV (Fig. 26).Ref. [180] assigned a charge transfer gap of 2.1 eV, and Ref. [54] a similar gap (2.2 eV), which probably corresponds to the direct gap.The QSG Ŵ result for Im ǫ xx shows sharper peaks than the experiment, as was shown for CuO ( §III F 5).For the same reason explained there, dynamical effects will smooth out Im ǫ xx .The conventional picture, originating from Verwey, is that the tetrahedral sites are Fe 3+ and octahedral sites consist of equal numbers of Fe 2+ and Fe 3+ .Below the Verwey temperature magnetite is a narrow-gap insulator with a bandgap 0.14-0.3eV.It was traditionally believed that above the Verwey temperature magnetite becomes a half-metal, in part because the conductivity increases by ∼100-fold across the transition, and evidence from interpretations of PES and STS experiments suggested a finite density-of-states at E F .However, a more recent high-resolution PES experiment [203] found that the band gap persists above the Verwey temperature, reduced by ∼50 meV.That a gap of order 0.2 eV persists was confirmed by a subsequent STS measurement [204].For a more detailed summary of the experimental status, see the work of Liu and Di Valentin [205].These authors applied various one-body techniques (LDA+U, hybrid functionals) to study magnetite and concluded that the traditional picture of magnetic order yields a metallic ground state.They argued that a gap appears because the octahedral sites disproportionate into two kinds of atoms, one with a high-spin (moment 4µ B ) and one with a low-spin (moment ∼3.5µ B ).We performed QSGW and QSG Ŵ calculations with symmetry suppressed, so every spin could assume an independent value.Spins on the tetrahedral sites converged to a common value and those on the tetrahedral sites converged to two values, one high spin and one low spin (Table XI).The calculation was extremely difficult to stabilize and a fully converged solution was never found, even after 100 iterations.Local moments were stable as iterations proceeded, but the bandgap fluctuated; for that reason error bars are given for the gaps shown in Table XI.Here we denote high-spin and low-spin sites as "o1" and "o2." The Fe o2 t ↓ 2g manifold splits off one band (actually two, because the unit cell consists of two Fe 3 O 4 formula units which weakly couple), and this split-off band forms the valence band maximum (Fig. 28).On the other hand, the Fe o1 t ↓ 2g manifold does not split in the same way, and it forms the conduction band minimum.Splitting of Fe o2 t ↓ 2g is the main way in which cubic symmetry is broken and a gap is formed.
Ladder diagrams have two major effects: first, they sharply narrow all of the octrahedral Fe-d bandwidths, e.g. the t ↓ 2g manifold forming the conduction band minimum is narrowed by ∼30%.Second, ladders cause band centers to shift in a highly orbital-dependent manner.Shifts on the occupied states are modest, but for the unoccupied states they can be quite large.Note for example the band center of the t ↓ 2g level around 1.5 eV is pushed down by ∼1 eV, while the t ↓ 2g and e ↓ g levels in the 5-6 eV range shift by ∼2 eV.Unfortunately, no experiments are available for benchmarking.
Ordering of the Fe levels is qualitatively in agreement with the picture of Ref [205]; see their Figure 5, and the QSG Ŵ spin moments are similar to their Table 1; we thus affirm their description of magnetite.

IV. CONCLUSIONS
We presented an extension of the QSGW approximation in which we include vertex corrections to W , calculated at the level of the BSE (QSG Ŵ ).The primary aim of this work was to establish to what extent QSG Ŵ rectifies the most severe errors in QSGW , with the ultimate aim to develop a high-fidelity, universally applicable theory.If low-order diagrams are sufficient to yield high-fidelity one-and two-particle properties, Green's function methods offer an enormous opportunity to be both high-fidelity and relatively efficient.QSG Ŵ supplies excitonic effects which QSGW omits.They are known to be important, as already extensively discussed in the literature (see e.g.Ref. 28), so it is of interest to see to what extent it captures the quasiparticle spectrum in insulators, and whether systematic errors can be discerned.
While the QSGW approximation has long been known to overestimate bandgaps, the discrepancies with experiments are much more systematic than more commonly used G 0 W 0 approaches with G and W constructed, e.g. from the LDA.Its systematic character is a consequence of self-consistency, in part because it does not rely on density-functional theory.We presented a few contexts where a non self-consistent approach is not important (e.g.Bi 2 Te 3 ), and others (e.g.TiSe 2 , La 2 CuO 4 ) where it is fundamentally problematic.Even when such an approach yields a good bandgap, it may occur for adventitious reasons.The relatively unsystematic nature of the errors in one-shot approaches make it difficult to assess what diagrams are essential to realize the goal of a universally applicable, high-fidelity theory. 12Thus, self-consistency is crucial for the aims of this work.
The present work surveyed a wide range of insulators, including tetrahedrally coordinated semiconductors where experimental information is reliable and abundant, and also a variety of other sp systems, d 0 oxides, and polar compounds, and a family of 3d transition metal antiferromagnetic oxides.Each materials system had a distinct set of characteristics, but apart from some important exceptions critically examined in this work, QSG Ŵ predicts with fairly high fidelity both one-particle and optical properties for all of the systems we studied.The exceptions are important and formed a major focus of this study.Two shortcomings clearly identified were the omission of electron-phonon interaction, which causes gaps to be too large in widegap systems, and the omission of the vertex in the exact self-energy.This vertex pushes down nearly dispersionless core-like states, and when they form the valence band maximum the bandgap is consistently underestimated.By constructing hybrid self-energies, we could account for both of these shortcomings in an approximate way, and draw the following conclusions: 1.At the QSG Ŵ level, there is a very close connection between the fidelity of the fundamental gap E G and the dielectric constant ǫ ∞ .When one is well described, so is the other, and vice-versa.This provides a much more robust benchmark of a theory than benchmarking one-particle properties alone.
2. If we take the first point as an ansatz for a general principle, it can be used in cases where experimental data is unavailable or inconsistent.We presented evidence for several systems (CeO 2 , SrTiO 3 , TiO 2 , ScN, CuAlO 2 , FeO) where the calculated results inform the experimental observations and indicate that accepted values of the one-particle properties need adjustment.For FeO, the revision is rather dramatic.
3. A low-order diagrammatic theory appears to describe the dielectric response with high fidelity for all the systems in this study, to the extent we are able to reliably extract experimental data.Ladder diagrams appear to be sufficient to capture well the main part of the optical response functions and one-particle Green's functions in most insulators, even strongly correlated ones.While such an assertion is likely not universally true [46], it appears to be the case for broad classes of materials.

4.
Ladders not only shift the bandgap but further narrow the d bandwidth in some systems (NiO, La 2 CuO 4 , Fe 3 O 4 ).It may be that the addition of a low-order GW-like theory accounting for spin fluctuations, such as the dual-trilex formulation of Stepanov et.al [59], may adequately account for spin and charge response functions even in strongly correlated materials.
These last two observations suggest the tantalizing possibility that, with some modest extensions that may be added hierarchically, a broadly applicable, high-fidelity ab initio approach to solving one-and two-particle properties of the many-body problem is within reach.
G −1 simplifies to a linear function of ω and thus reduces to a linear algebraic eigenvalue problem.The bar over the Z factor indicates that is not equivalent to Eq. 13, since it is defined at zero frequency ω=0: 1 − 1/ Zj = Σ ′ (k, 0) Evidently ǫ − µ + Σ(k, 0) is the eigenvalue of a hamiltonian defined as the one-body part of G −1 , but including the static part of Σ.The (linearized) energy-dependence of Σ modifies this eigenvalue to read E is identical to the QSGW quasiparticle energy if Σ is a linear function of ω.Now let us retain the quadratic term in Σ and determine the shift in E to estimate the difference between LQSGW and QSGW .Let us denote the LQSGW eigenvalue E − µ as E 0 .Expanding G −1 to second order we obtain, to lowest order in Σ ′′ (k, 0) : The lowest-order difference between LQSGW and QSGW QP levels is the second term in parenthesis. .Gcut(ψ) and Gcut(M ) are planewave cutoffs for the interstitial part of the one-particle and two-particle basis sets, in units of 2π/a; Σcut the energy cutoff for above which Σ is restricted to a diagonal part, as described in the text.N k is the number of divisions along each reciprocal lattice vector defining the k mesh.When two numbers appear, the c axis is assigned a mesh different than the basal plane.The latter number is selected to make the spacing between k points as similar as possible along the three directions.N1p is the total number of basis functions in the unit cell.Nv and Nc are the number of occupied and unoccupied eigenstates included in the construction of the vertex.N flt is the number of points in the interstitial where envelope functions were added to increase the basis completeness.ϕz lists the local orbitals (LO) for each element: n − s, n − p, n − d, where n is the principal quantum number of the LO.n without an overbar indicates the LO covers a core-like state, well below the linearization energy with a principal quantum number one less than that of the valence.n indicates the LO energy is far above the linearization energy, and is included to better treat unoccupied states well above the Fermi energy.Both kinds of LO are discussed in Ref. [38].Partial waves marked as p replace the l=1 partial wave with the corresponding p 1/2 partial wave computed from the Dirac equation, as discussed in Ref. [38].This is a small effect but it improves the matrix elements for spin-orbit coupling.
small value V TF , typically 2×10 −5 Ry, though sometimes somewhat larger values, up to 2×10 −4 Ry were used.The dielectric constant, ǫ ∞ , can vary by a few percent over this range.For that reason ǫ ∞ was calculated for several values of +V TF , e.g.1×10 −5 , 1×10 −5 , and 3×10 −5 Ry, and the reported value is the result when extrapolated to zero.
Frequency mesh: to construct the self-energy, an energy integration on the real frequency axis is taken.A regular quadratic mesh of the form ω i = dw×i + dw 2 i 2 /(2ω c ) is used, with i spanning ω i =0 and the largest eigenstate.Points are linearly spaced for dw ≪ ω c , but the spacing increases for dw ω c .It has been found empirically that results are essentially independent of mesh for dw<0.08Ry and ω c 0.1 Ry.In practice we use dw=0.02Ry and ω c =0.2 Ry to obviate the need for checking convergence.To pick up the poles of G and W to make Σ, the contour is deformed to include an integration on the imaginary axis of ω (I, §2F).In all the calculations used here, we used 6 points on a Legendre quadrature.A few checks showed that the result hardly depended on the number of points in the quadrature.
Manual vs auto-generated input: Questaal has an automatic generator, blm, to construct input files from structural data.Most input parameters are automatically generated by blm, such as the MT radii r MT , the product basis cutoffs, and the plane wave cutoffs, the Gaussian smoothing radius defining the envelope functions, and the placements for floating orbitals, when they are sought.Also for the vast majority of parameters, the code uses default values if inputs are not explicitly specified.For a few parameters, manual intervention is needed to monitor convergence, especially the number of k points and the plane wave cutoffs G cut (ψ) and G cut (M ).Hankel function energies E must be manually set, but usually fixed values as noted above are sufficient.Occasionally interpolation continues to be an issue and can be stabilized by making E deeper, e.g.E= − 0.6 Ry was needed to stabilize SrTiO 3 .Results are largely insensive to the choice of E, provided it is not pushed too deep.

startFIG. 1 .
FIG.1.Flow chart of the QSGW cycle.Non-interacting eigenfunctions and energies (ψ, ǫ) are calculated self consistently (blue).These are used to construct the non-interacting Green's function G 0 , Coulomb interaction v, RPA polarization P 0 = − iG 0 G 0 , and W =(1 − vP 0 ) −1 v. W is used to make a vertex and better P via Eq.14, which gives the improved Ŵ .One cycle makes the static self-energy Σ 0 , that is passed to H 0 (green) and the cycle repeated to self-consistency.

F
. Divergence of the macroscopic dielectric function at q = 0

FIG. 4 .
FIG. 4. (top)Fundamental bandgap for selected materials calculated within the LDA (red squares), QSGW (blue circles), and QSG Ŵ (green diamonds).The black crosses add to QSG Ŵ an estimate for the gap correction from the electron-phonon interaction when it exceeds 0.19 eV.Where available, this was taken from Ref.[41]; otherwise it was estimated from the Frölich expression, Eq. (24).(bottom) ε RPA

FIG. 5 .
FIG. 5. Middle panel: Energy bands in NiO, within the QSG Ŵ approximation (colored bands), with valence band maximum at 0. Green projects onto Ni d character, blue onto Ni sp character and red onto O p character.Shown for comparison are the QSGW results (light gray bands).Left panel: QSG Ŵ (green) and QSGW (gray) DOS, compared against PES data (E<EF ) and BIS data (E>EF ) from Ref. 76 (circles).Right: experimental dielectric function Im ǫ(ω) from Ref.[77], compared to results calculated at three levels of approximation: RPA@QSGW , BSE@QSGW , and BSE@QSG Ŵ .As is typical, the shoulder of RPA@QSGW is blue-shifted, by roughly 2 eV in this case.Adding ladders (BSE@QSGW ) shifts the shoulder towards the experiment, but it is still ∼1 eV too high, as a consequence of overestimate of the QSGW fundamental gap (TableIX).BSE@QSG Ŵ describes Im ǫ(ω) rather well, including peaks around 6 eV, 13.5 eV and 17 eV.However the shoulder around 3.5 eV is slightly red-shifted compared to experiment, indicating that the fundamental gap is underestimated.ε∞ is also overestimated.See §III F 1 for more details.

FIG. 12 .FIG. 13 .
FIG.12.Imaginary part of the macroscopic dielectric function for LiF.The experimental data (blue squares)[131] is compared with the results from the QSG Ŵ 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 (∼22 eV) is 1.25 eV.
FIG. 14. (Left) energy bands for ScN.Red and green depict projection onto Sc and N orbitals, respectively.(Right) real and imaginary parts of the BSE dielectric function.The grey dashed line depicts the RPA dielectric function computed from the QSG Ŵ self-energy.The dotted line is a guide to the eye, extrapolating the second shoulder in Im ǫ to zero.Table below: Experimental and QSG Ŵ gaps and values for ǫ∞.

FIG. 15 .
FIG.15.Left: energy band structure in CeO2 computed by QSG Ŵ .Red, green, and blue correspond to Ce-f, Ce-d, and O-p character, respectively.Light gray dashed lines show corresponding QSGW bands.Right: Corresponding dielectric function.Circles are absorpance data digitized from Ref.[142].Red curve is the BSE absorption α generated from the QSG Ŵ hamiltonian (α=4πn/λ, n 2 =dielectric function).The energy axis for the calculated function is redshifted by 0.5 eV as described in the text.

FIG. 16
FIG. 16. (Left) energy bands for TiO2.Colored bands are taken from QSG Ŵ calculations, with red and blue showing projections onto Ti d and O, respectively.Dashed gray lines show corresponding QSGW results.(Right) dielectric function for polarization parallel and perpendicular to the z axis.Circles connected by dotted red lines are data digitized from Ref. 105.BSE Im ǫ ⊥ (ω) and Im ǫ || (ω), computed from QSG Ŵ , are shown in green.Light dashed lines compare the RPA results generated from the same QSG Ŵ hamiltonian.

3 FIG. 17 .
FIG. 17. (Left) energy bands for SrTiO3.Colored bands taken from QSG Ŵ calculations, with red and blue showing projections onto Ti 4d and O 2p, respectively.Gray lines show corresponding QSGW results.(Right) dielectric function.Circles connected by dotted red lines are data digitized from Ref. 105.The black line shows the Im ǫ BSE (ω) generated from a QSG Ŵ hamiltonian on a 10×10×10 k-mesh (the waviness is an artifact of incomplete k convergence).The grey line shows Im ǫ RPA (ω) generated from the same hamiltonian.Im ǫ RPA (ω) vanishes at the fundamental direct gap.
FIG. 18. (Left) energy bands for CuAlO2.Colored bands are taken from QSG Ŵ calculations, with red, blue, and green showing projections onto Cu d, O, and Cu sp, respectively.Dashed gray lines show corresponding QSGW results.Points S, F, and L correspond to (2/3,1/3,0), (1/2,1/2,0), and (0,0,-1/2), respectively, as multiples of the reciprocal lattice vectors.(Right) dielectric function (average of x and z directions) as function of frequency ω (eV).Solid line shows Im ǫ BSE (ω), generated from QSG Ŵ .The green dashed line shows Im ǫ RPA also generated from the QSG Ŵ hamiltonian.Note Im ǫ BSE and Im ǫ RPA approach 0 near the fundamental direct gap at 4 eV; however, Im ǫ BSE has additional subgap peaks from strongly bound excitons.The blue dotted line shows Im ǫ RPA (ω) generated from QSGW .The RPA functions look similar except for a 1 eV shift, because of the difference in bandgaps.

6 Z
FIG.19.Left: energy band structure of VO2 in the M1 phase within the QSGW approximation.Red, green, and blue project onto V d(m=−2,−1,0), V d(m=1,2) and O orbitals, respectively.The QSG Ŵ bands (not shown) are essentially the same but with a 0.3 eV smaller bandgap (see text).Right: corresponding bands of the M2 phase.In this figure red and green project onto the dimerized and undimerized V d orbitals, respectively.

3 FIG. 25 .
FIG. 25. (Left) Density-of-states generated from G compared to that generated from the QSG Ŵ G0.Zero energy corresponds to the valence band maximum.(Right) spectral function from interacting G generated by QSG Ŵ .Yellow lines below −2 eV trace out the QSG Ŵ bands, and are equivalent to those in the left panel of Fig. 24.

FIG. 26 .
FIG. 26.Energy band structure of La2CuO4.(a) and (b) panels apply to QSG Ŵ and QSGW approximations, respectively.Colors depicts the following orbital character: red: Cu d x 2 −y 2 ; green: Cu d 3z 2 −1 ; blue: O 2p; cyan: La 5d.Black valence bands are the remaining Cu 3d orbitals.Panel (c) is the LDA band structure, with the same color scheme.The narrow window of black bands at 3-4 eV are of La 4f character.Panel (d) is G LDA W LDA result, with Z=1 and including the off-diagonal parts of Σ.

3 FIG. 27 .FIG. 28 .
FIG. 27. (Left) BSE dielectric functionIm ǫxx in La2CuO4, computed from QSG Ŵ , compared to reflectivity data from Ref.[180] (Expt1) and conductivity data from Ref.[54] (Expt2).Data in Ref.[180] is anomalously small.They also report ǫ∞∼5, which does not seem to be compatible with their scale for Im ǫxx, from the Kramers Kronig relation.Expt 1 shown in the figure scales data taken from Ref.[180] by a factor of 5 to bring it approximately in line with Ref.[54].(Right) Spectral function from the interacting G generated by QSG Ŵ .Yellow lines are energy bands from the QSG Ŵ G0.
[92]in the G LDA W LDA context[92]may be a factor.It cannot be ruled out that the nearly perfect agreement for E1 and E2 transitions in zincblende semiconductors, where experimental data is available.Red squares, blue circles, and green diamonds correspond to LDA, QSGW and QSG Ŵ matching

TABLE III .
[91]while those at L are not.E0 in narrow gap zincblende semiconductors is overestimated by 0.1 eV.The tendency does not hold for narrow-gap semiconductors that form in other structures, shown as the entries in the second half of the Table.An estimate for low-temperature bandgap for Ti2Se2 (24-atom P 3c1 CDW structure) is taken from Ref.[91].The table below shows the mean difference with the photoemission data, and the RMS fluctations about the mean, for the different levels of approximation.

TABLE IV .
Bandgaps in systems with valence band maximum formed from a core-like d state.The bandgap of FeS2 is not well known, but the QSG Ŵ gap lies below the most likely value of 0.9 eV.

TABLE VI .
Birefringence in selected insulators.

Table below :
Experimental and QSG Ŵ gaps and values for ǫ∞.
[140].[140].b Based on analogy with InN, as described in the text

TABLE IX .
Bandgap; ε∞; local magnetic moment; band-and k-averaged Z factor in selected antiferromagnetic insulators.Strongly orbital dependent: Z∼0.65 for the (mostly O) valence bands, and ∼0.45 for the (mostly Cu) conduction bands c Strongly state-dependent: Z∼0.65 for the highest valence bands, and ∼0.4 for lowest conduction band d Refs. b

TABLE X .
(21)cal properties in FeO as a function of hybridization parameter β, Eq.(21).EG is the fundamental gap in eV; EG (Γ→Γ) is the direct gap at Γ. Labels in the second column are used in Fig.23.

TABLE XI .
Spin moments and bandgap of Fe3O4 calculated with the QSGW and QSG Ŵ approximations.Fet, Feo1, Feo2 indicate Fe on tetrahedral sites, and the two inequivalent octahedral sites.

TABLE XII .
a(a.u.)Gcut(ψ, M ) Σcut(Ry) N k N1p Nv Nc N flt ϕz Materials parameters: a is lattice constant (quantity in parenthesis is c/a where applicable)