**4. Method**

In the following, we label the full four-component "perturbed wave function", associated with absorption of one photon with angular frequency Ω and a hole in orbital *b*, by *ρ*Ω,*<sup>b</sup>* , including both radial and spin-angular parts implicitly. As in [41–43], we use here the RPAE-approximation for the many-body response to the absorption of an XUV-photon, albeit within a relativistic framework.

#### *4.1. The Form of the Light–Matter Interaction*

The standard expression for light–matter interaction,

$$h\_I = \text{ec}\mathfrak{a} \cdot \mathbf{A}(\mathbf{r}, t), \tag{42}$$

comes from applying minimal coupling: **p** → **p** + *e***A** to the Dirac Hamilonian in Equation (1). Within the dipole approximation, the vector potential is assumed to be space-independent: **<sup>A</sup>**(**<sup>r</sup>**, **t**) → **<sup>A</sup>**(**t**). This is often referred to as the "velocity gauge" expression for light—matter interaction:

$$h\_I^{\text{velocity}} = \text{eca} \cdot \mathbf{A}(t). \tag{43}$$

A unitary transformation can be made to recast the interaction in the alternative "length gauge" form

$$h\_I^{\text{length}} = \boldsymbol{\epsilon} \mathbf{r} \cdot \mathbf{E}.\tag{44}$$

For details see e.g., ref. [87]. Because our interest here is low-energy photoelectrons, we will stick to the dipole approximation. It is well known that the two gauge forms give identical results when evaluated by an exact wave function, but also for approximations that employ a local potential to describe electron–electron interaction. The non-local exchange potential in the Hartree–Fock approximation can lead to different results in the two gauges when static orbitals are assumed [65,66]. As was shown in the 1970s, the gauge invariance for one-photon processes is restored by the RPAE class of many-body effects [88]. Recently, this was discussed in connection with the calculation of two-photon processes, as needed for the calculation of atomic delays [44], and it was shown that gauge invariance required a full two-photon RPAE treatment. Because ref. [44] also showed that the length gauge results are completely dominated by the time-order where the XUV photon is absorbed first and much less sensitive to final state interactions (after absorption of two photons) than velocity gauge, only the length form will be used here.

With linearly polarized light, we may now write the lowest order approximation of the transition matrix elements from Equation (8) as

$$M^{(1)}(q,\Omega,b) = \langle q \mid ez \mid b \rangle E\_{\Omega \prime} \tag{45}$$

and similarly the two-photon matrix element in Equation (37) as

$$M^{(2)}(q,\omega,\Omega,b) = \lim\_{\xi \to 0^{+}} \oint\_{p} \frac{\langle q \mid ez \mid p \rangle \langle p \mid ez \mid b \rangle}{\varepsilon\_{b} + \hbar \Omega - \varepsilon\_{p} + i\xi} E\_{\omega} E\_{\Omega \prime} \tag{46}$$

where intermediate states, *p*, are to be summed and integrated over for the bound and continuum part of the spectrum respectively. An important difference compared to the onephoton matrix element is that the two-photon matrix element is intrinsically complex for the above threshold ionization, i.e., when *h*¯ Ω exceeds the binding energy, even if correlation effects are neglected.

## *4.2. Diagrammatic Perturbation Theory*

The approximation is illustrated by the diagrams in Figure 1, and a detailed derivation can be found in ref. [44]. The solution of the RPAE equations is done iteratively as indicated in Figure 1 and includes the linear response to the interaction with the XUV photon,

$$\left(\varepsilon\_{b}\pm\hbar\Omega-h\_{\mathbf{x}}^{DF}\right)\mid\rho\_{\Omega,b}^{\pm}\rangle=\sum\_{p}^{\mathrm{exc}}\mid p\rangle\langle p\mid\left(d\_{\Omega\_{j}}+\delta u\_{\Omega}^{\pm}\right)\mid b\rangle,\tag{47}$$

