Relativistic correction scheme for core-level binding energies from GW Open Access

Core-level binding energies (BEs), measured by x-ray photoemission spectroscopy (XPS), are element-specific but also depend 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–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.⁴ 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⁷ and is the focus of this work.

The most common approach to calculating core-level BEs is the delta self-consistent field (ΔSCF) method,⁸ wherein one computes the total energy difference between the ground and core-ionized states 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 the experiment for small molecules.^9,10 A similar accuracy has been obtained with high-level wavefunction-based delta coupled cluster methods,^11–13 albeit at much higher computational cost. The introduction of occupation constraints and the explicit generation of a charged system in Δ-based approaches lead to a plethora of problems.^14,15 Most importantly, the application to periodic systems requires 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–20 One particularly promising response approach is the GW approximation²¹ to many-body perturbation theory. GW offers access to the core-state quasiparticle energies, which can be directly related to the core-level BEs. The GW approach is the standard approach to compute valence excitations (band structures) in solid-state physics.²² With the advent of efficient implementations with localized basis sets,^23–28 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,29–31 GW has also been successfully applied to deep valence and semi-core states with excitation energies up to 80 eV.^32–34 However, the application of GW for deep core states with BEs larger than 100 eV, which are the ones relevant for the chemical analysis, has rarely been attempted. The first GW studies^19,35–37 for deep core states reported partly large deviations of several electronvolts from the experiment. Recently, we have shown that GW can also be successfully applied for relative and absolute 1s molecular core-level BEs when utilizing highly accurate techniques for the integration of the self-energy³⁸ and eigenvalue self-consistency in the Green’s function.²⁰ 

The application of GW (or any other method) to absolute core-level binding energies requires an accurate treatment of relativistic effects. Already for the 1s core levels of the p-block second 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, the CORE65 benchmark set), obtaining a mean deviation from the experiment of 0.3 eV.²⁰ 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 enable the description of noncollinear electronic states. GW equations with an explicit spin-dependent electron–electron interaction were only developed 12 years ago^39,40 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 two-component Dirac–KS-DFT calculations,^41,42 KS-DFT with second-variation SOC⁴³ or KS-DFT with fully relativistic pseudopotentials.^44,45 These implementations are in the following referred to as fully relativistic GW approaches.

Fully relativistic GW calculations 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,^46,47 transition metal chalcogenides^43,44 and dichalcogenides,^45,48,49 perovskites,^44,50 and bismuth-based topological insulators.^51–56 Compared to non-relativistic GW, the fully relativistic approach is computationally at least four times more expensive.⁴² Alternatively, the SOC can be added as a perturbative correction to the quasiparticle band structure, which was common in earlier GW studies.⁵⁷ Since this approach is computationally less expensive, it continues to be employed.^58–60 However, it has been shown that the SOC post-correction scheme can fail, for example, to describe the band inversion in topological insulators.⁵³ 

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,²² fully relativistic GW calculations are rare and have been mostly conducted for diatomic molecules,^41,44 but recently also for transition metal halide complexes.⁴² 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).³¹ 

Relativistic considerations for the innermost core-levels differ from those for valence excitations. SOC leads, for core states with the 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 fourth-period elements to several tens of eV.⁶¹ SOC affected core states require, thus, a noncollinear treatment, and we are only aware of one GW study⁴¹ 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 eV to 700 eV are measured in XPS. While scalar relativistic effects are heightened in the proximity of the poorly screened nuclear charge, SOC is suppressed for 1s core states. In contrast to valence excitations, the most common scalar-relativistic approximation, ZORA (or, in the variant used in this work and explained in Sec. II D 2, atomic ZORA^62,63), performs poorly for the absolute eigenvalues associated with the innermost core levels.^64,65 In order to avoid the computational expense of a fully Dirac-like, two-component GW calculation, which is not essential for the physics of 1s core excitations, we derive here an element-specific relativistic corrective term for non-relativistic and scalar-relativistic reference states, which we add in a post-GW perturbative step.

