An equation-of-motion coupled cluster singles and doubles approach to Auger–Meitner spectra based on the one-center approximation

Abstract

A robust and computationally efficient methodology to compute Auger decay rates is presented that combines equation-of-motion coupled cluster singles and doubles two-particle Auger density matrices (also known as two-particle Dyson matrices) with precalculated bound-continuum integrals from atomic calculations, known as the one-center approximation. Illustrative applications include KLL Auger electron spectra (AES) of several small and medium-sized molecules.

---

1 Introduction

Absorption of X-ray radiation by an atom or molecule leads to the formation of highly energetic and short-lived core-hole states. Such states can relax through two primary mechanisms: radiative decay (X-ray fluorescence) or non-radiative Auger electron emission.^1–3 For molecules containing atoms from the first three rows of the periodic table, Auger decay is the predominant relaxation channel. In this process, a valence electron fills the core-hole, while another electron is emitted into the continuum, carrying away the excess energy. Auger Electron Spectroscopy (AES) measures the intensity and energy distribution of these emitted electrons, providing a map of the electronic states of the system that are accessible during the decay process.⁴ Due to the atomic specificity of X-ray absorption – where selective excitation of a core electron at a particular atomic site is possible – AES has been used for decades as a powerful tool for probing the chemical environment of an atom within a molecule in gas phase and in materials science.^5–10

More recently, growing interest in the properties of core-hole states and the applications of AES has been driven by the rapid development of novel laser and accelerator-based light sources, such as X-ray free electron lasers.^11–13 These advancements enable an entirely new class of experimental studies, including real-time tracking of the evolution of the molecular system in electronically excited states.^14–21

The most commonly applied theoretical description of Auger decay was originally provided by Wentzel,²² who used perturbation theory and Fermi's golden rule to calculate transition rates into the continuum for a system with two active electrons. A major assumption in this treatment is that Auger decay is a two-step process, where the emission of the Auger electron is decoupled from the preceding core-hole creation by X-ray photoionization or photoexcitation. As a result, the initial state in the Auger decay is considered to be a metastable excited state, which undergoes relaxation through autoionization that can be described in a time-independent manner.²³

The two major difficulties in the theoretical modeling of AES arise from the need to include the state of the emitted electron in the calculations, and the proper inclusion of electronic correlation.^24–27 The description of the emitted electron presents a significant challenge, as the vast majority of quantum chemical methods are designed to describe electrons in molecules within the bound domain, rather than in the continuum. Therefore, standard methods of electronic structure theory, which rely solely on L²-integrable functions, cannot be directly applied to describe Auger decay.

Similarly, accurate modeling of AES spectral features critically depends on accounting for electronic correlation, since the high energies of core-hole states can produce hundreds of distinct decay channels with stable products. These stable products often correspond to states involving multiple electronic excitations with several closely lying electronic configurations, which may require the use of highly-correlated or multi-reference methods for an accurate description.

Over the years, multiple approaches have been proposed to incorporate the state of the unbound electron in the modeling of Auger decay. The most rigorous approach to describing the emitted electron involves solving the scattering problem by imposing proper asymptotic boundary conditions through the Schrödinger equation or the Lippmann–Schwinger equation.²⁸ This has been demonstrated using both numerical integration of the continuum wavefunction,^29,30 and expansion in a finite set of Gaussian-type orbitals.³¹ However, a significant drawback of this approach is the very slow convergence of the continuum orbital expansion, which significantly increases the computational cost and makes it feasible only for atoms or the smallest molecules.

In principle, it is also possible to model the coupling between bound and continuum states in Auger decay using purely L²-integrable functions. This can be achieved through the Stieltjes imaging procedure, which provides a systematic approach for approximating the continuum by a set of discrete states obtained through standard diagonalization of the Hamiltonian in a finite basis set.^32–34 Alternatively, techniques from non-Hermitian quantum mechanics, such as complex scaling or the Complex Basis Function method, can be employed to compute Auger decay widths using only square-integrable functions.^35,36 However, all such approaches, where the continuum states are treated implicitly, require the use of specially crafted orbital basis sets, which may significantly limit their applicability. Furthermore, since these methods are not directly rooted in scattering theory, additional modeling is necessary to extract partial Auger decay widths to obtain direct observables for AES.³⁷

Since the early ab initio studies of Auger decay, it has been recognized that core orbitals are highly localized on specific atomic centers. This implies that only the portions of the electronic densities from valence orbitals that are situated in close proximity to the atom bearing the core hole contribute significantly to the Auger process. This observation forms the basic justification behind the one-center approximation (OCA) method for modeling Auger decay in polyatomic molecules.^38–43 In the OCA approach, the two-electron Coulomb matrix elements, which determine the partial Auger decay rates, are calculated under the assumption that the process is strictly localized on a specific atom. Consequently, only orbitals centered on that atom contribute to the decay rate. This allows to utilize pre-tabulated Coulomb integrals corresponding to pure atomic Auger rates to estimate decay rates in complex molecules. Therefore, the OCA avoids the need to explicitly compute the electron scattering wavefunction for every molecular system. Effectively, the OCA method turns a multi-center problem into a pure single-center (atomic) problem, at the same time drastically reducing the computational cost.

The use of the OCA approximation for modeling Auger decay rates has a long history, in particular in connection with Green's function and configuration interaction approaches for the bound states.^38–44 Recently, novel implementations of the OCA have been presented, that incorporate highly correlated electronic structure methods within modern software for quantum-chemical calculations. Specifically, methods combining the OCA with correlated bound-domain electronic densities derived from multi-state restricted active space perturbation theory (MS-RASPT2),^45–48 multiconfiguration self-consistent field (MCSCF), or multiconfiguration pair-density functional theory (MP-DFT) have been reported.⁴⁹ These implementations demonstrate the effectiveness of the OCA approach in accurately reproducing various observables in Auger electron spectra for small to medium-sized organic molecules such as pyrimidine and oxazole.

In this study, we aim to further advance these developments by proposing a method that integrates OCA with bound-domain wave functions described by the equation-of-motion coupled-cluster singles and doubles (EOM-CCSD) approach.^50–52 Both CCSD and its excited-state counterpart, EOM-CCSD, offer distinct advantages when compared to other classes of excited-state electronic structure methods. The single-reference-based EOM-CCSD provides an accurate description of multi-configurational states while maintaining a simple computational framework and systematic and controllable treatment of electron correlation.

A few years ago, the EOM-CCSD approach was applied to model Auger spectra using simplified continuum orbital representations assumed as either plane waves or Coulomb waves.^53,54 Although results for certain molecules appeared promising, simplified continuum orbital models introduced significant inaccuracies in predicting Auger decay partial widths, especially concerning the ratio of singlet to triplet decay channels. In the present work, we demonstrate that a straightforward solution to this issue can be achieved by incorporating pre-calculated continuum atomic integrals within the OCA framework. Therefore, our newly introduced methodology, termed OCA-EOM-CCSD, combines the simplicity and effectiveness of both the OCA technique and the EOM-CCSD wave function.

There is certainly a need for new and robust computational methods capable of first-principle modeling Auger decay spectra in complex molecules. Typically studied polyatomic molecules can easily possess a large number of possible decay channels, which surpasses the capability of approaches relying either on direct calculation of continuum wave functions or the explicit definition of an active space.^29,30 The OCA-EOM-CCSD method circumvents both these limitations, providing a computationally efficient framework suitable for a large set of molecules.

The paper is structured as follows: In the next section, we review the theoretical frameworks of both the OCA approach and the EOM-CCSD wave function. Then, we give the details of the methodology and calculations. Subsequently, we present computational illustrations of the OCA-EOM-CCSD method in the application to KLL Auger spectra for a set of benchmark molecules. Finally, we summarize our findings and conclusions in the last section.

2. Theoretical background