where *δu*<sup>±</sup> Ω is the (linearized) corrections to the Dirac–Fock potential induced by the electromagnetic field (cf. Figure 1c–f,i–l). The Dirac–Fock potential, cf. Equation (4) is defined from its matrix element between orbitals *m*, *n* (occupied or unoccupied):

$$
\langle \langle m \mid u\_{DF} \mid n \rangle = \sum\_{c}^{cor} \langle \{mc\} \mid V\_{12} \mid \{nc\} \rangle,\tag{48}
$$

where curly brackets denote anti-symmetrization. *V*12 denotes here the Coulomb interaction. It is also possible to define a Hartree–Fock type potential for the Breit interaction [89,90], but this aspect of the electron–electron interaction is neglected here. In addition to the Dirac– Fock potential, we usually add a so-called projected potential, *uproj*, to the Hamiltonian in Equation (4). Aiming for a final state with a hole in one of the originally occupied orbitals, the projected potential cancels the removed electron's monopole interaction with

all unoccupied orbitals, without affecting the interaction between the electrons in the ground state. More details can be found in [44]. Through this extra potential, some of the contributions from Figure 1c, precisely those which ensure that the photoelectron feels the correct long-range potential, are accounted for already in lowest order. When converged, the iterative procedure gives the same results if the projected potential is used or not, but the convergence is often much improved in the latter case, especially close to ionization thresholds.

**Figure 1.** Goldstone diagrams illustrating RPAE for the many-body screening of the photon interaction. The set labelled with (**<sup>a</sup>**,**g**) are forward and backward propagation, respectively. Diagrams (**b**,**h**) are the lowest order contributions, while (**<sup>c</sup>**–**f**) and (**i**–**l**) give the many-body response. The sphere indicates the correlated interaction to infinite order. The wavy line indicates the photon interaction and the dashed line the Coulomb interaction. Downward lines (labelled with a, b) stand for holes created when electrons are removed from initially occupied orbitals, and upward lines (labelled with o, p) for initially unoccupied orbitals.

The calculations are performed with a basis set obtained through diagonalization of the radial one-particle Dirac–Fock Hamiltonians in a primitive basis of B-splines [91], defined on a knot sequence in a spherical box. B-splines are piecewise polynomials of a given order *k*. The radial components *f* and *g* of the relativistic wave function, (cf. Equation (5)) are expanded in B-splines of different orders: typically *k* = 7 and *k* = 8 respectively. It has been shown by Froese, Fischer, and Zatsarinny [92] that the use of different B-spline orders is a way to ge<sup>t</sup> rid of the so-called spurious states, which are known to appear in the numerical spectrum after discretization of the Dirac Hamiltonian. Details of the use of B-splines to solve the Dirac equation can be found in ref. [93].

We use further exterior complex scaling (ECS)

$$r \rightarrow \begin{cases} r, & 0 < r < R\_{\text{C}} \\ R\_{\text{C}} + (r - R\_{\text{C}})e^{i\varphi}, & r > R\_{\text{C}} \end{cases} \tag{49}$$

and thus the eigenenergies of the virtual orbitals are complex in general. As a consequence, the energy integration path avoids the pole in Equation (46) and thus the sum and integration over unoccupied states *p* can be represented by a finite sum [41,42].

With converged RPAE results the two-photon matrix elements in Equation (37) can be calculated for the absorption as well as the emission path to sideband *n*:

$$M\_{\mathbf{a}/\mathbf{e}}^{(2)} = \langle q \mid ez \mid \rho^{+}\_{(2n \mp 1)\omega, b} \rangle, \text{ where } \epsilon\_q = \epsilon\_b + 2n\omega,\tag{50}$$

where length gauge has been assumed. The integration in Equation (50) involves two continuum functions and will not converge on any finite interval. The integrand is therefore divided into two parts. The first is an inner region 0 ≤ *r* < *R* < *RC*, where the perturbed wave function and final state can be determined numerically on the B-spline basis. The second is an outer region *R* ≤ *r* < ∞ where the functions can be assumed to be solutions to the pure Coulomb problem, albeit with a possible phase shift. By using different breakpoints *R*, we can check that the result is independent on where the change from numerical to the analytical integration is done. This procedure was described in ref. [42] but has to be slightly changed for the relativistic case, as will be discussed in the next subsection.

