Accurate Absolute and Relative Core-Level Binding Energies from GW.

We present an accurate approach to compute X-ray photoelectron spectra based on the GW Green's function method that overcomes the shortcomings of common density functional theory approaches. GW has become a popular tool to compute valence excitations for a wide range of materials. However, core-level spectroscopy is thus far almost uncharted in GW. We show that single-shot perturbation calculations in the G0W0 approximation, which are routinely used for valence states, cannot be applied for core levels and suffer from an extreme, erroneous transfer of spectral weight to the satellite spectrum. The correct behavior can be restored by partial self-consistent GW schemes or by using hybrid functionals with almost 50% of exact exchange as a starting point for G0W0. We also include relativistic corrections and present a benchmark study for 65 molecular 1s excitations. Our absolute and relative GW core-level binding energies agree within 0.3 and 0.2 eV with experiment, respectively.

Core-level spectroscopy techniques, such as X-ray photoelectron spectroscopy (XPS), are important tools for chemical analysis and can be applied to a broad range of systems including crystalline[1,2] and amorphous materials,[3−5] liquids,[6−8] adsorbates at surfaces,[9] or 2D materials.[10,11] XPS measures core-level binding energies (BEs), which are element-specific but depend on the local chemical environment. For complex materials, the assignment of the experimental XPS signals to the specific atomic sites is notoriously difficult due to overlapping spectral features or the lack of well-defined reference data.[4] Accurate theoretical tools for the prediction of core excitations are therefore important to guide the experiment. Calculated relative BEs, that is, BE shifts with respect to a reference XPS signal, are particularly useful for the interpretation of experimental spectra. However, the prediction of accurate absolute core-level energies is equally important, in particular, when reference core-level energies are not available.

The most common approach to compute core-level BEs is the Delta self-consistent field (ΔSCF) method, which is based on Kohn–Sham density functional theory (KS-DFT). In ΔSCF, the core-level BEs are calculated as the total energy difference between the neutral and the ionized systems.[12] Relative core-level BEs from ΔSCF generally compare well to experiment. For small molecules, deviations typically lie in the range of 0.2 to 0.3 eV,[13] which is well within or close to the chemical resolution required for most elements. The dependence of the relative BEs on the exchange-correlation (XC) functional is almost negligible for small systems,[13] but can be more severe for complex materials.[11,14] Absolute ΔSCF BEs can differ by several electronvolts from the experimental data. This deviation is quite sensitive to the XC functional.[15] The best results for absolute core excitations have been obtained using the TPSS[16] and SCAN[17] meta-generalized gradient approximations. The reported mean absolute deviations from the experiment lie in the range of ∼0.2 eV for benchmark sets of small molecules. For medium-sized to large molecules, however, the accuracy of ΔSCF can quickly reduce by an order of magnitude for absolute BEs.[18] This behavior can be partly attributed to an insufficient localization of the core hole in the calculation for the ionized system. Constraining the core hole in a particular state can be difficult, and variational instabilities are not uncommon.[19]

Most importantly, ΔSCF cannot be applied without further approximations to periodic systems, such as surfaces, where the ionized calculation would lead to a Coulomb divergence.[20] Such divergences can be circumvented by using cluster models,[17,21] by neutralizing the unit or supercell with compensating background charges,[22] or by adding the compensating electrons to the conduction band.[23−25] However, these approximations can obscure the calculations and even lead to qualitatively wrong results, as recently demonstrated for oxide surfaces.[26]

Higher level theoretical methods such as Delta coupled-cluster (ΔCC) approaches yield highly accurate relative and absolute core ionization energies.[27−29] ΔCC also requires the computation of a core-ionized system, leading to the same conceptual problems as in ΔSCF. Response theories, for example, equation-of-motion coupled cluster, avoid these problems but deviate by several electronvolts from experiment and require at least triples contributions for quantitative agreement.[30] Good accuracy for deep states was reported for a recently introduced direct approach based on effective one-particle energies from the generalized KS random phase approximation.[31] However, the application of these higher level methods is restricted to small- or medium-sized systems due to unfavorable scaling with the system size and large computational prefactors.

