The GW Compendium: A Practical Guide to Theoretical Photoemission Spectroscopy

Abstract

The GW approximation in electronic structure theory has become a widespread tool for predicting electronic excitations in chemical compounds and materials. In the realm of theoretical spectroscopy, the GW method provides access to charged excitations as measured in direct or inverse photoemission spectroscopy. The number of GW calculations in the past two decades has exploded with increased computing power and modern codes. The success of GW can be attributed to many factors: favorable scaling with respect to system size, a formal interpretation for charged excitation energies, the importance of dynamical screening in real systems, and its practical combination with other theories. In this review, we provide an overview of these formal and practical considerations. We expand, in detail, on the choices presented to the scientist performing GW calculations for the first time. We also give an introduction to the many-body theory behind GW, a review of modern applications like molecules and surfaces, and a perspective on methods which go beyond conventional GW calculations. This review addresses chemists, physicists and material scientists with an interest in theoretical spectroscopy. It is intended for newcomers to GW calculations but can also serve as an alternative perspective for experts and an up-to-date source of computational techniques.

Keywords: GW approximation; Hedin's equations; band structure; electron affinity; ionization potential; quasiparticle; self-energy; theoretical spectroscopy.

Figure 1

Schematic of the photoemission (PES) and inverse photoemission (IPES) process. In PES (A) an electron is excited by an incoming photon from a previously occupied valence state (lower shaded region) into the continuum (gray shaded region, starting above the vacuum level E_vac). In IPES (B) an injected electron with kinetic energy E_kin undergoes a radiative transition into an unoccupied state (upper shaded region) thus emitting a photon in the process.

Figure 2

(A) Schematic representation of an ARPES experiment. By varying the angles θ and ϕ with respect to the crystallographic axes (a_i), the measured spectrum is direction, or k, dependent. In practice, the detector angle is usually varied with respect to a fixed beam. (B) A typical spectral function features a sharp peak attributed to the quasiparticle, an incoherent background, and satellites away from the single particle peak. (C) ARPES data of the upper valance bands of ZnO (Kobayashi et al., 2009). The corresponding G₀W₀ band structure of ZnO is shown in Figure 27.

Figure 3

Top: Depiction of the quasiparticle concept. (A) A crowd of people is analogous to the electronic ground state. A new person (that represents an additional electron) enters the crowd in (B). The new person begins to interact with other people who, in turn, interact back with the new person in (C) and form a polarization cloud. An effective, or renormalized, object, the quasi-person, moves through the crowd in (D). Even though it is an interacting system, the many-person state in (D) can still be connected to, or identified by, the single person added to the crowd. This connection allows us to identify the quasi-person. Bottom: Schematic representation of photoemission spectroscopy.

Figure 4

X-ray photoemission spectrum with 800 eV incident energy compared to two calculated spectra. The red line shows the evGW₀@LDA spectrum (see section 5), whereas the blue spectrum contains additional vertex corrections in form of a cumulant expansion (see section 11). The evGW₀@LDA+C^* spectrum contains the addition of the Shirley background (shown by the black dashed line) and loss effects of the outgoing photoelectron. Data retrieved from Guzzo et al. (2011), where the GW results are labeled as G₀W₀. However, self-consistency in the eigenvalues was in fact applied, which is in our notation evGW₀ (Private Communication).

Figure 5

The most basic pieces of diagrammatic perturbation theory are G₀ and v. From these, all other quantities can be built. The interaction v(1, 2) is instantaneous. Therefore, the dashed line is perpendicular to the time axis. The arrows in G₀ and G point in only one direction, but both time orderings are included.

Figure 6

The exact G contains amplitudes from all possible paths between 1 and 2. Amplitudes from all of these paths are represented by the rectangle placed between the field operators, the action of which is represented by * symbols. These terms can be calculated order-by-order with perturbation theory. At a given order n, we must connect n interaction lines at internal times in all possible − and allowed − ways. Concrete examples of diagrams are in Figure 7. All terms of the topology which can be inserted between two G₀ lines form the reducible self-energy.

Figure 7

At first order, n = 1, there are only two possible self-energy diagrams. These are the diagrams of the Hartree-Fock approximation, the direct electrostatic interaction (left) and exchange (right). Two possible n = 2 diagrams are also shown (there are others). The bottom diagrams are forbidden because they do not have two G₀ lines at each end of the interaction lines. When drawing the diagrams, a certain degree of flexibility is allowed and they must be interpreted carefully. For example, the curved interaction lines above must still be treated as instantaneous in a calculation.

Figure 8

