Ab initio calculation of the electron capture spectrum of ¹⁶³Ho: Auger–Meitner decay into continuum states

If one looks at the spectral shape of the differential electron capture decay rate one finds that these spectra resemble x-ray core level photo electron spectra (XPES). At the edges one finds an enhanced electron capture decay rate as here the electronic states are at resonance. Above an edge this decay rate is higher than below the edge, because additional channels for electronic excitations are possible. Excess energy can be transferred to an emitted, free electron via the Auger–Meitner process. Hereby the Dy atom is left in an ionized state. These final states are similar to those reached in a photo electron emission event, thereby explaining the similarity between the XPES and EC line-shape. If one simplifies the free electron density of states as well as the Coulomb matrix elements coupling the bound-states to states with an Auger–Meitner electron, one can derive the resulting line-shape to be Mahan like [14, 27]. As we want a highly accurate description of the EC spectral shape, one should not make these approximations, but include the details of the decay channels in full complexity. The free electron density of states, as well as the Auger–Meitner Coulomb integrals can be calculated such that one can derive the EC line-shape without these approximations. This leads to equation (5), where the broadening now is given by a state dependent self energy as defined in equation (4).

In order to understand how this results in the line-shape given in figure 1 we here focus on a single decay channel. We define ψ_b to be a single bound excited state of Dy, reached after an electron capture event. This state couples to the continuum via equation (4), defining the self energy Σ_bb (ω). The resulting differential electron capture decay rate for state ψ_b becomes:

$$\begin{equation}\frac{\mathrm{d}{{\Gamma}}_{b}}{\mathrm{d}\omega }\propto \mathrm{Im}\left(\frac{1}{\omega -{E}_{b}-{{\Sigma}}_{bb}\left(\omega \right)}-\frac{1}{\omega +{E}_{b}+{{\Sigma}}_{bb}\left(\omega \right)}\right)\end{equation} \tag{ 6 }$$

with E_b the excitation energy of the state ψ_b with respect to the Ho ground-state. According to Fermi's golden rule, the life-time of resonance ψ_b is proportional to the inverse transition rate at the resonance energy ${\tau }_{b}^{-1}\propto -\mathrm{Im}\left[{{\Sigma}}_{bb}\left({E}_{b}\right)\right]$. Hence, life-time and energy-dependent broadening are encoded in the diagonal elements of the self-energy. Furthermore, following equation (4) the diagonal parts of the self energy depend on the operator U_A which includes the Coulomb matrix elements between free Auger–Meitner electrons and electrons in orbitals describing the bound-state ψ_b .

To study the impact of those Coulomb matrix elements we focus on a specific case. As an example we take a bound state that belongs to the 4s edge of the EC-spectrum. We denote this state by $\vert \left[\mathrm{H}\mathrm{o}\right]\underline{4s}\rangle$ where the underscore implies that there is a core hole in the 4s orbital due to electron capture. Note that this is a multi-slater-determinant state, but we will be concerned with this technicality later. This state can decay via many different channels. We focus on the channel where a 4p electron fills the hole and transfers the energy-difference to a 4d electron which leaves the atom with a kinetic energy and angular momentum l = 3

$$\begin{equation}\left[\mathrm{H}\mathrm{o}\right]\underline{4s}\to \left[\mathrm{H}\mathrm{o}\right]\underline{4p}\underline{4d}+{\mathrm{e}}_{{\epsilon}f}^{-}.\end{equation} \tag{ 7 }$$

