Relativistic correction scheme for core-level binding energies from $GW$

We present a relativistic correction scheme to improve the accuracy of 1s core-level binding energies calculated from Green's function theory in the $GW$ approximation, which does not add computational overhead. An element-specific corrective term is derived as the difference between the 1s eigenvalues obtained from the self-consistent solutions to the non- or scalar-relativistic Kohn-Sham equations and the four-component Dirac-Kohn-Sham equations for a free neutral atom. We examine the dependence of this corrective term on the molecular environment and on the amount of exact exchange in hybrid exchange-correlation functionals. This corrective term is then added as a perturbation to the quasiparticle energies from partially self-consistent and single-shot $GW$ calculations. We show that this element-specific relativistic correction, when applied to a previously reported benchmark set of 65 core-state excitations [J. Phys. Chem. Lett. 11, 1840 (2020)], reduces the mean absolute error (MAE) with respect to experiment from 0.55 to 0.30 eV and eliminates the species dependence of the MAE, which otherwise increases with the atomic number. The relativistic corrections also reduce the species dependence for the optimal amount of exact exchange in the hybrid functional used as starting point for the single-shot $G_0W_0$ calculations. Our correction scheme can be transferred to other methods, which we demonstrate for the Delta self-consistent field ($\Delta$SCF) approach based on density functional theory.


I. INTRODUCTION
Core-level binding energies (BEs), measured by X-ray photoemission spectroscopy (XPS), are element-specific, but depend also on the local chemical environment and thus afford access to information about the chemical bonding, oxidation state, and coordination of a given element in a sample. [1][2][3] The energetic differences (chemical shifts) between atomic species of the same type can be smaller than 0.5 eV for second-row elements and can be as low as 0.1 eV for carbon 1s excitations. 4 The interpretation of an XPS spectrum can be very difficult due to overlapping features or the lack of well-defined reference data. 5,6 Highly accurate theoretical tools for the prediction of relative and absolute BEs are therefore necessary to guide the experiment and its interpretation. The reliable computation of absolute core-level energies is generally more challenging than the calculation of energy shifts 7 and is the focus of this work.
The most common approach to calculating core-level BEs is the Delta self-consistent field (∆SCF) method, 8 wherein one computes the total energy difference between the ground and core-ionized state using Kohn-Sham density functional theory (KS-DFT). The best absolute core-level BEs have been obtained with meta-generalized gradient approximation (meta-GGA) functionals, yielding mean deviations of ≈0.2 eV with respect to experiment for small molecules. 9,10 Similar accuracy has been obtained with high-level wavefunction-based Delta coupled cluster methods, [11][12][13] albeit at much higher computational cost. The introduction of occupation constraints and the explicit generation of a charged system in ∆-based approaches leads to a plethora of problems. 14,15 Most importantly, the application to periodic systems rea) Electronic mail: Levi.Keller@aalto.fi. quires further approximations, e.g., neutralizing the unit or supercell. 16,17 These problems do not occur in response theories, which avoid the explicit introduction of a core-ionized system and recently emerged as viable alternatives to ∆-based schemes for core-level calculations. [18][19][20] One particularly promising response approach is the GW approximation 21 to many-body perturbation theory. GW offers access to the core-state quasiparticle energies, which can be directly related to the corelevel BEs. The GW approach is the standard approach to compute valence excitations (band structures) in solid-state physics. 22 With the advent of efficient implementations with localized basis-sets, [23][24][25][26][27] GW has also become the method of choice for calculating the BEs of frontier orbitals in molecules and is nowadays routinely applied to large molecular structures with several hundred atoms. [27][28][29] GW has also been successfully applied to deep valence and semi-core states with excitation energies up to 80 eV. [30][31][32] However, the application of GW for deep core states with BEs larger than 100 eV, which are the ones relevant for chemical analysis, has rarely been attempted. The first GW studies 19,[33][34][35] for deep core states reported partly large deviations of several electronvolts from experiment. Recently, we have shown that GW can also be successfully applied for 1s molecular core states when utilizing highly accurate techniques for the integration of the self-energy 36 and eigenvalue self-consistency in the Green's function. 20 The application of GW (or any other method) to absolute core-levels binding energies requires an accurate treatment of relativistic effects. Already for the 1s core levels of the pblock 2nd period elements, the magnitude of relativistic effects begins to exceed the size of the chemical shifts. Recently, we have applied the GW method to a benchmark set of 65 core-level BEs of the second-period elements C, N, O and F in small and medium-sized molecules (henceforth CORE65 benchmark set), obtaining a mean deviation from experiment of 0.3 eV. 20 Therein, we applied a simple relativistic correction to the GW computed core-level energies. The purpose of the present article is to describe this correction.
Relativistic effects enter the GW formalism via the underlying reference calculation. A fully-relativistic one-particle reference is described by a four-component Dirac spinor, or approximately by two-component spinors. Both types of spinors can describe noncollinear electronic states and thus spin-dependent electron-electron interactions. Explicitly spin-dependent GW equations using spinors as input were only developed 12 years ago, 37,38 and provide a framework to properly describe spin-orbit coupling (SOC) effects in GW . Several GW codes emerged over the past years that implement these equations, employing spinors from all-electron 2-component Dirac-KS-DFT calculations, 39,40 from KS-DFT with second-variation SOC 41 or from KS-DFT with fullyrelativistic pseudopotentials. 42,43 These implementations are in the following referred to as fully-relativistic GW approaches.
Fully-relativistic GW calculation have been primarily used to compute valence and conduction bands of solids with strong SOC effects. Fully-relativistic results have been reported for actinide, metals 44,45 transition metal chalcogenides 41,42 and dichalcogenides, 43,46,47 perovskites 42,48 as well as bismuth-based topological insulators. [49][50][51][52][53][54] Compared to non-relativistic GW , the fullyrelativistic approach is computationally at least four times more expensive. 40 Alternatively, the SOC can be added as a perturbative correction to the quasiparticle band structure, which was common in earlier GW studies. 55 Since this approach is computationally less expensive, it continues to be employed. [56][57][58] However, it has been shown that the SOC post-correction scheme can fail, for example, to describe the band inversion in topological insulators. 51 An explicit treatment of the SOC is typically not necessary for valence excitations of molecules, unless very heavy atoms are involved. Since molecular GW studies have mainly focused on organic semiconductors, 22 fully-relativistic GW calculations are rare and have been mostly conducted for diatomic molecules, 39,42 but recently also for transition metal halide complexes. 40 If relativistic effects are considered in molecular GW calculations at all, they are more commonly included by employing one-component reference states that capture only scalar-relativistic effects. This approach has been employed in GW calculations with scalar-relativistic pseudopotentials 23,27 or the zeroth-order regular approximation (ZORA). 29 Relativistic considerations for the innermost core-levels differ from those for valence excitations. SOC leads, for core states with angular quantum number l > 0, to a splitting. Already for 2p states of third-row elements, e.g. sulfur, this splitting lies in the range of 1 eV and increases for 4th-period elements to several tens of eV. 59 SOC affected core states require thus a noncollinear treatment and we are only aware of one GW study 39 that reports results for deep, spin-orbit-coupled p states. In this study we focus on organic molecules and small inorganic molecules with elements C, N, O and F, for which typically the 1s excitations in the energy range from 250 to 700 eV are measured in XPS. While scalar relativistic effects are heightened in the proximity of the poorly-screened nuclear charge, the SOC operator does not directly affect core states. In contrast to valence excitations, the most common scalar relativistic approximation, ZORA, performs poorly for the absolute eigenvalue associated with the innermost core levels. 60,61 In order to avoid the computational expense of a fully Dirac-like, two-component GW calculation, which is not essential to the physics of 1s core excitations, we derive here an element-specific relativistic corrective term for nonrelativistic and scalar-relativistic reference states, which we add in a post-GW perturbative step.
The remainder of this article describes this correction. In section II we describe the GW formalism, highlighting aspects that are particularly relevant for core-level calculations, and follow this up with an overview of the aspects of relativistic theory relevant to this work. We then describe the methods employed in section III. We present and discuss the results of our correction schemes, which we apply to GW and ∆SCF computed core-level BEs of the CORE65 benchmark set in section IV and finally draw conclusion in section V.