Below we outline the theoretical framework of the OCA-EOM-CCSD methodology using normal (non-resonant) Auger decay as an example. However, this framework is fully applicable to resonant Auger emission as well.

Auger decay is an autoionization process, which can be represented as

Ψ^N_I → Ψ^N−1_F + e⁻ (1)

where Ψ^N_I denotes the N-electron initial state with a core vacancy, and Ψ^N−1_F represents the (N − 1) -electron final state (channel) of the stable product. The emitted Auger electron is characterized by its momentum k⃑ and energy E_k = E_I − E_F, as required by energy conservation. Following the perturbative approach of Wentzel,²² the total probability rate for decay from the initial state Ψ^N_I into the final state Ψ^N−1_F is given by:where Ĥ_el is the total electronic Hamiltonian of the system, and is the wave function of the core ionized state with the incoming asymptotic boundary condition. Within the single-channel approximation, the wave function can be represented as an antisymmetrized product of the bound-state wave function Ψ^N−1_F and the continuum electron state:We do not explicitly consider spin adaptation of the wave function here. For the spin-free electronic Hamiltonian Ĥ_el and the spin-adapted Ψ^N−1_F, it is sufficient to include a spin-degeneracy factor g_α in the final expression for Γ_IF.

In this work, the bound-domain wave functions in the initial state Ψ^N_I and the final state Ψ^N−1_F are obtained using the EOM-CCSD method, fully analogous to the approach presented in ref. 53 and 54. As demonstrated in those references, the initial state wave function Ψ^N_I can be conveniently computed within the fc-CVS-EOM-CCSD framework,⁵⁵ employing the core-valence separation (CVS) technique^56,57 combined with a frozen-core approximation. Assuming that the continuum orbital ϕ_k⃑ is orthogonal to the bound-domain orbitals and the frozen core is applied, the only contribution from the Ĥ_el − E_I coupling to the decay widths Γ arises from the two-electron interaction 1/r_ij in Ĥ_el. Consequently, the transition probability for the Auger decay can be expressed through the two-particle Dyson matrix elements R_IF;pqr:⁵⁸

where 〈pq||kr〉 are mixed two-electron integrals between the molecular orbitals from the bound domain (p,q,r) and the continuum orbital ϕ_k⃑ describing the emitted Auger electron. The two-body Dyson matrices R_IF;pqr connect the bound components of the many-body wave function between the initial and final states:Since we intend to employ the EOM-CCSD method to obtain the tensor R_IF;pqr, we cannot assume that it is Hermitian conjugate. Therefore, both tensors R_IF;pqr and R_FI;pqr must be calculated separately. The explicit expressions for calculating the matrices R_IF;pqr using the fc-CVS-EOM-CCSD wave function can be found in ref. 53. In that study the continuum orbital ϕ_k⃑ was approximated as either a plane wave or a Coulomb wave. Consequently, all mixed integrals 〈pq||kr〉 had to be explicitly evaluated, followed by the angular integration over k̂. In this work, these costly computational steps will be entirely eliminated by applying the one-center approximation.

Within the OCA, it is assumed that the only significant contributions to the integrals 〈pq||kr〉 arise from orbitals localized on atom A, the site of the initial core hole. The molecular orbitals {ϕ_p} are projected onto an auxiliary basis set {χ_i}, known as the Minimal Basis Set (MBS), which simplifies the evaluation of two-electron integrals by including only atomic orbitals essential for describing core and valence electrons. Consequently, the problem of evaluating integrals 〈pq||kr〉 becomes strictly single-centered. Hence, it is convenient to transform the integration over angular coordinates k̂ into a summation over partial waves:

where ϕ_E,lm denotes the continuum orbital characterized by energy and angular momentum quantum numbers l and m.

To calculate 〈pq||E_lmr〉 within the OCA, we follow the procedure presented in ref. 45. In the first step, the molecular orbitals {ϕ_p} are projected onto the space spanned by the MBS:

The transformation matrix D_μp is defined as:^42,45

D = T⁻¹UC. (9)

Here, T represents the overlap matrix of the MBS, U is the overlap between the MBS and the original atomic basis set {f_κ}, and C is the molecular orbital coefficient matrix expressed in the original atomic basis set {f_κ}.

Within the OCA, the two-electron integrals are approximated as:⁴²

Here, the integrals 〈χ^A_μχ^A_ν||χ^A_E,lmχ^A_ρ〉 involve only occupied atomic orbitals (in our case 1s, 2s and 2p) localized on the atom A that bears the core hole.

A second assumption within the OCA is that the continuum orbital χ^A_E,lm remains unchanged by the molecular environment compared to an isolated atom A. Therefore, the integrals 〈χ^A_μχ^A_ν||χ^A_E,lmχ^A_ρ〉 can be computed once for each specific atom and then tabulated. Due to atomic spherical symmetry, only radial integrals need to be explicitly computed, which are subsequently combined with precomputed angular integration coefficients, as outlined in ref. 45. This approach removes the necessity for explicit evaluation of continuum orbital and corresponding hybrid integrals, greatly simplifying the computational procedure.

3. Methodological and computational details

We illustrate our implementation of the OCA-EOM-CCSD methodology by simulating normal Auger decay spectra corresponding to the KLL transitions in a set of benchmark molecules: H₂O, HNCO, CO₂, benzene, and oxazole. In normal Auger decay, the initial state Ψ^N_I is a singly charged cation generated by ionization from a core orbital, and the final states Ψ^N−1_F correspond to doubly charged cations. The initial states Ψ^N_I were obtained using the CVS-EOM-IP-CCSD variant of the equation-of-motion coupled-cluster theory with the frozen-core approximation,⁵⁵ which allows for single ionization from the core shell. The final states Ψ^N−1_F were computed using the double-ionization variant of EOM-CCSD (EOM-DIP-CCSD),^59,60 which allows the removal of two electrons from valence orbitals. The method bears similarities with other 2h-type propagator/Green's function approaches which have been used in the context of Auger spectroscopy, see for instance ref. 40, 41, 44 and 61.

We emphasize that in our approach, the initial and final states, obtained through the CVS-EOM-IP-CCSD and EOM-DIP-CCSD methods, are computed using a common set of canonical molecular orbitals. As a result, any relaxation effects in the target states are included implicitly through the higher-order EOM operators. The accuracy of this strategy for core ionization energies has been tested previously,^55,62–64 and its effectiveness stems from the fact that the studied core-ionized states are dominated by one-electron transitions, for which EOM-CCSD is sufficient to capture the essential electron correlation and relaxation effects.

The frozen core space was defined to include the necessary 1s orbitals of carbon, nitrogen, and oxygen atoms. For all atoms, we used the Dunning correlation-consistent basis set cc-pVTZ. The molecular geometries were taken from previously published studies,^46,54,65 and their Cartesian coordinates are provided in the SI.

Our implementation of the OCA method does not impose any constraints on the choice of the minimal basis set (MBS). Therefore, we tested two different choices of MBS: one based on the subset of the cc-pVTZ basis set which is contracted as to represent the atomic orbitals 1s, 2s and 2p (as done in ref. 45), and another one using the standard STO-3G basis set. We refer to ref. 42 for a quantitative criterion to assess the accuracy of representing the valence molecular orbital space with the selected minimal basis set, and to Section S3 in the SI for a comparison of such measures for the two choices of MBS applied here. We employed precalculated atomic integrals with continuum orbitals as reported in ref. 66. For a meaningful comparison with experimental data, the discrete stick spectra were convoluted using a uniform Gaussian broadening with full-width-at-half-maximum (FWHM) values chosen between 0.7 and 1.4 eV. This procedure approximately accounts for the average vibrational broadening observed in the measured spectra. The experimental Auger spectra shown in the figures below were generated with the help of software for digitization.⁶⁷ The ab initio energies of the core-ionized and final dicationic states were used without the application of any empirical shifts. In the following, we report the energetics of the AES in two ways: either as the energy of the emitted electron (see eqn (1)), which is the difference between the initial core-ionization energy and the energies of the doubly ionized states:

