**2. Theory**

## *2.1. The Dirac Equation*

The starting point for calculations in a relativistic framework is the Dirac equation. We aim here for calculations on many-electron systems, and as a first approximation we let the electron–electron interaction be approximated by an average potential: the relativistic version of the Hartee–Fock (HF) potential, usually called the Dirac–Fock (DF) potential. Each electron is then governed by the one-particle Hamiltonian:

$$
\hbar\_{DF} = c\mathbf{a} \cdot \mathbf{p} + \left(\mu\_{DF} - \frac{c^2}{4\pi\varepsilon\_0} \frac{Z}{r}\right) \mathbf{1}\_4 + mc^2 \beta\_\prime \tag{1}
$$

with eigenvalues labeled by *E*, and where *α* is expressed in Pauli matrices and *β* has the corresponding form

$$\mathfrak{a} = \begin{pmatrix} 0 & \sigma \\ \sigma & 0 \end{pmatrix}, \quad \mathfrak{F} = \begin{pmatrix} \mathbf{1}\_2 & 0 \\ 0 & -\mathbf{1}\_2 \end{pmatrix}. \tag{2}$$

For closed shell atoms, as the rare gases treated here, the Dirac–Fock potential is spherically symmetric, and the two-component radial part of the wave function can be separated out and determined by the radial Hamiltonian

$$
\varepsilon \left( h\_{\mathbf{x}}^{DF}(r) - mc^2 \right) \left( \begin{array}{c} f\_{\mathbf{x}}(r) \\ g\_{\mathbf{x}}(r) \end{array} \right) = \left( E - mc^2 \right) \left( \begin{array}{c} f\_{\mathbf{x}}(r) \\ g\_{\mathbf{x}}(r) \end{array} \right) = \varepsilon \left( \begin{array}{c} f\_{\mathbf{x}}(r) \\ g\_{\mathbf{x}}(r) \end{array} \right) \tag{3}
$$

with

$$\begin{pmatrix} h\_{\mathbf{x}}^{DF}(r) - mc^2 \end{pmatrix} = \begin{pmatrix} u\_{DF}(r) - \frac{\varepsilon^2}{4\pi\varepsilon\_0} \frac{Z}{r} & -c\hbar \left(\frac{d}{dr} - \frac{\kappa}{r}\right) \\\ c\hbar \left(\frac{d}{dr} + \frac{\kappa}{r}\right) & u\_{DF}(r) - \frac{\varepsilon^2}{4\pi\varepsilon\_0} \frac{Z}{r} - 2mc^2 \end{pmatrix},\tag{4}$$

where *fκ* is the upper, typically larger, component, and *gκ* the lower, typically smaller, component. The four-component eigenfunction to the one-particle Hamiltonian in Equation (1) can now be written as [78]

$$\Psi\_{n\ell j\text{m}}(r,\theta,\phi) = \begin{pmatrix} \frac{f\_{n\ell j}(r)}{r} \chi\_{\text{xm}}(\theta,\phi) \\ \frac{i\xi\_{n\ell j}(r)}{r} \chi\_{-\text{xm}}(\theta,\phi) \end{pmatrix}$$

$$= \begin{pmatrix} \frac{f\_{n\ell}(r)}{r} \sum\_{\nu,\mu} \langle \ell \mu \text{s} \nu \mid j m \rangle \xi\_{\nu} \chi\_{\ell,\mu}(\theta,\phi) \\ \frac{i\xi\_{n\ell j}(r)}{r} \sum\_{\nu,\mu} \langle (2j-\ell) \mu \text{s} \nu \mid j m \rangle \xi\_{\nu} \chi\_{(2j-\ell),\mu}(\theta,\phi) \end{pmatrix},\tag{5}$$

where *χκm*(*<sup>θ</sup>*, *φ*) is a vector coupled function of a spherical harmonic and a spin function *ξν*. The relativistic quantum number *κ* is defined by the eigenvalue equation (*σ* · *-* + <sup>1</sup>)*χκm* = <sup>−</sup>*κχκm* and takes the value *κ* = ( + 1) − *j*(*j* + 1) − 1/4. When *κ* is negative, (*j* = + 1/2), the spherical harmonic associated with the small component, will be one unit of orbital angular momenta larger than that for the large component, and vice verse for positive *κ* (*j* = − 1/2).