The GW method[32] expresses the response of the system to a core hole through perturbation theory and also avoids, like other response theories, an explicit optimization of the wave function of a core-ionized system. GW is thus a promising method to improve upon the limitations of traditional Δ approaches and has become a widespread tool for the accurate prediction of electron removal energies of valence states in molecular and solid-state systems.[33] It is routinely applied to systems with several hundred atoms[34−36] and recently even to system sizes with >1000 atoms.[37,38] However, the computation of deep core levels, as measured by XPS, has been rarely attempted with GW. The first exploratory studies for solid-state systems obtained partly promising results.[39,40] The few existing studies for molecular core excitations also give a mixed first impression, as anything between a 0.5[18] and 10 eV deviation from experiment has been reported.[31,41] In this work, we show how reliable and highly accurate core-level BEs can be obtained from GW and explain why large deviations from experiment were previously reported. We also present a GW benchmark set for 1s core states complementary to the popular GW100 benchmark set[42] for valence excitations.

First, we introduce the GW framework. The central object of GW is the self-energy Σ, which contains all quantum-mechanical exchange and correlation interactions of the hole created by the excitation process and its surrounding electrons. The self-energy is calculated from the Green’s function G and the perturbation expansion in the screened Coulomb interaction W, as formulated by Hedin in the 1960s.[32] The poles of G directly correspond to the excitation energies, as measured in photoemission spectroscopy.

In practice, GW is performed within the first-order perturbation theory (G0W0) and starts from a set of mean-field single-particle orbitals {ψ} and corresponding eigenvalues {ε}. These are usually obtained from a preceding KS-DFT or Hartree–Fock (HF) calculation. The GW quasiparticle (QP) energies are computed by iteratively solvingfor , where vxc is the XC potential from DFT and spin variables are omitted. In the following, we use the notation Σ = ⟨ψ|Σ|ψ⟩ and vxc = ⟨ψ|vxc|ψ⟩ for the (n, n) diagonal matrix elements of the self-energy and XC potential. The QP energies are related to the BE of state n by , and the self-energy Σ is given bywhere η is a positive infinitesimal. The self-energy is typically split into a correlation Σ and an exchange part Σ, Σ = Σ + Σ, where Σ is computed from W0 = W0 – v and Σ is computed from the bare Coulomb interaction v. The mean-field Green’s function G0 is given bywhere εF denotes the Fermi energy. W0 in eq is the screened Coulomb interaction in the random-phase approximation (RPA) and is computed from the dielectric function, as described in ref (33).

We now discuss the application of GW to core-level spectroscopy. The basic requirement to obtain computational XPS data from GW is an explicit description of the core electrons. We treat the latter efficiently by working in a local all-electron basis of numeric atom-centered orbitals (NAOs). Furthermore, we showed that highly accurate frequency integration techniques for the computation of the self-energy (eq ) are required for core states.[18] Unlike for valence states, the self-energy has a complicated structure with many poles in the core region, as displayed in Figure a. For such complex pole structures, the analytic continuation that is frequently employed in GW calculations for valence states to continue Σ from the imaginary to the real frequency axis fails completely.[18] We showed that the contour deformation (CD) technique, in which a full-frequency integration on the real frequencies axis is performed, yields the required accuracy. Results from CD exactly match the computationally demanding fully analytic solution of eq .[18] Our CD-GW implementation is computationally efficient, enabling the computation of system sizes exceeding 100 atoms; see ref (18) for details of our GW core-level implementation in the all-electron code FHI-aims.[43] Numerically stable and precise algorithms for the computation of the self-energy are only the first step toward reliable core-level excitations from GW. In the following, the failure of standard G0W0 schemes for core states is explored.

Figure 1

Core excitation for a single water molecule from G0W0 and evGW0. (a) Real part of the correlation self-energy Σ(ω) using the PBE functional as a starting point. Diagonal matrix elements for the oxygen 1s orbital. (b) Self-energy from G0W0@PBE shifted by Δev relative to the evGW0 result. The intersection with the red dashed line is the graphical solution of the QP equation (eq ). (c) Spectral function A(ω) as defined in eq from G0W0@PBE and evG0W0@PBE. (d) Self-energy and (e) spectral function using PBEh(α = 0.45) as a starting point. (f) Spectral function obtained from G0W0@PBEh for different amounts of exact exchange α. The vertical gray dashed line indicates the QP solution from evGW0@PBE. The self-energy is slightly broadened for better visualization, and each sharp peak actually corresponds to a pole.