The coupling strength of this decay is determined by its corresponding Coulomb–Slater integral. This integral depends on the energy of the free electron (). The size of this integral for the specific case of our example (4p and 4d scatter to 4s and _f ) is shown in the left panel of figure 2 as a function of the Auger–Meitner electron's kinetic energy. One can observe a very distinct energy dependence, with a clear maximum at 322  eV and a zero at 2037  eV. We can understand this energy dependence by looking at the wave-functions involved. The Coulomb–Slater integrals include the product of the bound radial wave-function (4d orbital) and the Auger–Meitner electron wave function. The right panel of figure 2 shows the radial wave functions for the bound 4d_5/2 orbital with dashed black lines and for unbound states with three different kinetic energies going from small (top panel) to large (bottom panel). The kinetic energies of the radial wave functions used in the right panels is indicated in the left panel. For Auger–Meitner electrons with a low kinetic energy the unbound electron has a very small amplitude (red line, top panel) in the vicinity of bound 4d electrons. This leads to a small matrix element for Auger–Meitner decay at low kinetic energy. With increasing energy and decreasing wavelength of the free electron the amplitude of its wave-function in the region of the bound orbitals increases (blue line, middle panel). With this also the size of the Coulomb matrix elements increases. At a kinetic energy of the free electron of 322 eV, the channel indicated in equation (7) has a maximal scattering amplitude. For even higher kinetic energies of the free electron and corresponding smaller wave-lengths (orange line, bottom panel), oscillations of the free electron wave-function in the spatial region of the bound electron reduce the Coulomb–Slater integrals. At 2037  eV this integral becomes zero. Beyond this point the oscillatory behavior of the Auger–Meitner electron's wave-function increases the Slater integral's value again, until the next minimum at even higher kinetic energy and smaller wave-lengths is reached (not shown).

Zoom In Zoom Out Reset image size

Once the Auger–Meitner energy dependent Coulomb–Slater matrix elements are determined, one can calculate the self energy of state $\vert {\psi }_{b}\rangle =\vert \left[\mathrm{H}\mathrm{o}\right]\underline{4s}\rangle$ for the single channel used in our example. The operator U_A in (4) couples state ψ_b to states with one additional Auger–Meitner electron which can have all possible kinetic energies. As these states are independent of each other, one can sum them separately and obtain the full self energy; i.e. the cross section of Auger–Meitner scattering starting from the state $\left[\mathrm{H}\mathrm{o}\right]\underline{4s}$ as shown in figure 3(a). At energies below 400  eV the self energy is small, as there is not enough energy to generate free electron states combined with a 4p core hole. Above 400  eV the 4s resonance couples to a continuum of Auger–Meitner-states due to Coulomb interaction. The self energy has a maximum around 750  eV, slowly decays and becomes zero around 2500  eV. This is directly related to the Slater integral given in figure 2 which shows, as a function of the kinetic energy of the free electron state, a similar behavior, albeit shifted by the energy of the 4p core hole.

Zoom In Zoom Out Reset image size

With the self energy of this specific state for the given single decay channel one can now calculate the differential decay rate. Figure 3(b) shows in red the resulting spectrum. In grey the spectrum represented by a single Lorentzian is shown. One can clearly see the asymmetric line-shape of the resonance, which is a consequence of the energy-dependent Coulomb–Slater matrix elements and the resulting energy dependent self energy. The specific energy-dependence of the self energy in figure 3(a) induces not only the asymmetric broadening, but also the small bump to the right of the peak, as well as the increased spectral weight in the resonance's wing. Even the second rise in self energy above 2500  eV modifies the spectrum, although this effect is suppressed by the neutrino phase space factor.

In order to calculate the spectral shape as shown in figure 1, one needs to include all possible bound-states ψ_b and all possible decay channels where one bound electron scatters into the created core hole and one bound electron scatters into a free electron state. While the energy dependence of these processes can be understood qualitatively by the overlap between single-particle wave-functions before and after the Auger–Meitner decay, a quantitative treatment needs to include the multi-configurational nature of initial and final states.

In the previous subsection we demonstrated how energy-dependent Coulomb matrix elements affect the line broadening. We focused on the decay of a single state only, which is encoded in the corresponding diagonal entry of the self-energy as defined in equation (4). Including multiple scattering channels leads to mixing between bound-state resonances, which in turn is reflected in the off-diagonal elements of the self energy. These describe how different resonances are coupled to the same final states via different scattering channels. For example one can reach a state with a hole in the 4p and 4d orbital and a free electron by first absorbing a 4s electron into the nucleus followed by an Auger–Meitner decay, scattering a 4p electron into the 4s orbital and emitting a 4d electron. One can reach the same state by first absorbing a 3s electron into the nucleus followed by an Auger–Meitner decay, scattering a 4p electron into the 3s orbital and emitting a 4d electron. Both pathways lead to the same final state and each pathway comes with its own phase and amplitude. As these are incommensurate, the different channels interfere, leading to Fano-like line-shapes.

