Multi-reference quantum chemistry protocol for simulating autoionization spectra: Test of ionization continuum models for the neon atom

TL;DR

This paper introduces a quantum chemistry protocol combining electronic structure and continuum models to accurately simulate Auger decay spectra, validated on neon atom data.

Contribution

The work presents a novel combined approach for calculating Auger decay rates using restricted active space methods and continuum orbitals, with systematic testing against reference data.

Findings

Continuum orbital accuracy critically affects decay rate predictions.

Auger energies are well reproduced with small active spaces.

Direct Coulomb potential inclusion improves continuum wave functions.

Abstract

In this contribution we present a protocol to evaluate partial and total Auger decay rates combining the restricted active space self-consistent field electronic structure method for the bound part of the spectrum and numerically obtained continuum orbitals in the single-channel scattering theory framework. On top of that, the two-step picture is employed to evaluate the partial rates. The performance of the method is exemplified for the prototypical Auger decay of the neon $1s^{-1}3p$ resonance. Different approximations to obtain the continuum orbitals, the partial rate matrix elements, and the electronic structure of the bound part are tested against theoretical and experimental reference data. It is demonstrated that the partial and total rates are most sensitive to the accuracy of the continuum orbitals. For instance, it is necessary to account for the direct Coulomb potential of the ion for the determination of the continuum wave functions. The Auger energies can be reproduced quite well already with a rather small active space. Finally, perspectives of the application of the proposed protocol to molecular systems are discussed.

Tables

Table 1. Table 1: QC setups and number of states used in the state-averaged RASSCF calculations

			States
QC	Basis Set	Active Space (AS)	Ne/Ne+	X2C
I	[7s6p3d2f]	RAS$(8;1,1)$	41/132	no
II	[22s6p3d2f]-rcc	RAS$(8;1,1)$	41/132	yes
III	[9s8p5d4f]	RAS$(26;1,1)$	131/420	no

Table 2. Table 2: Maximum number of configurations for each AS

RAS	RAS$(8;1,1)$	RAS$(8;26,1)$	RAS$(8;33,1)$
Ne	41	131	166
Ne+	197	629	797

Equations

$${\Gamma }_{i\alpha }=2\pi \text{absolutevalue}{\text{matrixelement}{\Psi }_{\alpha }\mathcal{H}-{\mathcal{E}}_i{\Psi }_i}^2.$$

$${\color[rgb]{0,0,0}\ket{\Psi}=\sum_{j}C_{j}\cdot a^{\dagger}_{j_{1},\sigma_{1}}\cdots a^{\dagger}_{j_{N},\sigma_{N}}\ket{0}\,.}$$

$${\psi }_{\alpha ,\sigma }(r,\vartheta ,\varphi )=\frac{1}{r}{w}_l^{fk}(r)Y_l^m(\vartheta ,\varphi )\zeta (\sigma ),$$

$$\left(\frac{d^2}{d{r}^2}+2\left(\frac{k^2}{2}-V_f(r)\right)-\frac{l(l+1)}{r^2}\right)w_l^{fk}(r)=0.$$

$$|{\Psi }_{\alpha }⟩=\sum _{M^{+}=-S^{+}}^{S^{+}}\sum _{\sigma =-\frac{1}{2},\frac{1}{2}}{C}_{S^{+},M^{+};\sigma }^{S,M}|{{\rm Y}}_{\alpha }^{M^{+},\sigma }⟩,$$

$$|{{\rm Y}}_{\alpha }^{\sigma M^{+}}⟩=a_{\alpha ,\sigma }^{†}|{\Psi }_{f,M^{+}}^{+}⟩.$$

$$V^{\text{free}}$$

$$V^{\text{eff}}(r)$$

$$V_f^{\text{scr}}(r)$$

$$V_f^{\text{J}}(r)$$

$$V_f^{\text{JX}}(r)$$

$$\psi^{\text{free}}_{\alpha,\sigma}(r,\vartheta,\phi)={\color[rgb]{0,0,0}\sqrt{\frac{2}{\pi}k}}\cdot j_{l}(kr)Y_{l}^{m}(\vartheta,\phi)\zeta(\sigma)$$

$$\innerproduct{\psi_{\alpha,\sigma}}{\psi_{\alpha^{\prime},\sigma^{\prime}}}=\delta_{ll^{\prime}}\delta_{m_{l}m_{l}^{\prime}}\delta_{\sigma\sigma^{\prime}}{\color[rgb]{0,0,0}\delta(\varepsilon_{\alpha}-\varepsilon_{\alpha^{\prime}})}\,,$$

$$\psi^{\text{eff}}_{\alpha,\sigma}(r,\vartheta,\phi)={\color[rgb]{0,0,0}\sqrt{\frac{2}{\pi k}}}\cdot\frac{F_{l}(\eta,kr)}{r}Y_{l}^{m}(\vartheta,\phi)\zeta(\sigma)\,,$$

$$w_l^{fk}(r\to 0)$$

$$w_l^{fk}(r\to \infty )$$

$$+\sin {\delta }_l^f(k)G_l(\eta ,kr)).$$

$$\displaystyle\Gamma_{i\alpha}=2\pi\Biggl{|}$$

$$\displaystyle\left.+\sum_{u<v}\matrixelement{\Upsilon_{\alpha}^{\sigma M^{+}}}{\frac{1}{r_{uv}}}{\Psi_{i}}-\mathcal{E}_{i}\innerproduct{\Upsilon_{\alpha}^{\sigma M^{+}}}{\Psi_{i}}\right]\Biggr{|}^{2}\,.$$

$$\text{innerproduct}{{\rm Y}}_{\alpha }^{\sigma M^{+}}{\Psi }_i=\text{innerproduct}{\psi }_{\alpha ,\sigma }{\Phi }_{i\alpha }^{M^{+}}.$$

$$|{\Phi }_{i\alpha }^{M^{+}}⟩=\sum _{s=1}^{N_{\text{orb}}}\underset{{\varphi }_{\alpha ,s}^{M^{+}}}{\underbrace{\text{matrixelement}{\Psi }_{f,M^{+}}^{+}{a}_s{\Psi }_i}}|{\phi }_s⟩.$$

$$\sum _u\text{matrixelement}{{\rm Y}}_{\alpha }^{M^{+},\sigma }{h}_u{\Psi }_i$$

$$\sum _{u<v}\text{matrixelement}{{\rm Y}}_{\alpha }^{\sigma M^{+}}\frac{1}{r_{uv}}{\Psi }_i$$

$$+\text{innerproduct}{\psi }_{\alpha ,\sigma }{\tilde{\Phi }}_{i\alpha }^{2,M^{+}}.$$

$$|{\Xi }_{i\alpha }^{M^{+},q}⟩=\sum _{s_1\ne s_2}^{N_{\mathrm{orb}}}\text{matrixelement}{\Psi }_{f,M^{+}}^{+}({\hat{a}}_q^{+}{)}^{†}{\hat{a}}_{s_1}{\hat{a}}_{s_2}{\Psi }_i|{\phi }_{s_1}{\phi }_{s_2}⟩,$$

$$\text{innerproduct}\alpha abc=\int {\psi }_{\alpha }^{\ast }({\mathbf{r}}_1){\chi }_b({\mathbf{r}}_1)\underset{f_{ac}({\mathbf{r}}_1)}{\underbrace{\left[\int \frac{{\chi }_a^{\ast }({\mathbf{r}}_2){\chi }_c({\mathbf{r}}_2)}{r_{12}})d{\mathbf{r}}_2^3\right]}}d{\mathbf{r}}_1^3.$$

$$\text{innerproduct}\alpha abc$$

$$\times \underset{F_{abc}(r)}{\underbrace{\left[{\int }_{\Omega }{Y}_l^m(\Omega ){\chi }_b(\mathbf{r}(r,\Omega ))f_{ac}(\mathbf{r}(r,\Omega ))d\Omega \right]}}dr.$$

$${\mathrm{AES}}_{1s^{-1}3p}(\epsilon )=\sum _{\alpha }{\Gamma }_{1s^{-1}3p\alpha }G(\epsilon -{\epsilon }_{\alpha },\gamma )$$