AE_IF = E_I(CVS-EOM-IP) − E_F(EOM-DIP), (11)

or directly using the binding energies of the final doubly ionized states:

BE_F = E_F(EOM-DIP). (12)

While the first version is advantageous for direct comparison with experimental data, the second choice is more convenient in the context of molecules with multiple K-edges, allowing one to clearly identify how different final states contribute to each given edge. The conversion between theoretical AEs and BEs involves core-ionization potentials, and the present CVS-EOM-IP values for each studied molecule are provided in the SI (Table S1).

All calculations were performed using a developer version of Q-Chem.⁶⁸

4. Results

4.1. Water

The normal Auger electron spectrum of the singly ionized water molecule has been extensively studied, both experimentally and theoretically.^{24,32,38,69–71} Despite its apparent simplicity, it has been shown that modeling the Auger spectrum of H₂O poses significant challenges for ab initio methods. The most sophisticated calculations, conducted by Inhester et al.,⁷⁰ revealed the importance of accounting for both electronic correlation and nuclear dynamics effects to correctly reproduce the main features of the experimental Auger spectrum.

Fig. 1 presents our simulated spectrum based on the OCA-EOM-CCSD methodology, along with the experimental spectrum and theoretical spectra obtained using two other representations of the continuum orbital combined with the EOM-CCSD wave functions. In all theoretical cases, the continuous spectra were generated by applying a uniform Gaussian broadening of 1 eV to the discrete stick spectra. The full OCA-EOM-CCSD spectrum was constructed based on a total of 188 dicationic states (75 singlet and 113 triplet states), obtained through the EOM-DIP-CCSD method. The exact Auger decay widths for all the 2h dominated channels are provided in the SI (Table S2).

Fig. 1 H₂O. Experimental⁶⁹ and computed normal Auger spectra of a singly ionized molecule. Theoretical curves were obtained from the stick spectra by applying a constant Gaussian broadening of 1.0 eV (FWHM). Gray sticks correspond to singlet decay channels, and blue sticks to triplet channels. The spectra labeled as EOM-CCSD/PW and EOM-CCSD/CW are based on plane wave and Coulomb wave models for the continuum electron, respectively, as described in ref. 54. The two OCA-EOM-CCSD results differ in the choice of minimal basis set, using either STO-3G or a subset of the cc-pVTZ basis set.

The experimental spectrum shows significant broadening due to nuclear motion and can be roughly divided into three segments:

1. From 500 eV to 485 eV, where two electrons are removed from the outer valence orbitals (3a₁, 1b₁, 1b₂);

2. From 485 eV to 460 eV, where one electron is removed from the outer valence shell and one from the inner valence orbital (2a₁);

3. From 460 eV to 450 eV, where both electrons are removed from the inner valence shell.

The highest intensity is observed in the first energy region.

Two variants of OCA-EOM-CCSD calculations were performed: one based on the STO-3G minimal basis set, and the other on a subset of the cc-pVTZ basis set (labelled “subVTZ”). Inspection of Fig. 1 clearly shows that the present OCA-EOM-CCSD model successfully reproduces the main features of the experimental spectrum, despite the neglect of explicit nuclear dynamics and the use of a fixed broadening of 1 eV. As apparent from Fig. 1, the choice of MBS does not visibly affect the spectral profile.

A direct comparison should be made with the results presented in ref. 54, where the same bound-domain electronic transition densities (i.e., two-body Dyson functions) and channel energies (based on EOM-CCSD calculations) were used, but the continuum electron was modelled explicitly using either a plane wave (PW) or a Coulomb wave (CW). Fig. 1 shows that the current OCA-EOM-CCSD model outperforms both the PW and CW models and is carried out at a significantly lower computational cost. The PW model, in particular, tends to overestimate contributions from triplet decay channels. Although this faulty behaviour is greatly improved by the CW model, the predicted intensities for different singlet channels in the CW model remain somewhat unsatisfactory, in particular for the ¹A₁ (2a₁⁻¹3a₁⁻¹) and ¹A₁ (2a₁⁻²) channels. For these two channel, the CW modes yields higher intensities than for the reference channel ¹A₁ (1b₁⁻²).

A more quantitative comparison of these observations is presented in Table 1, where we report the relative partial widths, normalized to the dominant ¹A₁ (2b₁⁻²) channel. We compare our OCA-EOM-CCSD results with those obtained using the PW and CW models, as well as with OCA calculations based on RASSCF/RASPT2 wave functions from the literature.⁴⁵ The results of Inhester et al.,⁷⁰ obtained using MRCI electronic structure calculations combined with a numerically solved coupled radial Schrödinger equations for the continuum electron, are given in the last column. The representation of the continuum electron of Inhester et al. can be considered as nearly exact and thus serves as a benchmark here.

Table 1 Relative partial widths Γ_i (with respect to the first ¹A₁ channel) for the main channels in the AES of H₂O (1s⁻¹), obtained using different ab initio methods. PW and CW refer to the plane wave and Coulomb wave models, respectively, as described in ref. 54. OCA(subVTZ) and OCA(STO-3G) represent the present results obtained with the OCA-EOM-CCSD method using different choices of the minimal basis set (MBS). The last column contains benchmark results based on a numerical continuum orbital and MRCI bound-domain electronic structure calculations.⁷⁰ The total Auger widths Γ_tot are provided in the last row. Numbers in parentheses are the kinetic energies of the Auger electron (in eV) for the given channel

Channel	EOM-CCSD				RASSCF/PT2 OCA(subVTZ)^45	MRCI numerical^70
	PW	CW	OCA(subVTZ)	OCA(STO-3G)
^3B1 (3a1^−11b1^−1)	19	4	3	4 (499.54)	1 (499.93)	2 (500.67)
^1A1 (1b1^−2)	100	100	100	100 (498.54)	100 (498.65)	100 (499.39)
^1B1 (3a1^−11b1^−1)	105	95	98	94 (497.09)	98 (497.33)	95 (497.98)
^3A2 (1b1^−11b2^−1)	0	0	0	0 (495.57)	0 (495.64)	0 (496.60)
^1A1 (3a1^−2)	63	67	70	63 (494.00)	70 (493.86)	69 (494.64)
^1A2 (1b1^−11b2^−1)	94	80	76	68 (493.79)	86 (493.95)	80 (494.68)
^3B2 (3a1^−11b2^−1)	14	3	2	2 (493.70)	1 (493.82)	1 (494.63)
^1B2 (3a1^−11b2^−1)	77	71	67	58 (491.60)	73 (491.61)	69 (492.36)
^1A1 (1b2^−2)	39	53	41	34 (486.87)	52 (486.54)	52 (487.45)
^3B1 (2a1^−11b1^−1)	140	31	12	12 (481.00)	14 (481.78)	16 (482.30)
^3A1 (2a1^−13a1^−1)	132	29	13	12 (479.25)	10 (480.73)	14 (480.58)
^3B2 (2a1^−11b2^−1)	103	22	8	8 (475.37)	6 (477.14)	8 (476.82)
^1B1 (2a1^−11b1^−1)	9	71	56	56 (474.24)	15 (476.37)	52 (475.76)
^1A1 (2a1^−13a1^−1)	20	102	63	66 (472.98)	37 (473.54)	58 (473.27)
^1B2 (2a1^−11b2^−1)	8	47	33	29 (468.10)	11 (469.26)	35 (468.75)
^1A1 (2a1^−2)	53	115	34	34 (454.49)	16 (456.21)	21 (457.19)
Γtot/meV	175.1	121.7	199.3	178.8	180.4	145.6

---