Figure a shows the G0W0 self-energy matrix elements for the O1s state of an isolated water molecule using the Perdew–Burke–Ernzerhof (PBE)[44] functional for the underlying DFT calculation (G0W0@PBE). Instead of iterating eq , we can obtain its solution graphically by finding the intersections of the straight line ω – ε + vxc – Σ with the self-energy matrix elements Σ. As apparent from Figure a, a clear single solution is missing. Many intersections are observed, which are all valid solutions of eq .

To further investigate this multisolution behavior, we calculate the spectral function A(ω)[18,33]where m runs over all occupied and virtual states. In eq , we include also the imaginary part of the complex self-energy and use, unlike in fully self-consistent GW,[45,46] only the diagonal matrix elements of Σ. The spectral function for the oxygen 1s excitation of an isolated water molecule is reported in Figure c. We observe many peaks with similar spectral weights. No distinct peak can be assigned to the QP excitation. In other words, G0W0@PBE does not provide a unique QP solution for the 1s excitation. This is in sharp contrast with the valence case, where G0W0@PBE is routinely applied to molecules and a clear single solution has been reported in the vast majority of cases.[42] Multiple solutions have been reported for a couple of systems.[42,47−49] However, typically not more than two possible solutions were found, whereas a multitude of peaks are observed in the spectral function of molecular 1s states at the G0W0@PBE level without exception.

Figure c illustrates that for core states, the QP energy and the satellite spectrum have merged. Satellites are, for example, due to multielectron excitations such as shake-up processes[50,51] and have typically much smaller spectral weights than the QP peak. The fact that we observe the opposite for G0W0@PBE implies that almost all spectral weight has been transferred from the QP peak to the satellites. We will next investigate the origin of this behavior and provide a solution using a strategy previously also employed for 3d states in transition-metal oxides[52,53] or semicore states.[54]

We start by updating the KS eigenvalues {ε} in the Green’s function with the G0W0 QP energies, re-evaluate eq , and iterate until G is self-consistent in the eigenvalues. For most valence and virtual states, a unique QP solution exists at the G0W0 level, whereas for core states, we initialize the iteration in G with an approximation of the QP energy. This procedure yields a partially eigenvalue self-consistent scheme denoted as evGW0@PBE, where W is kept fixed at the W0 level. Iterating the eigenvalues in G shifts the onset of the pole structure of the self-energy to lower energies; see Figure a. The pole structure of the self-energy looks similar in G0W0@PBE and evGW0@PBE but is shifted by a constant amount. Figure b shows that the G0W0@PBE self-energy is indeed almost identical to that of evGW0@PBE when shifted by Δev = −28.7 eV. The effect of this shift is that the graphical solution now produces a clear QP solution and a satellite spectrum with much lower intensity, as displayed in Figure c. In other words, eigenvalue self-consistency in G achieves a separation of the QP peak and the satellite spectrum for deep core states, which can be understood as follows.

Satellites occur in frequency regions, where the real part of Σ has poles and its imaginary part has complementary peaks, which is shown in detail in our recent GW review article.[33] As obvious from eq , large imaginary parts correlate with low spectral weights, that is, satellite character. Rewriting the self-energy into analytic form reveals its pole structurewhere Ω are charge-neutral excitations and P are transition amplitudes.[33] The G0W0 self-energy therefore has poles at ε – Ω and ε + Ω, where i indicates occupied states and a indicates virtual states. Each of these poles gives rise to satellite features and can be understood as an electron or hole excitation coupled to a neutral excitation.

For G0W0@PBE, ε are PBE eigenvalues and Ω are close to PBE eigenvalue differences between occupied and virtual states. The neutral excitations Ω are typically underestimated at the PBE level, whereas the eigenvalues are overestimated by several electronvolts in the valence region of the spectrum and by 20–30 eV for the 1s core states.[18,41] The ε1s – Ω poles in the self-energy are therefore considerably too high in energy and start to energetically overlap with the QP energy of the core state. This explains why the satellites have such high spectral weight in G0W0@PBE and why no distinct QP peak can be found.