In order to determine the interference between different resonances we calculated 226 states ψ_b that represent the bound-state spectrum within the energy range between 0 to 2838 eV with an experimental resolution of 8  eV well. Note that these are not necessarily eigen-states of the Hamiltonian, but linear combinations of states in an energy interval such that the spectrum is sufficiently reproduced. See appendix D for more information. In figure 4 we show in blue the spectrum one obtains with these 226 resonances each convoluted by a single Lorentzian line-shape for later comparison. We also calculated the self-energy for these states due to the Auger–Meitner decay as given in equation (4). This produces a 226 by 226 matrix. With the use of equation (5), these 226 states ψ_b and the self-energy one can calculate the full spectrum including Auger–Meitner decay and interference between different decay channels. The resulting spectrum is given by the orange line in figure 4.

Zoom In Zoom Out Reset image size

In order to understand the influence of the Fano-effect, i.e. interference between different channels that reach the same final state, we can remove the interference from the full calculation. We then can compare the calculations with and without Fano effect. Removing interference between locally bound states ψ_b that can decay to the same state with an additional Auger–Meitner electron, is done in practice by neglecting the off diagonal elements of the self energy Σ_bb' → δ_bb'Σ_bb . The resulting spectrum is given by the green line in figure 4.

If we compare (figure 4) the calculation with (orange) and without (green) interference between the resonances, we find that the former leads to a more pronounced asymmetry in line-shape. Due to interference, the intensity is reduced left to the resonance and increased right to it. Similar to what was demonstrated in Fano's original paper [28], one observes destructive interference on the low energy tails of resonances and constructive interference on the high energy tails. From figure 4 it also becomes clear that, with or without Fano interference between the resonances, the line-shape is always rather asymmetric. The orange and green line both look very asymmetric compared to the blue line. The most dominant impact on the energy-dependent, asymmetric line-broadening comes from the diagonal elements of the self-energy. The energy dependence of the Auger–Meitner matrix elements, combined with the density of states of the free electron, plays a major role in the determination of the final line-shape.

Due to the finite number of active orbitals used for the calculation of bound states (see [13] for details) one should expect a substantial error in the total electronic binding energy. As long as the errors for ground-state ψ_Ho and final states ψ_b are the same, this influences the spectrum only marginally. Adding orbitals to the active space in the calculation will lower the energies of both ground-state and excited states in almost equal measure. Thereby the excitation energy is changed only marginally.

The inclusion of unbound electrons massively increases the active Hilbert space for the final states. This leads to an energy reduction for all of them and shifts the spectrum to lower energies. These energy-shifts of excited states are of similar magnitude but varying for different resonances. In the calculation energy-shifts are contained in the real part of the self-energy (equation (5)), which renormalises the energies of the bound states. As the final states are now treated more accurately than the ground-state, all resonances are shifted to energies which are too low compared to experiment. In order to correct for this, one needs to add the renormalisation of the ground-state energy due to the unbound orbitals. This leads to a shift in ground-state energy by 15.76  eV and improves the agreement between calculated and experimental energies of the resonances.

The calculation of the Coulomb matrix elements (B.2) which couple bound and unbound electrons requires single particle wave-functions of the Auger–Meitner electrons. As the wave-functions of unbound states are not square integrable we use a countable, orthonormal basis of square integrable wave-packets (WP). We follow the methods developed by Rubtsova, Kukulin and Pomerantsev [31]. These WPs are constructed from the spinors ς_qjlm (r, θ, φ) which are eigen-functions of the free Dirac-Hamiltonian in spherical coordinates. We discretized the continuous energy spectrum by dividing the spectrum of radial momentum q into intervals. We take the values for momenta such that we obtain a set of electrons with equidistant energy spacing $\left\{{{\epsilon}}_{n}=c\sqrt{{q}_{n}^{2}+{m}^{2}{c}^{2}}\enspace \enspace \vert \enspace \enspace n=1\dots \infty \right\}$. We take an energy spacing of 2 eV between successive _n . This yields a reasonable energy resolution for the self energy. The set of free electron basis functions is truncated at an energy of 4 keV. This yields 2000 wave-packets per angular momentum to describe the Auger–Meitner electrons. The cut-off scale at 4 keV is chosen to be larger than the Q-value in order to avoid artificial side-effects due to a truncated basis set. The wave-packets are obtained as

$$\begin{equation}{\varphi }_{{{\epsilon}}_{n}jlm}\left(r,\theta ,\varphi \right)={\int }_{{q}_{n}}^{{q}_{n}+\delta q} {\varsigma }_{qjlm}\left(r,\theta ,\varphi \right)\mathrm{d}q.\end{equation} \tag{ A.1 }$$

