Many-body calculations of two-photon, two color matrix elements for attosecond delays

TL;DR

This paper presents advanced many-body calculations of atomic delays in noble gas photoionization, demonstrating gauge invariance and examining the validity of common delay interpretations, with findings relevant for attosecond physics.

Contribution

It introduces a comprehensive two-color two-photon Random-Phase Approximation with Exchange approach for calculating atomic delays, highlighting gauge invariance and convergence properties.

Findings

Gauge invariance of atomic delays is demonstrated.

Universal continuum--continuum contribution is generally valid.

Special cases like the 3s Cooper minimum show deviations.

Abstract

We present calculations for attosecond atomic delays in photoionization of noble gas atoms based on full two-color two-photon Random-Phase Approximation with Exchange in both length and velocity gauge. Gauge invariant atomic delays are demonstrated for the complete set of diagrams. The results are used to investigate the validity of the common assumption that the measured atomic delays can be interpreted as a one-photon Wigner delay and a universal continuum--continuum contribution that depends only on the kinetic energy of the photoelectron, the laser frequency and the charge of the remaining ion, but not on the specific atom or the orbital from which the electron is ionized. Here we find that although effects beyond the universal IR--photoelectron continuum--continuum transitions are rare, they do occur in special cases such as around the $3s$ Cooper minimum in argon. We conclude also that in general the convergence in terms of many-body diagrams is considerably faster in length gauge than in velocity gauge.

Equations

$$u_{q,\Omega ,a}^{(1)}\left(r\right)\approx -\pi M^{(1)}\left(q,\Omega ,a\right)\sqrt{\frac{2m}{\pi k{\hslash }^2}}$$

$$\times e^{i\left(kr+\frac{Z}{k{a}_0}\ln 2kr-\ell \frac{\pi }{2}+{\sigma }_{Z,k,\ell }+{\delta }_{k,\ell }\right)},$$

$${\sigma }_{Z,k,\ell }=\arg \left[\Gamma \left(\ell +1-i\frac{Z}{k{a}_0}\right)\right].$$

$$u_{q,\omega ,\Omega ,a}^{(2)}\left(r\right)\approx -\pi M^{(2)}(q,\omega ,\Omega ,a)\sqrt{\frac{2m}{\pi k{\hslash }^2}}$$

$$e^{i\left(kr+\frac{Z}{k{a}_0}\ln kr-\ell \frac{\pi }{2}+{\sigma }_{\ell ,Z}\left(k\right)+{\delta }_{\ell }\left(k\right)\right)},$$

$$h=\frac{1}{2m}{\mathbf{p}}^2+V+\frac{e}{m}\mathbf{p}\cdot \mathbf{A}+\frac{e^2}{2m}\cdot {\mathbf{A}}^2$$

$$U=\exp \left[\frac{i}{\hslash }\left(e^2{\int }^{t}d{t}^{\prime }{A}^2(t^{\prime })\right)\right],$$

$$Uh{U}^{†}U\Psi =i\hslash U\frac{\partial }{\partial t}{U}^{†}U\Psi .$$

$$h_I^{\mathrm{velocity}}=\frac{e}{m}\mathbf{p}\cdot \mathbf{A},$$

$$U=\exp \left[\frac{i}{\hslash }\left(e\mathbf{r}\cdot \mathbf{A}(t)+e^2{\int }^{t}d{t}^{\prime }{A}^2(t^{\prime })\right)\right]$$

$$h_I^{\mathrm{length}}=e\mathbf{r}\cdot \mathbf{E},$$

$$M^{(1)}\left(q,\Omega ,a\right)=⟨q\mid ez\mid a⟩E_{\Omega },$$

$$M^{(1)}\left(q,\Omega ,a\right)=⟨q\mid \frac{e}{m}{p}_z\mid a⟩A_{\Omega }.$$

$$M^{(2)}(q,\omega ,\Omega ,a)=$$

$$\underset{\xi \to 0^{+}}{\lim }\sum _p\int \frac{⟨q\mid ez\mid p⟩⟨p\mid ez\mid a⟩}{{ϵ}_a+\hslash \Omega -{ϵ}_p+i\xi }{E}_{\omega }{E}_{\Omega }.$$

$${\varphi }_{\mathrm{a}}$$

$${\varphi }_{\mathrm{e}}$$

$$M_{\mathrm{a}}(\ell )$$

$$M_{\mathrm{e}}(\ell )$$

$${\eta }_{Z,k,\ell }$$

$${\tau }_A=\frac{{\varphi }_{\mathrm{e}}-{\varphi }_{\mathrm{a}}}{2\omega }.$$

$${\varphi }_{>}$$

$${\varphi }_{<}$$

$${\tau }_W=\frac{{\varphi }_{>}-{\varphi }_{<}}{2\omega }.$$

$$h_0^{\ell }\left(r\right)=-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial r^2}+\frac{{\hslash }^2}{2m}\frac{\ell \left(\ell +1\right)}{r^2}$$

$$-\frac{e^2}{4\pi {ϵ}_0}\frac{Z}{r}+u_{\mathrm{HF}}+u_{\mathrm{proj}}.$$

$$h_I=$$

$$=$$

$$\mid {\psi }_a\left(t\right)⟩=\mid a⟩e^{-i{ϵ}_{a}t/\hslash }+\sum _j\sum _{\sigma =\pm 1}\mid {\rho }_{{\Omega }_j,a}^{(1:\sigma )}⟩e^{-i\left(\sigma {\Omega }_j+{ϵ}_a/\hslash \right)t}$$

$$+\sum _{j,j^{\prime }}\sum _{\sigma =\pm 1}\sum _{{\sigma }^{\prime }=\pm 1}\mid {\rho }_{{\Omega }_{j^{\prime }},{\Omega }_j,a}^{(2:{\sigma }^{\prime },\sigma )}⟩e^{-i\left({\sigma }^{\prime }{\Omega }_{j^{\prime }}+\sigma {\Omega }_j+{ϵ}_a/\hslash \right)t}+\dots ,$$

$$\left(i\hslash \frac{\partial }{\partial t}-h_0^{\ell }\right)\mid {\psi }_a\left(t\right)⟩=h_I\mid {\psi }_a\left(t\right)⟩,$$

$$\left(\sigma \hslash {\Omega }_j+{ϵ}_a-h_0^{\ell }\right)\mid {\rho }_{{\Omega }_j,a}^{(1:\sigma )}⟩=d_{{\Omega }_j}\mid a⟩.$$