II. THEORY
A. G 0 W 0 quasiparticle energies In practice, GW is often performed as one-shot perturbative approach (G 0 W 0 ) on top of an underlying mean field theory calculation. Possible mean field theories are Hartree-Fock (HF), KS-DFT or hybrid DFT which yield the molecular orbitals (MOs) {φ n } and eigenvalues { n } used as input (starting point) for the G 0 W 0 calculation. The G 0 W 0 quasiparticle (QP) energies G 0 W 0 n are obtained by solving the QP equation and can be related to the BE of state n by BE n = − The nth diagonal elements of the KS exchange-correlation (XC) potential and self-energy operator Σ are denoted by v XC n = φ n | v XC |φ n and Σ n = φ n | Σ |φ n , respectively. Note that we have omitted spin indices. The self-energy operator is given by (2) where G 0 is the non-interacting KS Green's function, W 0 is the screened Coulomb interaction, and η is a positive infinitesimal. The KS Green's function is obtained from the KS orbitals and eigenvalues by where F is the Fermi energy, and the sum runs over both occupied and virtual KS orbitals. Within the random phase approximation, the screened Coulomb interaction is given by where the dielectric function ε is ε(r, r , ω) = δ (r, r ) − dr χ 0 (r , r , ω) |r − r| (5) and the irreducible polarizability χ 0 is given in the real-space Adler-Wiser representation 62,63 as The index i runs over occupied orbitals, and a runs over virtual orbitals.
Equation (1) is non-linear, and can be solved iteratively or approximately by linearization using a Taylor expansion to first order around n . 22 As we pointed out in our previous work, 20,36 the linearization error increases rapidly with increasing BE and can already amount to 0.5 eV for deeper valence states. 22 The magnitude of this error is in the range of the chemical shifts expected for 1s excitations. In addition, core-level BEs are an order of magnitude larger than deep valence BEs, potentially leading to even larger linearization errors. We therefore always solve the QP equation iteratively.

