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.

C ore-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 corelevel 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 coupledcluster (Δ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 oneparticle 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 quantummechanical 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 (G 0 W 0 ) and starts from a set of meanfield single-particle orbitals {ψ n } and corresponding eigenvalues {ε n }. These are usually obtained from a preceding KS-DFT or Hartree−Fock (HF) calculation. The GW quasiparticle (QP) energies ε n G W 0 0 are computed by iteratively solving where v xc is the XC potential from DFT and spin variables are omitted. In the following, we use the notation Σ n = ⟨ψ n |Σ|ψ n ⟩ and v n xc = ⟨ψ n |v xc |ψ n ⟩ 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 ε = − BE n n G W 0 0 , and the selfenergy Σ is given by  where η is a positive infinitesimal. The self-energy is typically split into a correlation Σ c and an exchange part where Σ c is computed from W 0 c = W 0 − v and Σ x is computed from the bare Coulomb interaction v. The mean-field Green's function G 0 is given by where ε F denotes the Fermi energy. W 0 in eq 2 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 2) 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 1a. For such complex pole structures, the analytic continuation that is frequently employed in GW calculations for valence states to continue Σ c 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 2. 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 corelevel 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 corelevel excitations from GW. In the following, the failure of standard G 0 W 0 schemes for core states is explored. Figure 1a shows the G 0 W 0 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 (G 0 W 0 @PBE). Instead of iterating eq 1, we can obtain its solution graphically by finding the intersections of the straight line ω − ε n + v n xc − Σ n x with the self-energy matrix elements Σ n c . As apparent from Figure 1a, a clear single solution is missing. Many intersections are observed, which are all valid solutions of eq 1.
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 4, 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  1c. We observe many peaks with similar spectral weights. No distinct peak can be assigned to the QP excitation. In other words, G 0 W 0 @PBE does not provide a unique QP solution for the 1s excitation. This is in sharp contrast with the valence case, where G 0 W 0 @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 G 0 W 0 @PBE level without exception. Figure 1c 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 G 0 W 0 @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 {ε m } in the Green's function with the G 0 W 0 QP energies, re-evaluate eq 1, and iterate until G is self-consistent in the eigenvalues. For most valence and virtual states, a unique QP solution exists at the G 0 W 0 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 evGW 0 @PBE, where W is kept fixed at the W 0 level. Iterating the eigenvalues in G shifts the onset of the pole structure of the self-energy to lower energies; see Figure 1a. The pole structure of the self-energy looks similar in G 0 W 0 @PBE and evGW 0 @ PBE but is shifted by a constant amount. Figure 1b shows that the G 0 W 0 @PBE self-energy is indeed almost identical to that of evGW 0 @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 1c. In other words, eigenvalue selfconsistency 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 Σ n c 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 4, large imaginary parts correlate with low spectral weights, that is, satellite character. Rewriting the selfenergy into analytic form reveals its pole structure where Ω s are charge-neutral excitations and P s are transition amplitudes. 33 The G 0 W 0 self-energy therefore has poles at ε i − Ω s and ε a + Ω s , 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 G 0 W 0 @PBE, ε i are PBE eigenvalues and Ω s are close to PBE eigenvalue differences between occupied and virtual states. The neutral excitations Ω s 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 − Ω s poles in the self-energy are therefore considerably too high in energy and start to energetically overlap with the QP energy of the The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter core state. This explains why the satellites have such high spectral weight in G 0 W 0 @PBE and why no distinct QP peak can be found.
In evGW 0 , we replace the KS-DFT eigenvalues {ε m } in eq 3 by ε m + Δε m , where Δε m is the GW correction. For a PBE starting point, Δε m is negative for occupied states, and the poles of Σ n c shift to lower energies, away from the QP energy. The poles in the core region are now located at ε 1s + Δε 1s − Ω s 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 evGW 0 leads to consistent improvements over G 0 W 0 . 53, 57 We therefore iterate the eigenvalues only in G.
The effect of self-consistency in G can be mimicked in a G 0 W 0 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 evGW 0 @ PBE. For α = 0.45, we obtain approximately the same shift of the pole structure as in evGW 0 @PBE and observe a distinct QP peak at the same frequency; see Figure 1d,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 1f, G 0 W 0 @PBE0 still suffers from a large transfer of spectral weight to the satellites.
Previous GW core-level studies for small molecules 31,41 reported G 0 W 0 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 1c. Linearizing the QP equation by a Taylor expansion to first-order around ε n , 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 1 should always be solved iteratively for core states.
We now assess the accuracy of evGW 0 @PBE and G 0 W 0 @ 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 2a. 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.
For evGW 0 @PBE, the mean absolute error (MAE) of the absolute BEs with respect to experiment is reported in Figure  2a 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 2b shows the MAE at the G 0 W 0 @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 The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter agrees nicely with our analysis of the self-energy in Figure 1d, where we found that α ≈ 0.45 reproduces the evGW 0 selfenergy 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 2b. Figure 3 compares the absolute BEs obtained from experiment to the theoretical BEs computed at the ΔSCF, evGW 0 @PBE, and G 0 W 0 @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.
The evGW 0 @PBE approach consistently yields excellent agreement of the absolute BEs with experiment for all data points, as shown in Figure 3. With an overall MAE of 0.3 eV, the accuracy of evGW 0 @PBE is well within the chemical resolution required for the interpretation of most XPS spectra. G 0 W 0 @PBEh(α = 0.45) yields a similar overall MAE, which, however, depends to some extent on the species; see Table 1. As shown in Figure 3d, F1s removal energies are systematically underestimated with G 0 W 0 @PBEh. Results for this element could, in principle, be improved by using an element-specific optimized α value based on the analysis in Figure 2b.
Relative BEs are very well reproduced with all three theoretical methods, as shown in Table 1 and in more detail in Figure S1. With ΔSCF and evGW 0 @PBE, we obtain MAEs <0.2 eV and slightly larger errors between 0.2 and 0.3 eV with G 0 W 0 @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 G 0 W 0 setups routinely used for valence excitations cannot be employed. For core states, G 0 W 0 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 evGW 0 can be reproduced in G 0 W 0 , 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   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 evGW 0 , 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 selfconsistent evGW 0 calculations, we use the PBE functional as a starting point, and for G 0 W 0 @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 R 2 is in most cases >0.9; see Tables S2 and S3. Alternatively, the extrapolation can be performed with respect to C n −3 , where C n 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 G 0 W 0 @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.