In evGW0, we replace the KS-DFT eigenvalues {ε} in eq by ε + Δε, where Δε is the GW correction. For a PBE starting point, Δε is negative for occupied states, and the poles of Σ shift to lower energies, away from the QP energy. The poles in the core region are now located at ε1s + Δε1s – Ω and the corresponding satellite peaks are separated from the QP peak and reduced in spectral weight. Note that for a quantitative comparison of the satellite positions to experiment the cumulant expansion is needed, which includes vertex corrections beyond the GW self-energy.[54−56] We focus in this work on the QP energies.

Inserting the QP energies also in W (evGW) leads to an underscreening, if not compensated by missing vertex corrections. This is in line with the observation that band gaps are considerably overestimated[57] in evGW and that evGW spectra are overly stretched,[58] whereas evGW0 leads to consistent improvements over G0W0.[53,57] We therefore iterate the eigenvalues only in G.

The effect of self-consistency in G can be mimicked in a G0W0 calculation, which is computationally less demanding, by including exact exchange in the DFT functional. We employ the PBE-based hybrid (PBEh) functional family,[59] which is characterized by an adjustable fraction α of HF exchange and corresponds for α = 0.25 to the PBE0[60,61] functional. Here we optimize α to reproduce the self-energy structure of evGW0@PBE. For α = 0.45, we obtain approximately the same shift of the pole structure as in evGW0@PBE and observe a distinct QP peak at the same frequency; see Figure d,e. Increasing the amount of exact exchange, the QP peak in the spectral function moves to lower energies, which is in agreement with the starting point optimization studies conducted for valence excitations.[33,59,62] However, distinct QP peaks are obtained only for α > 0.3. As shown Figure f, G0W0@PBE0 still suffers from a large transfer of spectral weight to the satellites.

Previous GW core-level studies for small molecules[31,41] reported G0W0 calculations performed on top of generalized gradient approximation (GGA) or hybrid functionals with a low amount of exact exchange. Our analysis presented in this Letter demonstrates that those studies cannot have found the QP solution because their spectral function would look like the yellow spectrum in Figure c. Linearizing the QP equation by a Taylor expansion to first-order around ε, as done for such a spectral function in refs (31), (41), and (63), leads to uncontrollable results, which partly explains the large deviation of the reported results from experiment.[31,41] Furthermore, the linearization error rapidly increases with increasing BE and may already amount to 0.5 eV for deeper valence states, as shown in ref (33). As already pointed out in our previous work,[18]eq should always be solved iteratively for core states.

We now assess the accuracy of evGW0@PBE and G0W0@PBEh with respect to experiment for a benchmark set of 65 1s BEs of gas-phase molecules, denoted in the following as CORE65. This benchmark set contains 30 C1s, 21 O1s, 11 N1s, and 3 F1s excitations from 32 small, inorganic and organic molecules up to 14 atoms; see Table S1 for details. The CORE65 benchmark covers a variety of different chemical environments and bonding types and the most common functional groups. As with all correlated electronic structure methods, GW converges slowly with respect to the basis set size.[33] Even at the quadruple-ζ level, the BEs deviate by 0.2 to 0.4 eV from the complete basis set limit. (See Tables S2 and S3.) All GW results are thus extrapolated to the complete basis set limit using the Dunning basis set family cc-pVnZ (n = 3–6).[64,65]

Because we expect relativistic effects to become important for heavier elements, we add relativistic corrections for the 1s excitations as a postprocessing step to the GW calculation. Our relativistic corrections have been obtained by solving the radial KS and four-component Dirac-KS equations self-consistently for a free neutral atom at the PBE level and evaluating the difference between their 1s eigenvalues; see Figure a. Details of our relativistic correction scheme, which is similar to the one reported in ref (16), and its comparison to other relativistic methods will be described in a forthcoming paper.

Figure 2

Mean absolute error (MAE) of the absolute BEs for the CORE65 benchmark set with respect to experiment. (a) evGW0@PBE with and without relativistic correction. (b) MAE for G0W0@PBEh dependent on the amount of exact exchange α in the PBEh functional. Relativistic corrections are included.