Inspection of Table 1 shows that the OCA-EOM-CCSD results agree very well with the benchmark values, particularly for the triplet channels, where the relative widths match almost exactly. A slight disagreement is observed for the singlet channel involving Auger electron emission from the ¹A₁ (2a₁⁻²) state, which appears overestimated in our OCA-EOM-CCSD calculations. This discrepancy is likely attributed to electron correlation effects in the bound states calculations rather than inaccuracies in the continuum electron description.

Comparison between the current OCA-EOM-CCSD results and the earlier OCA-RASSCF-based results also shows strong consistency. The relative widths agree very well, except for singlet channels that involve the removal of one outer-valence electron and one inner-valence electron, where the OCA-EOM-CCSD relative widths are visibly higher. This region of the spectrum is, however, particularly challenging for ab initio modeling, as it is influenced by shake-up and shake-off transitions⁷² as well as by interchannel coupling effects in the continuum,^73,74 which are neglected in the present independent-channel model.

Table 1 also includes the kinetic energies of the emitted electron for each channel AE_IF, calculated at different levels of theory: EOM-CCSD, RASPT2, and MRCI. For the most dominant channels, the computed energies differ by less than 1 eV, which is comparable to the estimated vibrational broadening in the experimental spectra. However, there are also cases where the differences reach up to 2.7 eV, such as for the ¹A₁ channel involving double ionization from the inner-valence 2a₁ orbital. In this case, the MRCI energy appears to most closely match the experimental peak in that region, while the EOM-CCSD energy tends to visibly underestimate the value (i.e., overestimating the EOM-DIP energy for that channel). This discrepancy may be attributed to the multi-reference character of the ¹A₁ (2a₁⁻²) state, which is not fully captured at the EOM-CCSD level and would require going beyond the double-excitation manifold in the EOM/coupled-cluster operator. On the other hand, multi-reference methods such as RASPT2 and MRCI—with properly selected active spaces—describe the nature of such states more accurately. Nevertheless, it should be noted that selecting appropriate active spaces in multi-reference calculations becomes increasingly challenging for larger molecules.

Finally, we also compared the total Auger decay widths obtained with the different methodologies. The total widths are reported at the end of Table 1. The best experimental estimate for the total width of singly core-ionized water is (160 ± 5) meV.⁷⁵ Our current OCA-EOM-CCSD results predict total widths of 199 meV and 179 meV, depending on the choice of MBS. These values are consistent with the OCA-RASSCF result (180 meV), but notable larger than both the best experimental estimates and the benchmark values from ref. 70. The observed 20 meV difference in the total Auger decay widths between the two minimal basis set choices suggests that, although the relative partial widths remain robust, the absolute total widths predicted by OCA-based methods may vary quite significantly. This effect likely originates from differences in the overlap of the minimal basis set with the full orbital basis set used for the bound-domain calculations. Nonetheless, the relative partial widths remain largely unaffected by the choice of minimal basis set, as demonstrated in Table 1.

4.2. Isocyanic acid

Isocyanic acid (HNCO) is a particularly appealing molecule for benchmarking Auger spectra theoretical models because it is the simplest molecule containing the three most relevant elements for organic molecules: carbon, nitrogen, and oxygen. Moreover, highly accurate experimental data for both normal and resonant Auger processes have recently been reported.⁷⁶ We applied our OCA-EOM-CCSD methodology, combined with the STO-3G minimal basis set, to simulate Auger spectra following ionization at three different K-edges. Fig. 2 and 3 present the calculated Auger spectra alongside the corresponding experimental ones. The theoretical spectral profiles were obtained by applying a uniform Gaussian broadening of 0.8 eV. For ease of comparison, the spectra are plotted with two scales: the upper scale corresponds to the Auger electron kinetic energy (AE), and the lower scale corresponds to the binding energy of the resulting dicationic final states (BE). Fig. S2 in the SI shows a comparison of the OCA-EOM-CCSD spectra at the carbon K-edge using either STO-3G or subVTZ as MBS.

Fig. 2 HNCO. Experimental⁷⁶ and computed normal Auger spectra for the carbon K-edge (top) and nitrogen K-edge (bottom) obtained using the OCA(STO-3G)-EOM-CCSD method. The theoretical curves were generated from the stick spectra by applying a constant Gaussian broadening of 0.8 eV (FWHM). Gray sticks correspond to singlet decay channels, and blue sticks to triplet decay channels.

Fig. 3 HNCO. Experimental⁷⁶ and computed normal Auger spectrum for the oxygen K-edge, obtained using the OCA(STO-3G)-EOM-CCSD method. The theoretical curve was generated from the stick spectrum by applying a constant Gaussian broadening of 0.8 eV (FWHM). Gray sticks correspond to singlet decay channels, and blue sticks to triplet decay channels.

The total spectra were generated by including 450 final dicationic states (290 singlet states and 160 triplet states), covering a binding energy range of approximately 40 eV (from about 34 eV to 74 eV). Detailed energies and widths for the most dominant decay channels are provided in the SI (Table S3). Similarly to the case of water, our calculations predict that the Auger spectra at all edges are dominated by singlet decay channels. The calculated total widths arising from singlet channels are 34.6 meV, 102.6 meV, and 128.4 meV for the C, N, and O K-edges, respectively. For the triplet channels, the corresponding values are 2.1 meV, 5.7 meV, and 9.8 meV, respectively.

The Auger spectra of HNCO display a rich structure with multiple distinct peaks. Although it can be compared with isoelectronic CO₂, the lower symmetry of HNCO (C_s point group) results in a greater number of distinguishable transitions. When comparing the spectra for the three distinct K-edges, it is clear that their shapes vary significantly depending on which atomic site is core ionized. Notably, the nitrogen K-edge spectrum exhibits its maximum peak at the lowest binding energies (corresponding to the highest Auger electron energies), with the highest intensities arising from the first two HNCO²⁺ states: ¹A′ [(2a′′)⁻²] and ¹A′′ [(9a′)⁻¹ (2a′′)⁻¹)]. In contrast, the carbon and oxygen K-edge spectra show their maximum intensities at intermediate binding energies, between 45 and 50 eV, with the dominant contributions originating from channels with main configurations ¹A′ [(8a′)⁻² +(1a′′)⁻²] and ¹A′′ [(8a′)⁻¹ (1a′′)⁻¹)]. This behavior reflects the localized nature of the Auger decay process and the distribution of valence orbitals across the constituent atoms.

When evaluating the performance of our OCA-EOM-CCSD method, it must be noted that in all three cases the agreement with experimental spectra is satisfactory. The carbon K-edge spectrum is particularly well reproduced, and the relative intensities of all peaks are reasonably well captured in the nitrogen K-edge spectrum as well. For the oxygen K-edge, even though the major features are reproduced, a noticeable discrepancy appears in the binding energy range 40–45 eV. Overall, our results, obtained without introducing any empirical energy scaling, are in fair agreement with the more accurate results by CASCI calculations in ref. 76, which explicitly account for vibrational broadening of the spectra through a model based on moment theory.

The oxygen K-edge spectrum of HNCO has also recently been calculated using the OCA-RASSCF/RASPT2 method.⁴⁷ While the RASSCF/RASPT2 approach offers visibly better accuracy compared to our OCA-EOM-CCSD results for the oxygen edge, it was limited to a smaller portion of the spectrum and it was not demonstrated how the methodology would perform for the carbon or nitrogen K-edges.

4.3. Carbon dioxide

CO₂ is another molecule frequently used as a benchmark for testing different methodologies to calculate Auger electron spectra.^26,77,78 In contrast to HNCO, CO₂ is linear and possesses higher symmetry, which results in some of its electronic states being doubly degenerate. This makes it a computationally more tractable system for theoretical modeling as one can take advantage of its higher symmetry.