#### *4.3. The Continuum–Continuum Transition*

To evaluate Equation (50), we need the final continuum state, *q*, for a photoelectron of energy  *<sup>q</sup>*, obtained in a relativistic framework and with the phase shift it gets from the many-body environment. A good approximation is found as the solution of

$$
\hbar \psi\_{\emptyset} = \epsilon\_{\emptyset} \psi\_{\emptyset} \tag{51}
$$

where *h* = *hDF κ* + *uproj*. By expanding the radial functions *f* and *g* (cf. Equation (5)), in B-splines, *fq*(*r*) = ∑ *ciBi*(*r*), and vice versa for *g*, we can reformulate Equation (51) to a system of linear equations for the coefficients *ci*. Exclusion of the first B-spline yields a regular solution, that is zero at the origin. This determines *ψq* up to a normalization constant. After normalization, which will be discussed below, *ψq* is used for the first part of the integration in Equation (50), i.e., from zero to *R*. We note in passing that in practice it is enough to obtain the large component of the relativistic wave function for a specific energy, because in the region dominated by the Coulomb potential, the Dirac equation gives the small component directly from the large one:

$$\mu\_{\mathcal{K}}(r) = \frac{c\hbar\left(\frac{d}{dr} + \frac{\kappa}{r}\right)\mu\_f(r)}{\epsilon + 2mc^2 + \frac{c^2 Z\_{eff}}{4\pi\epsilon\_0}\frac{1}{r}},\tag{52}$$

where *Zeff* is the effective Coulomb potential felt by the escaping electron.

For the second part of the integration, from *R* to infinity, we need to extract information from the numerical representations of *q* and *ρ* to perform the rest of the integral in Equation (50) analytically as was described in ref. [42]. The final state *q*, is a phase shifted regular solution to the Coulomb problem, which should be correctly normalized, and the perturbed wave function *ρ*, is a phase-shifted outgoing solution with an amplitude determined by the photoionization process.

The outgoing solution *ρ* well outside the ionic core can easily be compared with the pure Coulomb solutions, Equations (23) and (24), combined as in Equation (31), to determine the phase shift, ˜ *δ* in Equation (25). It can easily be checked that the obtained phase shift is independent of *r*, and then Equations (23) and (24) can be used again to construct the solution at any large *r*.

The final state phase-shifted regular solution from Equation (51) can be complimented by its irregular counterpart through Equation (20), evaluated with the relativistic forms of ¯ *f* and *g*¯, and then again the phase shift can be determined from comparision with Equations (23) and (24), combined as in Equation (31), and finally Equation (23) can be used to construct the final state at any *r*.

An additional advantage with the possibility to complement a regular solution with its corresponding irregular solution, and be able to construct the outgoing function, is that it is easy to normalize. The probability flux through the surface of a sphere of radius *R* is

$$\mathcal{J}(\mathbb{R}) = i c \Big( u^f(r)^\* u^g(r) - u^g(r)^\* u^f(r) \Big)\_{r=R} \tag{53}$$

and is constant for any large value of *R*, far outside the core. Because the asymptotic expressions for the large and small components are simple oscillating waves and their relation is *ζ* (cf. Equation (8)), the rate should be <sup>2</sup>*cζ*|*A*|<sup>2</sup> and from that we can determine the amplitude *A*. From the expression for *ζ* in Equation (9), we note the close resemblance with the non-relativistic rate ¯*hk*|*A*|2/*<sup>m</sup>*, just slightly adjusted for relativistic effects.

The last part of the integral, from *R* to *r* → <sup>∞</sup>, in Equation (50) can now be calculated as was described in ref. [42], but now with continuum solutions obtained from Equations (23) and (24).