The RRPA method, which is also known as the linear response within the timedependent Dirac–Fock (TDDF) formalism, will be used to describe the atomic response to electromagnetic radiation. It accounts for the interaction with the electromagnetic field in lowest order, including also corrections to the static Dirac–Fock potential by field-perturbed orbitals [47,79]. The method is discussed further in Section 4. In the next section, we will discuss expressions for the radial continuum wave functions at large, but not infinite distances from the ion.

#### *2.2. The Scattering Phase of the Photoelectron*

Although the total photoionization cross section is determined by the amplitude of the outgoing electron wave packet, its phase is crucial for its angular dependence as well as its delayed appearance in the continuum. In the following, we discuss the difference of the scattering phase in a relativistic formulation compared to the non-relativistic one.

We consider first an *N*-electron atom that absorbs a photon and subsequently ejects a photoelectron from orbital *b*. The radial photoelectron wave function will in the nonrelativistic case be described by an outgoing phase-shifted Coulomb wave that asymptotically has the form

$$
\mu^{(1)}\_{q,\Omega,b}(r) \approx -\pi M^{(1)}\_{\text{rrel}}(q,\Omega,b) \sqrt{\frac{2m}{\pi k \hbar^2}} \, e^{i\left(kr + \frac{Z}{Ea\_0} \ln 2kr - \ell \frac{\pi}{2} - \sigma\_{Z,k,\ell} + \delta\_{k,\ell}\right)} . \tag{6}$$

Here energy normalization is assumed, and *M*(1) nrel is the non-relativistic electric dipole transition matrix element to the final continuum state *q* with momenta *k*, , and *m*. Although *M*(1) nrel can be chosen to be real in a one-electron context it will be complex when correlation effects are considered. The Coulomb phase is

$$\sigma\_{Z,k,\ell} = \arg\left[\Gamma\left(\ell + 1 + \frac{iZ}{ka\_0}\right)\right],\tag{7}$$

for a photoelectron in the field from a point charge of *Ze*. Note that in Equations (6) and (7), we use the negative Coulomb phase convention, rather than the equivalent positive sign convention that is more commonly used: cf. Equations (1) and (2) in Ref. [44], in order to easily relate the phase expressions to existing relativistic theory in the literature [80]. The additional phase shift *<sup>δ</sup>k*, comes from the short range many-body potential of the final state. The Bohr radius is here denoted with *a*0. In the relativistic case, the asymptotic radial wave function will have an upper and a lower component, cf. Equation (5), which will have the form [80]

$$u\_{q,\Omega,b}^{(f,1)}(r) \approx -\pi M^{(1)}(q,\Omega,b)\sqrt{\frac{2m}{\pi k\hbar^2}\left(1+\frac{\varepsilon}{2mc^2}\right)} \times e^{i\left(kr+\eta\ln 2kr - \gamma\frac{\pi}{2} - \bar{\upsilon}\_{\mathcal{Z},k,\gamma} + \upsilon + \bar{\mathcal{Z}}\_{\mathcal{Z},\gamma}\right)},$$

$$u\_{q,\Omega,b}^{(g,1)}(r) \approx -i\zeta\pi M^{(1)}(q,\Omega,b)\sqrt{\frac{2m}{\pi k\hbar^2}\left(1+\frac{\varepsilon}{2mc^2}\right)} \times e^{i\left(kr+\eta\ln 2kr - \gamma\frac{\pi}{2} - \bar{\upsilon}\_{\mathcal{Z},k,\gamma} + \upsilon + \bar{\mathcal{Z}}\_{\mathcal{Z},\lambda,\gamma}\right)},\tag{8}$$

where the superscripts *f* and *g* indicate the large (upper) and small (lower) components respectively and

$$\zeta = \sqrt{\frac{E - mc^2}{E + mc^2}} = \frac{k\hbar}{2mc} \frac{1}{\left(1 + \frac{\varepsilon}{2mc^2}\right)}\tag{9}$$

is the relation between the large and small component at infinity. This asymptotic relation is given directly by Equation (4), with  = *E* − *mc*<sup>2</sup> being the kinetic energy at infinity. The form of the components in Equation (8) is indeed the same as in the non-relativistic case, but the parameters have slightly changed definition: *M*(1) is now the relativistic matrix element, and *k* is calculated from the relativistic kinetic energy as

$$k = \frac{\sqrt{E^2 - m^2 c^4}}{\hbar c} = \frac{\sqrt{2\varepsilon m}}{\hbar} \sqrt{1 + \frac{\varepsilon}{2mc^2}}.\tag{10}$$