$${\rho }_f(r)=\frac{1}{4\pi }{\int }_0^{4\pi }{\rho }_f(\mathbf{r})d\Omega .$$

$$Z_f(r)=Z-{\int }_0^r{\rho }_f(r^{\prime })r^{\prime 2}d{r}^{\prime }.$$

$$J_f(r)=J_f(0)-\frac{4\pi }{r}{\int }_0^{r}d{r}^{\prime }{\int }_0^{r^{\prime }}{r}^{\prime \prime }{\rho }_f(r^{\prime \prime })d{r}^{\prime \prime }.$$

$$\frac{J_f(0)}{4\pi }=\underset{r^{\prime }\to \infty }{\lim }{\int }_0^{r^{\prime }}{r}^{\prime \prime }{\rho }_f(r^{\prime \prime })d{r}^{\prime \prime }.$$

$$X_f^{\text{S}}(r)=-3{\left(\frac{3}{8\pi }{\rho }_f(r)\right)}^{\frac{1}{3}},$$

$$A_{i\alpha }^{\sigma M^{+}}=\underset{\text{NO}}{\underbrace{\underset{\text{SO}}{\underbrace{\text{matrixelement}{\psi }_{\alpha ,\sigma }h{\Phi }_{i\alpha }^{M^{+}}+\stackrel{r^{-1}\text{ coupling}}{\overbrace{\sum _{q=1}^{N_{\mathrm{orb}}}\text{matrixelement}{\psi }_{\alpha ,\sigma }{\phi }_q^{+}\frac{1}{r_{12}}{\Xi }_{i\alpha }^{M^{+},q}}}}}+\text{innerproduct}{\psi }_{\alpha ,\sigma }{\tilde{\Phi }}_{i\alpha }^{1,M^{+}}+\stackrel{r^{-1}\text{ coupling}}{\overbrace{\text{innerproduct}{\psi }_{\alpha ,\sigma }{\tilde{\Phi }}_{i\alpha }^{2,M^{+}}}}-{\mathcal{E}}_i\text{innerproduct}{\psi }_{\alpha ,\sigma }{\Phi }_{i\alpha }^{M^{+}}}}.$$

Full text

SCF Self-Consistent Field PT2 Second Order Perturbation Theory RAS Restricted Active Space CAS Complete Active Space AS Active Space MCSCF Multi-configurational Self-Consistent Field (SCF) RASSCF Restricted Active Space (RAS) SCF

CASSCF Complete Active Space (CAS) SCF

RASPT2 Restricted Active Space Second-order Perturbation Theory RASSI

RAS State Interaction HF Hartree-Fock CI Configuration Interaction MCDF Multi-configurational Dirac-Fock ADC Algebraic Diagrammatic Construction Fano-ADC Fano-Stieltjes Algebraic Diagrammatic Construction DKH Douglas-Kroll-Hess X2C Exact Two Component Decoupling ANO Atomic Natural Orbital QC Quantum Chemistry GTO Gaussian Type Orbital ICD Interatomic Coulombic Decay ETMD Electron Transfer Mediated Decay PES Photoelectron Spectroscopy AES Auger Electron Spectrum AIS Autoionization Spectroscopy RPES Resonant Photoelectron Spectroscopy (PES) SE Schrödinger Equation SD Slater Determinant SO Strong Orthogonality GTO Gaussian Type Orbital DO Dyson Orbital MO Molecular orbital GS Gram-Schmidt

I Introduction

Ionization triggered by photon absorption occurs along two pathways. In direct photoionization, the energy is transferred to an ejected electron. Alternatively, the system can be first put into a metastable state by a resonant excitation and afterwards decay via an autoionization mechanism. Autoionization can be approximately understood as a two-step process Wentzel (1927), in which the decay can be considered independently from the excitation process and interferences between direct and autoionization are neglected. For example, let us consider an atomic species, such as a neon atom that is prepared in a highly excited state $\ket{\Psi_{i}}$ above the continuum threshold at $E=0$ eV, Fig. 1. This state spontaneously decays into the continuum state $\ket{\Psi_{\alpha}}$ comprising the discrete state $\ket{\Psi_{f}^{+}}$ of the ion and the emitted electron, $\ket{\psi_{\alpha}}$, carrying the excess energy $\varepsilon_{\alpha}=\mathcal{E}_{i}-\mathcal{E}_{f}$.

The system’s electronic structure is thus encoded into the kinetic energy spectrum of the ionized electrons. PES and Autoionization Spectroscopy (AIS) map bound states to the continuum which makes them less sensitive to selection rule suppression and more informative than spectroscopies involving optical transitions between bound states Van der Heide (2012); Hofmann (2013); Hüfner (2003).

Autoionization processes, predominantly Auger decay Meitner (1922), but also Interatomic Coulombic Decay (ICD) Cederbaum et al. (1997) and Electron Transfer Mediated Decay (ETMD) Zobeley et al. (2001) are particularly interesting on their own. Due to their correlated nature, they not only probe but also initiate or compete with intricate ultrafast electronic and nuclear dynamics e.g. Slavíček et al. (2016, 2014); Brown et al. (2009); Stumpf et al. (2016); Wang et al. (2018). Additionally, they provide the main channel for the decay of core vacancies de Groot and Kotani (2008) and play a key role in biological radiation damage, creating highly charged cations while a cascade of highly reactive low energy electrons is emitted Howell (2008); Alizadeh et al. (2015); Stumpf et al. (2016); Martin and Feinendegen (2016); Yokoya and Ito (2017). Further, in free electron laser experiments operating with ultrashort intense X-ray pulses, autoionization after multiple photoionization induces the Coulomb explosion of the target, which limits the achievable spectroscopic and temporal resolution Neutze et al. (2000). Due to this wealth of applications, autoionization and especially the local Auger effect have been studied extensively both theoretically and experimentally since its discovery by Meitner Meitner (1922) and description by Wentzel Wentzel (1927).

Remarkably, AIS simulations of molecular systems remain challenging until today, although the fundamental theory is known for decades Åberg and Howat (1982); Åberg (1980); Fano (1961). For atoms, methods combining highly accurate four-component Multi-configurational Dirac-Fock (MCDF) calculations with multichannel scattering theory are publicly available Fritzsche (2012), whereas no such general purpose code exists for molecules. The main complication of the molecular case lies in the construction of molecular continuum states $\ket{\Psi_{\alpha}}$. The approaches to the simulation of AIS published during the last decades can be classified into two families – those that circumvent the continuum orbital problem and those that treat the continuum orbital explicitly.

The first family comprises the following flavors: The simplest method that allows to assign experimental AIS is to evaluate the energetic peak positions Ågren et al. (1975); Ågren (1981); Hillier and Kendrick (1976). On top of that simple estimates for the partial decay rates can be obtained based on an electron population analysis Mitani et al. (2003). More advanced approaches rely on an implicit continuum representation with Stieltjes imaging Carravetta et al. (2000), a Green’s operator Schimmelpfennig et al. (1992, 1995), or a propagator Oddershede (1987); Liegener (1982) formalism. From this group, the Fano-Stieltjes Algebraic Diagrammatic Construction (Fano-ADC) method Averbukh and Cederbaum (2005) has been used to evaluate Auger, ICD and ETMD decay rates of van der Waals clusters Kolorenč et al. (2008), first row hydrides Kolorenč and Averbukh (2011) and the $[\text{Mg}(\text{H}_{2}\text{O})_{6}]^{2+}$ cluster Stumpf et al. (2016).

Therein, the continuum is approximately represented with spatially confined basis functions. However, the description of Auger electrons with kinetic energies of several hundreds of eV requires large basis sets leading to computationally demanding simulations.

The second family, relying on an explicit representation of the continuum wave function, consists of the following approaches: The one-center approximation which uses atomic continuum functions centered at the vacancy-bearing atom to describe the outgoing electron in the evaluation of partial decay rates Siegbahn et al. (1975); Jennison (1980); Larkins et al. (1990); Fink (1995). This approximation can be applied on top of high-level electronic structure methods Travnikova et al. (2009). It is also applied in the XMOLECULE package Inhester et al. (2016); Hao et al. (2015) which is based on very cost efficient electronic structure calculations and atomic continuum functions. This allows for the evaluation of ionization cascades but may limit the applicability for strongly correlated systems, e.g., possessing a multi-configurational character.