These functions span the space of unbound single-particle states, but are not orthogonal to the bound-state wave-functions. The final basis $\left\{{\phi }_{{{\epsilon}}_{n}jlm}\left(\mathbf{r}\right)\right\}$ of unbound states is constructed by orthonormalizing the ${\varphi }_{{{\epsilon}}_{n}jlm}\left(\mathbf{r}\right)$ with respect to the bound-state wavefunctions ${\psi }_{njlm}\left(\mathbf{r}\right)$. The latter are obtained from the density functional theory program FPLO [32–34] as discussed in [13].

This basis allows for a construction of the second quantized operators described in the next section.

With the discretization of the continuous energy spectrum from appendix A and the corresponding WPs as basis functions, we can construct the kinetic energy operator for Auger–Meitner electrons as

$$\begin{equation}{K}_{\mathrm{A}}=\sum _{{{\epsilon}}_{n}}\sum _{l=0}^{5}\sum _{j=l-\frac{1}{2}}^{l+\frac{1}{2}}\sum _{m=-j}^{j}{{\epsilon}}_{n}{c}_{{{\epsilon}}_{n}ljm}^{{\dagger}}{c}_{{{\epsilon}}_{n}ljm}\end{equation} \tag{ B.1 }$$

where ${c}_{{{\epsilon}}_{n}ljm}^{{\dagger}}/{c}_{{{\epsilon}}_{n}ljm}$ creates/annihilates an electron in the unbound state ${\phi }_{{{\epsilon}}_{n}jlm}\left(\mathbf{r}\right)$ with kinetic energy _n . In the following we will denote the four quantum-numbers in the above equation collectively by ${\epsilon}\equiv \left\{{{\epsilon}}_{n},l,j,m\right\}$ and a sum over is understood as four sums over all quantum-numbers.

The knowledge of the Auger–Meitner electron wave-functions ϕ (r) enables us further to calculate the matrix elements for Coulomb scattering of bound electrons into unbound states via

$$\begin{equation}{U}_{{\epsilon}abc}=-\int \frac{\left({\phi }_{{\epsilon}}^{{\ast}}\left({\mathbf{r}}_{1}\right)\cdot {\psi }_{b}\left({\mathbf{r}}_{1}\right)\right)\left({\psi }_{a}^{{\ast}}\left({\mathbf{r}}_{2}\right)\cdot {\psi }_{c}\left({\mathbf{r}}_{2}\right)\right)}{\vert {\mathbf{r}}_{1}-{\mathbf{r}}_{2}\vert }\mathrm{d}{\mathbf{r}}_{1}\enspace \mathrm{d}{\mathbf{r}}_{2}\end{equation} \tag{ B.2 }$$

where Roman indices denote quantum numbers $\left\{n,j,l,m\right\}$ of bound state wavefunctions ψ_njlm (r). The operator that governs all Auger–Meitner scatterings is given by

$$\begin{equation}{U}_{\mathrm{A}}=\sum _{{\epsilon}}\sum _{ijk}{U}_{{\epsilon}ijk}{c}_{{\epsilon}}^{{\dagger}}{c}_{i}^{{\dagger}}{c}_{k}{c}_{j}+\mathrm{h}.\mathrm{c}.\equiv \sum _{{\epsilon}}{U}_{{\epsilon}}.\end{equation} \tag{ B.3 }$$

The full dynamics of the remaining electrons after EC is described by the Hamiltonian

$$\begin{equation}H={H}_{\mathrm{D}\mathrm{y}}+{K}_{\mathrm{A}}+{U}_{\mathrm{A}}\equiv {H}_{0}+\sum _{{\epsilon}}{U}_{{\epsilon}}.\end{equation} \tag{ B.4 }$$

Here, interactions between Auger–Meitner electrons are neglected, since we assume that in most of the relaxation processes there is only a single Auger–Meitner electron present at a time.

With the wave-functions for Auger–Meitner electrons in appendix A and their second quantized operators in appendix B we now derive (4) for the self-energy.