The first factor on the right-hand side of Equation (10) is identical to the non-relativistic expression for *k*, which is thus only slightly adjusted as long as the kinetic energy of the released electron is modest: *ε mc*2. The constant *η* is given by

$$\eta = Z a \mathbb{E} \sqrt{\frac{1}{E^2 - m^2 c^4}} = \frac{Z}{a\_0 k} \left(\frac{\varepsilon}{mc^2} + 1\right) \tag{11}$$

where *α* is the fine structure constant, *α* = *h*¯ /(*<sup>a</sup>*0*mc*). In the non-relativistic limit *η* will thus tend to *<sup>Z</sup>*/(*<sup>a</sup>*0*k*) as expected by comparison with Equation (6). The relativistic Coulomb phase is

$$\mathfrak{v}\_{Z,k,\gamma} = \arg\left[\Gamma(\gamma + i\eta)\right] \tag{12}$$

with

$$
\gamma = \sqrt{\kappa^2 - \mathfrak{a}^2 Z^2} \tag{13}
$$

and

$$\nu = \frac{1}{2} \arg \left[ \frac{-\varkappa + \frac{iZ}{ka\_0}}{\gamma + i\eta} \right]. \tag{14}$$

The phase induced by the short-range part of the many-body potential for the final state is denoted with ˜ *<sup>δ</sup>Z*,*k*,*κ*.

#### *2.3. Phase-Shifted Relativistic Coulomb Functions at Large Distances*

Calculations on many-body systems have to be done numerically. While the wave function for the escaping photoelectron will differ from the analytically known Coulombic ones at short distances, it will approach a combination of a phase-shifted known regular and irregular Coulomb function outside the core of the remaining ion. Because transition matrix elements between continuum states do not converge on a finite grid, it is convenient to have access to continuum solutions, with a possible phase shift *δ*, that can be used to continue the integration to infinity. We are here interested to find expressions for the relativistic case, but it is illustrative to compare with the more studied non-relativistic formulation.

The solutions to the hydrogen-like Schrödinger equation with positive energy is given by the Coulomb functions (see e.g., [81]). The regular Coulomb function is in particular

$$F\_{\ell}(\eta\_{\text{nrel}},kr) = \frac{1}{2}e^{\frac{\pi}{2}\eta\_{\text{nrel}}}\frac{|\Gamma(\ell+1+i\eta\_{\text{nrel}})|}{(2\ell+1)!}e^{-i\text{k}r}(2kr)^{\ell+1}M(\ell+1+i\eta\_{\text{nrel}},2\ell+2,2ikr),\tag{15}$$

where *M* is the confluent hypergeometric function, *σ* is defined in Equation (7) and *η*nrel = *<sup>Z</sup>*/(*<sup>a</sup>*0*k*). Non-relativistic Coulomb functions expressions, valid for large *kr*, are provided in Ref. [82]:

$$F\_{\ell} = \oint \cos \Delta\_{\text{nrel}} + \tilde{f} \sin \Delta\_{\text{nrel}} \tag{16}$$

$$G\_{\ell} = f \cos \Delta\_{\text{nrel}} - \xi \sin \Delta\_{\text{nrel}} \tag{17}$$

for the regular, *F*, and irregular, *G*, Coulomb functions respectively, where

$$
\Delta\_{\text{nrel}} \equiv kr + \frac{Z}{ka\_0} \ln 2kr - \frac{\pi}{2} \ell - \sigma\_{Z,k,\ell} + \delta \tag{18}
$$

and ¯ *f* and *g*¯, which depend on *Z*, *r*, *k* , and , can be obtained through simple recursive formulas given in ref. [82]. When *r* → <sup>∞</sup>, *g*¯ → 0 and ¯ *f* → 1 and thus the regular function approaches a sin-function, and the irregular a cos-function, both with amplitude one. The combination

$$F\_{\ell}(\eta\_{\text{nrel}}, kr) - iG\_{\ell}(\eta\_{\text{nrel}}, kr) \tag{19}$$