The remainder of this article describes this correction. In Sec. 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 Sec. 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 Sec. IV and finally draw conclusion in Sec. V.

In practice, GW is often performed as a one-shot perturbative approach (G₀W₀) 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 an input (starting point) for the G₀W₀ calculation. The G₀W₀ quasiparticle (QP) energies $ϵ n G 0 W 0$ are obtained by solving the QP equation,

and can be related to the BE of state n by $B E n σ =− ϵ n σ G 0 W 0$⁠, where σ denotes the spin index. The nth diagonal elements of the KS exchange–correlation (XC) potential and self-energy operator Σ are denoted by $v n σ x c = ϕ n σ | v x c | ϕ n σ$ and $Σ n σ = ϕ n σ | Σ σ | ϕ n σ$⁠, respectively. The self-energy operator is given by

where $G 0 σ$ is the non-interacting KS Green’s function, W₀ 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 the occupied and virtual KS orbitals. Within the random phase approximation, the screened Coulomb interaction is given by

with the bare Coulomb interaction v(r, r′) = 1/|r − r′| and the dielectric function, ϵ,

where the irreducible polarizability χ₀ is given in the real-space Adler–Wiser representation^66,67 as

The index i runs over occupied orbitals, and a runs over virtual orbitals. Note that we drop the spin index σ in the following for simplicity.

We split the self-energy in a correlation part Σ_c and an exchange part Σ_x, where Σ = Σ_c + Σ_x. The exchange term is frequency-independent and calculated as the HF-like exchange with the KS orbitals, as described in Ref. 38. Details on the calculation of the correlation part are given in Sec. II B.

Equation (1) is non-linear and can be solved iteratively or approximately by linearization using a Taylor expansion to the first order around ϵ_n.²² As we pointed out in our previous work,^20,38 the linearization error increases rapidly with an increase in BE and can already amount to 0.5 eV for deeper valence states.²² 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.