Further, the influence of the molecular field may be taken into account perturbatively Faegri and Kelly (1979); Higashi et al. (1982) or in a complete manner with, e.g., the single-center approach, where the whole molecular problem is projected onto a single-centered basis Zähringer et al. (1992); Demekhin et al. (2007, 2009, 2011); Inhester et al. (2012); Banks et al. (2017).

Finally, multi-channel scattering theory methods that combine finite multi-centered basis sets with the appropriate boundary conditions to represent the molecular continuum have been developed Colle and Simonucci (1989) and applied to the AIS of a variety of small systems Colle and Simonucci (1989, 1993); Colle et al. (2004). For instance, the recently developed XCHEM approach Marante et al. (2017a) has been applied to simulate PES and AIS of atoms Marante et al. (2017b) and small molecular systems Klinker et al. (2018). These techniques represent the most general and accurate quantum-mechanical treatment of the problem, thus potentially serving as a high-level reference, although connected to substantial computational effort.

Summarizing, most of the mentioned methods have been applied only to simple diatomics, first row hydrides, halogen hydrides, and small molecules consisting of not more than two heavy atoms. Studies of larger molecular systems, such as tetrahedral molecules, small aldehydes, and amides Ortiz (1984); Correia et al. (1991), solvated metal ions Stumpf et al. (2016) and polymers Endo et al. (2016) are very scarce. In fact, the Fano-ADC Averbukh and Cederbaum (2005); Kolorenč et al. (2008) and the XMOLECULE Hao et al. (2015); Inhester et al. (2016) approaches are the only publicly available tools that allow to simulate AIS for a variety of systems without restricting the molecular geometry. Further, both methods are not suited to treat systems possessing multi-configurational wave functions. This puts studies of some chemically interesting systems having near-degeneracies, for example, transition metal compounds, or of photodynamics in the excited electronic states, e.g., near conical intersections, out of reach. To keep up with the experimental advancements, the development of a general purpose framework to evaluate autoionization decay rates (Auger, ICD, and ETMD) for molecular systems is warrant. Such a framework should be kept accessible, transferable and easy to use, i.e. it should be based on widespread robust and versatile Quantum Chemistry (QC) methods.

Here, we present a protocol that combines multi-configurational (RASSCF) bound state wave functions with single-centered numerical continuum orbitals in the single-channel scattering theory framework Burke (2011). We have chosen the RASSCF approach, since it is known to yield reliable results for core-excited states Bokarev and Kühn (2019), needed in the simulation of X-ray absorption Bokarev et al. (2013); Pinjari et al. (2016), resonant inelastic scattering Josefsson et al. (2012); Bokarev et al. (2015),and photoemission spectra Grell et al. (2015); Golnak et al. (2016); Norell et al. (2018) suggesting its application to AIS.

Although the ultimate goal is to investigate molecules, this proof-of-concept contribution focuses on the simulation of the prototypical neon $1s^{-1}3p$ Auger Electron Spectrum (AES) to calibrate the approach, since highly accurate reference data are available from both theory Stock et al. (2017) and experiment  Kivimäki et al. (2001); Yoshida et al. (2000); Ueda et al. (2003); Tamenori and Suzuki (2014). Special attention is paid to the representation of the radial continuum waves, which is investigated herein by a thorough test of different approximations. Note that our implementation allows to calculate molecular AIS as well as PES which will be presented elsewhere.

We commence this article with an introduction to the underlying theory and further give important details of our implementation. It continues with the computational details and the benchmark of our results against theoretical and experimental references. Finally, we conclude the discussion and present perspectives for the molecular application.

II Theory

The approach for the calculation of partial autoionization rates comprises the following approximations:

(i) The two-step model Wentzel (1927) is employed, i.e., excitation and decay processes are assumed to be decoupled and interference effects between photoionization and autoionization are neglected, Fig. 1. Within this approximation, the partial autoionization rate for the decay $i\rightarrow\alpha$ reads Åberg and Howat (1982)

[TABLE]

Atomic units are used everywhere, unless explicitly stated otherwise.

(ii) We require that all bound state wave functions have the form of Configuration Interaction (CI) expansions in terms of $N$-electron Slater determinants:

[TABLE]

The $a^{\dagger}_{i,\sigma_{i}}$ are the usual fermionic creation operators in the spin-orbital basis $\{\varphi_{i,\sigma_{i}}\}$. In this work, multi-configurational bound state wave functions for the unionized and ionized states $\ket{\Psi_{i}}$ and $\ket{\Psi^{+}_{f}}$ are obtained with the RASSCF or Restricted Active Space Second-order Perturbation Theory (RASPT2) method. However, the presented protocol can employ any CI-like QC method.

(iii) The limit of weak relativistic effects is assumed, thus the total spins $S,S^{+}$ and their projections onto the quantization axis $M,M^{+}$ of the bound unionized and ionized system are good quantum numbers. Further, $S$ and $M$ of the unionized species are conserved during the process.

(iv) The single-channel scattering theory framework is employed, disregarding interchannel coupling, as well as correlation effects between the bound and outgoing electrons.

(v) The continuum orbitals are treated as spherical waves, subject to the spherically averaged potential $V_{f}(r)$ of the ionic state $\ket{\Psi_{f}^{+}}$. Hence, the continuum orbitals have the form

[TABLE]

with spherical harmonics $Y_{l}^{m}(\vartheta,\phi)$ being the angular part and $\zeta(\sigma)$ the spin function. For brevity, $\zeta(\sigma)$ is generally omitted and present only when needed. The composite channel index $\alpha=(f,l,m,k)$ contains the index of the ionized state $f$, the orbital and magnetic quantum numbers $l$, $m$ and the wave number $k=\sqrt{2\varepsilon_{\alpha}}$ of the continuum orbital. This notation uniquely identifies the total energy $\mathcal{E}_{\alpha}=\mathcal{E}_{f}+\varepsilon_{\alpha}$, the continuum orbital and bound state for each channel $\ket{\Psi_{\alpha}}$. Generally, indices $i$ and $f$ always refer to bound states of the unionized and ionized species and $\alpha$ denotes decay channels.

The radial part $w_{l}^{fk}(r)/r$ is determined by solving the radial Schrödinger equation

[TABLE]

With the assumptions (i)-(v), the $N$-electron ionized continuum states with conserved total spin and projection of the unionized states, $S$ and $M$, can be written as:

[TABLE]

where $\sigma$ is the spin projection of the outgoing electron. The Clebsch-Gordan coefficients $C_{S^{+},M^{+};\sigma}^{S,M}=\innerproduct{S,M}{S^{+},M^{+};\frac{1}{2},\sigma}$ couple the channel functions $\ket{\Upsilon_{\alpha}^{\sigma M^{+}}}$. These are constructed by inserting an additional electron with the continuum orbital $\ket{\psi_{\alpha,\sigma}}$ into the bound ionic state with spin projection $M^{+}$, retaining the anti-symmetry:

[TABLE]

Note that in contrast to the continuum state $\ket{\Psi_{\alpha}}$ which is an eigenstate of $\mathbf{S}^{2}$ and $S_{z}$, the $\ket{\Upsilon_{\alpha}^{\sigma M^{+}}}$ are not eigenfunctions of $\mathbf{S}^{2}$, but only of $S_{z}$.

III Details of the implementation

III.1 Model potentials

We commence with the discussion of different approximations to the potential $V_{f}(r)$ in Eq. (4) which are central for the quality of the free electron function. In this work, we have used the following models:

[TABLE]

In words, we use (a) no potential; (b) an effective Coulomb potential with a fixed charge $Z_{\text{eff}}$; (c) the screened Coulomb potential of an ionic state $f$ with the charge $Z_{f}(r)$ varying with distance; (d) the nuclear, $-Z/r$, and spherically averaged direct potential of the ionic state, $J_{f}(r)$; (e) the potential of (d) augmented with Slater’s exchange term $X^{\text{S}}_{f}(r)$ Slater (1951). Details on the definition of $Z_{f}(r)$, $J_{f}(r)$ and $X^{\text{S}}_{f}(r)$ are given in Appendix A.