will thus asymptotically approach an outgoing wave, with modulus square equal to unity. Energy normalized continuum functions are obtained by multiplications with 2*m*/*πkh*¯ 2 .

It is interesting to note that Equations (16) and (17) imply that the irregular (regular) function can readily be obtained when the regular (irregular) one is at hand. In the former case, the irregular solution is found as

$$\mathbf{G}\_{\ell} = \frac{\left(\frac{dF\_{\ell}}{dr} - \frac{F\_{\ell}}{\xi^2 + f^2} \left(\frac{d\mathcal{g}}{dr}\mathcal{g} + \frac{df}{dr}\mathcal{f}\right)\right)}{k + \eta/r + \frac{1}{\xi^2 + f^2} \left(\frac{d\mathcal{g}}{dr}\mathcal{f} - \frac{df}{dr}\mathcal{g}\right)}. \tag{20}$$

Turning to the relativistic Coulomb problem, we set out to find the relativistic counterparts to Equations (16) and (17), which to the best of our knowledge, are not available in the

literature. The exact two-component relativistic regular, *F* ˜ *γ*, and irregular, *G* ˜ *γ*, solutions are given in a pioneering article by Johnson and Cheng [80]. In particular the regular solution is

$$F\_{\gamma}(\eta, kr) = \sqrt{\frac{E + mc^2}{2E}} \frac{1}{2} e^{\frac{\pi}{2}\eta} \frac{|\Gamma(\gamma + i\eta)|}{\Gamma(2\gamma + 1)} (-2ikr)^{\gamma} e^{ikr}$$

$$\begin{pmatrix} \left(-\kappa + \frac{iZ}{ka\_0}\right) M\_{\gamma} + (\gamma - i\eta) M\_{\gamma + 1} \\\ -i\zeta \left(\left(-\kappa + \frac{iZ}{ka\_0}\right) M\_{\gamma} - (\gamma - i\eta) M\_{\gamma + 1}\right) \end{pmatrix} \tag{21}$$

with *ζ*, *γ*, *η* and *k* given in Equations (9)–(11) and (13), and the short-hand notation

$$\begin{array}{ll} M\_{\gamma} = & M(\gamma - i\eta, 2\gamma + 1, -2iz) \\ M\_{\gamma + 1} = & M(\gamma + 1 - i\eta, 2\gamma + 1, -2iz) \end{array} \tag{22}$$

has been used for the confluent hypergeometric functions.

An asymptotic expansion of the confluent hypergeometric function, *M* can be found in ref. [83], which indeed can be used to obtain asymptotic expansions for *F* ˜ *γ* and *G* ˜ *γ* on forms similar to Equations (16) and (17):

$$\bar{F}\_{\gamma} = \sqrt{\frac{E + mc^2}{2E}} \begin{pmatrix} \bar{f}\_{\text{large}} \cos \Delta - \bar{\mathcal{g}}\_{\text{large}} \sin \Delta\\ -\zeta \left( \bar{\mathcal{g}}\_{\text{small}} \cos \Delta + \bar{f}\_{\text{small}} \sin \Delta \right) \end{pmatrix} \tag{23}$$

and

$$\mathbf{G}\_{\gamma} = \sqrt{\frac{\mathbf{E} + m\mathbf{c}^2}{2E}} \begin{pmatrix} -\left(\ddot{\mathbf{g}}\_{\text{large}}\cos\Delta + \ddot{f}\_{\text{large}}\sin\Delta\\ -\ddot{\zeta}\left(\ddot{f}\_{\text{small}}\cos\Delta - \ddot{g}\_{\text{small}}\sin\Delta\right) \end{pmatrix} \end{pmatrix} \tag{24}$$

with

$$
\Delta = kr + \eta \ln 2kr - \pi \gamma / 2 - \vartheta\_{\mathbb{Z}, k, \gamma} + \nu + \delta \tag{25}
$$

with *σ*˜ and *ν* given in Equations (12) and (14). The possible extra phase shift is denoted by ˜ *δ*. In the non-relativistic limit Δ → Δnrel ± *π*/2, for *κ* > 0 and *κ* < 0 respectively, and thus the sin/cos—functions in Equations (16) and (17) are replaced with ∓ cos / ± sin in the upper components of Equation (23) and (24). The relativistic ¯ *f* , *g*¯ functions are obtained as

$$f\_{\text{large}/\text{small}} = \text{Re}(\aleph\_{\pm})\tag{26}$$

$$
\delta\_{\text{large}/\text{small}} = \text{Im}(\aleph\_{\pm}) \tag{27}
$$

from