The exact vertex Γ, shown in (A), can be replaced with approximations to simplify the calculation. The approximation in (B) is referred to as a “single spacetime point” because the vertex has no internal structure. In contrast, the vertex in (C) has internal structure. The diagram shown here is only an example to demonstrate the role of Γ and does not correspond to the exact self-energy or the GW self-energy.

Figure 9

Diagrammatic representation of Equations (13–17). The GW approximation reduces the self-energy to a product of G with W. The first equation (Dyson's equation) has a G line on the left- and right-hand sides. This equation can be iterated, inserting G₀ + G₀ΣG in place of each G on the RHS, forming the Dyson series. The same iterative procedure for W forms its own Dyson series.

Figure 10

Flowchart for a G₀W₀ calculation starting from a KS-DFT calculation. The KS energies $\left\{{ϵ}_s^{\text{KS}}\right\}$ and orbitals $\left\{{\varphi }_s^{\text{KS}}\right\}$ are used as input for the G₀W₀ calculation. For the full expressions of χ₀, ε and $W_0^c$ see Equations (25–28) and (30). The spin has been omitted for simplicity.

Figure 11

(A) Real and (B) imaginary part of the self-energy Σ^c(ω). Displayed is the diagonal matrix element ${\Sigma }_s^c=\langle s\left|{\Sigma }^c(\omega )\right|s\rangle$ for the HOMO of the water molecule. The gray-dashed line at ≈ − 12.0 eV indicates the QP solution ϵ_s. (C) Spectral function A_ss(ω) computed from Equation (37). The PBE functional is used as starting point in combination with the cc-pV4Z basis set. Further computational details are given in Appendix C.

Figure 12

Error introduced by linearizing the QP equation, ${\Delta }_{\text{Zshot}}=|{ϵ}_s^{\text{iter}}-{ϵ}_s^{\text{Zshot}}|$, where ${ϵ}_s^{\text{iter}}$ has been obtained from the iterative procedure and ${ϵ}_s^{\text{Zshot}}$ from Equation (34). “HOMO-x” indicates deeper valence states. The PBE functional is used as starting point in combination with the cc-pV4Z basis set. Further computational details are given in Appendix C.

Figure 13

Contour deformation technique: Integration paths in the complex plane to evaluate Σ^c(ω). Γ⁺ and Γ⁻ are the integration contours, which are chosen such that the poles of G₀, but not the poles of W₀ are enclosed. Γ⁺ encircles the upper right and Γ⁻ the lower left part of the complex plane. ω′ denotes frequencies of the integration grid and ω the frequency at which Σ^c is calculated.

Figure 14

G₀W₀@PBE self-energy matrix elements for the HOMO of benzene obtained with different frequency integration techniques: contour deformation (CD) and analytic continuation (AC) using the Padé model with 128 parameters and the 2-pole model. See Appendix C for further computational details.

Figure 15

Schematic representation of the projector augmented wave (PAW) scheme. The all-electron wave function ϕ is constructed from the smooth auxiliary function $\tilde{\varphi }$ and corrections from the hard and smooth atom-centered auxiliary wave functions ϕ^a and $\widetilde{{\varphi }^a}$, respectively.

Figure 16

Basis set convergence for G₀W₀ calculations. (A) Convergence for a plane wave basis set. Bandgap of wurtzite ZnO dependent on the number of bands and on the corresponding cutoff energy (data from SI of Yan et al., 2012). (B) Convergence and extrapolation procedure for a localized basis set. Ionization potential (IP) for the HOMO of benzene plotted with respect to the inverse of the number of basis functions N_func using the cc-pVnZ basis set series. Further computational details are given in Appendix C.

Figure 17

Non-self-consistent Sternheimer approach for obtaining W₀ without empty states. Δϕ⁰ and v have a parametric dependence on the real space point r. Δϕ⁰ depends additionally on the frequency ω.

Figure 18

Starting point dependence of G₀W₀: the left side shows the G₀W₀ HOMO energy of the water molecule for hybrid functional starting points with different amounts of exact exchange. The HOMO energy in self-consistent GW (scGW) is shown on the right. The dashed line marks the experimental value of 12.62 eV (Page et al., ; Lias and Liebman, 2003). All GW values are extrapolated to the exact basis set limit using the cc-pVnZ (n = 3–5) basis sets. Further computational details are given in Appendix C.

Figure 19

CSP scheme representative for a small molecule. $\Delta v_s^c$ (Equation (70)) is plotted with respect to $\Delta v_s^x$ (Equation (69)) for a set of occupied states s. The HOMO and HOMO-1 states are indicated. The new α value is obtained from the slope of the straight line fitted through the red symbols. Data retrieved from Körzdörfer and Marom (2012).

Figure 20

Schematic representation of the straight line condition for total energies E (left) and derivatives ∂E/∂f (right). f is the occupation number. The DSLE is shown for three different cases: convex (blue), concave (red), and mixed curvature (green). Reprinted with permission from Dauth et al. (2016). Copyright (2016) by the American Physical Society.

Figure 21

Scaling of state-of-the-art G₀W₀ implementations with respect to system size using graphene nanoribbons as a benchmark system. (A) Smallest graphene nanoribbon unit with 114 atoms. (B) Comparison of the scaling of the canonical G₀W₀ (Wilhelm et al., 2016) and the low-scaling implementation (Wilhelm et al., 2018). The latter requires operations of at most $\mathcal{O}(N^3)$ complexity (red diamonds). Dashed lines represent least-square fits of exponent and prefactor. Data retrieved from Wilhelm et al. (2018). Both algorithms are implemented in the CP2K program package.

Figure 22

G₀W₀ and scGW in terms of Feynman diagrams. In G₀W₀, the irreducible self-energy is constructed from G₀ and W₀. The Green's function updated with the lowest order self-energy, G₁ (shown as the bold Green's function line), contains an infinite series of Σ insertions. In fully self-consistent GW, the starting point dependence is removed and all quantities in the diagrammatic expansion are fully dressed. Here, we assume a true G₀ starting point instead of a mean-field G₀ so that subtraction of v^MF is not necessary to include in the diagrams.