If one includes Auger–Meitner scattering, calculating the EC spectrum (2) involves the inversion of the Hamiltonian (B.4). Here, the interaction U (B.3) creates or annihilates one Auger–Meitner electron with energy . The operator only couples sub-spaces of Fock-space that differ in the number of Auger–Meitner electrons with energy by one. After the EC event, there are only bound electrons and one hole. Final states become more unlikely to influence the spectrum the higher the number of Auger–Meitner electrons they contain. To make use of this property, we divide the full Fock-space $\mathcal{F}$ into subspaces with fixed numbers of Auger–Meitner electrons by introducing the following orthogonal projections. P⁽⁰⁾ projects onto the sub-space of all configurations that do not involve Auger–Meitner electrons. We denote this subspace by ${P}^{\left(0\right)}\mathcal{F}\equiv \left\{{P}^{\left(0\right)}\psi \enspace \vert \enspace \psi \in \mathcal{F}\right\}$. For every Auger–Meitner energy the operator ${P}_{{\epsilon}}^{\left(1\right)}$ projects onto the sub-space of configurations that involve exactly one Auger–Meitner electron with energy . Onto the sub-space with two Auger–Meitner electrons with energies ₁ and ₂ we can project with ${P}_{{{\epsilon}}_{1},{{\epsilon}}_{2}}^{\left(2\right)}$ and so on. By construction the decomposition of $\mathcal{F}$ and the completeness relation read

$$\begin{equation}\mathcal{F}=\bigcup _{n=0}^{N}\left(\bigcup _{{{\epsilon}}_{1},\dots ,{{\epsilon}}_{n}}{P}_{{{\epsilon}}_{1},\dots ,{{\epsilon}}_{n}}^{\left(n\right)}\mathcal{F}\right) \mathbb{1}=\sum _{n=0}^{N}\sum _{{{\epsilon}}_{1},\dots ,{{\epsilon}}_{n}}{P}_{{{\epsilon}}_{1},\dots ,{{\epsilon}}_{n}}^{\left(n\right)}\end{equation} \tag{ C.1 }$$

where N is the total number of electrons in the system. Since Auger–Meitner energies have been discretized and cut-off at 4 keV, the Fock-space is finite dimensional and we can express the Hamiltonian as a block-matrix by using the projection operators. We restrict us on the subspace that includes configurations with at most one Auger–Meitner electron

$$\begin{equation}H{\vert }_{\left({P}^{\left(0\right)}\mathcal{F}\right)\bigcup \enspace \left(\bigcup _{{\epsilon}}\underset{{\epsilon}}{\overset{\left(1\right)}{P}}\mathcal{F}\right)}=\left(\begin{matrix}\hfill {P}^{\left(0\right)}{H}_{0}{P}^{\left(0\right)}\hfill & \hfill {P}^{\left(0\right)}{U}_{{{\epsilon}}_{1}}{P}_{{{\epsilon}}_{1}}^{\left(1\right)}\hfill & \hfill {P}^{\left(0\right)}{U}_{{{\epsilon}}_{2}}{P}_{{{\epsilon}}_{2}}^{\left(1\right)}\hfill & \hfill \cdots \hfill & \hfill \cdots \hfill & \hfill {P}^{\left(0\right)}{U}_{{{\epsilon}}_{n}}{P}_{{{\epsilon}}_{n}}^{\left(1\right)}\hfill \\ \hfill {P}_{{{\epsilon}}_{1}}^{\left(1\right)}{U}_{{{\epsilon}}_{1}}{P}^{\left(0\right)}\hfill & \hfill {P}_{{{\epsilon}}_{1}}^{\left(1\right)}{H}_{0}{P}_{{{\epsilon}}_{1}}^{\left(1\right)}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill {P}_{{{\epsilon}}_{2}}^{\left(1\right)}{U}_{{{\epsilon}}_{2}}{P}^{\left(0\right)}\hfill & \hfill 0\hfill & \hfill {P}_{{{\epsilon}}_{2}}^{\left(1\right)}{H}_{0}{P}_{{{\epsilon}}_{2}}^{\left(1\right)}\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill {\vdots}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \ddots \hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill {\vdots}\hfill \\ \hfill {P}_{{{\epsilon}}_{n}}^{\left(1\right)}{U}_{{{\epsilon}}_{n}}{P}^{\left(0\right)}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill {P}_{{{\epsilon}}_{n}}^{\left(1\right)}{H}_{0}{P}_{{{\epsilon}}_{n}}^{\left(1\right)}\hfill \end{matrix}\right).\end{equation} \tag{ C.2 }$$