B. Frequency treatment for core-states
The accurate frequency integration of the self-energy (Equation (2)) is one of the major challenges for the calculation of deep core states. A common approach for valence states is to evaluate the self-energy for imaginary frequencies and analytically continue it to the real frequency axis by fitting the self-energy matrix elements to a multipole model. Analytic continuation is employed in many state-of-the-art GW implementations 24,27,28,64 and yields accurate results for valence states. 65 The structure of Σ n for a valence state is typically smooth in the frequency region where the QP solution is expected, and is well reproduced by analytic continuation. For core states, the self-energy has a complicated structure with many poles. We showed that analytic continuation becomes numerically unstable in the core region and completely fails to reproduces the self-energy structure. 36 For core states, a more accurate evaluation of the selfenergy on the real frequency axis is required. We employ the contour deformation (CD) technique, 22,23,36,66,67 where the numerically unstable integration along the real frequency axis is avoided by extending the integrand to the complex plane. The contours are chosen such that only the poles of G 0 are enclosed in the contours and the contour integral is evaluated using the residue theorem; see Refs. 22,36 for details. The integral along the real frequency axis is then evaluated as where θ is the Heaviside step function and i refers to occupied and a to unoccupied orbitals. The CD technique reproduces the self-energy structure for core excitations exactly and matches the results from the computationally more expensive fully-analytic solution of Equation (2). 36 Recently, combining the CD approach with analytic continuation of W 0 has been proposed as alternative approach, 68 but has not been tested for core excitations.
C. Restoration of the core-level quasiparticle peak In this section we briefly present the GW variants that we used in this work. By now many different GW flavours have emerged in practical calculations. 22 The most common flavor, G 0 W 0 based on a semi-local DFT starting point, breaks down for core states, as we will detail in the following. We therefore need to go beyond the most common approach.
The problem with the conventional GW approach is related to a loss of spectral weight in the quasiparticle peak in the spectral function A(ω) = 1/π ∑ m ImG m , where m runs over all occupied and virtual states and G = G 0 + G 0 ΣG. For molecular valence states, G 0 W 0 performed on top of a DFT calculation with the Perdew-Burke-Ernzerhof (PBE) 69 functional (G 0 W 0 @PBE) yields a clearly identifiable QP peak. This peak corresponds to a distinct solution of Equation (1). Multiple solutions, that would indicate spectral weight transfer to other peaks in the spectral function, have only been observed for a few systems for frontier orbitals. 35,65,[70][71][72] These are rare cases and usually not more than two possible solutions are observed.
The situation is dramatically different for deep core states. The analysis of the spectral functions in our recent work showed that a unique QP solution is not obtained with G 0 W 0 @PBE for 1s states. 20 The spectral function shows a multitude of peaks with similar spectral weight, but no distinct QP excitation. A spurious transfer of spectral weight from the QP peak to plasmon satellites has been previously observed for deep valence states of transition metal oxides 73,74 and semi-core excitations of sodium. 31 However, for deep core states the transfer of spectral weight is far more extreme resulting in a spectral function, where satellite spectrum and quasiparticle have completely merged. Such a spectral function contradicts the expected physics. Photoemission spectra of molecular 1s excitations show strong QP peaks accompanied by satellites features due to multi-electron excitations such as shake-up processes, which are orders of magnitudes smaller than the main excitation. 75,76 We showed that the 1s QP peak can be correctly restored by including eigenvalue self-consistency in G, while keep-ing PBE as starting point. 20 This scheme is referred to as evGW 0 @PBE. In evGW 0 , the screened Coulomb interaction is kept fixed at the W 0 level and the Green's function is recomputed replacing the mean-field eigenvalues with the QP energies from Equation (1). We enforce eigenvalue selfconsistency only in G. Inserting the QP energies also in W 0 (evGW ) reduces the screening, which is not advantageous because the overscreening in W at the PBE level compensates the underscreening due to missing vertex corrections. It has been shown that evGW 0 yields band gaps in good agreement with experiment, 77 while underscreening errors in the evGW scheme lead to too large band gaps 77 and overly stretched spectra. 78 .
Higher-level self-consistency schemes, such as fullyselfconsistent GW (scGW ) 79,80 or quasiparticle self-consistent GW (QSGW ) 81 , are expected to restore the 1s QP peak as well, but might not yield better agreement with experiment than evGW 0 . It has been shown that scGW overestimates molecular HOMO excitations 82 and band gaps in solids 83 . Similar underscreening effects are also expected for core states. A first exploratory study seems to confirm this assumption for QSGW , reporting an overestimation of 2 eV for 1s core states of small molecules. 35 evGW 0 is computationally more demanding than a G 0 W 0 calculation because the QP equation is not only solved for the 1s core states of interest, but repeatedly for all occupied and virtual states until convergence in G is reached. We showed that the core-level QP peak can be also restored in a G 0 W 0 calculation by using a XC functional with a high fraction of exact exchange as starting point. 20 We employ the PBEh(α)functional family with an adjustable amount α of HF exact exchange. 84 The XC energy E xc is given by where E EX x denotes the HF exchange energy. E PBE x and E PBE c are the PBE exchange and correlation energy, respectively. In this work, we followed Ref. 20 and used both evGW 0 and G 0 W 0 @PBEh(α). We analyze how the relativistic corrections we devised affect the two schemes.