III.2 Continuum orbitals

The solutions to the radial Schrödinger equation, Eq. (4), using the potentials (7a)–(7e) can be understood as follows: In (a), assuming a free particle, we completely neglect any influence of the ion onto the outgoing electron. Here, the radial part of $\psi_{\alpha,\sigma}(r,\vartheta,\varphi)$ corresponds to spherical Bessel functions $j_{l}(kr)$ Olver et al. (2010):

[TABLE]

The prefactor $\sqrt{\frac{2}{\pi}k}$ ensures the correct normalization,

[TABLE]

which is exact only for the free particle approach, and approximate for all other potentials, due to the long range Coulomb distortion.

An effective Coulomb form (b) is the simplest approach to approximately account for the ionic potential. It has the advantage that the solutions to Eq. (4) are still analytically available in the form of regular Coulomb functions $F_{l}$ Olver et al. (2010):

[TABLE]

with $\eta=-Z_{\text{eff}}/2k$.

The approaches (c-e) employ numerically obtained potentials according to Eqs. (7c)–(7e) and thus require a numerical solution of the radial Schrödinger equation (4).

It is done using Numerov’s method and the numerically obtained radial waves are scaled such as to satisfy the asymptotic boundary conditions:

[TABLE]

The scaling factors and phase shifts $\delta^{f}_{l}(k)$ are obtained by matching the numerical solutions and their first derivatives with a linear combination of the regular and irregular Coulomb functions $F_{l}$ and $G_{l}$, see Ref. 74. This matching is carried out in the asymptotic region where the potential is well approximated by the Coulomb potential of the ions net charge $V_{f}(r)\approx-Z_{\text{net}}/r$.

Since the screened Coulomb potential model, Eq. (7c) is an ad-hoc assumption based on the simple idea that the nuclear charge is screened by the integrated electron density, the applicability of this model needs to be tested. In turn, (d) and (e), employing the spherically averaged direct Coulomb and exchange terms according to Eqs. (7d) and (7e), in a sense correspond to the Hartree and Slater-$X_{\alpha}$ levels of accuracy, respectively.

The continuum orbitals obtained by any of these methods are not orthogonal to the bound orbitals, which is in contrast to the behavior that an exact continuum orbital would possess.

III.3 Continuum matrix elements

Using the decomposition of the continuum states in terms of channel functions in Eq. (5) and separating the Hamiltonian into its one and two electron parts $\mathcal{H}=\sum_{u}h_{u}+\sum_{u<v}1/r_{uv}$, the partial decay rate becomes:

[TABLE]

Here, $u$ and $v$ are electron indices, $h_{u}$ contains the electronic kinetic energy and electron-nuclei attraction terms and $1/r_{uv}$ is the electron repulsion. The expressions for the overlap, one- and two-electron matrix elements in Eq. (12) are obtained by using Löwdins Slater determinant calculus Löwdin (1955). Here one has to take into account that the unionized and ionized bound states are obtained in separate SCF calculations. Consequently, they have different sets of $N_{\rm orb}$ spin-orbitals $\{\varphi_{i}\}$ and $\{\varphi^{+}_{i}\}$ that are not mutually orthogonal. The spin coordinates are implicitly assumed to be assigned as introduced in Eq. (2). Then, the respective creation and annilation operators are $a_{i}^{\dagger},\,a_{i}$ and $({a^{+}_{i}})^{\dagger},\,a^{+}_{i}$.

With this, the overlap integral in Eq. (12) can be rearranged into the overlap of the continuum orbital and the Dyson Orbital (DO) $\ket{\Phi_{i\alpha}^{M^{+}}}$,

[TABLE]

The DO is generally defined as the $N-1$ particle integral over the transition density of the unionized and ionized states $\ket{\Psi_{i}}$ and $\ket{\Psi_{f,M^{+}}^{+}}$ that are associated with the channel $\alpha$ and can be expressed in second quantization form as a linear combination of the spin-orbitals $\{\varphi_{s}\}$ of the unionized species, with the coefficients $\phi^{M^{+}}_{\alpha,s}$

[TABLE]

Staying on this route, the one- and two-electron transition matrix elements in Eq. (12) can be expressed as:

[TABLE]

and

[TABLE]

Here, $\ket{\Xi_{i\alpha}^{M^{+},q}}$ is the two-electron reduced transition density

[TABLE]

and $\ket{\tilde{\Phi}_{i\alpha}^{n,M^{+}}}$ are the one- and two-electron conjugated Dyson orbitals for $n=1$ and 2, respectively, see Eqs. (28) and (29). A simpler formulation of the NO terms proposed in Ref. 76 is not used herein, because in practice it is not always strictly equivalent to our approach, check Supplement: Section II for details.

The Strong Orthogonality (SO) approximation, i.e., the assumption that the overlaps of the continuum and all unionized orbitals are zero, $\innerproduct{\psi_{\alpha,\sigma}}{\varphi_{i,\sigma_{i}}}=0$, implies that the overlap integrals between the continuum and the ordinary and conjugated Dyson orbitals in Eqs. (13), (15), and (16) vanish. Since the evaluation of the conjugate DOs can be quite involved, SO approximation substantially simplifies the computation.

A similar effect is achieved by using the Gram-Schmidt (GS) orthogonalization to project the orbitals of the ionized system out of the continuum functions and enforce $\innerproduct{\psi^{\text{GS}}_{\alpha,\sigma}}{\varphi^{+}_{i,\sigma_{i}}}=0$. This approach has been tried before, e.g., in Ref. 44 and will be tested herein as well.

In the following, we will use the labels SO, NO and GS to indicate that the overlap terms in Eqs. (13), (15), and (16) have been neglected (SO), fully included (NO), and that GS-orthogonalized continuum functions have been used to evaluate the partial rates including the NO terms as well. Further, we also compare the results of the “full” Hamiltonian coupling, against the popular choice to account only for the two-electron terms in Eq. (12), denoting them as $\mathcal{H}$ and $r^{-1}$ coupling, respectively. The combination of matrix elements corresponding to each approach is detailed in Appendix B.

To sum up, the present article reports on the influence of different combinations of the introduced approximations to the potentials, transition matrix elements, and continuum orbitals on the partial Auger decay rates of the exemplary neon $1s^{-1}3p$ resonance. As a shorthand notation for the combination of the different approximations we will use the $coupling$ $\cdot$ $potential$ $\cdot$ $nonorthogonality$ notation, were applicable. For instance, $\mathcal{H}\cdot V_{f}^{\text{JX}}(r)\cdot\text{NO}$ denotes partial rates obtained using the $\mathcal{H}$ coupling, with radial waves corresponding to the $V_{f}^{\text{JX}}(r)$ potential including the overlap terms.

IV Computational Details

IV.1 Quantum chemistry for bound states

The wave functions of the bound neutral and ionic states have been obtained in separate state-averaged RASSCF calculations with a locally modified version of MOLCAS 8.0 Aquilante et al. (2016). To prevent mixing of different angular momentum basis functions into one orbital, the ”atom” keyword has been employed. The QC schemes used to evaluate the bound states $\ket{\Psi_{i}}$ and $\ket{\Psi_{f}^{+}}$ are presented in Table 1.

The RAS formalism is a flexible means to select electronic configurations. Therein, the AS is subdivided into three subspaces RAS1, RAS2 and RAS3. The RAS notation that is used throughout the paper is to be understood as follows: In all spaces, the $1s$ orbital forms the RAS1 subspace and the $2s$ and $2p$ orbitals build up the RAS2 one. The occupation of the $2s$ and $2p$ orbitals in RAS2 is not restricted, while only $h$ electrons may be removed from RAS1. Finally the RAS3 subspace contains $v$ virtual orbitals that can be occupied by at most $p$ electrons. Thus, we can herein uniquely specify each AS as RAS$(v;h,p)$. The ASs used in this study are: RAS$(8;1,1)$, containing $3s$, $3p$, $4s$, and $4p$ orbitals in the RAS3; RAS$(26;1,1)$, enlarging RAS3 by the $3d$, $4d$, $5s$, $5p$, $6s$, and $6p$ orbitals and RAS$(33;1,1)$, adding the $4f$ orbitals. The number of configurations possible with each AS is shown in Table 2.