From the Lehmann representation of the spectrum (2) we see that the calculation only involves P⁽⁰⁾(z − H)⁻¹ P⁽⁰⁾ which we can directly read off from (C.2)

$$\begin{equation}{P}^{\left(0\right)}{\left(z-H\right)}^{-1}{P}^{\left(0\right)}={\left({P}^{\left(0\right)}z-{P}^{\left(0\right)}{H}_{0}{P}^{\left(0\right)}-{\Sigma}\left(z\right)\right)}^{-1}\end{equation} \tag{ C.3 }$$

where

$$\begin{equation}{\Sigma}\left(z\right)=\sum _{{\epsilon}}\left[{P}^{\left(0\right)}{U}_{{\epsilon}}{P}_{{\epsilon}}^{\left(1\right)}{\left(z-{H}_{0}\right)}^{-1}{P}_{{\epsilon}}^{\left(1\right)}{U}_{{\epsilon}}{P}^{\left(0\right)}\right].\end{equation} \tag{ C.4 }$$

This expression can be further simplified by recognizing that ${P}^{\left(0\right)}{U}_{{\epsilon}}{P}_{{{\epsilon}}^{\prime }}^{\left(1\right)}=0$ if ≠ ', since U annihilates (creates) an Auger–Meitner electron with energy . Thus, using the completeness relation (C.1) one gets ${P}^{\left(0\right)}{U}_{{\epsilon}}\mathbb{1}={P}^{\left(0\right)}{U}_{{\epsilon}}{P}_{{\epsilon}}^{\left(1\right)}$ such that the self-energy is reduced to

$$\begin{equation}{\Sigma}\left(z\right)={P}^{\left(0\right)}\sum _{{\epsilon}}\left[{U}_{{\epsilon}}{\left(z-{H}_{0}\right)}^{-1}{U}_{{\epsilon}}\right]{P}^{\left(0\right)}.\end{equation} \tag{ C.5 }$$

The great advantage of this equation is that it involves only the projection onto the space without Auger–Meitner electrons which is a large reduction of active space. If we use the Krylov-basis of bound multi-configuration states $\left\{{\psi }_{b}\right\}$ introduced in section 2 to express the projection ${P}^{\left(0\right)}={\sum }_{{\psi }_{b}}\vert {\psi }_{b}\rangle \langle {\psi }_{b}\vert$, the self-energy Σ(z) becomes equivalent to (4). Taking the expectation value of (C.3) with respect to T|ψ_Ho⟩ one arrives at the final expression for the EC spectrum (5).

In this derivation two approximations have been made. First we neglected the interaction between Auger–Meitner electrons in the construction of the Hamiltonian (B.4). Second we neglected subspaces with more than one Auger–Meitner electron. The first one is a fundamental simplification, since it determines the form of the Hamiltonian as matrix in (C.2). This form is essential for the whole derivation. Further interactions between unbound electrons would change this form, drastically increase the size of Fock-space and couple multiple subspaces making the calculation no longer feasible. The other approximation is not this essential but yields an important performance boost in the calculation. One could extend it to include states with two or more unbound electrons, which would in analogy yield a self-energy for the states with a single unbound electron. However, this already shows that the effect is expected to be small as it involves two Auger–Meitner-decays and hence is of order $O\left({U}_{\text{A}}^{2}\right)$.

In this section we describe how to calculate the expression for the self-energy Σ_bb'(ω) derived in appendix C numerically. As a first step we describe how to obtain an optimum basis set ψ_b to represent the spectrum for bound states only. The second step describes the calculation of the self energy on this basis.

In order to calculate the set given by the states ψ_b we start from the many-body Ho ground-state wave-function. We then annihilate one of the ns or np_1/2 electrons creating 20 new many-electron states. Starting from these states we generate a Krylov basis using a block Lanczos routine. These states represent the spectral function well, but are not necessarily eigen-states of the Hamiltonian. In total we include 2000 states in our Krylov basis. Most of these states are outside the spectral window of interest. In order to reduce the dimension of the self energy we diagonalize the Hamiltonian on a basis of the Krylov states and only keep the states that have an energy below the Q value. In total we kept 226 states. For more detail on the generation of the Krylov basis and the block Lanczos procedure see [13].