D. Relativistic methods
It has long been recognized that relativistic effects play a large role in the chemistry of heavy elements. 85,86 In this work, we treat the core states of light elements carbon through fluorine, whose BEs are in the range of 250 -700 eV, and for which relativistic effects are usually smaller than 1 eV. However, the accuracy required to resolve XPS spectra of 1s excitations of 2nd row elements is in the range of some tenths of an electronvolt and therefore on the same order of magnitude as the relativistic effects.
Our relativistic correction scheme for GW is based on two different relativistic KS-DFT methods. The first is the 4-component Dirac Kohn-Sham (4c-DKS) approach, further also referred to as fully-relativistic scheme. The second uses the scalar-relativistic ZORA.

Fully relativistic Dirac approach
The relativistic description of a non-interacting electron in an external potential V is given by the Dirac equation 87 where c the speed of light, σ is the Pauli vector and p is the momentum operator. For electron-like states in the nonrelativistic limit (c → ∞), the large component φ reduces to the wave function of the Schrödinger equation, while the small component vanishes. For the case with N interacting electrons, the electronic relativistic Hamiltonian is where g(i, j) is the electron-electron interaction. In the nonrelativistic case, g(i, j) corresponds to the Coulomb operator, where the interaction between two electrons is instantaneous. When including relativity, this cannot be correct because the Coulomb interaction between electrons involves the exchange of photons traveling at the speed of light. The relativistic electron-electron operator is much more complicated than the non-relativistic one and cannot be written in closed form. Its perturbation expansion in terms of 1/c 2 yields where g 0 is the instantaneous Coulomb interaction and the first order correction g 1 is the Breit term, 88,89 which introduces magnetic and retardation effects. In our 4c-DKS calculations, we include only the g 0 term. The contributions from the Breit term are believed to be small [90][91][92][93] and they are therefore neglected in most relativistic calculations.
In relativistic KS-DFT, [94][95][96][97][98] the XC functional should in principle also include relativistic effects and should be formulated in terms of the four-current density. 99 The latter is the basic density variable in relativistic KS-DFT. A relativistic generalization of the local density approximation (RLDA) 96 has been proposed, as well as a semi-empirical gradient corrected variant (RGGA). 95,100 Common practice, however, is to use a non-relativistic XC functional in conjunction with the Dirac kinetic energy, 93,99 which is the procedure we follow in this work. We use here a 4c-DKS approach with non-relativistic GGA and hybrid GGA functionals.

Scalar relativistic ZORA approach
The computational cost for a fully relativistic 4c-DKS approach is significantly higher than for the non-relativistic Schrödinger Kohn-Sham (SKS). The scalar relativistic ZORA approximation retains the computational effort of an SKS calculation, and has been shown to capture relativistic effects in good agreement with other scalar-relativistic all-electron schemes 101 .
The ZORA scheme is derived by solving one of the two coupled equations in Equation (9) for the small component χ and inserting it into the other equation, which yields the following (still exact) expression for the unnormalized large component Expanding the parenthetical term as a geometric series yields the regular approximation. 98 Retaining only the scalar part of the zeroth order term transforms Equation (13) to where the ZORA kinetic energy Hamiltonian T ZORA is defined as 102 Since the potential enters non-linearly in the denominator of Equation (15) it is clear that ZORA is gauge dependent, i.e., a constant shift of the electrostatic potential does not lead to a constant shift in the energy. Different methods have been proposed to restore gaugeinvariance. One of them is the popular scaled ZORA approximation, 103 where the eigenvalues are rescaled after self-consistency is reached almost, but not completely restoring gauge-invariance. Full gauge-invariance is achieved in the atomic ZORA scheme (aZORA), which we use in this work. In aZORA the potential in the denominator of Equation (15) is replaced with the onsite free-atom potential v at ( j) near the nucleus, on which the localized basis function ϕ j is centered. The aZORA Hamiltonian depends therefore explicitly on the atom index of the basis function ϕ j it acts upon. We employ the aZORA as defined in Refs. 104, 105 and benchmarked in Ref. 101. Since the kinetic term T aZORA depends on ϕ j , the matrix elements need to be symmetrized to restore Hermiticity, which finally gives While the absolute values of the scaled ZORA eigenvalues are closer to the 4c-DKS reference, 103 we expect the relative shifts with respect to 4c-DKS, which are relevant for the proposed correction scheme, to be more consistent with aZORA. The reason is that the latter, unlike scaled ZORA, restores the gauge invariance completely.

E. Atomic relativistic corrections for GW
For GW , we have developed three simple correction schemes to account for relativistic effects: I) Atomic relativistic corrections are added to the QP energies. II) The aZORA Hamiltonian is used for the underlying DFT calculation and the obtained KS eigenvalues and MOs are used as a starting point for GW . III) aZORA is used as in II and atomic relativistic corrections are added to the QP energies. The atomic corrections are always added as a post-processing step to the converged QP energies and have been obtained as follows.
For scheme I, the atomic relativistic corrections ∆ SKS 1s,at are computed as difference between the non-relativistic SKS 1s eigenvalues ( SKS 1s,at ) and the fully relativistic 4c-DKS 1s eigenvalues ( 4c-DKS 1s,at ), The are computed selfconsistently at the PBE level by solving the radial SKS, aZORA and 4c-DKS equations respectively.