The carbon K-edge Auger spectrum calculated using the present OCA(STO-3G)-EOM-CCSD method is shown in Fig. 4. We refer to Fig. S3 in the SI for a comparison with the OCA(subVTZ)-EOM-CCSD results. The experimental spectrum closely resembles that of HNCO at the carbon site. However, it can be observed that the main peaks in CO₂ are somewhat sharper, with less vibrational broadening. Our calculated spectrum reproduces the experimental features very well, similar to the case of HNCO. The theoretical spectrum was obtained by including contributions from 260 dicationic final states in total (120 singlet and 140 triplet states). The continuous curve was generated by applying a constant Gaussian broadening of 0.7 eV to the stick spectrum. The dominant contributions arise from dicationic states with binding energies in the range of 45–50 eV, primarily involving vacancies in the valence orbitals 1π_u and 3σ_u. The two most intense peaks in the calculated spectrum correspond to the states ¹Δ_g(1π_u⁻²) and ¹Σ_g⁺(1π_u⁻²), with binding energies of 48.02 and 48.82 eV, respectively.

Fig. 4 CO₂. Experimental⁶⁹ and computed Auger spectra at the carbon K-edge. The bottom panel shows our results based on the OCA(STO-3G)-EOM-CCSD method. The two middle panels display spectra computed using explicit representations of the continuum orbital – either a plane wave (PW) or a Coulomb wave (CW) – combined with the same EOM-CCSD electronic structure method, as shown in ref. 54. All theoretical spectral profiles were generated from stick spectra by applying a constant Gaussian broadening of 0.7 eV (FWHM). Gray sticks correspond to singlet decay channels, and blue sticks to triplet decay channels. See Fig. S3 in the SI for a comparison of the OCA(STO-3G) and OCA(subVTZ) results.

Similar to the case of H₂O, we also compare our results with those obtained using explicit continuum orbital representations, namely the plane wave (PW) and Coulomb wave (CW) models, in combination with the same EOM-CCSD densities in the bound domain. As with water, the OCA-EOM-CCSD method outperforms both the PW and CW models. While both simplified models can qualitatively reproduce the main spectral features, the PW model tends to significantly overestimate contributions from triplet decay channels. The CW model improves upon this, but shows exaggerated intensities from deeper bound states at binding energies around 65 eV.

Let us now consider the oxygen K-edge Auger spectrum. The present OCA-EOM-CCSD results, compared with the digitized experimental spectrum,⁶⁹ are shown in Fig. 5. This case is different in that it involves equivalent core sites, as the two oxygen atoms in CO₂ are symmetry-equivalent. As a result, in standard calculations, the molecular orbitals describing O(1s) electrons are delocalized over both centers. Within the OCA framework, one must account for this symmetry in order to obtain a meaningful physical result. In ref. 45, this was done by applying a Cholesky localization procedure to localize the two core orbitals, effectively reducing the symmetry of the molecule (from D_2h to C_2v). Here, we adopt a different strategy to address the presence of equivalent atoms. Our calculations were carried out using the full Abelian point group symmetry, and the final Auger width for each distinct decay channel Ψ^N−1_F was obtained by summing the contributions from all symmetry-equivalent core-ionized initial states:

where the summation over I_A includes all core-ionized initial states Ψ^N_I corresponding to ionization from symmetry-equivalent atomic sites. In the case of CO₂, this summation contains two terms, as there are two symmetry-equivalent core-ionized states at the oxygen atoms. The OCA atomic integrals in this approach are associated with one of the equivalent atoms A_i – the specific choice does not affect the result. We verified numerically that this symmetry-preserving procedure yields results equivalent to those obtained using orbital localization combined with symmetry reduction.

Fig. 5 CO₂. Experimental⁶⁹ and computed Auger spectra at the oxygen K-edge. The bottom panel shows the present calculation based on the OCA(STO-3G)-EOM-CCSD method. Theoretical curve was generated from stick spectra by applying a constant Gaussian broadening of 0.7 eV (FWHM). Gray sticks correspond to singlet decay channels, while blue sticks correspond to triplet decay channels.

The calculated OCA-EOM-CCSD spectrum shown in Fig. 5 shows good agreement with the experimental data. Once again, the spectrum is dominated by singlet decay channels (total width of 146.8 meV versus 10 meV for triplet dicationic states). The most intense peaks arise from states with binding energies around 50 eV. A noticeable discrepancy is observed at the high-binding-energy tail of the spectrum, where the position of experimental peak at Auger kinetic energy of approximately 472 eV is understimated by about 4–5 eV in our calculations and is more intense than observed. Similarly to H₂O, this requires more accurate treatment of the dicationic states in that energy region.

4.4. Benzene

Core-level spectroscopy of the benzene molecule has been the subject of numerous experimental and theoretical studies (see, e.g., ref. 41, 65, 79 and 80 for examples related to Auger) because of its unique importance in organic chemistry and the wide range of derivatives with multiple applications. Recently, the normal Auger spectrum of benzene was revisited in a thorough theoretical study by Jayadev et al.⁶⁵ Benzene is a particularly interesting system due to its high symmetry and the presence of six equivalent carbon atoms. This means that the carbon K-edge Auger spectrum can originate from six different core-ionized initial states, with symmetries ²A_1g, ²E_1u, ²E_2g, ²B_1u (with the two E states being doubly degenerate).

As in the case of the oxygen K-edge spectrum for CO₂, in our OCA-EOM-CCSD calculations, we used the full Abelian point group symmetry of C₆H₆ (D_2h), and performed the summation for each individual EOM-DIP channel in eqn (13) over all six possible initial states. Our calculated spectrum, together with the most recent experimental data, is shown in Fig. 6. The continuous theoretical spectrum was obtained by applying a constant Gaussian broadening of 0.8 eV and includes contributions from a total of 520 final dicationic states (230 singlet and 290 triplet states).

Fig. 6 Benzene. Experimental⁸⁰ and computed Auger spectra. The bottom panel shows the present OCA(STO-3G)-EOM-CCSD results. We refer to Fig. S4 in the SI for a comparison with the OCA(subVTZ)-EOM-CCSD results. Theoretical results digitized from ref. 65, obtained using EOM-CCSD combined with the plane-wave model and the complex basis function approach (CBF-ΔCCSD) are shown in the middle panel. Gray sticks correspond to singlet channels, and blue sticks to triplet channels. The theoretical spectral profile (red curve) was obtained from the stick spectra by applying a constant Gaussian broadening of 0.8 eV (FWHM).

In Fig. 6, we also include the theoretical spectrum reported by Jayadev et al.,⁶⁵ obtained using the EOM-CCSD method combined with a PW model for the continuum electron, as well as the complex basis function approach based on the CCSD wave function. A visual comparison of Fig. 6 shows that, while all three theoretical methods reproduce the general shape of the spectrum reasonably well, the present OCA-EOM-CCSD spectrum agrees remarkably well with the experimental data and clearly outperforms the results from ref. 65. This is particularly evident in the region corresponding to the highest Auger electron energies, which arises from the population of the most weakly bound dicationic states, predominantly involving vacancies in the HOMO (1e_1g), HOMO−1 (3e_2g), and HOMO−2 (1a_2u) orbitals of benzene. In the binding energy range between 25 and 32 eV, three well-pronounced peaks are observed in the spectrum. The first peak is dominated by transitions to the ¹E_2g (1e_1g⁻²) and ¹A_1g (1e_1g⁻²) states, with binding energies of 25.67 and 26.12 eV, respectively. The second peak arises mainly from a cluster of states with symmetries ¹B_1g, ¹E_1g, and ¹B_2g, and dominant configurations (3e_2g⁻¹1e_1g⁻¹), with a total binding energy around 28.3 eV. The third peak is dominated by a single transition to the ¹E_1u (1a_2u⁻¹1e_1g⁻¹) state at a binding energy of 30.85 eV.