Figure 23

Two of the diagrams of the screened Coulomb interaction in scGW that are not present in G₀W₀.

Figure 24

Schematic of Hedin's full set of equations (A) and Hedin's GW approximation (B–D). In (A), all five quantities are iterated to self-consistency. In (B), self-consistent GW (scGW), Γ is set to a single spacetime point and the remaining four quantities are determined self-consistently. Eigenvalue self-consistent GW shown in (C), evGW, updates only the quasiparticle energies while leaving the wave functions unchanged. In the scGW₀ or evGW₀ procedures shown in (D), one iterates G to self-consistency in Dyson's equation but does not update χ₀ or W.

Figure 25

Self-energy matrix elements for the HOMO of a single water molecule obtained with G₀W₀, evGW₀ and level-aligned G₀W₀. In all three cases PBE is used as starting point. The inlet shows the graphical solutions of the QP equation. See Appendix C for further computational details.

Figure 26

(A) Band gaps of semiconductors and insulators computed with PBE, G₀W₀, and evGW₀ in the all-electron, linearized augmented plane wave (LAPW) framework. Data taken from Jiang and Blaha (2016). (B) Band gaps of semiconductors and insulators computed with PBE, G₀W₀, QSGW, and scGW in the projector-augmented-wave (PAW) framework. Data taken from Grumet et al. (2018).

Figure 27

G₀W₀ band structure of ZnO superimposed on experimental ARPES data (Yan et al., 2012). The experimentally measured lifetimes of the states are indicated by the shading, with white shading indicating long lifetime. The G₀W₀ calculations are based on the optimized effective potential approach for exact exchange mentioned in section 5.3 that includes LDA correlation [OEPx(cLDA)]. Reprinted with permission from Yan et al. (2011). Copyright (2011) by IOP Publishing Ltd.

Figure 28

G₀W₀ band structure of K₂Sn₃O₇ (McAuliffe et al., 2017). The main panel illustrates the difference between the PBE (dark gray lines) and the G₀W₀@PBE (red lines) band structure. The unoccupied states of the PBE band structure have been shifted up in energy for better visibility so that the bottom of the conduction bands coincide in both band structures. The right panel shows the G₀W₀@PBE density of states (DOS) resolved into s, p and d angular momentum channels.

Figure 29

For systems with mid-gap defect levels, computing the band gap alone is not enough to test the material for potential applications. The position of the defect level can also be computed with GW. In these cases, the absolute position of the VBM, CBM, and defect level are important.

Figure 30

Illustration of the image effect. (B) Shows the image charge and image potential induced by an additional electron (e.g., anionic charge on a molecule) outside a surface. (A) Provides a graphic illustration how the image potential of a germanium (Ge) surface could be probed with a carbon monoxide (CO) test molecule. By adding thicker and thicker sodium chloride (NaCl) layers between CO and Ge, the CO molecule moves along the Ge image potential. The resulting CO gap will then depend on the NaCl layer thickness, which is indeed the case as (C) illustrates. Subfigure (A) adapted from Freysoldt et al. (2009) under the terms of the Creative Commons Attribution 3.0 License. Data for (C) obtained from Freysoldt et al. (2009).