III. COMPUTATIONAL DETAILS
All GW and ∆SCF calculations are performed with the allelectron FHI-aims program package, 24,104,106 which is based on numerically tabulated atom-centered orbitals (NAOs). Core-level BEs from G 0 W 0 , evGW 0 and ∆SCF calculations are calculated for the CORE65 benchmark set introduced in Ref. 20, which contains 65 1s binding energies of secondrow elements (C, N, O and F) for small organic and inorganic molecules. The settings for G 0 W 0 , evGW 0 and ∆SCF are the same as in our previous work 20 and are summarized in the following.
The ∆SCF calculations are performed with the PBE0 107,108 hybrid functional employing def2 quadruple-ζ valence plus polarization (def2-QZVP) 109 basis sets. The all-electron def2-QZVP Gaussian basis sets are treated numerically in FHIaims for compliance with the NAO scheme. We decontract the def2-QZVP basis sets to enable a full relaxation of the other electrons in the presence of a core-hole; see Ref. 20 for further details and an explanation of the basis-set choice.
For the GW calculations, the QP equation (Equation (1)) is always solved iteratively. For the partially self-consistent evGW 0 scheme, we iterate the eigenvalues additionally in G. We use the PBE functional 69 as a starting point for the evGW 0 calculations. For G 0 W 0 , we employ the PBEh(α) hybrid functionals 84 for the underlying DFT calculation, where α indicates the fraction of HF exchange in the functional. The core-level BEs are extrapolated to the complete basis set limit to account for the slow convergence of the GW QP energies with respect to basis set size. 22,27,65,110,111 The extrapolation is performed by a linear regression with respect to the inverse of the total number of basis functions using the Dunning basis set family cc-pVnZ (n=3-6). 112,113 Details are given in the Supporting Information (SI) in Table S3 and comprehensive convergence studies are presented in Figure S1. Furthermore, we use the CD technique 36 to compute the GW self-energy. The integral over the imaginary frequency axis in Equation (7) is computed using modified Gauss-Legendre grids 24 with 200 grid points.
Relativistic effects for GW are included in three different ways as described in Section II E. For ∆SCF, we account for relativistic effects self-consistently using the aZORA approximation. 104 We also apply the atomic relativistic schemes introduced in Section II E to ∆SCF for comparison. To obtain the atomic relativistic corrections for equations (17) and (18), the radial DKS, SKS and aZORA-KS equations are solved self-consistently on numerical real-space grids with the DFTATOM code 114 incorporated in FHI-aims.
We investigate the dependence of the relativistic eigenvalue corrections on the molecular environment, XC functional and basis set using the DIRAC program, 115,116 which features a 4c-DKS DFT implementation for the 3D electronic wave function, enabling also molecular calculations. Similar to Equation (17), we define the molecular corrections as where 4c-DKS 1s,mol are molecular 1s eigenvalues of the 4c-DKS Hamiltonian. The corresponding non-relativistic eigenvalues SKS 1s,mol are here obtained from a 4c-DKS calculation, resetting the speed of light to the non-relativistic limit (c → ∞).
The DIRAC calculations are performed for the molecular structures of the CORE65 benchmark set, excluding the spinpolarized O 2 case, using all-electron Dyall basis sets 117 of triple-zeta quality and the PBE functional. We define the difference ∆MOL between molecular and atomic eigenvalue correction as The functional dependence of the atomic corrections is assessed for the PBEh(α) hybrid family. We also study the basis set dependence for the Dyall series 117 with reference to the fully converged radial solution from DFTATOM In pursuit of open materials science, 118 we made the results of all relevant calculations available on the Novel Materials Discovery (NOMAD) repository. 119