For evGW0@PBE, the mean absolute error (MAE) of the absolute BEs with respect to experiment is reported in Figure a for relativistic and nonrelativistic calculations. We find that, already for second-row elements, relativistic effects start to dominate the error in the QP energies. In the nonrelativistic case, the core-level BEs are generally underestimated (see Table S2), and the MAE increases with the atomic number. Accounting for relativistic effects, the species dependence in the MAE is largely eliminated.

Figure b shows the MAE at the G0W0@PBEh(α) level with respect to the amount of exact exchange α in the PBEh functional, including relativistic corrections. These α-dependent calculations are performed for a subset of 43 excitations of the CORE65 benchmark set, for which the mapping between the core state and the atom is trivial and requires no analysis of, for example, molecular orbital coefficients. The smallest MAE is obtained for α values around 0.45. This observation agrees nicely with our analysis of the self-energy in Figure d, where we found that α ≈ 0.45 reproduces the evGW0 self-energy best. For smaller α values, the BEs are underestimated, and for larger values, the BEs are increasingly overestimated. The species dependence of the optimal α values are mostly reduced when taking relativistic effects into account. The optimal α values increase only slightly with the atomic number ranging from 0.44 to 0.49; see Figure b.

Figure compares the absolute BEs obtained from experiment to the theoretical BEs computed at the ΔSCF, evGW0@PBE, and G0W0@PBEh(α = 0.45) levels. The ΔSCF calculations are performed with the PBE0[60,61] functional and are carefully converged by adding additional tight basis functions to standard Gaussian basis sets for the core-hole calculation.[66] Following a recently proposed ΔSCF simulation protocol,[17] we include scalar relativistic effects self-consistently via the zeroth-order regular approximation (ZORA).[67] Our BEs obtained from ΔSCF-PBE0, with an overall MAE of 0.3 eV, agree much better with experiment than those reported in previous studies (0.7 eV),[15] which must be attributed to incomplete basis sets and the neglect of relativistic effects.

Figure 3

Absolute C1s (a), N1s (b), O1s (c), and F1s (d) binding energies (BEs) for the CORE65 benchmark set comparing calculated values at the ΔSCF, evGW0@PBE, and G0W0@PBEh(α = 0.45) levels to experiment. The respective computational method underestimates the BE when the data point is below the black line and overestimates the BE when the data point is above the black line.

The evGW0@PBE approach consistently yields excellent agreement of the absolute BEs with experiment for all data points, as shown in Figure . With an overall MAE of 0.3 eV, the accuracy of evGW0@PBE is well within the chemical resolution required for the interpretation of most XPS spectra. G0W0@PBEh(α = 0.45) yields a similar overall MAE, which, however, depends to some extent on the species; see Table . As shown in Figure d, F1s removal energies are systematically underestimated with G0W0@PBEh. Results for this element could, in principle, be improved by using an element-specific optimized α value based on the analysis in Figure b.

Table 1

Mean Absolute Error (MAE) in eV with Respect to Experiment for the CORE65 Benchmark Seta

	ΔSCF		evGW0@PBE		G0W0@PBEh(α = 0.45)
core-level	absolute BEs	relative BEs	absolute BEs	relative BEs	absolute BEs	relative BEs
all	0.33	0.14	0.30	0.18	0.33	0.26
C1s	0.36	0.10	0.27	0.18	0.24	0.29
N1s	0.32	0.08	0.30	0.14	0.16	0.23
O1s	0.32	0.22	0.32	0.22	0.48	0.25
F1s	0.12	0.13	0.44	0.05	0.83	0.11

MAE for absolute and relative core-level BEs, where the latter is the shift of the BE with respect to a reference molecule. CH4, NH3, H2O, and CH3F have been used as reference molecules for C1s, N1s, O1s, and F1s, respectively. Relativistic effects are accounted for in all three methods.

Relative BEs are very well reproduced with all three theoretical methods, as shown in Table and in more detail in Figure S1. With ΔSCF and evGW0@PBE, we obtain MAEs <0.2 eV and slightly larger errors between 0.2 and 0.3 eV with G0W0@PBEh(α = 0.45). Results for F1s are reported for the sake of completeness. However, note that we have only two data points for the relative BEs, and the experimental uncertainties are generally larger for fluorine than for the lighter elements. Except for F1s BEs, the MAEs are not species-dependent.