Figure 31

The image potential of a repeated slab system (B) differs from that of an isolated surface (A). The dashed lines in (C) mark the difference that can be computed with a suitable correction scheme (Freysoldt et al., 2008). As the charge moves across the interface, the ratio of dielectric constants for the “charged” and “uncharged” regions changes. As a result, the image potential changes sign.

Figure 32

(A) Graphene has two hexagonal sublattices (A and B) in its honeycomb structure with translation vectors a₁ and a₂. (B) The Brillouin zone is hexagonal with two symmetry inequivalent corners labeled K and K′. (C) Near the Dirac points at K and K′, the dispersion is linear. The band structure is computed at the PBE level and taken from the Computational 2D Materials Database (Haastrup et al., 2018) with Fermi energy set to zero.

Figure 33

Top and side views of (A) hexagonal boron nitride, (B) hydrogenated silicene (silicane), (C) phosphorene, (D) and 2H-MoS₂. Structures taken from the Computational 2D Materials Database (Haastrup et al., 2018).

Figure 34

Band structure of MoS₂ in the 2H phase at the PBE (black) and G₀W₀ (red) levels. The G₀W₀ bands include spin-orbit coupling but the PBE bands do not. The Fermi energies for each case are indicated by horizontal dotted lines. Data taken from the Computational 2D Materials Database (Haastrup et al., 2018).

Figure 35

Ionization spectrum of pyridine. G₀W₀@PBE0 QP energies compared to the experimental photoemission spectrum (Liu et al., 2011). The calculated spectrum has been artificially broadened; the position of the QP energies is indicated with vertical bars. All QP energies are extrapolated using the cc-pVnZ (n=3–6) basis sets, see Appendix C for further computational details. The QPs of the first valence states are colored in red, green and blue.

Figure 36

GW100 benchmark comparing IP_HOMO energies computed at the G₀W₀@PBE level. FHI-aims is set as reference: ΔIP_HOMO = IP_HOMO(FHI-aims)−IP_HOMO(X). (A) Comparison of extrapolated/converged results for VASP (Maggio et al., 2017), WEST (Govoni and Galli, 2018), BerkeleyGW (van Setten et al., 2015). Shown are the results from full-frequency treatments and iterative solutions of the QP equation. (B) Comparison of localized basis set codes using the Gaussian basis set def2-QZVP (Weigend and Ahlrichs, 2005) for Turbomole (no-RI) (van Setten et al., 2015) and the N³ implementation in CP2K (Wilhelm et al., 2018). Note that BN, O₃, MgO, BeO, and CuCN are excluded for WEST, VASP and CP2K and that the BerkeleyGW and Turbomole data contain only a subset of 19 and 70 molecules, respectively. Box plot: Outliers represented by dots; boxes indicate the “interquartile range” measuring where the bulk of the data are.

Figure 37

Fundamental gaps of gas-phase benzene and band gap of the benzene crystal (space group Pbca). PBE was used as starting point for the G₀W₀ calculations. Data retrieved from Refaely-Abramson et al. (2013).

Figure 38

Density difference for the CO molecule between Hartree-Fock (HF) and PBE (left), coupled cluster singles-doubles (CCSD) and self-consistent GW (right). Charge depletion in the three methods is encoded by blue and charge accumulation by red colors. The same computational settings as in Caruso et al. (2013a) have been used.

Figure 39

Total energy of atoms computed with three different GW variants for atoms and small molecules plotted as a difference to the essentially exact Configuration Interaction (CI) results. Data retrieved from Caruso et al. (2012a).

Figure 40

Diagrammatic representation of Hedin's equations. All 5 quantities are coupled to all others. Here, we omit the Hartree potential from the G diagram, though it must also be included when translating the diagrams to the equations in Appendix B.

Figure 41

(A) Schematic representation of optical absorption. (B) Diagrammatic representation of GW/BSE. The electron and hole are represented by G lines, computed in the GW approximation, and their direct interaction is through the screened Coulomb interaction.

Figure 42

Schematic representation of diagrams included with (A) GW, (B) particle-hole T-matrix, and (C) particle-particle T-matrix. In each case, channels going into the box are correlated further with additional interactions.

Figure 43

Spectral function A(ω) of graphene on SiC at the Dirac point for electron doping density of n = −5.9 × 10¹³ cm⁻². Data taken from Lischner et al. (2013).