IV. RESULTS AND DISCUSSION
We first present the atomic relativistic corrections and discuss their dependence on technical and convergence parameter, the XC functional and the molecular environment. We  (17)) and aZORA Hamiltonian (Equation (18)) with respect to the atomic number.  (17)) and aZORA Hamiltonian (Equation (18) proceed with a discussion of non-relativistic results for the CORE65 benchmark set and demonstrate how our simple correction schemes, based on these atomic corrections, improve the agreement of the computed absolute 1s BEs to experiment.
A. Atomic relativistic corrections Figure 1(a) shows a sketch of the 1s eigenvalues from the non-relativistic SKS, 4c-DKS, and scalar relativistic aZORA calculations. The SKS eigenvalues are generally overestimated with respect to the 4c-DKS reference, while aZORA underestimates the 1s eigenvalues by nearly as much. The atomic eigenvalue corrections ∆ 1s,at for SKS (Equation (17)) and aZORA (Equation (17)) are given in Table I (21) between the Dyall all-electron basis set at the double, triple and quadruple-ζ levels 117 and the radial solution on numeric real-space grids. b) Atomic eigenvalue correction ∆ SKS 1s,at (Equation (17)) computed with the PBEh(α) functional with 3 different values of α. c) Difference ∆MOL as defined in Equation (20) between the molecular eigenvalue correction and atomic value for C1s, N1s, O1s and F1s excitations for the CORE65 benchmark set. Bars contain the 2nd and 3rd quartiles, whiskers extend to encompass 95% of the results, and outliers are shown as dots.
atomic number, ranging from −4 meV for Li to −13.4 eV for Ar. The aZORA corrections are positive and increase from 4 meV (Li) to 12.4 eV (Ar).
The atomic corrections given in Table I are visualized in Figure 1(b). We observe that the magnitude of the atomic corrections for both SKS and aZORA, depends on the fourth power of atomic number Z. This dependence is known from relativistic correction to the exact energy of the hydrogenic orbital, whose leading order term in a perturbative expansion scales as the fourth power of Z. 87,[120][121][122] As the 1s orbitals are poorly-screened by the outer orbitals, the magnitude of the relativistic correction trends similarly to that of the unscreened hydrogenic orbitals.
The results reported in Table I are obtained from the selfconsistent solutions of the radial SKS and 4c-DKS equations.
The radial equations enforce a spherical symmetry of the solution. However, most atoms have ground states with nonspherical symmetry. For the second-row, this applies to B, C, O and F. For these elements, the spherical solutions are too high in total energy by several tenths of eV and assume fractional occupation numbers. The 1s eigenvalue corrections ∆ SKS 1s,at obtained from the radial SKS and 4c-DKS equation are therefore an approximation. To estimate the error introduced by this approximation, we solved the 3D SKS equations for the free neutral atom and compared the 1s eigenvalues of the spherical and non-spherical solution. The spherical solution is also obtained with the 3D equations and is identical to the radial one, if we do not break the symmetry and enforce integer occupations. The non-spherical 3D solution is obtained by employing occupation constraints. We find that the difference in the absolute 1s eigenvalues between spherical and non-spherical solution is less than 50 meV, see Table S1 (SI), which is an order of magnitude smaller than the relativistic corrections themselves. The error in the relative values, ∆ SKS 1s,at , is expected to be even smaller and we conclude thus that the radial approximation is sufficient.
The radial calculations are performed on a numeric realspace grid, which can be easily converged, whereas the 3D calculations rely on relativistic all-electron Gaussian basis sets, potentially introducing a basis set incompleteness error. Figure 2(a) shows the basis set convergence of the Dyall series with respect to the radial solution. At the double-ζ level, the error is within a few meV, and for the relevant 1s states, we reach convergence already at the triple-ζ level, see Figure 2(a). We use the quadruple-ζ basis set for the calculations shown in Figure 2(b) and the triple-ζ basis set for the calculations in Figure 2(c).
In Figure 2(b), we examine the dependence of the atomic eigenvalue correction ∆ SKS 1s,at on the fraction α of exact HF exchange in the PBEh(α) functional for α = 0 (PBE), α = 0.25 (PBE0) and α = 0.5. The magnitude of the eigenvalue correction shows a slight dependence on α, increasing by an amount that is proportional to the fraction of exact exchange. At first glance, the α dependence seems more pronounced for heavier elements. However, this is only true for the absolute values: Setting the PBE functional as reference and comparing to PBEh(α = 0.5), the magnitude of the atomic correction increases by 12 meV for carbon 1s, which corresponds to 10.2%, and by 40 meV for fluorine 1s, which, however, corresponds only to 5.7%. In fact, the α dependence seems to decrease with the atomic number when comparing relative deviations; see Table S5 (SI) for the tabulated values. For all elements listed in Table I, we find that the α dependence is an order of magnitude smaller than the relativistic correction itself. We thus neglect it when applying our relativistic correction schemes to the 1s QP energies from G 0 W 0 @PBEh.
In Figure 2(c), we compare the atomic eigenvalue correction (Equation (17)) to the molecular eigenvalue correction (Equation (19)). For most of the excitations considered, the atomic eigenvalue correction slightly underestimates the molecular correction, but the difference between the two is under 5 meV for 49 of the 63 excitations considered, with a maximum deviation of 12.6 meV. This distribution of these differences is similar for the core excitations of different elements, and is small enough in comparison with the magnitude of the eigenvalue correction to justify the use of the atomic values irrespective of the chemical environment.
The atomic SKS corrections reported in Table I are Figure 2(b) suggest that these differences must be partly attributed to the exchange treatment. The remaining differences might be due to usage of non-relativistic basis sets in combination with the 4c-DHF Hamiltonian in Ref. 9. The atomic SKS corrections are also surprisingly similar to the corrections derived for second-period elements in an early work from the 1960s based on Pauli perturbation theory of charged 2-electron atoms. 123,124 Pauli perturbation theory is based on the first order in the expansion of Equation (13) in terms of 1/c 2 . It is highly singular in the deep-core region 125 and has been largely replaced by the regular approximation, which expands Equation (13) in terms of 2c 2 −V . The correspondence worsens when valence electrons are included. 126