In summary, we showed that GW is a reliable and accurate method to calculate 1s core excitations. However, standard G0W0 setups routinely used for valence excitations cannot be employed. For core states, G0W0 calculations starting from GGA or standard hybrid functionals experience a huge weight transfer from the QP peak to the satellites. In fact, this weight transfer is so extreme that a unique QP solution does not exist for the molecules we have investigated. We demonstrated for a PBE starting point that eigenvalue self-consistency in G is mandatory to achieve a proper separation between QP and satellite peaks in the GW calculation. The effects of evGW0 can be reproduced in G0W0, which is computationally less expensive, by using a hybrid functional with a high fraction of exact exchange as a starting point. We found that 45% of HF exchange is optimal. Furthermore, the inclusion of relativistic effects and a proper extrapolation to the complete basis limit are crucial to obtain accurate core-level BEs. Our work is an important stepping stone for the accurate calculation of XPS spectra of condensed systems, where Δ-based approaches face conceptual limitations. GW can be applied without restrictions to systems with periodic boundary conditions and is also for large molecular structures a reliable and numerically robust method. Furthermore, this work is fundamental for the calculation of X-ray absorption spectra (XAS) from the Bethe–Salpeter equation,[68] which uses the GW results as input.

Computational Details

All calculations are performed with the FHI-aims program package,[43,69,70] where the all-electron KS equations are solved in the NAO scheme. The structures of the CORE65 molecules have been optimized at the DFT level using NAOs of tier 2 quality[43] to represent core and valence electrons. The PBE functional[44] is used to model exchange and correlation in combination with the atomic ZORA[43,67] kinetic energy operator. Van der Waals interactions are accounted for by employing the Tkatchenko–Scheffler dispersion correction.[71]

Core-level BEs from ΔSCF are calculated using the PBE0[60,61] hybrid functional, (atomic) ZORA, and def2 quadruple-ζ valence plus polarization (def2-QZVP)[72] basis sets. The def2-QZVP basis sets are all-electron basis sets of contracted Gaussian orbitals, which are optimized to yield accurate total energies.[72] Gaussian basis sets can be considered as a special case of an NAO and are treated numerically in FHI-aims. To guarantee the full relaxation of other electrons in the presence of the core hole, we decontracted the def2-QZVP basis sets in the ΔSCF calculation to add tighter core functions. To properly localize the core hole at a specific atom, a Boys localization[73] has been performed at the end of the SCF cycle of the charge-neutral calculation. The localized wave function is used as an initial guess for the charged system, and the core electron is removed from the localized eigenstates.

For the GW calculations, we use the CD technique[18,33,34,74] to evaluate the frequency integral of the self-energy and employ a modified Gauss–Legendre grid[70] with 200 grid points for the imaginary frequency integral. The QP equation is always solved iteratively. For evGW0, we additionally iterate the QP energies in G including explicitly all occupied states and the first five virtual states in the iteration. Scissor shifts are employed for the remaining virtual states. For the partially self-consistent evGW0 calculations, we use the PBE functional as a starting point, and for G0W0@PBEh(α) calculations, the PBEh(α) hybrid functionals are used.[59] The core-level BEs are extrapolated to the complete basis set limit using the Dunning basis set family cc-pVnZ (n = 3–6),[64,65] which is the standard choice for correlated electronic-structure methods. The extrapolation has been performed with four points by a linear regression against the inverse of the total number of basis functions. The standard error of the extrapolation is <0.1 eV, and the correlation coefficient R2 is in most cases >0.9; see Tables S2 and S3. Alternatively, the extrapolation can be performed with respect to C–3, where C is the cardinal number of the basis set. We found that the difference between both extrapolation schemes is very small; for example, the average absolute deviation is only 0.04 eV for the CORE65 G0W0@PBEh(α = 0.45) data. Self-energy matrix elements and spectral functions are calculated at the cc-pV4Z level. The relativistic corrections for the GW energies are computed at the PBE level from free neutral atom calculations on numerical real space grids.[75] Note that the GW approximation is based on a propagator approach, avoiding the explicit generation of a core hole. A localization procedure as in ΔSCF is therefore not required.