$$\begin{split} \aleph\_{\pm} &= \sum\_{n=0} \frac{(\gamma - i\eta)\_n (-\gamma - i\eta)\_n}{n!} (2ikr)^{-n} \pm \\ \sum\_{n=0} \frac{(\gamma + 1 + i\eta)\_n (-\gamma + 1 + i\eta)\_n}{n!} (-2ikr)^{-n} \end{split} \tag{28}$$

where (*a*)*n* = *a*(*a* + <sup>1</sup>)(*a* + <sup>2</sup>)...(*<sup>a</sup>* + *n* − <sup>1</sup>), (*a*)0 = 1. Similarly to the non-relativistic case *g* ¯ large/small → 0 and ¯ *f*large/small → 1 when *r* → ∞. Thus the upper regular, and the lower irregular, approach cos Δ, whereas the lower regular and the upper irregular tend to sin Δ. The asymptotic expressions are thus

$$F\_{\gamma}(\eta, kr) \rightarrow \sqrt{\frac{E + mc^2}{2E}} \begin{pmatrix} \cos \Delta \\ -\zeta \sin \Delta \end{pmatrix} \tag{29}$$

$$\mathbf{G}\_{\gamma}(\eta, kr) \to \sqrt{\frac{E + mc^2}{2E}} \begin{pmatrix} -\sin \Delta \\ -\zeta \cos \Delta \end{pmatrix},\tag{30}$$

when *kr* → <sup>∞</sup>, and the combination

$$\mathcal{F}\_{\gamma}(\eta, kr) - i\mathcal{G}\_{\gamma}(\eta, kr) \to \sqrt{\frac{E + mc^2}{2E}} \begin{pmatrix} 1 \\ i\zeta \end{pmatrix} e^{i\Lambda} \tag{31}$$