For all QC schemes (see Table 1) all states are included in the RASSCF procedure for the Ne wave functions, while the core excited states have been excluded for the calculations of Ne+.

Atomic Natural Orbital (ANO) type basis sets have been employed using a $(22s17p12d11f)$ primitive set. It was constructed by supplementing the ANO exponents for neon Widmark et al. (1990) in each angular momentum with eight Rydberg exponents generated according to the scheme proposed by Kaufmann et al. Kaufmann et al. (1989). The contractions were then obtained with the GENANO module Almlöf and Taylor (1987) of OpenMolcas Galván et al. (2019) using density matrices from state averaged RASSCF calculations for Ne and Ne+ with the RAS$(33;1,1)$ AS. All possible 166 states for Ne have been taken into account but for Ne+, the core excited manifold has been excluded, leading to 532 states. To get the final basis set, the sets obtained for Ne and Ne+ have been evenly averaged. In this work, we use the contractions: [7s6p3d2f], obtained by the procedure described above; [9s8p5d4f], corresponding to [7s6p3d2f] supplemented with Rydberg contractions that resulted from the GENANO procedure using the ”rydberg” keyword; [22s6p3d2f]-rcc, similiar to [7s6p3d2f] but with uncontracted $s$ functions and scalar relativistic corrections according to the Exact Two Component Decoupling (X2C) scheme Peng and Reiher (2012), see Supplement: Sec. III.

All energies have been corrected using the single state RASPT2 method Malmqvist et al. (2008) with an imaginary shift of $0.01$ a.u. To be consistent with the basis set generation, scalar relativistic effects have only been included for QC II. The State Interaction (RASSI) Malmqvist et al. (2002) module of MOLCAS was used to compute the biorthonormally transformed orbital and CI coefficients Malmqvist (1986) for the atomic and ionic states: $\{\varphi_{i}\},\{C_{j}\}\rightarrow\{\tilde{\varphi}_{i}\},\{\tilde{C}_{j}\}$ such that $\innerproduct{\tilde{\varphi}^{+}_{i}}{\tilde{\varphi}_{j}}=\delta_{ij}$, while the total wave functions remain unchanged. This biorthonormal basis is used in the evaluation of the Dyson orbitals, Eqs. (14), (28), (29), and two-electron reduced transition densities, Eq. (17).

IV.2 Matrix elements in the atomic basis

The matrix elements between bound orbitals occuring in Eqs.(12), (15), and (16) are evaluated by transforming the orbitals to the atomic basis and calculating the atomic basis integrals with the libcint library Sun (2015). The most time consuming part in the computation of the partial decay rates using Eq. (12) is the estimation of the two-electron continuum-bound integrals. Transforming the two-electron reduced transition densities and the orbitals to the atomic basis $\{\chi_{i}\}$ and neglecting the spin integration, the two-electron continuum-bound integrals used in the practical evaluation of Eq. (16) read as:

[TABLE]

The function $f_{ac}(\mathbf{r}_{1})$ is similar to an atomic nuclear attraction integral, and is evaluated as a function of $\mathbf{r}_{1}$ using the libcint library Sun (2015). An important point is to exploit the fact that the kinetic energy of the continuum electron is only encoded in its radial part. Transforming to spherical coordinates $\mathbf{r}\rightarrow(r,\Omega)$ centered at the origin of the outgoing electron allows to separate off the radial integration:

[TABLE]

The angular integral in $F_{abc}(r)$ is determined numerically by using the adaptive two-dimensional integration routine cuhre from the Cuba 4.2 library Hahn (2005). To reduce the number of points at which $F_{abc}(r)$ is evaluated, an adaptive spline interpolation is used, which was developed by us and implemented in our code. Therein, the grid spacing is adjusted such that the absolute error estimate of the interpolation is kept lower than $10^{-6}$ a.u. on each region with a different spacing. Note that the $F_{abc}(r)$ are determined only once, while the final radial integration in Eq. (IV.2) needs to be evaluated for every transition $i\rightarrow\alpha$. The radial integration in Eq. (IV.2) is carried out using the Simpson rule. For the one-electron continuum-bound integrals that are needed in the evaluation of the one-electron matrix elements in Eq. (12), an analogous approach has been implemented. This protocol has been developed aiming for the application to molecular systems. Hence, it does not exploit that the atomic orbitals are eigenstates of angular momentum operator, which would be usual in a purely atomic approach.

V Results and Discussion

Here, we perform a thorough benchmark of the approaches to evaluate AES presented in Section II on the exemplary Auger decay of the neon $1s^{-1}3p$ resonance. First, in Section V.1, we compare the AES modelled with our protocol against data obtained from an atomic MCDF calculation Stock et al. (2017), which serves as a high level theoretical reference. Second, in Section V.2, we undertake the comparison to experimental results. The comparison against both theory and experiment is needed since no uniform and highly resolved experimental data covering the full energy range discussed herein has been published to date.

V.1 Benchmark of theoretical models

Panel (a) of Fig. 2 shows neon AESs, resulting from the autoionization of the $1s^{-1}3p$ states, obtained using bound state wave functions from QC scheme I (cf. Table 1) with radial waves corresponding to the spherically averaged direct-exchange potential $V^{\text{JX}}_{f}(r)$ Eq. (7e). The partial rates have been evaluated using the full $\mathcal{H}$ coupling as well as the approximate $r^{-1}$ coupling in Eq. (12). Further, the nonorthogonality of the continuum and bound orbitals was accounted for by including the overlap terms, NO, see Eqs. (13), (15), and (16). A spectrum employing $r^{-1}$ coupling at the MCDF level, obtained by Stock et al. Stock et al. (2017) with the RATIP package Fritzsche (2012), serves as a theoretical benchmark. Therein, the atomic structure was obtained with a configuration space including single electron excitations from the $1s$, $2s$, and $2p$ to the $np$ up to $n=7$ and $3d$ orbitals. The four-component continuum orbitals have been obtained as distorted waves within the potential of the respective ionized atomic state.

All spectra have been constructed by assigning a Gaussian lineshape with an FWHM of $\gamma=0.1$ eV to each channel.

[TABLE]

Where $G(\varepsilon,\gamma)=\sqrt{\ln{2}/(\pi\gamma)}\,\,\mathrm{exp}(-{4\ln{2}}\varepsilon^{2}/{\gamma^{2}})$, $\varepsilon$ is the kinetic energy of the emitted electrons and $\varepsilon_{\alpha}=\mathcal{E}_{1s^{-1}3p}-\mathcal{E}_{f}$ is the Auger electron energy of the channel $\ket{\Psi_{\alpha}}$. The spectra were normalized to the peak height at $811.5$ eV, and shifted globally by $-5.35$ eV (QC I) and $-2.45$ eV (MCDF) such that the peak at $811.5$ eV is aligned to the experimental data taken from Kivimäki et al. Kivimäki et al. (2001) (see Fig. 4). In addition, the dominant continuum orbital angular momentum contribution to the intensity is indicated for each peak. This shows that the regions $743-765$ eV, $765-800$ eV and $800-825$ eV correspond almost exclusively to the emission of $s$, $p$ and $d$ waves, respectively, with the exception that the peaks at $803$ eV and $808$ eV lie in the $d$ region but are due to $s$ wave emission. Hence, we will refer to this regions as $s$, $p$, $d$ rather than using the energies in what follows.

It is evident from Fig. 2 (a) that the normalized AES obtained for the $\mathcal{H}$ and $r^{-1}$ couplings are almost indistinguishable in the $s$ and $d$ regions. In contrast, the features in the $p$ region at $771$ eV and $776$ eV are enhanced by $10\%$, while the small bands in the range $783-788$ eV, are reduced by about $40\%$, if the $\mathcal{H}$ coupling terms are included. The overall effect on the spectrum, however, is still rather small and neglecting the contributions beyond the $r^{-1}$ coupling in the partial rate evaluation seems to be justified in this case.