Once we have a basis for the bound states that represent the spectral function well, we need to calculate the Auger–Meitner decay of these states. Firstly we need to determine the coupling between these bound states and states with one electron less and an additional Auger–Meitner electron. Given the energy discretization of the unbound electrons from appendix A, the operator U_A contains about 10⁷ terms. The number of involved Slater determinants for each of the wave-functions ψ_b is of similar magnitude. Hence, it is essential to rewrite the self-energy such that the calculation can be easily parallelized on multiple CPUs and such that the amount of repeated operations is reduced to a minimum. To achieve this, we insert U from (B.3) into (C.5) and make use of ${K}_{\mathrm{A}}{c}_{{\epsilon}}^{{\dagger}}\vert {\psi }_{b}\rangle ={\epsilon}{c}_{{\epsilon}}^{{\dagger}}\vert {\psi }_{b}\rangle$ yielding

$$\begin{align}\hfill {{\Sigma}}_{b{b}^{\prime }}\left(\omega \right)& =\sum _{{\epsilon}}\sum _{ijk}\sum _{{i}^{\prime }{j}^{\prime }{k}^{\prime }}{U}_{ijk{\epsilon}}^{{\ast}}\langle {\psi }_{b}\vert {c}_{i}^{{\dagger}}{c}_{j}^{{\dagger}}{c}_{{\epsilon}}{c}_{k}\enspace {\left(\omega +\text{i}{0}^{+}-{H}_{\mathrm{D}\mathrm{y}}-{K}_{\mathrm{A}}\right)}^{-1}{c}_{{\epsilon}}^{{\dagger}}{c}_{{i}^{\prime }}^{{\dagger}}{c}_{{k}^{\prime }}{c}_{{j}^{\prime }}\vert {\psi }_{{b}^{\prime }}\rangle {U}_{{\epsilon}{i}^{\prime }{j}^{\prime }{k}^{\prime }}\hfill \\ \hfill & =\sum _{{\epsilon}}\sum _{ijk}\sum _{{i}^{\prime }{j}^{\prime }{k}^{\prime }}{U}_{ijk{\epsilon}}^{{\ast}}\langle {\psi }_{b}\vert {c}_{i}^{{\dagger}}{c}_{j}^{{\dagger}}{c}_{k}\enspace {\left(\omega +\text{i}{0}^{+}-{H}_{\mathrm{D}\mathrm{y}}-{\epsilon}\right)}^{-1}{c}_{{i}^{\prime }}^{{\dagger}}{c}_{{k}^{\prime }}{c}_{{j}^{\prime }}\vert {\psi }_{{b}^{\prime }}\rangle {U}_{{\epsilon}{i}^{\prime }{j}^{\prime }{k}^{\prime }}\hfill \\ \hfill & \equiv \sum _{{\epsilon}}\sum _{\tau {\tau }^{\prime }}{U}_{\tau {\epsilon}}^{{\ast}}\langle \tau \vert {R}_{{\mathit{\text{H}}}_{\mathrm{D}\mathrm{y}}}\left(\omega +\text{i}{0}^{+}-{\epsilon}\right)\vert {\tau }^{\prime }\rangle {U}_{{\epsilon}{\tau }^{\prime }} \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \vert \tau =\left(i,j,k\right)\rangle ={c}_{k}^{{\dagger}}{c}_{j}{c}_{i}\vert {\psi }_{b}\rangle .\hfill \end{align} \tag{ D.1 }$$

In the last step we introduced the resolvent of the Dy Hamiltonian

$$\begin{equation}{R}_{{\mathit{\text{H}}}_{\mathrm{D}\mathrm{y}}}\left(z\right)={\left(z-{H}_{\mathrm{D}\mathrm{y}}\right)}^{-1}.\end{equation} \tag{ D.2 }$$

The resolvent matrix elements $\langle \tau \vert {R}_{{H}_{\mathrm{D}\mathrm{y}}}\vert {\tau }^{\prime }\rangle$ do not depend on the created Auger–Meitner electron and its kinetic energy. Therefore, the resolvent only needs to be calculated once—independent of the unbound electron's energy. The energy dependent line-broadening enters via the Slater integrals U_ijk which can be calculated fully in parallel. (D.1) yields an efficient form to evaluate the self-energy. The most resource consuming task is the evaluation of the matrix elements of the resolvent on the sub-space spanned by the vectors |τ⟩. This evaluation can be done by a conventional block Lanczos routine as implemented in Quanty [23, 29, 30].