will, in close analogy with the non-relativistic expression in Equation (19), asymptotically approach an outgoing wave, with modulus square unity. The energy normalized functions are again obtained by multiplication with <sup>2</sup>*m*/*πkh*¯ 2. We note finally that Equation (20) holds also in a relativistic framework. It provides the irregular solution, *G* ˜ *γ* from *F* ˜ *γ*, if ¯ *f* and *g*¯ are just replaced with ¯ *f*large and *g*¯large or ¯ *f*small and *g*¯small for the upper and lower components respectively.

## **3. Delay in Photoionization**

We will here briefly discuss the calculation of delays in laser-assisted photoionization, emphasizing the differences compared to the non-relativistic description. A detailed account of the latter can be found in refs. [42,44].

#### *3.1. The Wigner Delay*

The concept of delay was introduced by Wigner [35], Smith [36] and Eisenbud [37] as the derivative of the scattering phase with respect to energy. With a finite difference approximation of the derivative Δ*ω* = 2*ω*, the Wigner contribution to the atomic delay measured in a RABITT experiemnt is

$$
\pi\_W = \frac{\phi\_{>} - \phi\_{<}}{2\omega},
\tag{32}
$$

where *φ*>/< refer to the phases acquired in the XUV absorption step in the two paths where either the higher or the lower harmonic is absorbed. Non-relativistically, and for detection of the photoelectron in the **zˆ** direction, these phases are

$$\boldsymbol{\phi}\_{>}^{\text{nrel}} = \arg \left( \sum\_{\ell} M\_{>}^{\text{nrel}}(\ell) e^{i \left( -\ell \frac{\pi}{2} - \sigma\_{\mathbb{Z}, k > \ell} + \delta\_{k > \ell} \right)} Y\_{\ell, 0}(\mathbf{z}) \right)$$

$$\boldsymbol{\phi}\_{<}^{\text{nrel}} = \arg \left( \sum\_{\ell} M\_{<}^{\text{nrel}}(\ell) e^{i \left( -\ell \frac{\pi}{2} - \sigma\_{\mathbb{Z}, k < \ell} + \delta\_{k < \ell} \right)} Y\_{\ell, 0}(\mathbf{z}) \right) \tag{33}$$

where the short-hand notation for the one-photon matrix elements, *<sup>M</sup>*>/<() ≡ *<sup>M</sup>*(1)(*<sup>q</sup>*>/<, Ω>/<, *b*), with final photoelectron wave number *k*>/< and angular momentum , after absorption of a photon with angular frequency Ω>/<, is used. Relativistically the corresponding amplitudes have two components and it is more appropriate to define the Wigner delay as

$$\tau\_{W} = \frac{1}{2\omega} \arg\left(\sum\_{m=\pm\frac{1}{2}}$$

$$\left(\sum\_{\mathbf{x}} M\_{<} \left(\begin{array}{c} \chi\_{\text{x}m}(\mathbf{\hat{z}})\\ i\zeta\_{\text{X}-\text{x}m}(\mathbf{\hat{z}}) \end{array}\right) \mathbf{e}^{i\left(-\gamma\frac{\pi}{2} - \hat{\sigma}\_{\text{Z},k,\gamma} + \upsilon + \delta\_{\text{Z},k,\gamma}\right)}\right)^{\dagger}$$

$$\left(\sum\_{\mathbf{x}'} M\_{>} \left(\begin{array}{c} \chi\_{\text{x}'m}(\mathbf{\hat{z}})\\ i\zeta\_{\text{X}-\text{x}'cm}(\mathbf{\hat{z}}) \end{array}\right) \mathbf{e}^{i\left(-\gamma'\frac{\pi}{2} - \hat{\sigma}\_{\text{Z},k,\gamma'} + \upsilon' + \delta\_{\text{Z},k,\gamma'}\right)}\right)\right),\tag{34}$$

where the calculation of the delay of electrons emitted along the z-axis requires an incoherent sum over *m* = ±1/2. The two incoherent contributions to the Wigner delay are due to unresolved photoelectron spin in the final state.

## *3.2. The Atomic Delay*

We now consider measurements that employ the RABBIT technique [2], where an XUV comb of odd-order harmonics of a fundamental laser field with angular frequency *ω*, is combined with a synchronized, weak laser field with the same angular frequency. In RABBIT, the one-photon ionization process is assisted by an IR photon that is either absorbed or emitted. The same final state is reached when both an XUV harmonic with energy *h*¯ Ω< = (<sup>2</sup>*n* − <sup>1</sup>)*h*¯ *ω* and an IR photon is absorbed, as when the next XUV harmonic, *h*¯ Ω> = (<sup>2</sup>*n* + <sup>1</sup>)*h*¯ *ω*, is absorbed while an IR photon is emitted. This gives rise to modulated sidebands in the photoelectron spectrum at energies corresponding to the absorption of an even number of IR photons. Schematically the intensity of such a sideband can be written as [25]

$$S = |\left|A\_{\mathfrak{a}} + A\_{\mathfrak{e}}\right|^2 = |\left|A\_{\mathfrak{a}}\right|^2 + |\left|A\_{\mathfrak{e}}\right|^2 + A\_{\mathfrak{e}}^\* A\_{\mathfrak{a}} + A\_{\mathfrak{e}} A\_{\mathfrak{a}}^\*$$

$$\left|\left|A\_{\mathfrak{a}}\right|^2 + \left|\left|A\_{\mathfrak{e}}\right|\right|^2 + 2\left|A\_{\mathfrak{e}}\right|\left|\left|A\_{\mathfrak{a}}\right|\cos\left[\arg\left(A\_{\mathfrak{e}}\right) - \arg\left(A\_{\mathfrak{a}}\right)\right]\right|\tag{35}$$

where *A*a/e are the complex quantum amplitudes for the two-photon processes involving absorption (a) or emission (e) of an IR photon, and leading to the same final energy. The modulation arises from the last term in Equation (35) and can be shown to be governed by the delay between the IR and XUV pulses, *τ*, the group delay of the attosecond pulses in the train, *τXUV*, and by a contribution from the atomic system which is due the phase difference between the emission and the absorption paths in the atom:

$$\cos[\arg(A\_{\mathfrak{e}}) - \arg(A\_{\mathfrak{a}})] = \cos[2\omega(\tau - \tau\_{XUV}) + \phi\_{\mathfrak{e}} - \phi\_{\mathfrak{a}}].\tag{36}$$

The atomic contribution can be interpreted as an atomic delay: *τA* = (*φ*e − *φ*a)/2*<sup>ω</sup>*. Because the delay between the two light fields is controlled in the experiments and the pulse train group delay can be canceled through relative measurements, the atomic contribution can be extracted. A recent review of the experimental method can be found in [84]. In the following, we discuss the determination of *φ*a and *φ*e.

The outgoing radial wave function for the large component, after interaction with two photons, will, in accordance with the one-photon situation in Equation (8), have the asymptotic form

$$u^{(f2)}\_{\
q,\omega,\Omega,b}(r) \approx -\pi M^{(2)}(q,\omega,\Omega,b)\sqrt{\frac{2m}{\pi k\hbar^2}\left(1+\frac{\varepsilon}{2mc^2}\right)}e^{i\left(kr+\eta\ln 2kr-\gamma\frac{\omega}{2}-\vartheta\_{Z,k\gamma}+\upsilon+\delta\_{Z,k\gamma}\right)},\tag{37}$$

where the important difference compared to the one-photon case lies in the presence of the two-photon transition element *<sup>M</sup>*(2), which connects the initial state *b* to the continuum state *q* through all dipole-allowed intermediate states. The small component follows as in Equation (8). The phases acquired in the absorption and emission paths (cf. Equation (36)), are given by the corresponding two-photon matrix element and the phase of the photoelectron. In the non-relativistic case, and for photoelectrons with momentum along the common polarization axis of the fields, **zˆ**, they are given as

$$\begin{split} \boldsymbol{\phi}\_{\text{a}}^{\text{nrel}} &= \arg \left( \sum\_{\ell} M\_{\text{a,nrel}}(\ell) e^{i \left( -\ell \frac{\pi}{2} - \sigma\_{\text{Z},k,\ell} + \delta\_{\text{k},\ell} \right)} Y\_{\ell,0}(\mathbf{z}) \right) \\ \boldsymbol{\phi}\_{\text{c}}^{\text{nrel}} &= \arg \left( \sum\_{\ell} M\_{\text{e,nrel}}(\ell) e^{i \left( -\ell \frac{\pi}{2} - \sigma\_{\text{Z},k,\ell} + \delta\_{\text{k},\ell} \right)} Y\_{\ell,0}(\mathbf{z}) \right), \end{split} \tag{38}$$

where the short-hand notation

$$M\_{\rm a,nrel}(\ell) = M\_{\rm nrel}^{(2)}(q, \omega, \Omega\_{<\prime}, b),$$

$$M\_{\rm e,nrel}(\ell) = M\_{\rm nrel}^{(2)}(q, -\omega, \Omega\_{>\prime}, b) \tag{39}$$

has been used and the subscripts a and e stand for IR absorption and emission, respectively. For photoelectron emission along the **zˆ**-direction, i.e., *θ* = 0, the spherical harmonic is non-zero only for azimutal quantum number *m* = 0. The atomic delay, defined as the phase difference divided by 2 *ω*, can subsequently be calculated as

$$
\tau\_A = \frac{\phi\_\text{v} - \phi\_\text{a}}{2\omega}.\tag{40}
$$

In the Dirac case, there are two distinct differences. First, a sum over *m* = ±1/2 is required, because both spin-directions contribute to the emission along the **zˆ**-direction. Second, due to the multi-component wave function the Dirac case the expression gets more involved:

$$\tau\_{A} = \frac{1}{2\omega} \arg\left(\sum\_{m=\pm\frac{1}{2}}$$

$$\left(\sum\_{\mathbf{x}} M\_{\mathbf{3}}\left(\begin{array}{c} \chi\_{\mathbf{x}m}(\mathbf{\hat{z}})\\ i\!\!\!\!\!\chi\text{-}\_{\mathbf{x}\mathbf{m}}(\mathbf{\hat{z}}) \end{array}\right) \mathbf{e}^{i\left(-\gamma\frac{\mathbf{x}}{2} - \upvarphi\_{\mathbf{Z},k,\gamma} + \upsilon + \delta\_{\mathbf{Z},k,\gamma}\right)}\right)^{\dagger}$$

$$\times \left(\sum\_{\mathbf{x'}} M\_{\mathbf{e}}\left(\begin{array}{c} \chi\_{\mathbf{x'}m}(\mathbf{\hat{z}})\\ i\!\!\!\!\!\chi\text{-}\_{\mathbf{x'}}m(\mathbf{\hat{z}}) \end{array}\right) \mathbf{e}^{i\left(-\gamma'\frac{\mathbf{x}}{2} - \upvarphi\_{\mathbf{Z},k,\gamma'} + \upsilon + \delta\_{\mathbf{Z},k,\gamma'}\right)}\right)\right),\tag{41}$$

where the sum over *m* = ±1/2 is done incoherently (see e.g., the discussion in ref. [85]) , whereas the sum over *κ* is done coherently. The two incoherent contributions to the atomic delay are due to unresolved photoelectron spin in the final state. The expression for the Wigner and atomic delay for electrons detected along an arbitrary direction have been discussed in ref. [86]