The comparison with the MCDF data in the $d$ region shows that the Auger energies and relative intensities are overall quite well reproduced by our approach. Only the intensities of two small satellite peaks at $803$ eV and $809$ eV are overestimated and three tiny features around $805.5-807.5$ eV are not present when using our approach. The latter deficiency can be safely attributed to the smaller configuration space employed in the QC I scheme as compared to the one used to obtain the MCDF results in Ref. 68. Looking at the $s$ and $p$ regions, the Auger energies from the QC model I become slightly but increasingly blue shifted with respect to the MCDF ones at the lower energy flank of the spectrum. The intensities in turn are considerably overestimated by factors of about five and two for the $s$ and $p$ regions, respectively. Note that very similiar spectra are obtained if the RASPT2 energy correction is not used (Supplement: Fig. S1)

Panels (b) and (c) of Fig. 2 show the integrated decay rates for the $s$, $p$ and $d$ regions evaluated using QC model I with the $r^{-1}$ (b) and $\mathcal{H}$ couplings (c), comparing different approaches to compute the partial decay rates. The respective spectra are given in Supplement: Figs. S2-S7. The rates were obtained for all combinations of potential models in Eqs. (7a) - (7e) with the different nonorthogonality approaches SO, NO, and GS for the continuum orbitals. Since no absolute rates for the MCDF results are available, these have been scaled to the total decay rate obtained using the $r^{-1}\cdot V_{f}^{\text{JX}}(r)\cdot\text{NO}$ treatment. The data in panels (b) and (c) show that the decay rates corresponding to the $s$, $p$ and $d$ spectral regions converge for both couplings as the quality of the potential is increased from the free particle approximation $V^{\text{free}}$ to the spherically averaged direct-exchange potential $V_{f}^{\text{JX}}(r)$. The $s:p:d$ - ratio obtained from the MCDF spectrum, however, is not matched. Our approaches systematically overestimate the decay rates due to the $s$ and $p$ channels, in line with the mismatch of intensities unveiled in panel (a).

Prominently, using the approximate $r^{-1}$ coupling together with the NO overlap terms leads to an overestimation of the $s$ and $p$ regions with respect to the SO results by an order of magnitude for $V^{\text{free}}$ and $-1/r$, and by factors of $4$ and $3$ for $-6/r$. In contrast, employing the SO approximation or GS orthogonalized continuum orbitals results in more realistic decay rates. Here, the GS rates reproduce the SO ones for the $d$ region but are considerably smaller in the $s$ and $p$ regions. However, when either of the more physically sound potentials $V_{f}^{\text{scr}}(r)$, $V_{f}^{\text{J}}(r)$, or $V_{f}^{\text{JX}}(r)$ is used, the choice of the nonorthogonality treatment has only a weak influence on the integrated decay rates. This characteristic is due to the continuum-bound orbital overlaps that decrease if more realistic potentials are chosen, causing the NO overlap terms in Eqs. (13), (15), and (16), as well as the effect of the GS-orthogonalization to become negligible. Adding to this, the general insensitivity of the $d$ region to the nonorthogonality treatment is explained by the fact that the overlap of the $d$ continuum waves with the bound orbitals vanishes independently of the potential model. The overlap integrals are in Supplement: Table S1.

Generally, the inclusion of the $\mathcal{H}$ coupling, Fig. 2 (c), yields mostly smaller total decay rates and less pronounced differences between the SO, NO and GS results. The SO decay rates are usually the largest when $\mathcal{H}$ coupling is applied, while the NO and GS ones are smaller, underestimating the $s$ and $p$ regions for $V^{\text{free}}$, $-1/r$, and $-6/r$, with respect to SO. Noteworthy, the huge overestimation due to the combination of $r^{-1}$ coupling and NO terms is mitigated completely. The reason is that the NO contribution to the matrix element is determined by cancellation effects between different one- and two-electron NO terms rather than their individual magnitude. Neglecting the one-electron contributions in $r^{-1}$ coupling prohibits this error cancellation, leading to the immense overestimation observed for $V^{\text{free}}$, $-1/r$, and $-6/r$. Similar, the agreement between the GS and NO results for these potentials is not because the effect of the GS-orthogonalization is minute, but again due to the cancellation.

We conclude the discussion of this figure with the observation that the total rates obtained with the $V_{f}^{\text{J}}(r)$, and $V_{f}^{\text{JX}}(r)$ potentials reproduce the experimentally observed decay rate Avaldi et al. (1995). Independent of the chosen approach, the best match is obtained with the $V_{f}^{\text{J}}(r)$ potential, whereas the $V_{f}^{\text{JX}}(r)$ leads to slightly underestimated total rates, which can be attributed to the inclusion of the attractive exchange term, Eq. (25). Similarly, we assign the considerable underestimation of the total decay rates by the $V_{f}^{\text{scr}}(r)$ potential to the fact that it is more attractive in the core region than $V_{f}^{\text{JX}}(r)$, see Fig. 3.

The $-6/r$ potential yields total rates and spectra that are in approximate agreement with those obtained for $V_{f}^{\text{J}}(r)$ and $V_{f}^{\text{JX}}(r)$ apart from the case when $r^{-1}$ coupling and the NO terms are used (spectra in Supplement: Figs. S2-S7). Notably, the total rates obtained in $\mathcal{H}$ coupling with the effective potentials $V^{\text{free}}$ and $-1/r$ agree well with the experimental reference. However, despite the good agreement in the total rates the spectra for these approaches deviate in the $p$-region considerably from the ones obtained for the more accurate potentials (Supplement: Figs. S5-S7), indicating that this agreement is rather accidental.

To shed light on the influence of different potentials on the total decay rates and spectra obtained with QC model I, $\mathcal{H}$ coupling, and the NO approach, potentials, radial continuum functions and spectra are presented in Fig. 3.

Panel (a) contains the normalized AESs obtained using the potentials $-1/r$, $-6/r$, $V_{f}^{\text{scr}}(r)$, and $V_{f}^{\text{JX}}(r)$. The spectra have been shifted vertically for the sake of clarity. Shifts and broadening parameters are the same as in Fig. 2 (a). The $d$ region is represented very well with all potential models, with the exclusion that the satellite bands at $803$ eV and $808$ eV are barely present when using the $-1/r$ potential. In contrast, the relative intensities of the $s$ and $p$ regions of the spectra are strongly affected by the choice of the potential. Here, the $p$ region comprises two peak groups with different character.

Using the $V_{f}^{\text{JX}}(r)$ results as a reference for this discussion, the peaks around $771$ eV and $776$ eV are overestimated by $30\%$ with $V_{f}^{\text{scr}}(r)$, almost reproduced with $-6/r$ and underestimated by a factor of four when using $-1/r$. In contrast, the bands between $783$ eV - $788$ eV behave in an inverse manner. Namely, they are underestimated by $50\%$ with $V_{f}^{\text{scr}}(r)$ and overestimated by $50\%$ and a factor of four with $-6/r$ and $-1/r$, respectively. This behavior of the screened and the effective Coulomb potentials, is caused by the fact that $V_{f}^{\text{scr}}(r)$ is more attractive and $-1/r$ as well as $-6/r$ – less attractive than $V_{f}^{\text{JX}}(r)$, as depicted in panels (b)-(d). This trend is solely dictated by the potential but not the coupling or nonorthogonality treatment; note the exception of $r^{-1}\cdot\text{NO}$ combination (Supplement: Figs. S2-S7). Finally in the $s$ region, the tighter potential results in larger intensities: $-1/r$ and $-6/r$ lead to an underestimation by an order of magnitude and $35\%$, respectively, whereas $V_{f}^{\text{scr}}(r)$ yields an overestimation by $25\%$.

For each region, the radial waves and respective potentials, including the angular momentum term, are shown in the panels (b)-(d) corresponding to the characteristic peaks denoted in panel (a).

Note that the $V^{\text{free}}$ and $V_{f}^{\text{J}}(r)$ cases are very similar to the $-1/r$ and $V_{f}^{\text{JX}}(r)$ ones, respectively, and have been excluded from the discussion of this figure. The respective spectra are in Supplement: Fig. S7.