B. Non-relativistic quasiparticle energies
In our previous work 20 we briefly discussed the effect of relativistic corrections, comparing non-relativistic and relativistic 1s BEs from evGW 0 to experiment. We will now analyze non-relativistic evGW 0 in more detail and additionally include the non-relativistic G 0 W 0 @PBEh and ∆SCF results in the discussion. Figure 3 displays the mean absolute error (MAE) of the absolute 1s BEs with respect to experiment for the CORE65 benchmark set. The MAEs obtained from nonrelativistic evGW 0 calculations increase with atomic number (Figure 3(c)) and the magnitude of this increase is within the range of the atomic relativistic corrections given in Table I. The distribution of these errors is shown in Figure 4(a), where the grouping in species is evident. evGW 0 systematically underestimates the 1s BEs for all 65 excitations. The non-relativistic ∆SCF calculations underestimate the 1s BEs as well and the MAEs show a very similar trend with respect to atomic number, see Figure 3(c). Comparing the overall MAE for the non-relativistic calculations, we find that the MAE for ∆SCF is with 0.71 eV slightly larger than the 0.55 eV MAE for evGW 0 .
Relativistic effects are also apparent when considering the optimal α for use in a G 0 W 0 @PBEh(α) scheme. Figure 3(a) shows the MAE for non-relativistic G 0 W 0 @PBEh(α) calculations with respect to the fraction of exact exchange α. These calculations have been carried for a subset of 43 excitations of the CORE65 benchmark set, for which the mapping between core state and atom does not require analysis of, e.g., MO coefficients. In our previous work 20 , we reported the α dependence of the MAE including relativistic effects (Figure 3(b)) and found that the smallest MAE is obtained for α values around 0.45. For MAEs smaller than 0.45, the BEs are underestimated and for larger α values increasingly overestimated. For the non-relativistic results we observe a much stronger species dependence of the optimal α value. As the non-relativistic Hamiltonian underestimates the core-level BE, increasing the exact exchange reduces the screening, resulting in a larger BE. An increase in α can thus offset the relativistic error. Comparing the MAE from non-relativistic G 0 W 0 @PBEh(α) calculations (Figure 3(a)), we find that the optimal α indeed increases with atomic number, from 0.44 for C1s excitations to about 0.55 for F1s.

C. Atomic and scalar relativistic correction schemes
We investigate three simple schemes to account for relativistic corrections (RC) in 1s core-level BEs from GW . Additionally, we discuss the application of these three schemes to ∆SCF. The first approach is to add the atomic corrections ∆ SKS 1s,at defined in Equation (17) to the QP energies and corresponds to the scheme we employed in our latest work. 20 We label scheme I "method+RC". The second is to use aZORA for the underlying DFT calculations, and use the aZORA eigenvalues and MOs as the starting point for the GW calculation. We refer to scheme II as "method + aZORA". In the third scheme, we use scheme II to obtain the QP energies and add the atomic corrections ∆ aZORA 1s,at defined in Equation (18) afterwards. We label scheme III "method + aZORA + RC".
For evGW 0 we explored also a variant of scheme I, where we added the atomic corrections to the DFT eigenvalues instead to the QP energies. These corrected eigenvalues were then used as starting point for the evGW 0 calculation. This pre-correction variant yields with a mean absolute difference of 16 meV BEs that are extremely similar to the ones from evGW 0 +RC. Adding the atomic correction as post-processing step is transferable to non-relativistic GW results obtained from any code and we thus disregard the pre-correction variant in the following.
Compared with the non-relativistic energies, the evGW 0 +RC scheme reduces the error with respect to experiment, as shown in Figure 4(b). The errors are more tightly distributed, and the clustering by species is no longer evident. Generally, the BEs are still underestimated. However, the overall MAE is reduced from 0.55 to 0.3 eV and is now well within the accuracy required for chemical analysis. Furthermore, the species-dependence in the MAE is largely eliminated; see Figure 3(c) and Table II. Solely the MAE for the F1s excitations is with 0.44 eV slightly larger than for the other elements. This might be attributable to poor statistics since our benchmark set contains only 3 F1s excitation.
Scheme I has also been successfully employed for G 0 W 0 @PBEh. The range of optimal α is reduced by a factor of two for the G 0 W 0 @PBEh(α)+RC scheme vis-a-vis the non-relativistic one, see Figure 3(a,b). With the relativistic correction, the value of α that minimizes the MAE ranges  Table I. from 0.44 for C1s excitations to 0.49 for F1s excitations. This shows also in a slight species-dependent of the MAE value we reported for the G 0 W 0 @PBEh(α=0.45) results with RC earlier. 20 Judging by MAE alone (Table II), the evGW 0 +aZORA results are an improvement over the evGW 0 +RC scheme. The overall MAE is 0.18 eV. Their distribution (Figure 4(c)) is more centered. A slight clustering by species is observed, although this is not as obvious as for the non-relativistic results shown in Figure 4(a). In contrast to the non-relativistic values, the 1s excitations of the lighter elements, such as carbon, tend to be more underestimated than the 1s BEs of oxygen.
The evGW 0 +aZORA+RC scheme performs worse than evGW 0 +aZORA approach, with an error distribution similar to evGW 0 +RC, see Figure 4(d). The MAE for the individual species are in the same range as for evGW 0 +RC, as is the overall MAE with 0.32 eV, see Table II. The evGW 0 +aZORA+RC scheme is the most sophisticated among the three relativistic corrections discussed here: scalar-relativistic effects are included in the MOs used as starting point for the evGW 0 calculation and the QP energies are corrected with respect to the fully-relativistic atomic reference. It is thus surprising that it performs worse than evGW 0 +aZORA. The similar performance of evGW 0 +RC and evGW 0 +aZORA+RC rather implies that the effect of including relativistic effects in the MOs is minimal.
To further investigate this surprising behavior, we visualized the mean errors (MEs) for C1s, N1s, O1s and F1s in Figure 5(b). The error is defined as BE experiment − BE theory . Negative MEs indicate thus a systematic underestimation of the BEs with respect to experiment. For the non-relativistic results, the ME corresponds directly to the MAE and is increasingly negative with atomic number. The MEs for evGW 0 +RC and evGW 0 +aZORA+RC are negative and show almost no species dependence, which is in agreement with our previous analysis of the MAE and the error distribution in Figure 4. For evGW 0 +aZORA, however, we observe a trend that is reverse to the non-relativistic results. The ME increases with atomic number and becomes even positive for F1s.
Comparing non-relativistic with the relativistic BEs, we find that the size of the relativistic correction is 2-3 times larger with evGW 0 +aZORA than with the other two schemes, see Figure 5(c). In combination with the upwards trend observed for the ME, this implies that aZORA is overestimating the relativistic correction for 1s states. This reflects the well-known tendency of aZORA to overestimate the relativistic correction to core-state eigenvalue . 60,61 For the 2nd period elements, where the relativistic error ranges from 0.2−0.7 eV, the aZORA overcorrection compensates in part a chronic underestimation of the BEs. While this error cancellation may seem fortuitous for 2nd row elements, the rapid growth of the relativistic correction with the atomic number implies that evGW 0 +aZORA might lead to large errors for 1s BEs of heavier elements. To illustrate this, we analyze a small number of phosphorus and sulfur containing small molecules alongside the CORE65 benchmark set: H 2 S, SO 2 , PH 3 , PF 3 and PF 5 (see Figure 5c). The relativistic corrections obtained with the evGW 0 +aZORA scheme are more than 10 eV larger than with evGW 0 +RC and evGW 0 +aZORA+RC. Also the difference between the evGW 0 +aZORA+RC and evGW 0 +RC, which is negligible for 2nd period elements, becomes more significant: 2.1 eV for the P1s excitations, and 2.6 eV for the S1s. This suggests that the use of a scalar-relativistic reference for the underlying DFT calculation becomes more rele-vant as the magnitude of relativistic effects increase. However the effect of the relativistic reference is only about one third the magnitude of the relativistic correction for these states.
We applied the three correction schemes also to 1s BEs obtained from ∆SCF and plotted the MEs by species in Figure 5(a). We observe the same trends as for evGW 0 . Note the similarity between Figure 5(a) and (b). ∆SCF+RC or ∆SCF+aZORA+RC largely eliminate the species dependence of the ME. Although it is not as marked as for evGW 0 +aZORA, the MEs increase also slightly with atomic number for ∆SCF+aZORA. The size of the relativistic correction is two times larger with ∆SCF+aZORA than with ∆SCF+RC or ∆SCF+aZORA+RC, which suggests that the relativistic 1s corrections are also overestimated at the ∆SCF+aZORA level. This is important detail to consider when comparing the performance of XC functionals for ∆SCF calculations of 1s excitations since the relativistic treatment makes a difference. For example, a recent study 10 with the SCAN functional uses the ∆SCF approach in combination with scaled ZORA, while an atomic correction scheme was used for a similar benchmark study 9 with the TPSS functional. Both studies report very good agreement with experiment. However, it is difficult to judge, which functional performs better, since the relativistic effects are not treated on equal footing.
Both schemes that consistently improve the agreement with experiment, evGW 0 +RC and evGW 0 +aZORA+RC, chronically underestimate the 1s BEs. This might be attributed to the broadening of the experimental spectra due to vibration effects, while GW yields vertical excitation energies. It has been demonstrated for GW -computed excitations of frontier orbitals that the deviation to experiment can be reduced to 0.1 eV when fully resolving the vibrational structure based on Franck-Condon multimode analysis, performed as postprocessing step to the GW calculation. 127 This level of accuracy is often not required to resolve most XPS spectra, but could be in principle reached by applying the same approach.