The accurate frequency integration of the self-energy [Eq. (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,29,68 and yields accurate results for valence states.⁶⁹ 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 reproduce the self-energy structure.³⁸ 

For core states, a more accurate evaluation of the self-energy on the real frequency axis is required. We employ the contour deformation (CD) technique to evaluate the correlation self-energy Σ_c,^{22,23,38,70–72} 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₀ are enclosed in the contours, and the contour integral is evaluated using the residue theorem; see Refs. 22 and 38 for details. The integral along the real frequency axis is then evaluated as

where θ is the Heaviside step function, $W 0 c ( r , r ′ , ω ) = W 0 ( r , r ′ , ω ) −v ( r , r ′ )$⁠, 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 Eq. (2).³⁸ Recently, combining the CD approach with analytic continuation of W₀ has been proposed⁷³ and might be an alternative approach for the accurate calculation of the self-energy for core states.⁷⁴ 

In this section, we briefly present the GW variants that we used in this work. By now, many different GW flavors have emerged in practical calculations.²² The most common flavor, G₀W₀ 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 QP peak in the spectral function $A ( ω ) =1/π ∑ m I m G m$⁠, where m runs over all occupied and virtual states and G = G₀ + G₀ΣG. For molecular valence states, G₀W₀ performed on top of a DFT calculation with the Perdew–Burke–Ernzerhof (PBE)⁷⁵ functional (G₀W₀@PBE) yields a clearly identifiable QP peak. This peak corresponds to a distinct solution of Eq. (1). Multiple solutions, which would indicate spectral weight transfer to other peaks in the spectral function, have only been observed for a few systems for frontier orbitals.^69,76–78 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₀W₀@PBE for 1s states.²⁰ 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^79,80 and semi-core excitations of sodium.³³ However, for deep core states, the transfer of spectral weight is far more extreme resulting in a spectral function, where the 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 magnitude smaller than the main excitation.^81,82

We showed that the 1s QP peak can be correctly restored by including eigenvalue self-consistency in G, while keeping PBE as the starting point.²⁰ This scheme is referred to as evGW₀@PBE. In evGW₀, the screened Coulomb interaction is kept fixed at the W₀ level, and the Green’s function is recomputed replacing the mean field eigenvalues with the QP energies from Eq. (1). We enforce eigenvalue self-consistency only in G. Inserting the QP energies also in W₀ (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₀ yields bandgaps in good agreement with the experiment,⁸³ while underscreening errors in the evGW scheme lead to too large bandgaps⁸³ and overly stretched spectra.⁸⁴ 

Higher-level self-consistency schemes, such as fully self-consistent GW (scGW)^85–87 or quasiparticle self-consistent GW (QSGW),⁸⁸ are expected to restore the 1s QP peak as well but might not yield better agreement with the experiment than evGW₀. It has been shown that scGW overestimates molecular HOMO excitations⁸⁹ and bandgaps in solids.⁹⁰ 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.³⁷ 

evGW₀ is computationally more demanding than a G₀W₀ calculation because the QP equation is not only solved for the 1s core states of interest but also repeatedly for all occupied and virtual states until convergence in G is reached. We showed that the core-level QP peak can also be restored in a G₀W₀ calculation by using a XC functional with a high fraction of exact exchange as the starting point.²⁰ We employ the PBEh(α) functional family with an adjustable amount α of HF exact exchange.⁹¹ The XC energy E_xc is given by

where $E x E X$ denotes the HF exchange energy. $E x P B E$ and $E c P B E$ are the PBE exchange and correlation energy, respectively.

In this work, we followed Ref. 20 and used both evGW₀ and G₀W₀@PBEh(α). We analyze how the relativistic corrections we devised affect the two schemes.

It has long been recognized that relativistic effects play a large role in the chemistry of heavy elements.^92,93 In this work, we treat the core states of light elements carbon through fluorine, whose BEs are in the range of 250 eV–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 second 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 four-component Dirac Kohn–Sham (4c-DKS) approach, further also referred to as the fully relativistic scheme. The second uses the scalar-relativistic atomic ZORA.

The relativistic description of a non-interacting electron in an external potential V is given by the Dirac equation,⁹⁴ 

where $Ψ= ϕ χ$ is a four-component Dirac spinor, which is comprised of a large component ϕ and a small component χ, each of which has two components for the spin functions. The Dirac Hamiltonian shifted by the rest energy c² is given by

where c is the speed of light, σ is a vector of Pauli spin matrices, and p is the momentum operator. For electron-like states in the non-relativistic 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 non-relativistic 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 the closed form. Its perturbation expansion in terms of 1/c² yields

where g₀ is the instantaneous Coulomb interaction and the first order correction g₁ is the Breit term,^95,96 which introduces magnetic and retardation effects. In our 4c-DKS calculations, we include only the g₀ term. The contributions from the Breit term are believed to be small,^97–100 and they are, therefore, neglected in most relativistic calculations.

In relativistic KS-DFT,^101–105 the XC functional should, in principle, also include relativistic effects and should be formulated in terms of the four-current density.¹⁰⁶ The latter is the basic density variable in relativistic KS-DFT. A relativistic generalization of the local density approximation (RLDA)¹⁰³ has been proposed, as well as a semi-empirical gradient corrected variant (RGGA).^102,107 Common practice, however, is to use a non-relativistic XC functional in conjunction with the Dirac kinetic energy,^100,106 which is the procedure we follow in this work. Here, we use a 4c-DKS approach with non-relativistic GGA and hybrid GGA functionals.

The computational cost for a fully relativistic 4c-DKS approach is significantly higher than for non-relativistic Schrödinger Kohn–Sham (SKS). The scalar-relativistic atomic ZORA [see Eq. (18)] 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.^62,63,108

The original ZORA scheme is derived by solving one of the two coupled equations in Eq. (9) for the small component χ and inserting it into the other equation, which yields the following (still exact) expression for the un-normalized large component:

and can be rewritten as

Expanding the parenthetical term in Eq. (13) as a geometric series yields the regular approximation,¹⁰⁵ 

Retaining only the scalar part of the zeroth order term (k = 0) transforms Eq. (15) to

where the ZORA kinetic energy T^ZORA is defined as¹⁰⁹ 

The zeroth order is the only order in this expansion that leaves the Hamiltonian independent of ϵ_n, enabling a standard numerical treatment of the ensuing KS equations. From Eq. (14), it is apparent that the ZORA is equivalent to the supposition that ϵ_n ≪ 2c². This relation is correct for valence states, i.e., ZORA correctly reflects their scalar-relativistic behavior. In contrast, ϵ_n ≪ 2c² is a potentially severe approximation for the deep core states. However, the core states are relatively inert, i.e., the chemical shifts are small compared to the absolute value of ϵ_n. The error introduced for core states by neglecting the higher order terms $( ϵ n / ( V − 2 c 2 ) ) k$ in Eq. (15) is, therefore, largely an atomic property. For a given atomic species, orbital, and orbital occupation, the dependence of this error on different chemical environments is minor. Physical quantities that depend on a difference in energies are, therefore, not strongly affected.

Since the potential V enters in the denominator of Eq. (17), rather than as ϵ_n − V, it is clear that ZORA is gauge dependent, i.e., a constant shift in the electrostatic potential does not lead to a constant shift in the energy. In practice, this dependence on the choice of the zero-energy level of the potential in different chemical environments would be a severe problem. Different methods have been proposed to restore gauge-invariance. One of them is the scaled ZORA,¹¹⁰ where the eigenvalues are rescaled after self-consistency is reached. Scaled ZORA restores almost, but not completely gauge-invariance.

Full gauge-invariance is achieved in the atomic ZORA scheme (aZORA), which we use in this work. We employ the aZORA approach as defined in Eqs. (55) and (56) of Ref. 62 and benchmarked in Refs. 63 and 108. In aZORA, the potential in the denominator of T^ZORA is replaced with the potential v_at(j) of a radially symmetric free atom in empty space, which is calculated once at the outset of the calculation and, therefore, does not depend on the chemical environment in any way. The potential v_at(j) is evaluated using boundary conditions with zero potential at infinity (i.e., non-periodic boundary conditions). The dependence on index j indicates that v_at(j) corresponds to the potential of the atomic species at the center where the basis function φ_j is localized. For semilocal DFT calculations with, e.g., the PBE functional, we construct v_at(j) using the same functional. For our hybrid DFT calculations with the PBEh family, v_at(j) is, instead, constructed using PBE exchange–correlation for simplicity. This approximation is entirely minor for the same reasons as using v_at(j) instead of V in the first place. The dominant effect of either V or v_at(j) in Eq. (17) with respect to the 2c² term is the deep electron–nucleus potential, modified by the overall electrostatic screening of the surrounding core electrons. In contrast, differences associated with v^xc are negligible compared to the magnitude of 2c².

The atomic potential v_at(j) introduces a dependence on the localized basis functions φ_j in the kinetic term T^aZORA. Therefore, 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,¹¹⁰ 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.

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 $Δ ϵ 1 s , a t S K S$ are computed as the difference between the non-relativistic SKS 1s eigenvalues $( ϵ 1 s , a t S K S )$ and the fully relativistic 4c-DKS 1s eigenvalues $( ϵ 1 s , a t 4 c - D K S )$⁠,

The label “at” indicates that the calculations are performed for a free neutral atom. For scheme III, we use the atomic corrections $Δ ϵ 1 s , a t a Z O R A$⁠,

evaluating the difference to the aZORA 1s eigenvalues (⁠$ϵ 1 s , a t a Z O R A$⁠) instead to the SKS eigenvalues. The atomic 1s eigenvalues $ϵ 1 s , a t S K S$⁠, $ϵ 1 s , a t a Z O R A$⁠, and $ϵ 1 s , a t 4 c - D K S$ are computed self-consistently at the PBE level by solving the radial SKS, aZORA, and 4c-DKS equations, respectively.

All GW and ΔSCF calculations are performed with the all-electron FHI-aims program package,^24,62,111 which is based on numerically tabulated atom-centered orbitals (NAOs). Core-level BEs from G₀W₀, evGW₀, and ΔSCF calculations are calculated for the CORE65 benchmark set introduced in Ref. 20, which contains 65 1s binding energies of second row elements (C, N, O, and F) for 31 closed- and one open-shell molecule. The settings for G₀W₀, evGW₀, and ΔSCF are the same as in our previous work²⁰ and are summarized in the following.

The ΔSCF calculations are performed with the PBE0^112,113 hybrid functional employing def2 quadruple-ζ valence plus polarization (def2-QZVP)¹¹⁴ basis sets. The all-electron def2-QZVP Gaussian basis sets are treated numerically in FHI-aims 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 [Eq. (1)] is always solved iteratively. For the partially self-consistent evGW₀ scheme, we iterate the eigenvalues additionally in G. We use the PBE functional⁷⁵ as a starting point for the evGW₀ calculations. As discussed in Sec. II C, a distinct 1s QP solution does not exist at the G₀W₀@PBE level, i.e., for the initial evGW₀@PBE step. We, therefore, initialize the evGW₀ loop with approximate 1s QP solutions obtained, e.g., by linearizing the 1s QP equation or employing a larger value for the broadening parameter η. Note that this is only done for the first evGW₀@PBE steps. For G₀W₀, we employ the PBEh(α) hybrid functionals⁹¹ 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 the basis set size.^{22,27,69,115,116} All empty states that can be resolved by the given basis set are included in the calculation of the self-energy. 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).^117,118 Details are given in Table S3 of the supplementary material, and comprehensive convergence studies are presented in Fig. S1. Furthermore, we use the CD technique³⁸ to compute the GW self-energy. The integral over the imaginary frequency axis in Eq. (7) is computed using modified Gauss–Legendre grids²⁴ with 200 grid points.

Relativistic effects for GW are included in three different ways, as described in Section IIE. For ΔSCF, we account for relativistic effects self-consistently using the aZORA.⁶² We also apply the atomic relativistic schemes introduced in Sec. II E to ΔSCF for comparison. To obtain the atomic relativistic corrections for Eqs. (19) and (20), the radial DKS, SKS, and aZORA-KS equations are solved self-consistently on numerical real-space grids with the DFTATOM code¹¹⁹ incorporated into FHI-aims.

We investigate the dependence of the relativistic eigenvalue corrections on the molecular environment, XC functional, and basis set using the DIRAC program,^120,121 which features a 4c-DKS DFT implementation for the 3D electronic wave function, enabling also molecular calculations. Similar to Eq. (19), we define the molecular corrections as

where $ϵ 1 s , m o l 4 c - D K S$ are molecular 1s eigenvalues of the 4c-DKS Hamiltonian. The corresponding non-relativistic eigenvalues $ϵ 1 s , m o l S K S$ 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 spin-polarized O₂ case, using all-electron Dyall basis sets¹²² of triple-zeta quality and the PBE functional. We define the difference ΔMOL between the molecular and atomic eigenvalue corrections 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¹²² with reference to the fully converged radial solution from DFTATOM,

We first present the atomic relativistic corrections and discuss their dependence on technical and convergence parameter, the XC functional, and the molecular environment. We 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 the experiment.

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 [Eq. (19)] and aZORA [Eq. (20)] are given in Table I. The SKS corrections are negative and increase in magnitude with the 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 Fig. 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 for the Schrödinger equation is known from the 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,^94,123–125 and a similar result has been discussed for the ZORA Hamiltonian as well.⁶⁴ 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 self-consistent 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 non-spherical 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 $Δ ϵ 1 s , a t S K S$ obtained from the radial SKS and 4c-DKS equations 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 the spherical and non-spherical solutions is less than 50 meV, see Table S1 of the supplementary material, which is an order of magnitude smaller than the relativistic corrections themselves. The error in the relative values, $Δ ϵ 1 s , a t S K S$⁠, is expected to be even smaller, and thus, we conclude that the radial approximation is sufficient.

The radial calculations are performed on a numeric real-space 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 Fig. 2(a). We use the quadruple-ζ basis set for the calculations shown in Fig. 2(b) and the triple-ζ basis set for the calculations in Fig. 2(c).

In Fig. 2(b), we examine the dependence of the atomic eigenvalue correction $Δ ϵ 1 s , a t S K S$ 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 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 S6 of the supplementary material 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₀W₀@PBEh.

In Fig. 2(c), we compare the atomic eigenvalue correction [Eq. (19)] to the molecular eigenvalue correction [Eq. (21)]. For most of the excitations considered, the atomic eigenvalue correction slightly underestimates in magnitude the molecular correction, but the difference between the two is below 5 meV for 49 of the 63 excitations considered, with a maximum deviation of 12.6 meV. The distribution of these differences is similar for the core excitations of different elements and 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 very similar to the atomic corrections published in Ref. 9 for second row elements B → F. The atomic corrections in Ref. 9 are 10 meV–40 meV larger than ours and were computed by comparing the four-component Dirac–Hartree–Fock (DHF) energies with the non-relativistic HF energies. Our analysis of different PBEh functionals in Fig. 2(b) suggests that these differences must be partly attributed to the exchange treatment. The remaining differences might be due to the 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 two-electron atoms.^126,127 Pauli perturbation theory is based on the first order in the expansion of Eq. (13) in terms of 1/c². It is highly singular in the deep core region¹²⁸ and has been largely replaced by the regular approximation, which expands Eq. (13) in terms of $ϵ n V − 2 c 2$⁠. The correspondence worsens when valence electrons are included.¹²⁹ 

In our previous work,²⁰ we briefly discussed the effect of relativistic corrections, comparing non-relativistic and relativistic 1s BEs from evGW₀ to experiment. We will now analyze the non-relativistic evGW₀ calculations in more detail and additionally include the non-relativistic G₀W₀@PBEh and ΔSCF results in the discussion.

Figure 3 displays the mean absolute error (MAE) of the absolute 1s BEs with respect to the experiment for the CORE65 benchmark set. The MAEs obtained from non-relativistic evGW₀ calculations increase with the atomic number [Fig. 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 Fig. 4(a), where the grouping in species is evident. evGW₀ 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 the atomic number; see Fig. 3(c). Comparing the overall MAE for the non-relativistic calculations, we find that the MAE for ΔSCF is 0.83 eV, which is slightly larger than the 0.55 eV MAE for evGW₀.

Relativistic effects are also apparent when considering the optimal α for use in a G₀W₀@PBEh(α) scheme. Figure 3(a) shows the MAE for non-relativistic G₀W₀@PBEh(α) calculations with respect to the fraction of exact exchange α. These calculations have been carried out for a subset of 43 excitations of the CORE65 benchmark set, for which the mapping between the core state and the atom does not require the analysis of, e.g., MO coefficients. In our previous work,²⁰ we reported the α dependence of the MAE including relativistic effects [Fig. 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 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₀W₀@PBEh(α) calculations [Fig. 3(a)], we find that the optimal α indeed increases with the atomic number, from 0.44 for C1s excitations to about 0.53 for F1s.

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 $Δ ϵ 1 s , a t S K S$ defined in Eq. (19) to the QP energies and corresponds to the scheme we employed in our latest work.²⁰ 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 $Δ ϵ 1 s , a t a Z O R A$ defined in Eq. (20) afterward. We label scheme III “method + aZORA + RC.”

For evGW₀, we also explored a variant of schemes I and III, where we added the atomic corrections to the DFT eigenvalues instead to the QP energies. These corrected eigenvalues were then used as the starting point for the evGW₀ calculation. These pre-correction variants yield with a mean absolute difference of 16 meV (pre-correction variant of scheme I) and 2.2 meV (pre-correction variant of scheme III) BEs that are extremely similar to the ones from evGW₀ + RC and evGW₀ + aZORA + RC, respectively; see Table S5 and Fig. S2 of the supplementary material for details. Adding the atomic correction as a 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₀ + RC scheme reduces the error with respect to the experiment, as shown in Fig. 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 eV to 0.30 eV and is now well within the accuracy required for the chemical analysis. Furthermore, the species dependence in the MAE is largely eliminated; see Fig. 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 three F1s excitation.

Scheme I has also been successfully employed for G₀W₀@PBEh. The range of optimal α is reduced by a factor of two for the G₀W₀@PBEh(α) + RC scheme vis-á-vis the non-relativistic one; see Figs. 3(a) and 3(b). With the relativistic correction, the value of α that minimizes the MAE ranges from 0.44 for C1s excitations to 0.49 for F1s excitations. This also shows in a slight species dependence of the MAE values we reported for the G₀W₀@PBEh(α = 0.45) results with RC earlier.²⁰ 

Judging by the MAE alone (Table II), the evGW₀ + aZORA results are an improvement over the evGW₀ + RC scheme. The overall MAE is 0.18 eV. Their distribution [Fig. 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 Fig. 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₀ + aZORA + RC scheme performs worse than the evGW₀ + aZORA approach, with an error distribution similar to evGW₀ + RC; see Fig. 4(d). The MAE for the individual species is in the same range as for evGW₀ + RC, as is the overall MAE with 0.32 eV; see Table II. The evGW₀ + aZORA + RC scheme is the most sophisticated among the three relativistic corrections discussed here: scalar-relativistic effects are included in the MOs used as the starting point for the evGW₀ 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₀ + aZORA. The similar performance of evGW₀ + RC and evGW₀ + 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 Fig. 5(b). The error is defined as BE_theory − BE_experiment. Thus, negative MEs indicate a systematic underestimation of the BEs with respect to the experiment. For the non-relativistic results, the absolute value of the ME corresponds directly to the MAE, and the ME is increasingly negative with the atomic number. The MEs for evGW₀ + RC and evGW₀ + 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 Fig. 4. For evGW₀ + aZORA, however, we observe a trend that is reverse to the non-relativistic results. The ME increases with the 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₀ + aZORA than with the other two schemes; see Fig. 5(c). In combination with the upward 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 eigenvalues.^64,65 For the second period elements, where the relativistic error ranges from 0.2 eV to 0.7 eV, the aZORA overcorrection compensates, in part, a chronic underestimation of the BEs. While this error cancellation may seem fortuitous for second row elements, the rapid growth of the relativistic correction with the atomic number implies that evGW₀ + aZORA might lead to large errors for 1s BEs of heavier elements. To illustrate this, we analyze the relativistic corrections for a small number of phosphorus and sulfur containing small molecules alongside those for the CORE65 benchmark set: H₂S, SO₂, PH₃, PF₃, and PF₅ [see Fig. 5(c)]. The relativistic corrections obtained with the evGW₀ + aZORA scheme are more than 10 eV larger than with evGW₀ + RC and evGW₀ + aZORA + RC. Also the difference between evGW₀ + aZORA + RC and evGW₀ + RC, which is negligible for second 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 relevant as the magnitude of relativistic effects increases. However, the effect of the relativistic reference is only about one third of the magnitude of the relativistic correction for these states.

We also applied the three correction schemes to 1s BEs obtained from ΔSCF and plotted the MEs by species in Fig. 5(a). We observe the same trends as for evGW₀. Note the similarity between Figs. 5(a) and 5(b). ΔSCF + RC or ΔSCF + aZORA + RC largely eliminate the species dependence of the ME. Although it is not as marked as for evGW₀ + aZORA, the MEs also increase slightly with the 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. ΔSCF entails the difference in total energy between the ionized and neutral systems. However, the ZORA overcorrection of core states propagates to the final calculation of the BE in approximately the same manner as in the 1s levels. This can be seen from Slater transition state theory:^130,131 In ΔSCF, we calculate the energy difference ℏω = E(N − 1) − E(N), where N is the number of electrons and E(N − 1) and E(N) are the total energies of the core-ionized and neutral systems, respectively. This energy difference is approximately related to the 1s eigenvalue by $ℏω=∂E/∂ n 1 s +O ( ( δ n 1 s ) 3 )$ and ∂E/∂n_1s = ϵ_1s, where n_1s is the occupation number of the 1s state.¹³² This is an 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¹⁰ 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⁹ with the Tao-Perdew-Staroverov-Scuseria (TPSS) functional. Both studies report very good agreement with the 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 the experiment, evGW₀ + RC and evGW₀ + 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 the experiment can be reduced to 0.1 eV when fully resolving the vibrational structure based on the Franck–Condon multimode analysis, performed as a post-processing step to the GW calculation.¹³³ This level of accuracy is often not required to resolve most XPS spectra but could be, in principle, reached by applying the same approach.

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₀: A post-GW atomic correction (evGW₀ + RC) using a non-relativistic reference, employing an aZORA reference (evGW₀ + aZORA) and an aZORA reference along with a post-GW atomic correction (evGW₀ + aZORA + RC). All three schemes improve agreement with the experiment. The evGW₀ + RC and evGW₀ + RC + aZORA schemes reduce the mean absolute error to about 0.3 eV and eliminate the species dependence. The evGW₀ + 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₀ + RC and evGW₀ + aZORA + RC schemes indicates that the use of the scalar-relativistic reference has no significant effect on the result. Of the two, the evGW₀ + RC scheme offers the further advantage that it is readily applicable in codes that have not implemented the aZORA Hamiltonian, and this is our recommendation to the practitioner. The application of the scheme entails an a posteriori shift of the non-relativistic evGW₀ energies by the $Δ ϵ 1 s , a t S K S$ values in Table I. We have shown that the derived corrections for the non-relativistic reference also consistently improve core-level energies from G₀W₀ and ΔSCF, suggesting that they are generally applicable to core-level BEs from different theoretical methods.

The evGW₀ + RC and evGW₀ + 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 relativistically corrected GW also paves the way for the accurate calculation of XPS 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.^134,135

See the supplementary material for the comparison of 1s eigenvalues for spherical and non-spherical solutions with the KS equations for neutral atoms (Table S1); plot of basis set dependence of the extrapolation scheme (Fig. S1); convergence of KS 1s eigenvalues and total energies for cc-pVnZ and NAO-VCC-nZ basis set series (Table S2); scalar-relativistic evGW₀@PBE + aZORA results for basis set series cc-pVnZ (n = 3–6), extrapolated values, standard errors, and correlation coefficients (Table S3); results and experimental data for the CORE65 benchmark set for evGW₀ + aZORA and evGW₀ + aZORA + RC (Table S4); non-relativistic ΔSCF results and pre-corrected evGW₀ energies (Table S5); difference between the pre- and post-corrected schemes for the CORE65 benchmark set (Fig. S2); and tabulated values of Fig. 2(b) (Table S6).

In pursuit of open materials science,¹³⁶ we made the results of all relevant calculations available on the Novel Materials Discovery (NOMAD) repository.^137,138

The authors would like to thank the CSC–IT Center for Science for providing computational resources. D.G. acknowledges financial support from the Academy of Finland (Grant No. 316168).