One might understand the differences in the sensitivity of the $s$, $p$ and $d$ spectral regions to the potential by comparing the short range behavior of the radial waves $w_{l}^{fk}(r)/r$, Eq. (11a), with the electron density that is sharply peaked in the core region. This suggests that the matrix elements in Eq. (12) are very sensitive to the description of the continuum orbitals in this region. It is well known that only the $s$ waves have a considerable contribution at the core, since the radial functions tend to zero as $r^{l}$. In fact, the effective radial potential at the core is dominated by the angular momentum term for $l>0$, meaning that the influence of the present potential models on the total rates should decrease with increasing $l$. Regarding the relative changes of the integrated decay rates within each region, this is true for $\mathcal{H}$ and $r^{-1}$ coupling, when the NO overlap terms are included, irrespective of whether the GS orthogonalization is used in addition or not, Fig 2 (b) and (c). However, under the SO approximation for $r^{-1}$ couplings, the $p$ spectral region is more sensitive to variations in the potential model than the $s$ region.

In fact, for the $s$ waves it is seen that only the potentials $V^{\text{scr}}_{f}(r)$ and $V_{f}^{\text{JX}}(r)$ lead to similar radial waves, the slight differences being due to the fact that the screening model is too attractive in the core region. Notably, these slight deviations lead to an underestimation of the total decay rates obtained with the $V^{\text{scr}}_{f}(r)$ potential by about $25\%$ in comparison to those obtained for $V_{f}^{\text{JX}}(r)$ (see Fig. 2), underlining the sensitivity of the total decay rates to the choice of the model potential.

The effective Coulomb potential $-1/r$ provides a qualitatively wrong description in both, the core and outer regions, whereas the $-6/r$ potential leads to a sort of compromise in accuracy. It describes the core region much better than the $-1/r$ potential but as a trade off has a wrong asymptotic behavior in the valence region.

For all but the $r^{-1}\cdot\text{NO}$ approaches this suffices to predict spectra that are in qualitative agreement with the ones obtained for $V_{f}^{\text{JX}}(r)$, whereas using the $-1/r$ potential is only justified for the main features of the $d$ region, see Supplement: Figs. S2-S7.

If one is interested in the full spectrum, one should use a model taking into account the electronic potential of the ionized core, such as $V^{\rm scr}(r)$, $V_{f}^{\text{J}}(r)$ or $V_{f}^{\text{JX}}(r)$.

Summarizing the discussion until this point, it seems that satisfactory total decay rates and AESs are only obtained with the potentials $V_{f}^{\text{J}}(r)$ and $V_{f}^{\text{JX}}(r)$. Further, employing $r^{-1}$ coupling and the SO approximation seems to be very well justified for this models and we will use these in the comparison against experimental data below.

Generally, special caution has to be taken with the inclusion of the NO overlap terms for simple potential models. While the results obtained using $\mathcal{H}$ coupling greatly profit from error cancellation, this is not the case for $r^{-1}$ coupling, where the NO-terms strongly emphasize the deficiencies of an approximate potential. In addition, it remains to be clarified, whether this cancellation effects are a general feature, or a peculiarity of the neon $1s^{-1}3p$ AES.

V.2 Comparison to experimental data

In this section, the comparison of our theoretical results to the experimental data in the full spectral range is presented. In addition, to unravel the influence of the underlying QC onto the AES, we discuss spectra obtained with the QC schemes I-III as described in Table 1. The QC model II, which is more sophisticated than QC I, contains an uncontracted $s$ basis and accounts for scalar relativistic effects, whereas III employs an active space larger than in QC I .

To the best of our knowledge, no experimental Auger emission spectrum of the neon $1s^{-1}3p$ resonance that covers the full spectral range presented in Figs. 2 and 3 has been published to date. Hence, in Fig. 4 the spectra obtained with the QC models I-III and the $r^{-1}\cdot V_{f}^{\text{JX}}(r)\cdot\text{SO}$ approach are compared against experimental data taken from Kivimäki et al. Kivimäki et al. (2001), for the $d$ region (a), and Yoshida et al. Yoshida et al. (2000), for the $s$ and $p$ regions (b). Note that the $\mathcal{H}\cdot V_{f}^{\text{JX}}\cdot\text{NO}$ approach does not lead to a considerable improvement of the agreement with the experimental data, as shown in Supplement: Fig. S9. The spectra have been shifted by $-5.35$ eV, $-4.75$ eV, and $-5.18$ eV for QC models I, II, and III to align the peaks at $811.5$ eV. In panel (a), the spectra have been broadened using a Gaussian profile with an FWHM of $0.25$ eV, corresponding to the lineshape of the peak at $806.5$ eV in the experimental spectrum. Further, in panel (b) a Gaussian FWHM of $0.77$ eV was chosen to represent the shape of the high energy flank of the asymmetric peak at $778.5$ eV in the experimental spectrum. Additionally, the spectra have been normalized to the heights of the main peaks at $811.5$ eV (a) and $776$ eV (b). Finally, the experimental data have been recorded at angles of $54.7^{\circ}$ (a) and $56^{\circ}$ (b) to the polarization vector. Data for these angles has been chosen to rule out anisotropy effects which makes the comparison with our angle integrated spectra easier.

Panel (a) shows that the agreement between the experimental and the theoretical spectra is fairly good. In fact, the relative energetic positions and intensities of the main peaks are reproduced quite well. However, smaller satellite peaks at $803$ eV and $808.2$ eV are blue shifted by about $0.5$ eV. Generally, the agreement between theory and experiment is worse for the smaller features, although an unambiguous assignment is still possible. Using the uncontracted $s$ basis and including scalar relativistic effects in QC II does not influence the resulting spectrum. In turn, the larger active space incorporated in the QC scheme III leads to a better reproduction of some tiny features at $806.4$ eV, $807.5$ eV and $808.8$ eV. The blue shifts of the peaks already present with QC I and II are not affected, when employing QC III.

The comparison of the spectra covering the $s$ and $p$ regions with the experimental data is shown in panel (b). Here, the overall agreement is worse than in panel (a) both concerning the relative intensities and energetic positions of the peaks. Regarding the offset in energies, the low-energy part corresponds to transitions to the highest excited states of the ionic manifold, the energies of which are progressively overestimated. This is a typical situation when a smaller part of the electron correlation is recovered for the higher-lying states. Specifically, the positions of the peaks around $783$ eV and $788$ eV are reproduced well by all methods, while the intensity of the former, small peak is overestimated. Further, the peaks at $771$ eV and $776$ eV are red shifted by about $2$ eV, and the relative intensity of the former is overestimated by approximately $30\%$ (QC I and II) and $15\%$ (QC III). Finally, the peaks at $745$ eV and around $751$ eV, corresponding to the $s$ region are redshifted by $5$ eV (QC I-III) and $4$ eV (QC I, II), or $5$ eV (QC III), respectively. The intensities of these peaks are overestimated by about $30\%$ with respect to the experimental data. Hence, in the $s$ and $p$ regions, our results agree better with the experimental reference than with the MCDF spectrum. The total decay rates, however, are not visibly altered by the choice of either of the QC schemes I, II or III (cf. Supplement: Fig. S8).

To wrap up this discussion, we conclude that the RASSCF/RASPT2 electronic structure method combined with the $r^{-1}\cdot V_{f}^{\text{JX}}\cdot\text{SO}$ approach to construct continuum orbitals and evaluate the transition matrix elements provides Auger energies and intensities for the decay of the neon $1s^{-1}3p$ resonance of a similar quality as those obtained in Ref. 68 with the MCDF approach. In particular, a straightforward assignment of experimental results is possible. Further, the inclusion of scalar relativistic effects into the one-component electronic structure of the bound states has no notable influence on the spectra, while a large active space is necessary only to reproduce minor satellite features.

VI Conclusions and Outlook