V. CONCLUSION
Relativistic corrections for 1s core-level energies from GW have been derived for non-relativistic and scalar-relativistic starting points. We have investigated three schemes for 1s QP energies from evGW 0 : A post-GW atomic correction (evGW 0 +RC) using a non-relativistic reference, employ-ing an aZORA reference (evGW 0 +aZORA), and employing a aZORA reference along with a post-GW atomic correction (evGW 0 +aZORA+RC). All three schemes improve agreement with experiment. The evGW 0 +RC and evGW 0 +RC+aZORA schemes reduce the mean absolute error to about 0.3 eV and eliminate the species dependence. The evGW 0 +aZORA scheme further reduces the overall MAE to 0.2 eV, but does so inconsistently, due to a species-dependent overcorrection.
The similarity of the results for the evGW 0 +RC and evGW 0 +aZORA+RC schemes indicates that the use of the scalar-relativistic reference has no significant effect on the result. Of the two, the evGW 0 +RC scheme offers the further advantage that it is readily applicable in codes that have not implemented the aZORA Hamiltonian. We have shown that the derived corrections for the non-relativistic reference improve also consistently core-level energies from G 0 W 0 and ∆SCF, suggesting that they are generally applicable to core-level BEs from different theoretical methods.
The evGW 0 +RC and evGW 0 +aZORA+RC schemes correct the non-relativistic values in a consistent manner, and therefore form a solid foundation for further development and refinement. Further refinements may include the inclusion of vibrational effects to further improve the accuracy, in particular for C1s excitations, which are generally subject to very small chemical shifts. The development of relativisticallycorrected GW also paves the way for the accurate calculation of XPS the in condensed phase systems, where ∆-approaches are problematic. Relativistic corrections will also improve the accuracy of X-ray absorption spectra from the Bethe-Salpeter equation, which employs the quantities computed in the GW calculations. 128,129