$$\mid {\rho }_{0,{\Omega }_j,a}^{(1:\sigma )}⟩=\sum _p\frac{\mid p⟩⟨p\mid d_{{\Omega }_j}\mid a⟩}{\sigma \hslash {\Omega }_j+{ϵ}_a-{ϵ}_p},$$

$$\mid \Psi \left(t\right)⟩=\mid \left\{{\psi }_a\left(t\right),{\psi }_b\left(t\right),{\psi }_c\left(t\right),{\psi }_d\left(t\right),\dots \right\}⟩$$

$$\displaystyle r\rightarrow\Bigg{\{}\begin{array}[]{lr}r,&0<r<R_{C}\\ R_{C}+\left(r-R_{C}\right)e^{i\varphi},&r>R_{C},\end{array}$$

Full text

I Introduction

Techniques for probing ultrafast electronic dynamics, such as the Reconstruction of Attosecond Beating By Interference of Two-photon Transitions (RABBIT) Paul et al. (2001) or the attosecond streak camera Itatani et al. (2002), use delay-dependent modulations in photoelectron spectra to quantify the time it takes for an electron to escape an atomic potential  Schultze et al. (2010); Klünder et al. (2011); Guénot et al. (2012, 2014); Sabbar et al. (2015); Kotur et al. (2016); Gruson et al. (2016); Isinger et al. (2017a); Cirelli et al. (2018). These modulations arise since the interaction with the ionizing attosecond pulse (or pulse train) takes place in the presence of a laser field that is phase-locked to the attosecond light field(s). It has further been established Dahlström et al. (2012a); Dahlström et al. (2013); Pazourek et al. (2013) that it is meaningful to separate the measured atomic delay, $\tau_{A}\approx\tau_{W}+\tau_{cc}$, into a Wigner-like delay associated with the one-photon extreme ultraviolet (XUV) ionization process, $\tau_{W}$, and a contribution from the interaction with the infrared (IR) laser-field in the presence of the atomic potential, called the continuum–continuum delay, $\tau_{cc}$ (or Coulomb-Laser Coupling delay in the context of streaking). In this context $\tau_{cc}$ denotes the contribution from a single photoelectron in a Coulomb field that absorbs or emits an IR photon, as detailed in  Dahlström et al. (2012a); Dahlström et al. (2013). While the Wigner delay is known to be strongly dependent on the atomic origin of the electron Kheifets (2013), the contribution from the IR photon has been found to be more “universal” in the sense that it depends only on the kinetic energy of the photoelectron, the photon energy of the laser field and the charge of the remaining ion Dahlström et al. (2013). In the limit of weak fields, the physics can be described as interference effects between various two-photon processes. The validity of the $\tau_{cc}$ correction has been studied through many-body calculations of the two-photon process both for angular integrated measurements and for detection along the polarization axis Dahlström et al. (2012b); Dahlström and Lindroth (2014). Although there are exceptions, in particular close to ionization thresholds and at resonances Dahlström and Lindroth (2014); Sabbar et al. (2015); Kotur et al. (2016), the universality of the contribution from the laser photon has hitherto proved to be a good approximation, with the important practical consequence that Wigner delays, $\tau_{W}$, can be extracted from measured atomic delays, $\tau_{A}$. The many-body calculations themselves have been benchmarked, for example against the difference between $2s$ and $2p$ time delays in neon Isinger et al. (2017a), over a wide energy range, and against time delay differences between the outermost shells of rare gas atoms Guénot et al. (2014). Still, the procedure used so far has employed significant approximations. First, only the dominating time-order with the XUV photon being absorbed first, and the IR photon being exchanged subsequently in a continuum-continuum transition, has been studied in much detail Lindroth and Dahlström (2017). Second, a more careful account for many-body effects has only been done for the XUV photoionization process, while the interaction with the second photon has been calculated in the lowest-order/classical approximation Dahlström et al. (2012b); Feist et al. (2014). Although these approximations are reasonable they have important consequences: the results at this level of theory are expected to depend on whether the light–matter interaction is expressed in the length or velocity gauge. In addition, experimental results on the difference between $3s$ and $3p$ time delays in argon Klünder et al. (2011); Guénot et al. (2012) show a marked disagreement with theory in the region around the $3s$–Cooper minimum (at photon energies of $\sim 40$ eV). Therefore, it is important to push the study of atomic delays one step further. Here we have performed Random-Phase Approximation with Exchange (RPAE) type calculations for the complete two-photon process. The RPAE Approximation, which is identical to Time-dependent Hartree-Fock (TDHF), is known to account for the dominating many-body effects in one-photon ionization Amusia (1990). While the length and velocity form of the electric dipole interaction gives the same result for electrons in any local potential, the use of the Hartree-Fock exchange potential destroys this invariance. As was shown more than forty years ago Lin (1977), RPAE, which accounts fully for hole-particle excitations (including the effects usually called ground-state correlation, see below) is able to restore the gauge invariance. For two-photon processes L’Huillier et al. (1986); Pan et al. (1990) and beyond Sekino and Bartlett (1992) pioneering studies of many-body effects on the RPAE-level were done already in the nineteen eighties and nineties. The target at the time was absorption of equal energy photons for below-threshold ionization (one photon alone could not induce ionization). In contrast, our interest is in the interaction with two photons of very different energies, where one photon can initiate an above-threshold ionization process. We will demonstrate that, just as for one-photon ionization, gauge invariance is obtained when hole-particle excitations are fully accounted for including all time orders, i.e. in a complete two-photon RPAE calculation. We further show that the size of individual contributions is vastly different in the two gauges and that the common approximation to neglect the time-order where the IR photon is absorbed first leads to wrong results in velocity gauge.

In Sec. II we revisit the theory for atomic delays. The method of calculation is outlined in Sec. III and the results are presented in Sec. IV. In Sec. V we present a discussion of our findings and in Sec. VI we present our conclusions.

II Theory

Here we will briefly discuss the calculation of delays in laser-assisted photoionization. A more detailed account can be found in Ref. Dahlström and Lindroth (2014). We consider first an $N$-electron atom that absorbs one photon and subsequently ejects a photoelectron. The radial photoelectron wave function will asymptotically be described by an outgoing phase-shifted Coulomb wave

[TABLE]

where $M^{(1)}$ is the electric dipole transition matrix element to the final continuum state $q$ with momenta $k$, $\ell$, and $m_{a}$. When correlation effects are accounted for $M^{(1)}$ can contribute to the phase shift, and $\sigma_{Z,k,\ell}$ is the Coulomb phase for a photoelectron with wave number $k$ and angular momentum quantum number $\ell$ in the field from a charge of $Ze$:

[TABLE]

The phase $\delta_{k,\ell}$ in Eq. (LABEL:onephoton2) denotes the additional shift induced by the atomic potential at short range. In the following we label the full perturbed wave function associated with absorption of one photon with angular frequency $\Omega$ and a hole in orbital $a$, by $\left|\rho_{\Omega,a}\right>$, including both radial, angular and spin parts implicitly.

We will consider measurements that employ the RABBIT technique Paul et al. (2001), where an XUV comb of odd-order harmonics of a fundamental laser field, $\Omega=(2n+1)\omega$, is combined with a synchronized, weak laser field with angular frequency $\omega$. In RABBIT, the one-photon ionization process is assisted by an IR photon that is either absorbed or emitted. This gives rise to quantum beating of sidebands in the photoelectron spectrum at energies corresponding to the absorption of an even number of IR photons. The outgoing radial wave function after interaction with two photons (one XUV photon, $\Omega$, and one laser photon, $\omega$), will asymptotically have the form

[TABLE]

where the important difference compared to the one-photon case lies in the presence of the two-photon transition element $M^{(2)}$, which connects the initial state $\mid a\rangle$ to the continuum state $\mid q\rangle$ through all dipole-allowed intermediate states.

II.1 The form of the light–matter interaction

The standard expression for light matter interaction comes from minimal coupling ( $\mathbf{p}\rightarrow\mathbf{p}+e\mathbf{A}$ ), which gives the Hamiltonian

[TABLE]

If the spatial dependence of $\mathbf{A}$ can be neglected the diamagnetic term, $\sim\mathbf{A}^{2}$, can be removed through a unitary transformation of the wave function, $\Psi\rightarrow U\Psi$, with

[TABLE]

and the transformation of the time-dependent Schrödinger equation,

[TABLE]

This gives one remaining interaction term,

[TABLE]

which is usually referred to as the velocity gauge form. A different unitary transformation

[TABLE]

can be employed to find the alternative length gauge form

[TABLE]

for details see e.g. Ref. Selstø and Førre (2007). Here it is worth noting that in order to arrive at Eq. (7) from Eq. (5) it is necessary to assume that the potential term in the Hamiltonian in Eq. (4) commutes with $U$. This is obviously true for the Coulomb interaction with the nucleus, as well as between the electrons. However, due to the non-local nature of the Hartree-Fock exchange potential this is not the case within the Hartree-Fock approximation. Only by adding the RPAE class of many-body effects can the invariance between the two forms be restored Lin (1977). Close agreement between the two forms is often considered a quality mark for more elaborate calculations. Since the agreement is trivial for any local potential it is considered a necessary, albeit not sufficient, property.

With linearly polarized light we may now write the transition matrix elements from Eq. (LABEL:onephoton2) in length gauge as

[TABLE]

or in velocity gauge as

[TABLE]

These non-correlated transition matrix elements can be chosen to be real in Eq. (8) and imaginary in Eq. (9) by use of real radial wave functions.

Similarly the two-photon matrix element in Eq. (II) can be written as

[TABLE]

in length gauge (and similarly with the $ep_{z}/m$ operator and vector potentials $A_{\omega}A_{\Omega}$ in velocity gauge). An important difference compared to one-photon absorption is that the two-photon matrix element is intrinsically complex for above-threshold ionization, i.e. when the XUV photon energy exceeds the atomic binding energy, $\hbar\Omega>-\epsilon_{a}>0$.

The atomic contribution to the quantum beating of the side band at energy, $\epsilon_{q}=2n\hbar\omega+\epsilon_{a}$, in a RABBIT experiment is the phase difference between the quantum path where the XUV harmonic $\hbar\Omega_{>}=\left(2n+1\right)\hbar\omega$ is absorbed and an IR-photon is emitted and that where both an XUV harmonic, now of energy $\hbar\Omega_{<}=\left(2n-1\right)\hbar\omega$, and an IR-photon is absorbed. Eq. (II.1) shows the most important path, but contributions will also come from the reversed time-order where the IR photons are exchanged before absorption of any XUV photon. For this latter path there is in the general case no on-shell intermediate state that can contribute. It is thus assumed to be of less importance and is consequently often neglected. While this is a justified approximation for calculations in length gauge the situation is very different in velocity gauge as we will see below.

II.2 The time delay

Following the usual RABBIT formalism Dahlström and Lindroth (2014), we construct the phase shift of photoelectrons that take two different quantum paths leading to the same final state with momentum along the common polarization axis of the fields, $\mathbf{\hat{z}}$, as

[TABLE]

At this emission angle only the zero magnetic quantum number contributes to the ionization process, $m_{a}=0$. We use the following short-hand notation,

[TABLE]

where subscripts $\mathrm{a}$ and $\mathrm{e}$ stand for IR absorption and emission, respectively (do not confuse the subscript $\mathrm{a}$ with the quantum number label $a$ for the initial atomic state), where $M_{\mathrm{a/e}}(\ell)$ depend on angular momentum $\ell$ of the final $q$–state. The atomic delay for emission along $\mathbf{\hat{z}}$ can be calculated for sideband $2n$ as

[TABLE]

Similarly the one-photon phase shifts of the photoelectron in the $\mathbf{\hat{z}}$ direction are

[TABLE]

where we use short-hand notation for the one-photon matrix elements, $M_{>/<}(\ell)\equiv M^{(1)}(q_{>/<},\Omega_{>/<},a)$, with final photoelectron wave number $k_{>/<}$ and angular momentum $\ell$, after absorption of a photon with angular frequency $\Omega_{>/<}$. Eq. (II.2) can be used to compute the Wigner-like delay at sideband $2n$ along $\mathbf{\hat{z}}$ as

[TABLE]

We point out that the definition of Wigner delay using Eq. (14) breaks down at resonances that typically have large phase variations over the photon energy of the probe field Dahlström et al. (2012a). In the following we refer to the quantity $\tau_{A}-\tau_{W}$ as the delay difference induced by the laser field in RABBIT. We make a distinction between this “exact” delay difference and the approximate continuum-continuum delay that can be derived using asymptotic continuum functions, $\tau_{CC}\approx\tau_{A}-\tau_{W}$ Dahlström et al. (2013).

III Method

While the RPAE-approximation has been used to include electron correlation effects for the interaction with the ionizing XUV photon, the subsequent above-threshold interaction with the IR field has been limited to a static atomic interaction in our earlier studies Dahlström et al. (2012c); Dahlström and Lindroth (2014); Lindroth and Dahlström (2017). Here we discuss, in some detail, how this approximation can be lifted. The calculations are performed with a basis set obtained through diagonalization of effective one-particle Hamiltonians in a radial primitive basis of B-splines deBoor (1978), in a spherical box. For each angular momentum $\ell$ this one-particle Hamiltonian reads:

[TABLE]

It includes the (non-local) Hartree-Fock potential (HF), $u_{{\rm HF}}$, for the closed shell with $N$ electrons and a correction, $u_{{\rm proj}}$ (also non-local). The latter is called a projected potential (for the explicit form see Sec. III.1 below) and it ensures that any excited electron feels an approximate long-range potential with $N-1$ electrons remaining on the target. Since it is projected on virtual states it does not affect the occupied HF orbitals. The projected potential allows us to include some effects already in the basis set, that would otherwise be treated perturbatively through the RPAE-iterations. The eigenstates to $h_{0}^{\ell}$ form an orthonormal basis with eigenenergies $\epsilon_{i}$ that is used for the description of the occupied orbitals, but it is also used to span the virtual space of the photoelectron.

We start with writing the dipole interaction between one electron and the electromagnetic fields as

[TABLE]

where $d_{\Omega_{j}}$ is a time-independent operator that describe the coupling of the atom to the field with angular frequency $\Omega_{j}$ and $\xi$ is used to set the outgoing boundary condition for the interaction. With linearly polarized light the length gauge expression is $d_{\Omega_{j}}=ezE_{\Omega_{j}}$, where $ez$ is the dipole operator component along the field polarization. Consider now an electron in occupied orbital $\mid a\rangle$, i.e. in an eigenstate to the one-particle Hamiltonian in Eq. (III). When it absorbs ($\sigma=1$), or emits ($\sigma=-1$) photons it will acquire correction terms to its wave function of the type:

[TABLE]

where the superscripts and subscripts label sequences of interactions with photons (signs and angular frequencies) by joint increasing primes. Expressions for the correction terms can be found through the time-dependent Schrödinger equation

[TABLE]

by collecting the contributions that scale linearly with the field $E_{\Omega_{j}}$ and oscillate with $\exp[-i(\sigma\Omega_{j}+\epsilon_{a}/\hbar)t]$ as

[TABLE]

For a single electron case the desired one-photon correction to the wave function is simply obtained from Eq. (19), which we call the one-electron first-order perturbed wave function,

[TABLE]

where the sum over $p$ runs over all states (including also the continuum). For a many-electron system, however, there are more effects to consider. The starting point is then a Slater determinant $\mid\{ab\ldots n\}\rangle$ (where curly brackets denote anti-symmetrization, $\mid\{ab\ldots n\}\rangle\equiv(\mid ab\ldots n\rangle-\mid ba\ldots n\rangle\ldots)/\sqrt{n!}$ ). The field-corrected wave function will also be a Slater determinant, but now the orbitals are as given by Eq. III, i.e

[TABLE]

Since the interaction with the other electrons is accounted for by the Hartree-Fock potential the possible changes in it due to the interaction with the electromagnetic field have to be considered. We will return to this question below, but here we mark that the sum on $p$ in Eq. (20) is restricted to unoccupied states for the many-electron case (below the sums will be marked $exc$ to include only these “excited” states). This is simply what is expected from the Pauli exclusion principle. Alternatively, the restriction of $p$ to excited states can be understood from Eq. (21) where a $\sigma=+1$ excitation of orbital $a$ into $b$, as given in Eq. (20), will cancel the $\sigma=-1$ excitation of orbital $b$ into $a$ and vice verse.

We are interested in photoionization processes that happen when $\sigma\hbar\Omega+\epsilon_{a}>0$. This implies that there is a pole in the denominator of Eq. (20) that must be treated with the proper boundary condition and continuum integration. An efficient way to do this is to use exterior complex scaling (ECS) of the radial coordinate,

[TABLE]

which enforces the outgoing boundary condition for the unbound states. The eigenenergies of the orbitals are complex in general using ECS, and it has the advantage that the integration over the continuum is effectively performed by a sum over a discretized representation of all excited states, $p$, as written in Eq. (20).

The terms in Eq. (18) proportional to the product $E_{\Omega_{j}^{\prime}}E_{\Omega_{j}}$ that oscillates with $\exp[-i(\sigma^{\prime}\Omega_{j^{\prime}}+\sigma\Omega_{j}+\epsilon_{a}/\hbar)t]$ are

[TABLE]

for the case where the $\sigma^{\prime}$ interaction with frequency $\Omega_{j^{\prime}}$ happens after the $\sigma$ interaction with frequency $\Omega_{j}$. The second-order correction for a single electron are

[TABLE]

which simply builds on the first-order correction. Next, we need to define a notation for the corrections that includes all possible time orders by writing a square bracket around the signs, $[\sigma^{\prime},\sigma]$, and frequencies, $[\Omega_{j^{\prime}},\Omega_{j}]$, to include all joint permutations of primes on signs and frequencies. The second-order correction for a single electron with summed time orders is simply

[TABLE]

III.1 One-photon RPAE

The many-body response to the interaction with the photon is neglected in Eq. (20), but the bulk of these effects can be added through the RPAE method Amusia (1990), where certain sub-classes of many-body effects are included through the iterative solution of the equations for the coupled channels. Another name for RPAE is time-dependent Hartree-Fock Jamieson (1970), and we will here use that point of view to derive the expressions we need. With the HF-approximation each orbital is described as moving in an average potential from the other orbitals, and its matrix element between any orbitals $m,n$ (occupied or unoccupied), is:

[TABLE]

where the Coulomb interaction is given by

[TABLE]

Our starting point is a Slater determinant constructed from orbitals that are solutions to Eq. (III), with the Hartree-Fock potential defined as in Eq. (28). When the electrons interact with the field and acquire perturbations according to Eq. (III) the potential itself will change, $u_{HF}\rightarrow u_{HF}+\delta u$. This gives rise to additional paths for orbital $a$ to absorb or emit one photon with phase factor $\exp[-i\sigma\Omega_{j}t]$. In Eq. (28) we replace $b\rightarrow b+\rho_{\Omega_{j},b}^{(1:\sigma)}$, let the potential work on orbital $a$, and identify new terms to the excited states, $p$, that are linear in the electric field and oscillate with $\exp\left[-i(\sigma\Omega_{j}+\epsilon_{a}/\hbar)t\right]$ as

[TABLE]

In the case of absorption of a photon, $\sigma=1$, this implies that the second term in Eq. (30) is generated using a perturbed wave function that describes virtual emission of a photon, $\sigma=-1$. Adding Eq. (30) as an additional source term to the right-hand side of Eq. (19) leads to coupled equations for the correlated perturbed wave functions for absorption and emission of a photon,

[TABLE]

Use of Eqs. (20) and (30) leads to the final expression

[TABLE]

where the exchange interactions are written out explicitly. The upper part of Fig. 1 shows the Goldstone diagrams for $\mid\rho^{(1:\sigma)}_{\Omega_{j},a}\rangle$, where Fig. 1 (b) is the uncorrelated absorption of a photon $\Omega_{j}$, corresponding to the first term on the right-hand side of Eq. (III.1). Figs. 1 (c) and (d) account for the electron–hole interaction in forward propagation, corresponding to the second and third terms, while Figs. 1 (e) and (f) account for ground-state correlation effects, corresponding to the forth and fifth terms on the right-hand side of Eq. (III.1). The last term in Eq. (III.1) removes the projected potential, introduced in Eq. (III), which we take to be the monopole interaction with a given hole $c$,

[TABLE]

It will cancel the corresponding part of $\delta u$ of Fig. 1 (c) and (i) with $a=b=c$ and $K=0$. 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.

III.2 Two-photon RPAE

We now derive the interaction with two photons for the multi-electron case. The second interaction with the field can stimulate either the excited electron or the remaining hole from the first interaction. The latter effect arise when the staring point is a Slater determinant and the corrected wave function is of the form given in Eq. (21). The net result is a coupling of the wave functions associated different holes in Eq. (III), by the hole-hole dipole interaction in the source term of Eq. (18). Collecting the terms proportional to $E_{\Omega_{j^{\prime}}}E_{\Omega_{j}}$ from Eq. (18) that oscillate with $\exp[-i(\sigma^{\prime}\Omega_{j^{\prime}}+\sigma\Omega_{j}+\epsilon_{a}/\hbar)t]$, we write

[TABLE]

where the source terms on the right-hand side contain both time orders. In Eq. (III.2) the second line accounts for the interaction with the excited electron, Fig. 2 (a) and (b), while the third line accounts for hole transfer from another orbital, Fig. 2 (c) and (d). The minus on the third line comes from Wick’s theorem, which is evaluated using the Goldstone rules associated with the diagrams in Fig. 2 Lindgren and Morrison (1986).

The next step is to consider the many-body response. Second-order corrections to the Hartree-Fock potential can generate terms proportional to $E_{\Omega_{j^{\prime}}}E_{\Omega_{j}}$. By letting $b\rightarrow b+\rho_{\Omega_{j},b}^{(1:\sigma)}+\rho^{(2:[\sigma^{\prime},\sigma])}_{[\Omega_{j^{\prime}},\Omega_{j}],b}...$ in Eq.(28), and collecting the terms that oscillate with $\exp[-i(\sigma^{\prime}\Omega_{j^{\prime}}+\sigma\Omega_{j}+\epsilon_{a}/\hbar)t]$ we arrive at:

[TABLE]

The forward propagating ($\sigma=1$), direct contributions from lines two and three are depicted in Fig. 2 (k-l), and those from lines four and five in Fig. 2 (i) (only one of the two time-orders is shown). Another set of contributions, that will have the right oscillations, are the first order corrections from the Hartree-Fock potential when they, just as the dipole operator in Eq. (III.2), work on the corrected wave functions. This gives corrections

[TABLE]

for which the direct contributions are depicted in Fig 2 (e) and (g), and also

[TABLE]

where the direct contributions are depicted in Fig 2 (f) and (h). Note though that in both Eq. (III.2) and Eq. (III.2) the case when $j^{\prime}$ and $j$ are interchanged is to be added. Finally, there are second-order corrections that stem from the fact that the expression in Eq. (III) uses intermediate normalization, which means that the occupied orbitals, $\mid a\rangle$, are normalized and orthogonal to the corrections, $\mid\rho^{(1:\sigma)}_{\Omega_{j},a}\rangle,\,...$, while $\mid\psi_{a}\rangle$ is neither normalized nor orthogonal to $\mid\psi_{b}\rangle$. These corrections for the second-order interaction depend on the inner-product of the first-order corrections to the wave functions,

[TABLE]

Again the direct contribution is depicted in Fig. 2 (j) for one of the time-orders. The contributions from Eqs. (35 - 38) should now be added as source terms to Eq. (III.2) and we can write down the equation for the second order correction including the many-body response:

[TABLE]

The term with $-u_{\mathrm{proj}}$ compensates for the projected potential, which, as mentioned above, is important only for numerical convergence.

III.3 Calculation of two-photon matrix elements

For a RABBIT calculation with photoelectron energy $2n\hbar\omega+\epsilon_{a}$, we need two specific second-order correlated perturbed wavefunctions for orbital $a$ from Eq. (39),

[TABLE]

that include absorption of a smaller XUV photon and absorption, ($\mathrm{a}$), of a laser photon, as well as absorption of a larger XUV photon with emission, ($\mathrm{e}$), of a laser photon, denoted $\mid\rho_{\mathrm{a/e}}\rangle$ for brevity. Given $\mid\rho_{\mathrm{a/e}}\rangle$ we may directly extract the two-photon matrix elements needed for the calculation of the atomic delay, c.f. Sec.II.2. However, due to the on-shell above threshold contributions to the diagram in Fig. 2 (a), the construction of $\mid\rho_{\mathrm{a/e}}\rangle$ for the time-order where the XUV pulse is absorbed first involves an integration over a double pole and is not trivial. To circumvent this problem we first calculate the two-photon matrix element for the diagrams in Fig. 2 (a-d) and treat the additional corrections to $\mid\rho_{\mathrm{a/e}}\rangle$ separately. The different steps are detailed below.

The contributions from Fig. 2 (a-d) are calculated directly from the first order corrections $|\rho^{(1:+)}_{(2n\mp 1)\omega,a}\rangle$. In length gauge, the diagrams in Fig. 2 (a) and (c) amount to:

[TABLE]

where TO:1 stands for first time order, while the diagrams in Fig. 2 (b) and (d) amount to

[TABLE]

where TO:2 stands for the second time order. The final state $q$ is here an eigenstate to the effective one-particle Hamiltonian at the sideband kinetic energy $\epsilon=\epsilon_{a}+2n\hbar\omega$. As described in Ref. Dahlström and Lindroth (2014); Lindroth and Dahlström (2017) the numerical representation of the radial part of $\mid q\rangle$, denoted $P_{q}(r)$, is a solution of

[TABLE]

which can be reformulated as a system of linear equations for the unknown coefficients $c_{i}$ when expanded in B-splines

[TABLE]

For the case in Fig. 2 (a), where the first photon is of an XUV-wavelength causing ionization, and the second integral is between two continuum states, the integral in Eq. (41) will not converge for any finite interval on the real axis. The integration is instead performed numerically out to a distance far outside the atomic core, but within the unscaled region ($a_{0}\ll r<R_{C}$), while the final part of the integral is carried out using analytical Coulomb waves along the imaginary $r$-axis as described in Ref. Dahlström and Lindroth (2014). The numerical stability is monitored by comparison of different “break points” between the numerical and analytical descriptions. The integrals in Eq. (III.3), on the other hand, converge inside the numerical box since the IR-field can only induce a localized correction to the wave function.

The diagrams in Fig. 2 (e–j), and their exchange/switched time-order counterparts, can all be calculated by connecting converged first-order corrections with a single Coulomb interaction. Finally the diagrams in Fig. 2 (k–l) are found in an iterative procedure following that for the first order correction, Eq. (III.1). With the use of the projected potential, Eq. (33), all monopole terms are removed from the iterative procedure and the integral over the remaining Coulomb interaction does indeed converge on a finite interval. Therefore, it can be treated numerically inside the computational box. Separating the two-photon perturbed wave function in the lowest order contributions [Fig. 2 (a–d)] and the rest,

[TABLE]

and the the remaining contribution to the two-photon matrix element, $\delta M_{\mathrm{a/e}}$, can be deduced directly from $|\delta\rho_{\mathrm{a/e}}\rangle$, giving the final result:

[TABLE]

Of the four contributions in Eq. (41 - III.3) it is natural to assume that the first term in Eq. (41), Fig. 2 (a), is by far the dominating because it suffers from a zero in the denominator of the perturbed wavefunction. In contrast, Fig. 2 (b–d) are all connected with rather large denominators and should be small in general. The concept of cc-delays Dahlström et al. (2012a); Dahlström et al. (2013), where a photoelectron interact with a laser field after photoionization, derives from the assumption that the total two-photon process is well described by Fig. 2 (a) with use of a suitable long range potential, such as the projected potential in Eq. (33) Dahlström et al. (2012c); Dahlström and Lindroth (2014). Here we will show that this assumptions is close to the truth for calculations in length gauge, but wrong in velocity gauge.

III.4 Gauge-invariance

As discussed above the RPAE-approximation can be shown to produce gauge-invariant results Lin (1977) for the one-photon processes. This holds when the approximation is used consistently and without truncations. For example, the sum over core orbitals in Eq. (III.1) cannot be truncated and the orbital energies should be eigenvalues to the one-particle Hamiltonian used and cannot be replaced with experimental ionization energies. With these constraints we are here able to demonstrate gauge invariance also for the two-photon RPAE-approximation as will be seen below.

IV Results

Here we present calculations of atomic delays for photoelectrons emitted along the polarization axis $\mathbf{\hat{z}}$, as defined in Eqs. (II.2)–(14). In Figs 3–8 the horizontal axis labeled photon energy means the total photon energy absorbed by the photoelectron. In the case of atomic delays this implies the XUV-photon energy plus or minus the laser photon energy for laser absorption or emission, respectively.

IV.1 Neon 2p

The results for photoionization from the $2p$ orbital of neon are presented in Fig. 3 (length gauge) and Fig. 4 (velocity gauge). The RPAE iterations account for correlation effects from all three orbitals ($1s$, $2s$, $2p$), and the diagrams are evaluated using HF orbital energies. It is striking that the length gauge result is completely dominated by correlated XUV absorption followed by uncorreleted photoelectron–IR interaction, represented by the diagram in Fig. 2 (a). Only very small corrections, less than an attosecond, are found from the reversed time-order process (b), uncorrelated hole–field interations (c–d) and general correlated two-photon processes (e–l). This finding is supported by the comparison with experiment in Ref. (Isinger et al., 2017b), where good agreement was found over a large energy interval in length gauge using only the diagram in Fig. 2 (a) with experimental values substituted for the HF orbital energies. The results are more subtle in velocity gauge. The XUV first with uncorrelated photoelectron–IR interaction appears to be a reasonable approximation that deviates by a few attoseconds from the full calculation, but when the reversed time-order process is added (IR first) the deviation from the full calculation increases. Similarly, adding the hole–field interactions increases the deviation of the atomic delay further. Only the full two-photon RPAE calculation gives identical results in velocity and length gauge as seen in Fig. 4. The agreement between the gauges can be viewed as a validation test of the implementation.

IV.2 Argon 3p

The atomic delay for ionization from argon $3p$ is displayed in Fig. 5. The RPAE iterations account for effects from all five orbitals ($1s$, $2s$, $2p$, $3s$, $3p$), and the diagrams are evaluated using HF orbital energies. The delay is larger and changes more dramatically in argon as compared to neon. The velocity gauge result from Fig. 2 (a) alone underestimate the delay with around 40 as below the Cooper minimum and overestimate it by more than 50 as above. Including the full set of diagrams illustrated in Fig. 2 (a–l) leads to agreement between the length and velocity results within the numerical accuracy of the calculation.

In a truncated calculation, where the RPAE iterations account only for effects from the two outer shells $(3s,3p)$, there is a remaining difference between the two gauges, as shown in Fig. 6. The deviation from the full result is of the same order of magnitude for the two gauges, which implies that there is no clearly preferable gauge for the truncated two-photon RPAE calculation.

IV.3 Argon 3s

The atomic delay for ionization from argon $3s$ with photoelectrons emitted in the polarization direction is displayed in Fig. 7. The RPAE iterations account for effects from all five argon orbitals $(1s,2s,2p,3s,3p)$. While Koopman’s theorem states that the binding energy is equal to minus the HF orbital energy, which for $3s$ is $\sim 34.8$ eV, the true ionization energy is only $\sim 29.2$ eV. Therefore, we must substitute the HF orbital energies with experimental values for meaningful comparison with experiments. This simple procedure can be justified, since it corresponds to the inclusion of additional classes of diagrams Pan and Kelly (1989), but only at the price that the results again become gauge dependent. It is known from one-photon absorption experiments on argon that the cross section of $3s$ is affected by the strong $3p$ photoionization channel through electron correlation effects Amusia (1990). For two-photon processes a coupling from $3p$ to $3s$ can be directly stimulated by the second photon through the diagrams in Fig. 2 (c–d). The question now arises if such hole–field coupling effects can influence the atomic delay in argon?

The Cooper minimum in the $3s$ ionization cross section can be understood as a “replica” of the Cooper minimum in the $3p$ ionization channel. In more detail, the $3s$ minimum is caused by an interference effect between the direct path ($3s$) and correlated path ($3p\rightarrow 3s$), which results in very different ionization delays. While the $3p$ delays show a large negative peak (Fig. 5), the $3s$ delays show a large positive peak shown in Fig. 7. This conclusion is consistent with prior works based on RPAE (Dahlström et al., 2012b; Kheifets, 2013; Dahlström and Lindroth, 2014; Bray et al., 2018a) and Time-Dependent Local-Density-Approximation (TDLDA)  Magrakvelidze et al. (2015). Oddly, the large positive peak has not been observed in experiments  (Klünder et al., 2011; Guénot et al., 2012), while the negative peak has been reproduced experimentally using RABBIT Schoun et al. (2014); Palatchi et al. (2014). In our earlier studies of atomic delays, we have only accounted for Fig. 2 (a) and we have found that both the sign and position of the $3s$ delay peak is sensitive to correlation effects Dahlström and Lindroth (2014). Here, we find that the contributions from the remaining diagrams in Fig. 2, are not insignificant for $3s$ in argon, as seen in Fig. 7, and that the main additional contributions come from the hole–field coupling in Fig. 2 (c). However, the sum of all diagrams in our complete two-photon RPAE calculation does not resolve the discrepancy with argon $3s$ experiments at the Cooper minimum because the sign of our final delay peak remains positive and its position is not significantly altered (much less than an electron volt).

Very close to the Cooper minimum, where the one-photon matrix element goes through zero, it has been found harder to achieve good numerical accuracy. The scatter of break-points, see Sec. III.3, is indicated by error-bars in Fig. 7.

IV.4 Argon $3s-3p$

Finally, we show the difference in atomic delay for photoelectrons ionized from the two outer orbitals in argon, $\tau_{A}^{(3s)}-\tau_{A}^{(3p)}$, in Fig. 8. Here we display both the calculated result for electrons emitted in the polarization direction $\mathbf{\hat{z}}$, and for angular integrated detection, which is the configuration used in current RABBIT experiments Klünder et al. (2011); Guénot et al. (2012); Salières . We find that the atomic delay difference is not affected by the choice of detection in the region of the $3s$ Cooper minimum at $\sim 40$ eV. In contrast, the atomic delay difference is strongly altered close to the $3p$ Cooper minimum at $\sim 50$ eV due to the choice of detection, and the delay peak is reduced due to angular integration in agreement with the experimental angle-integrated results for argon $3p$ Palatchi et al. (2014) and calculations Dahlström and Lindroth (2016).

V Discussion

V.1 Gauge dependence

In general, our calculations show that much larger contributions arise from the reversed time-order, where the less energetic photon is absorbed first, in velocity gauge than in length gauge. This can be understood using a simple analytical calculation. For a Hamiltonian, $h$, with a local potential the length and velocity form of the dipole operator satisfy

[TABLE]

By assuming the dipole approximation, the vector potential for a given angular frequency and mode can be written as:

[TABLE]

where the two complex exponents can be physically interpreted as the drivers for emission and absorption of laser photon by the atom, respectively. Using the relation, $\mathbf{E}=-d\mathbf{A}/dt$, the expression for the electric field amplitude along the polarization axis $\mathbf{\hat{z}}$ is

[TABLE]

for emission and absorption, respectively. The gauge invariance of on-shell matrix elements follows from

[TABLE]

provided that $\left(\epsilon_{s}-\epsilon_{r}\right)=-\Omega$ for emission (photon creation) and $\left(\epsilon_{s}-\epsilon_{r}\right)=\Omega$ for absorption (photon annihilation).

Off-shell matrix elements are generally different. Consider the two-photon transition matrix element from initial state $\mid 0\rangle$ to a final state with $\epsilon_{f}=\epsilon_{0}+\Omega_{1}+\Omega_{2}$ via an intermediate state $\mid i\rangle$:

[TABLE]

the velocity gauge result is then a factor

[TABLE]

times the length gauge result. Eq. (51) has a maximum at $\epsilon_{i}=\epsilon_{0}+\left(\Omega_{1}+\Omega_{2}\right)/2$, and at this maximum it amounts to

[TABLE]

Therefore, we expect to find large differences for individual diagrams when, as in in typical RABBIT situation, the XUV photon has an energy of $20-40$ IR photons. We have indeed seen that the second time-order (TO2), which is always off-shell, is much more important in velocity gauge than in length gauge. Contributions for intermediate excited states and continuum states closely above the ionization threshold are likely to dominate, and for those the enhancement factor in Eq. (51) will be of quite some importance.

V.2 Universality of cc-delays in argon

In Fig. 9 we show the difference between atomic delay and Wigner delay, $\tau_{A}-\tau_{W}$, for argon from orbital $3p$ and $3s$. Panel (a) shows the low energy region with the $3s$ Cooper minimum, while panel (b) shows the high energy region with the $3p$ Cooper mininum. The approximate continuum–continuum delay, $\tau_{CC}$, is shown for comparison and it is calculated using the analytical expression of Eq. (100) from Ref. Dahlström et al. (2012a). The analytical cc-delay takes into account both long-range phase effects and long-range amplitude effects based on the Wentzel-Kramers-Brillouin (WKB) approximation, which gradually breaks down at low kinetic energies Dahlström et al. (2012a). Recently, excellent numerical agreement between $\tau_{A}-\tau_{W}$ for argon $3p$ and neon $2p$ was reported at very low kinetic energies using the diagram in Fig. 2 (a) Lindroth and Dahlström (2017). This suggested that the concept of “universality” goes beyond the analytical predictions of Ref. Dahlström et al. (2012a), down to much lower energies close to the ionization threshold, where the WKB approach is not applicable.

Surprisingly, the delay difference for argon $3s$ does not follow the universal curve [indicated by green dashed curve with data for $3p$ in Fig. 9 (a)], but instead shows irregular deviations at low electron energies close to the $3s$ Cooper minimum in Fig. 9 (a). The one-photon amplitude goes through zero at a photon energy of $\sim 40$ eV, which results in increased numerical uncertainty. The errorbars in Fig. 9 reflect the scatter between the different “break points”, c.f. Sec. III.3, and they signify that the observed deviations of $3s$ from $3p$ are real and that the universal trend is indeed broken, despite our limited numerical accuracy in this region. When only laser-stimulated continuum transitions are included in the calculation [Fig. 2 (a)], the deviation from the universal curve is not very large. The major part of the deviation comes from the remaining diagrams [Fig. 2 (b)–(l)], which suggests the importance of additional ways for the atom to interact with the fields when the single-photon XUV ionization process goes to zero. Because there are no resonances in the energy region shown for argon with RPAE, the irregular behaviour must be associated to the argon $3s$ Cooper minimum. Above and below the $3s$ Cooper minima, we find that the $3s$ delay difference agrees with the universal curve of $3p$, which indicates that correlation effects beyond the diagram in Fig. 2 (a), are significant only close to the exact photon energy region where the otherwise dominant one-photon correlated XUV ionization process vanishes.

What about the $3p$ Cooper minimum? This minimum is a “typical” Cooper minimum Cooper (1962) that arises due to a zero in the dipole transition $p\rightarrow d$ in XUV photoionization. The partial $3p$ cross-section does, however, not go to zero because the $p\rightarrow s$ dipole transition remains finite at all XUV energies. The deviation of $3p$ from the universal curve [indicated by the gray full curve in Fig. 9 (b)] is found to be small when laser stimulated continuum transitions are considered [Fig. 2 (a)] and very surprisingly even smaller when the full set of diagrams are included [Fig. 2 (a)–(l)].

Similar small deviations from the universal curve can be spotted in Fig. 5 of Ref. Dahlström et al. (2012b) for both argon $3s$ and $3p$, but because the effects are small compared to the associated atomic delays, they were not given much attention. More recently, irregular deviations from the “universal” curve of up to 20 as close to the argon $3p$ Cooper minimum was reported in Ref. Bray et al. (2018b). Our calculations show that such deviations are orders of magnitude too large and that they most likely arise due to an inconsistent description of combined correlation and field effects. The reduction of $3p$ deviation from the universal curve down to sub-attosecond precision in Fig. 8 (b) is most likely due to our improved description of the final state, where the effective spherical projected potential is substituted by self-consistent final state correlation effects [Fig. 2 (k)].

VI Conclusion

In this work we have shown that full two-photon RPAE calculations of atomic delays give gauge invariant results and that effects beyond the universal IR–photoelectron continuum–continuum transitions are rare, but do occur in special cases. In particular, we have found that the argon $3s$ Cooper minimum suffers from a non-universal delay because the correlated XUV dipole moment for photoionization vanishes, so that other processes, including XUV–hole interaction, may play an important role for the two-photon process. In contrast, we find that there are no such deviations from the universal delay curve for $3p$ in argon. Any $3p$ deviations that we find are on a sub-attosecond time scale, which disproves the strong $3p$ deviations recently proposed in Ref. Bray et al. (2018b) using a hybrid RPAE+TDSE approach.

Despite our best efforts, we have not been able to explain the discrepancy between theory and experiment for the argon $3s-3p$ atomic delays. This is because the full two-photon RPAE calculation still shows a positive peak in the atomic delay peak that is absent in experiments Klünder et al. (2011); Guénot et al. (2012); Salières . We note that recent simulations using Time-Dependent Density Functional Theory (TDDFT) (Sato et al., 2018) have generated results for the argon $3s-3p$ delay, in better agreement with experiments in this energy region. The authors of Ref. (Sato et al., 2018) attribute this success to their consistent treatment of the interaction with both light fields, as compared to the hybrid TDLDA+CC result in Ref. Magrakvelidze et al. (2015). However, the results of hybrid approaches, such as RPAE+CC Kheifets (2013) and TDLDA+CC Magrakvelidze et al. (2015) where the effect of the laser field is treated by a simple time shift given by analytical formulas Dahlström et al. (2012a); Pazourek et al. (2013), are mostly consistent with our new results. We cannot support the conclusion that an inconsistent description of the fields is the reason of the disagreement with experiments, because our present study does imply a consistent treatment of many-body effects for both fields. Still, it is hard to compare the many-body effects included with TDDFT (or TDLDA) with the present calculation and it is very well possible that the difference between the calculations lies here.

In closing, we wish to stress that XUV photoionization of argon is associated with strong satellite peaks that have not been considered in the present work, but have been studied in detail by Wijesundera and Kelly using many-body perturbation theory for the one-photon ionization process Wijesundera and Kelly (1987). A direct comparison between the partial cross-section for $3s^{-1}$ and the dominant satellite $3p^{-2}(^{1}D)3d$ from Ref. Wijesundera and Kelly (1987), shows that the satellite process does dominate in the photoelectron energy region of the $3s$ Cooper minimum with a cross-section of $\sim 0.02$ Mb as compared to our present value for the $3s$ partial cross-section of $\sim 0.003$ Mb located at a photon energy of $\sim 40$ eV using RPAE with all atomic orbitals and experimental energies in length gauge. Therefore, in order to better understand the discrepancy between experiments and theory it would be helpful to acquire atomic delays for larger energy ranges, but also to study the one-photon and two-photon partial cross-sections for $3s$ to be able to locate the exact position of the associated Cooper minima. Continued studies of shake-up processes in attosecond science, that go beyond the hybrid MCHF+CC approach of Ref. Feist et al. (2014), is desirable and maybe the right path to solve the long-standing argon delay puzzle.

Acknowledgments

The authors acknowledge support from the Knut and Alice Wallenberg Foundation and the Swedish Research Council, Grant No. 2014-3724, 2016-03789 and 2018-03845.