In this work, we have demonstrated an approach to the evaluation of autoionization rates on the example of the Auger decay from the neon $1s^{-1}3p$ resonance. The suggested protocol is based on the RASSCF/RASPT2 method to evaluate the bound state wave functions and energies, supplemented by a single-channel scattering model for the outgoing electron. Here, the single-center approximation is introduced to reduce the continuum orbital problem to the radial dimension by averaging over the angular structure of the ionized electron density. To model the true radial potential, six different models, $V^{\text{free}}$, $-1/r$, $-6/r$, $V^{\text{scr}}_{f}(r)$, $V^{\text{J}}_{f}(r)$, and $V^{\text{JX}}_{f}(r)$, have been discussed. Further, three different ways to account for the nonorthogonality of the continuum and bound orbitals, SO, NO, and GS, as well as the effect of using complete, $\mathcal{H}$, or approximate, $r^{-1}$, coupling in the partial rate evaluation has been investigated. All combinations of these sum up to 36 different variants to evaluate partial autoionization rates for an underlying bound state QC calculation and have been implemented in a standalone program.

Here we compared all these approaches with respect to their ability to reproduce the experimental Kivimäki et al. (2001); Yoshida et al. (2000) as well as theoretical AES obtained at the fully relativistic MCDF level Stock et al. (2017). The applied quantum chemistry protocol allows for a fairly good reproduction of the transition energies if compared to the theoretical reference and experiments. However, intensities are more difficult to reproduce. Here, the quality of the continuum orbital was shown to be the most important issue as it strongly influences the obtained AESs. Especially the core part of the ionic potential is of importance for AES. The effect of the potential is closely connected to the angular momentum of the outgoing electron as it governs the extent of the continuum orbital into the core region. For instance, we found that the $d$ region of the spectrum is rather insensitive to the choice of the model potential, whereas the $s$ and $p$ regions require to use one of the potentials $V^{\text{scr}}_{f}(r)$, $V_{f}^{\text{J}}(r)$ and $V^{\text{JX}}_{f}(r)$. Still, the MCDF intensities can only be reproduced in the $d$ region of the spectrum, while they are overestimated in the $s$ and $p$ parts. Further, the free particle model and the asymptotic Coulomb potential $-1/r$ fail to reproduce the complete spectrum.

Interestingly, inclusion of the NO terms in addition to using $r^{-1}$ coupling in the SO approximation does not in general lead to improved spectra, but rather strongly emphasizes the deficiencies of the $V^{\text{free}}$, $-1/r$ and $-6/r$ potentials. Due to pronounced error cancellation, however, this is diminished to large extent if the $\mathcal{H}$ coupling is used with the NO terms. Remarkably, the effective $-6/r$ potential yields qualitative agreement with the spectra obtained using the more accurate potentials for all approaches but the aforementioned combination of $r^{-1}$ coupling and NO terms.

In contrast, the spectra obtained with $V^{\text{scr}}_{f}(r)$, $V_{f}^{\text{J}}(r)$ and $V^{\text{JX}}_{f}(r)$, are weakly affected by the choice of both the coupling and nonorthogonality treatment.

Since it remains unclear whether the cancellation effects observed when using $\mathcal{H}$ coupling with the NO terms are a general feature, or a peculiarity of the neon AES, the NO terms should only be employed with caution. Generally, they should only be used, if the description of the potential and radial waves are sufficiently accurate. When dealing with approximate potentials however, our results suggest to employ the SO or GS approaches, that seem to be less sensitive to the quality of the continuum orbital. In addition, due to the computational simplicity, the SO approximation could be preferred.

The comparison with experimentally obtained spectra using the $r^{-1}\cdot V_{f}^{\text{JX}}(r)\cdot\text{SO}$ approach demonstrated the ability of the present method to accurately predict the neon $1s^{-1}3p$ AES, allowing a straightforward assignment of the experimental data. Interestingly, the general structure of the spectrum can be already reproduced quite well using a rather small active space, and is not sensitive to the inclusion of scalar relativistic effects. The best agreement with the experimental data is achieved by using a larger active space, including additional excitations to $3d$, $5s$, $5p$, $6s$, and $6p$ orbitals together with the spherically averaged direct exchange potential $V_{f}^{\text{JX}}(r)$. In addition our approach not only reproduces the experimentally measured AESs but the total decay rates of the neon $1s^{-1}3p$ resonance as well, when the potentials $-6/r$, $V_{f}^{\text{J}}(r)$, and $V^{\text{JX}}_{f}(r)$ are used. Since using the screened charge potential $V_{f}^{\text{scr}}(r)$ leads to notably underestimated absolute rates and is not computationally cheaper than using either $V_{f}^{\text{J}}(r)$ or $V^{\text{JX}}_{f}(r)$, providing a better accuracy, the latter two should be preferred.

Molecular systems can be treated with the presented method as well, however, in this case the molecular continuum has to be approximated by a single-centered spherically symmetric model. To keep the errors due to this approximation as small as possible, our findings suggest to use the $r^{-1}$ coupling together with the SO approximation to evaluate molecular AIS. The potentials should be modeled using either the direct $V_{f}^{\text{J}}(r)$ or direct-exchange $V_{f}^{\text{JX}}(r)$ variant. The applicability of these approximations has to be tested for the molecular case.

Currently our code is interfaced to the MOLCAS/openMolcas as well as to the Gaussian program packages, allowing to evaluate PES and AIS based on bound state calculations conducted with the RASSCF/RASPT2 Grell et al. (2015), as well as the linear-response time-dependent density functional theory method Möhle et al. (2018). A publication discussing the applicability of the present models to treat molecular systems will follow.

Appendix A Obtaining the model potentials

For the state-dependent models, Eqs. (7c)-(7e), the central quantity is the spherically-averaged electron density of the ionized state

[TABLE]

It determines the potentials in the following way: The screened charge $Z_{f}(r)$ in Eq. (7c) is evaluated as the difference between the nuclear charge $Z$ and the integrated number of electrons present in a sphere with radius $r$ around the atom:

[TABLE]

Further, the direct Coulomb potential of the ionized states in Eqs. (7d) and (7e) is the electrostatic potential of the spherically averaged electron density $\rho_{f}(r)$,

[TABLE]

$J_{f}(0)$ is defined by the asymptotic value of the integral over $r^{\prime\prime}$ in Eq. (23):

[TABLE]

Finally, a radial Slater type exchange Slater (1951),

[TABLE]

is employed.

Appendix B Coupling matrix elements

The total Auger transition matrix element $A_{i\alpha}^{\sigma M^{+}}=\matrixelement{\Upsilon_{\alpha}^{\sigma M^{+}}}{\mathcal{H}-\mathcal{E}_{i}}{\Psi_{i}}$ reads

[TABLE]

Here the contributions corresponding to the SO approximation, the NO terms and the $r^{-1}$ coupling have been indicated. If the GS approach is used, $\ket{\psi_{\alpha,\sigma}}$ is replaced by

[TABLE]

The conjugated one and two-electron Dyson orbitals are defined as:

[TABLE]

and

[TABLE]

Note that the sum of both terms, which occurs when the $\mathcal{H}$ coupling is used, takes the simple form

[TABLE]

as demonstrated in Manne and Ågren (1985), if the Hamiltonian eigenvalue equation $\mathcal{H}\ket{\Psi^{+}_{f,M^{+}}}=\mathcal{E}_{f}\ket{\Psi^{+}_{f,M^{+}}}$ is used. However, as detailed in the Supplement: Section II, both approaches are strictly equivalent only, if the occupied spin-orbitals of the bound ionized and unionized wave functions span the same space. This is generally also not fulfilled for RASSCF orbitals if separate SCF procedures are used to obtain ionized and unionized states since the respective transformation properties between orbitals from different subspaces are lost. Moreover, the formulation of Ref. 76 cannot be applied to the $r^{-1}$ coupling case whereas our formulation can.

Acknowledgements.

We would like to thank Prof. Dr. Stephan Fritzsche, from the Helmholtz Institute Jena, Randolf Beerwerth and Sebastian Stock from his group for fruitful discussions and for providing the unshifted results of the MCDF calculations from Ref. 68. Financial support from the Deutsche Forschungsgemeinschaft Grant No. BO 4915/1-1 is gratefully acknowledged.