Importantly, the experimentally observed intensity ratios between these three peaks are accurately reproduced only by the present OCA(STO-3G)-EOM-CCSD method and not by the methods tested in ref. 65.

Since benzene is the largest molecule in our benchmark set, it also serves as a good case for comparing the efficiency of the OCA approach with the PW and CW continuum orbitals models, as discussed above. To this end, we performed comparative calculations on a cluster node using 24 threads equipped with modern Intel Xeon Platinum processors. With the PW model (settings as reported in ref. 65), the bare calculation of the decay rate for one channel took 17 minutes. Using the CW model (with parameters of the Coulomb wave expansion as given in ref. 54) required 10 hours 57 minutes. In contrast, the OCA model needed only 12 seconds to calculate the same rate. These results clearly show that for the OCA-EOM-CCSD method the computational cost of the rate calculation itself can be considered as negligible compared to the bound-domain steps (i.e., CVS-IP and EOM-DIP calculations), which are common to all continuum models tested here.

Similarly to previous theoretical studies, our calculations also predict that the normal Auger spectrum of benzene is dominated by singlet channels, with calculated total widths of 39.6 meV for singlet versus 2.2 meV for triplet channels.

4.5. Oxazole

The last molecule in our study is 1,3-oxazole, C₃H₃NO. Oxazole is a particularly interesting molecule (structure and atom numbering shown in the SI, Fig. S1): it is planar, with C_s point group symmetry, and contains three, non equivalent, carbon atoms. As such, oxazole can serve as an intermediate case between benzene and HNCO, similar in size but with lower symmetry compared to benzene. Recently, core-level spectroscopy of oxazole, including normal Auger electron spectra at the C, N, and O K-edges, was measured and thoroughly analyzed using various theoretical methods by Schnack-Petersen et al.,⁴⁶ where also an excellent discussion of the electronic structure of oxazole was provided, which we will not repeat here.

We begin the discussion of our results with the carbon K-edge. Fig. 7 shows the Auger spectra calculated for the three distinct C(1s⁻¹) ionization states using our present OCA-EOM-CCSD method with STO-3G as the minimal basis set. The continuous spectrum was constructed by including transitions to 78 singlet states and 100 triplet dicationic states. Fig. 8 presents the cumulative theoretical spectrum alongside the experimental data. The carbon K-edge spectrum can be related to that of benzene; however, in contrast to benzene, the oxazole spectrum is broader, with less well-resolved peaks. This broadening likely results from the lower symmetry of the molecule and the higher density of dicationic states.

Fig. 7 1,3-Oxazole. Theoretical Auger electron spectra obtained from OCA(STO-3G)-EOM-CCSD calculations for the three distinct C(1s⁻¹) initial states. The continuous curves were generated from the stick spectra by applying a constant Gaussian broadening of 1.3 eV. Gray sticks correspond to singlet decay channels, and blue sticks to triplet decay channels.

Fig. 8 1,3-Oxazole. Experimental⁴⁶ and computed Auger spectra at the carbon K-edge. The theoretical spectrum was obtained from OCA(STO-3G)-EOM-CCSD calculations and represents the sum of contributions from the three distinct C(1s⁻¹) initial states. The continuous theoretical curve was generated from the stick spectrum by applying a constant Gaussian broadening of 1.3 eV. Gray sticks correspond to singlet decay channels, and blue sticks to triplet decay channels. Fig. S5 in the SI shows the OCA(subVTZ)-EOM-CCSD results.

The C2(1s⁻¹) and C5(1s⁻¹) Auger spectra are quite similar, with the highest intensity arising from transitions to the ¹A′ [(3a′′)⁻²] state at a binding energy of 27.81 eV, dominated by configurations with two vacancies in the HOMO orbital (3a′′). In contrast, the Auger spectrum at carbon C4(1s⁻¹) edge shows its strongest intensity from the ¹A′ [(2a′′)⁻¹ (3a′′)⁻¹] state, at a binding energy of 29.67 eV, which involves single vacancy in the HOMO and HOMO−1 orbitals. It should be noted that the overall shape of the OCA-EOM-CCSD spectrum closely resembles the one obtained with the EOM-CCSD method combined with the PW model,⁴⁶ with major differences only appearing in the high-binding-energy range (40–45 eV). In this region, the PW model predicts substantial contributions from triplet channels, whereas in our OCA-EOM-CCSD calculations, triplet contributions are very minor, as shown in Table 2. This is consistent with our earlier observations for benzene and HNCO.

Table 2 Cumulated Auger widths Γ in meV for different K-edges in 1,3-oxazole for the computed 200 (100 singlet and 100 triplet) final states

Atom	Singlet channels	Triplet channels	Total
C4	36.5	2.3	38.8
C5	38.1	2.1	40.2
C2	30.6	1.9	32.5
N	68.9	2.7	71.6
O	75.2	1.7	76.9

---

Fig. 9 shows our OCA(STO-3G)-EOM-CCSD Auger spectra compared to the experimental measurements for the nitrogen and oxygen K-edges. The theoretical spectra were obtained here by including 200 final states (100 singlet and 100 triplet states). Similarly to the carbon edge, the experimental spectra for the nitrogen and oxygen edges are broad, but exhibit shapes that differ significantly from the carbon K-edge spectrum. For the nitrogen K-edge, the highest intensity is observed in the binding energy range of 31–35 eV, arising from states with ¹A′ and ¹A′′ symmetry and dominant configurations such as (2a′′)⁻², (15a′)⁻², and [(15a′)⁻¹ (2a′′)⁻¹, involving vacancies in the HOMO−1 and HOMO−2 orbitals. For the oxygen K-edge, the maximum intensity occurs at even higher binding energies (lower Auger electron kinetic energy), around 46–48 eV, dominated by states with the (1a′′)⁻² configuration. Again, our calculations predict that triplet channel contributions are very minor in the overall Auger decay, as summarized in Table 2.

Fig. 9 1,3-Oxazole. Experimental⁴⁶ and computed Auger spectra at the nitrogen (left panel) and oxygen (right panel) K-edges. The theoretical spectra were obtained from OCA(STO-3G)-EOM-CCSD calculations. The continuous theoretical curves were generated from the stick spectra by applying constant Gaussian broadenings of 0.7 eV and 1.4 eV, respectively. Gray sticks correspond to singlet decay channels, and blue sticks to triplet decay channels.

Analysis of Fig. 8 and 9 shows that the OCA-EOM-CCSD method predicts the key features of the Auger spectra of oxazole very well. Interestingly, when comparing the performance of the OCA approach with the PW model in ref. 46, the overall agreement with the available experimental reference data is of comparable quality. The main difference is observed in the predicted contributions from triplet channels, which are noticeable in the PW model but almost negligible in the OCA-EOM-CCSD results.

5. Conclusions

Theoretical ab initio modeling of Auger electron spectra remains a significant challenge, even for the smallest molecular systems. One of the main limitations in the field is the lack of generally available software tools that can be broadly applied and that are based on modern electronic structure methods, unlike what is available for modeling purely bound-domain electronic spectroscopies.

In the case of normal Auger decay, one of the most promising and practical approaches is based on the one-center approximation, which offers simplicity and low computational cost compared to methods that explicitly describe the continuum electron. Here, we have presented a methodology that combines EOM-CCSD calculations for the bound-domain electronic structure of the initial core-ionized and final dicationic states with an OCA treatment of the continuum part. We tested the method on small to medium-sized molecules and obtained very good agreement with available experimental data. Importantly, the computational cost of this approach is substantially lower than those involving explicit representation of the continuum orbital, even when the continuum orbital is approximated by the simplest plane-wave function. Our findings show that the OCA-EOM-CCSD methodology can serve as an almost black-box tool. A key advantage of combining the OCA with the EOM-CCSD framework is the simplicity of the computational setup: there is no need to define active spaces or rely on specific density functional choices, as otherwise required in multi-reference methods or some DFT-based implementations.

Undoubtedly, the accuracy of presented Auger spectra computed using the EOM-CCSD method remains lower when compared to that achieved for typical X-ray absorption spectra (XAS).^55,63 This is attributed to the complexity of the Auger process, which is inherently a two-electron phenomenon. It requires an accurate treatment of the electronic continuum and effects due to electron correlation between the bound and continuum states. As a result, ab initio modeling of Auger spectra presents significantly greater challenges than modeling XAS spectra. Nevertheless, our results demonstrate that normal Auger decay spectra of molecules containing carbon, nitrogen, and oxygen atoms can now be computed with satisfactory accuracy on a routine basis. The main limitation we identified lies in the high-binding-energy region of the spectrum, where Auger decay transitions involve electrons from inner-valence orbitals and the EOM-CCSD description of electron correlation may require further refinement. Such refinement can be approached from two directions. One direction is to improve the description of electronic correlation in the bound domain by including the effect of triple excitations in the coupled-cluster operator.⁸¹ The other directions is to enhance the description of the electronic continuum, particularly by incorporating effects due to interchannel coupling.^73,74,82 In principle, the interchannel coupling could be included even within the OCA framework, although this has not yet been attempted.

Important open questions remain whether the OCA approach can be systematically improved, as discussed above, and whether it can be successfully extended to processes beyond KLL decay, such as Coster–Kronig transitions and decays from the L-edge. Work in these directions is currently ongoing.

Author contributions

W. Skomorowski: conceptualization (equal); methodology (equal); data curation (lead); formal analysis (lead); software (lead); writing original draft (lead); writing – review & editing (equal). B. N. C. Tenorio: conceptualization (equal); methodology (equal); formal analysis (supporting); software (supporting); writing – review & editing (supporting). S. Coriani: conceptualization (lead); methodology (equal); software (supporting); formal analysis (supporting); writing original draft (supporting); writing – review & editing (equal).

Conflicts of interest

There are no conflicts to declare.

Data availability

The main data supporting this article has been included as part of the Supplementary Information, or is available from the authors upon request. The approach will be available in an upcoming release of Q-Chem.

Supplementary information: fc-CVS-EOM-IP-CCSD ionization potentials, detailed partial decay widths, comparison of selected minimal basis sets, molecular geometries is available. See DOI: https://doi.org/10.1039/d5cp02277k

Acknowledgements

B. N. C. T and S. C. acknowledge invaluable discussions with Prof. Piero Decleva, University of Trieste. W. S. acknowledges financial support from the Polish National Agency for Academic Exchange within Polish Returns Programme, and from the National Science Centre, Poland (SONATA 2020/39/D/ST4/03048). Part of the work was carried out with support from the European Union's Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie Individual Fellowship (B.N.C.T., Grant Agreement No. 101027796). We acknowledge the Polish high-performance computing infrastructure PLGrid (HPC Center: ACK Cyfronet AGH) for providing computer facilities and support within computational grant no. PLG/2023/016563.

Notes and references

1. E. H. S. Burhop, The Auger effect and other radiationless transitions, Cambridge University Press, 2014.
2. F. De Groot, Chem. Rev., 2001, 101, 1779–1808.
3. J. A. Van Bokhoven and C. Lamberti, X-Ray Absorption and X-Ray Emission Spectroscopy, 2 Volume Set: Theory and Applications, John Wiley & Sons, 2016, vol. 1.
4. T. Carlson, Photoelectron and Auger spectroscopy, Springer Science & Business Media, 2013.
5. M. F. Chung and L. H. Jenkins, Surf. Sci., 1970, 22, 479–485.
6. R. Weissmann and K. Müller, Surf. Sci. Rep., 1981, 1, 251–309.
7. T. W. Haas, J. T. Grant and G. J. Dooley III, J. Appl. Phys., 1972, 43, 1853–1860.
8. L. A. Harris, J. Appl. Phys., 1968, 39, 1419–1427.
9. N. H. Turner and J. A. Schreifels, Analyt. Chem., 2000, 72, 99–110.
10. R. R. Rye, T. E. Madey, J. E. Houston and P. H. Holloway, J. Chem. Phys., 1978, 69, 1504–1512.
11. J. Feldhaus, J. Arthur and J. B. Hastings, J. Phys. B., 2005, 38, S799.
12. J. Ullrich, A. Rudenko and R. Moshammer, Annu. Rev. Phys. Chem., 2012, 63, 635–660.
13. J. Duris, S. Li, T. Driver, E. G. Champenois, J. P. MacArthur, A. A. Lutman, Z. Zhang, P. Rosenberger, J. W. Aldrich and R. Coffee, et al. , Nat. Photonics, 2020, 14, 30–36.
14. O. Gessner and M. Gühr, Acc. Chem. Res., 2016, 49, 138–145.
15. T. Driver, M. Mountney, J. Wang, L. Ortmann, A. Al-Haddad, N. Berrah, C. Bostedt, E. G. Champenois, L. F. DiMauro, J. Duris, D. Garratt, J. M. Glownia, Z. Guo, D. Haxton, E. Isele, I. Ivanov, J. Ji, A. Kamalov, S. Li, M.-F. Lin, J. P. Marangos, R. Obaid, J. T. O’Neal, P. Rosenberger, N. H. Shivaram, A. L. Wang, P. Walter, T. J. A. Wolf, H. J. Wörner, P. H. B. Zhen Zhang, M. F. Kling, A. S. Landsman, R. R. Lucchese, A. Emmanouilidou, A. Marinelli and J. P. Cryan, Nature, 2024, 632, 762–767.
16. A. Kamalov, A. L. Wang, P. H. Bucksbaum, D. J. Haxton and J. P. Cryan, Phys. Rev. A, 2020, 102, 023118.
17. V. Gruson, L. Barreau, Á. Jiménez-Galan, F. Risoud, J. Caillat, A. Maquet, B. Carré, F. Lepetit, J.-F. Hergott and T. Ruchon, et al. , Science, 2016, 354, 734–738.
18. S. Li, T. Driver, P. Rosenberger, E. G. Champenois, J. Duris, A. Al-Haddad, V. Averbukh, J. C. T. Barnard, N. Berrah and C. Bostedt, et al. , Science, 2022, 375, 285–290.
19. L. Cattaneo, J. Vos, R. Y. Bello, A. Palacios, S. Heuser, L. Pedrelli, M. Lucchini, C. Cirelli, F. Martn and U. Keller, Nat. Phys., 2018, 14, 733–738.
20. D. C. Haynes, M. Wurzer, A. Schletter, A. Al-Haddad, C. Blaga, C. Bostedt, J. Bozek, H. Bromberger, M. Bucher and A. Camper, et al. , Nat. Phys., 2021, 17, 512–518.
21. A. Barty, J. Küpper and H. N. Chapman, Annu. Rev. Phys. Chem., 2013, 64, 415–435.
22. G. Wentzel, Z. Phys., 1927, 43, 524–530.
23. T. Åberg and G. Howat, in Corpuscles and Radiation in Matter I, ed W. Mehlhorn, Springer-Verlag, Berlin, 1982, vol. 31 of Encyclopedia of Physics, pp. 469–619.
24. H. Ågren and H. Siegbahn, Chem. Phys. Lett., 1980, 69, 424–429.
25. G. Howat, T. Aberg and O. Goscinski, J. Phys. B, 1978, 11, 1575.
26. H. Ågren, J. Chem. Phys., 1981, 75, 1267–1283.
27. E.-L. Ohrendorf, F. Tarantelli and L. S. Cederbaum, J. Chem. Phys., 1990, 92, 2984–2999.
28. B. Schimmelpfennig, B. M. Nestmann and S. D. Peyerimhoff, J. Electron Spectrosc. Relat. Phenom., 1995, 74, 173–186.
29. G. Grell, O. Kühn and S. I. Bokarev, Phys. Rev. A, 2019, 100, 042512.
30. G. Grell and S. I. Bokarev, J. Chem. Phys., 2020, 152, 074108.
31. R. Colle and S. Simonucci, Phys. Rev. A, 1989, 39, 6247.
32. V. Carravetta and H. Ågren, Phys. Rev. A, 1987, 35, 1022.
33. V. Averbukh and L. S. Cederbaum, J. Chem. Phys., 2005, 123, 204107.
34. P. Kolorenč and V. Averbukh, J. Chem. Phys., 2020, 152, 214107.
35. F. Matz and T.-C. Jagau, J. Chem. Phys., 2022, 156, 114117.
36. A. Ghosh, N. Vaval and S. Pal, Chem. Phys., 2017, 482, 160–164.
37. F. Matz and T.-C. Jagau, Mol. Phys., 2023, 121, e2105270.
38. H. Siegbahn, L. Asplund and P. Kelfve, Chem. Phys. Lett., 1975, 35, 330–335.
39. C. Liegener, J. Chem. Phys., 1983, 79, 2924–2927.
40. H. Ågren, A. Cesar and C.-M. Liegener, Adv. Quant. Chem., 1992, 23, 1–82.
41. T. Porwol, G. Dömötör, H.-J. Freund, R. Dudde, C. M. Lingener and W. von Niessen, Phys. Scr., 1992, 1992, 197.
42. R. Fink, J. Electron Spectrosc. Relat. Phenom., 1995, 76, 295–300.
43. C. Liegener, J. Chem. Phys., 1996, 104, 2940–2945.
44. F. Tarantelli, A. Sgamellotti and L. S. Cederbaum, J. Chem. Phys., 1991, 94, 523–532.
45. B. N. C. Tenorio, T. A. Voß, S. I. Bokarev, P. Decleva and S. Coriani, J. Chem. Theory Comput., 2022, 18, 4387–4407.
46. A. K. Schnack-Petersen, B. N. C. Tenorio, S. Coriani, P. Decleva, J. Troß, K. Ramasesha, M. Coreno, R. Totani and A. Röder, J. Chem. Phys., 2022, 157, 214305.
47. B. N. C. Tenorio, K. B. Møller, P. Decleva and S. Coriani, Phys. Chem. Chem. Phys., 2022, 24, 28150–28163.
48. B. N. C. Tenorio, J. Pedersen, M. Barbatti, P. Decleva and S. Coriani, J. Phys. Chem. A, 2024, 128, 107–117.
49. A. E. A. Fouda, B. Jangid, E. Pelimanni, S. H. Southworth, P. J. Ho, L. Gagliardi and L. Young, J. Phys. Chem. A, 2025 DOI:10.1021/acs.jpca.5c01789.
50. A. I. Krylov, Annu. Rev. Phys. Chem., 2008, 59, 433–462.
51. J. F. Stanton and R. J. Bartlett, J. Chem. Phys., 1993, 98, 7029–7039.
52. R. J. Bartlett, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 126–138.
53. W. Skomorowski and A. I. Krylov, J. Chem. Phys., 2021, 154, 084124.
54. W. Skomorowski and A. I. Krylov, J. Chem. Phys., 2021, 154, 084125.
55. M. L. Vidal, X. Feng, E. Epifanovsky, A. I. Krylov and S. Coriani, J. Chem. Theory Comput., 2019, 15, 3117–3133.
56. L. S. Cederbaum, W. Domcke and J. Schirmer, Phys. Rev. A: At., Mol., Opt. Phys., 1980, 22, 206–222.
57. S. Coriani and H. Koch, J. Chem. Phys., 2015, 143, 181103.
58. R. Manne and H. Ågren, Chem. Phys., 1985, 93, 201–208.
59. M. Nooijen and R. J. Bartlett, J. Chem. Phys., 1997, 107, 6812–6830.
60. M. Wladyslawski and M. Nooijen, ACS Symp. Ser., 2002, 65–92.
61. C.-M. Liegener, Chem. Phys. Lett., 1982, 90, 188–192.
62. J. P. Carbone, L. Cheng, R. H. Myhre, D. Matthews, H. Koch and S. Coriani, Adv. Quantum Chem., 2019, 79, 241–261.
63. T. Fransson, I. E. Brumboiu, M. L. Vidal, P. Norman, S. Coriani and A. Dreuw, J. Chem. Theory Comp., 2021, 17, 1618–1637.
64. M. L. Vidal, P. Pokhilko, A. I. Krylov and S. Coriani, J. Phys. Chem. Lett., 2020, 11, 8314–8321.
65. N. K. Jayadev, A. Ferino-Pérez, F. Matz, A. I. Krylov and T.-C. Jagau, J. Chem. Phys., 2023, 158, 064109.
66. E. J. McGuire, Phys. Rev., 1969, 185, 1–6.
67. A. Rohatgi, WebPlotDigitizer, https://automeris.io.
68. E. Epifanovsky, A. T. B. Gilbert, X. Feng, J. Lee, Y. Mao, N. Mardirossian, P. Pokhilko, A. F. White, M. P. Coons and A. L. Dempwolff, et al. , J. Chem. Phys., 2021, 155, 084801.
69. W. E. Moddeman, T. A. Carlson, M. O. Krause, B. P. Pullen, W. E. Bull and G. K. Schweitzer, J. Chem. Phys., 1971, 55, 2317–2336.
70. L. Inhester, C. F. Burmeister, G. Groenhof and H. Grubmüller, J. Chem. Phys., 2012, 136, 144304.
71. A. Cesar, H. Ågren and V. Carravetta, Phys. Rev. A, 1989, 40, 187–206.
72. H. Ågren and V. Carravetta, J. Chem. Phys., 1987, 87, 370–380.
73. G. Howat, T. Aberg and O. Goscinski, J. Phys. B, 1978, 11, 1575.
74. R. Colle and S. Simonucci, Phys. Rev. A: At., Mol., Opt. Phys., 1990, 42, 3913–3925.
75. R. Sankari, M. Ehara, H. Nakatsuji, Y. Senba, K. Hosokawa, H. Yoshida, A. De Fanis, Y. Tamenori, S. Aksela and K. Ueda, Chem. Phys. Lett., 2003, 380, 647–653.
76. F. Holzmeier, T. J. A. Wolf, C. Gienger, I. Wagner, J. Bozek, S. Nandi, C. Nicolas, I. Fischer, M. Gühr and R. F. Fink, J. Chem. Phys., 2018, 149, 034308.
77. V. Feyer, P. Bolognesi, M. Coreno, K. C. Prince, L. Avaldi, L. Storchi and F. Tarantelli, J. Chem. Phys., 2005, 123, 224306.
78. J. Nordgren, L. Selander, L. Pettersson, C. Nordling, K. Siegbahn and H. Ågren, J. Chem. Phys., 1982, 76, 3928–3932.
79. F. Tarantelli, A. Sgamellotti, L. S. Cederbaum and J. Schirmer, J. Chem. Phys., 1987, 86, 2201–2206.
80. S. Carniato, P. Selles, A. Ferté, N. Berrah, A. H. Wuosmaa, M. Nakano, Y. Hikosaka, K. Ito, M. Žitnik, K. Bučar, K. Soejima, K. Jänkälä, D. Cubaynes, J.-M. Bizau, L. Andric, M. A. Khalal, J. Palaudoux, P. Lablanquie and F. Penent, J. Phys. B: At., Mol. Opt. Phys., 2020, 53, 244010.
81. D. A. Matthews, Mol. Phys., 2020, 118, e1771448.
82. T. Moitra, S. Coriani and P. Decleva, J. Chem. Theory Comput., 2021, 17, 5064–5079.

---