Photoelectron -- residual-ion entanglement in streaked shake-up ionization of helium
Hongyu Shi Uwe Thumm








Department of Physics, Kansas State University,
Manhattan, Kansas 66506, USA



§ ABSTRACT 

Streaked photoelectron emission spectra access the correlated dynamics of photoelectrons and residual target electrons with attosecond temporal resolution. We calculated ab initio single-ionization spectra for photoemission from helium atoms by co-linearly polarized ultrashort XUV and assisting few-femtosecond IR pulses.

Distinguishing direct and shake-up ionization resulting in ground-state and excited (n=2,3) residual ions, respectively, we examined the effects of the correlated photoemission dynamics on the photoelectron phase-accumulation as a function of the observable photoelectron detection direction and kinetic energy, and XUV – IR pulse delay. We tracked the dynamical evolution of the residual ion in relative streaked photoemission delays and found dominant contributions for shake-up emission from the residual ion – photoelectron interaction. 

These are in very good and fair agreement, respectively,  for n=2 and n=3 shake-up  photoemission along the pulse-polarization directions, with previous experimental and theoretical investigations [M. Ossiander et al, Nature Phys 13, 280–285 (2017)] and reveal a strong photoemission-direction dependence for shake-up ionization due to the coupling between the photoelectron and evolving residual-ion charge distribution in the IR-laser field. 



Photoelectron – residual-ion entanglement in streaked shake-up ionization of helium
    Hongyu Shi and Uwe Thumm
    September 9, 2024
===================================================================================





§ INTRODUCTION
 

Attosecond (1 as = 10^-18 s) streaking is an experimental technique for investigating the time-resolved ultrafast electronic dynamics in atoms <cit.>, molecules <cit.>, and solids <cit.>. Attosecond streaking experiments provided the first direct measurement of the oscillating electric field in intense light pulses and confirmed the generation of ultrashort sub-infrared (IR) cycle extreme ultraviolet (XUV) pulses with attosecond duration <cit.>. In a typical streaking experiment, an ultrashort XUV-pump pulse ionizes the target, emitting photoelectrons (PEs) into the field of a femtosecond infrared (IR) probe pulse  (art3D). This allows the recording of PE energy (or momentum) spectra for a set range of time delays τ between the two pulses. The τ-dependent spectra reveal information about the photoemission dynamics with attosecond resolution in time as "streaking traces", i.e., as an oscillating PE-energy (or momentum) dependence on τ that retraces the temporal profile of the IR streaking pulse with a characteristic absolute photoemission time delay relative to the IR-laser carrier electric field. While absolute streaking delays cannot be measured directly for photoemission from gaseous targets <cit.>, dynamical information can be accessed in streaked PE spectra from relative streaking delays between the centers of energy of streaking traces that are related to photoemission from two energetically distinct electronic levels <cit.>. 



Relative streaked photoemission time delays of 21± 5 as were measured between photoemission from the 2s and 2p orbitals of neon by Schultze et al. <cit.> in a milestone experiment that triggered an intense debate about the role of electronic correlation during the single ionization of an atomic target. This debate was fueled by several numerical models for streaked photoemission underestimating the measured relative photoemission time delay by approximately a factor of two, even though the calculations included electronic correlation at various levels of approximation <cit.>.

Accounting for electronic "shake-up" excitation during XUV photoemission from neon atoms, Isinger et al. <cit.> compared many-body perturbation theory calculations with measured energy-integrated sideband oscillations in interferometric photoemission spectra. Based on these RABITT (reconstruction of attosecond beating by interference of two-photon transitions) experiments, the authors suggested that the discrepancy between the relative streaking time delay of 21± 5 as obtained experimentally in Ref. <cit.> and in previous streaking calculations <cit.> might be due to shake-up processes. The authors of Ref. <cit.> worded this suggestion very carefully. They pointed out that two of their co-authors <cit.> could not explain the mismatch of streaking experiments and streaking calculations by numerically modeling RABITT spectra, employing a similar theoretical many-body-perturbation theory approach as Ref. <cit.> that includes electronic correlation. The authors provide good evidence, based on RABITT spectra, for the importance at the attosecond timescale of shake-up excitation during the photoemission of neon. However to the best of our knowledge and despite early theoretical testimony for the significance of dynamical electronic correlation on streaked photoemission <cit.>, the direct theoretical confirmation through streaking calculations of the 21± 5 as delay in Ref. <cit.> is still outstanding.

While the attosecond correlated ionization dynamics for multi-electron targets remains a subject of debate, remarkably, for prototypical helium, Ossiander et al. <cit.> found measured and calculated relative streaking time delays to agree at the (one) attosecond timescale. In their investigation, ground-state helium is singly ionized, and the dynamical response of the residual ion is analyzed by comparing direct ionization, resulting in ground-state (n=1) and shake-up ionization to a coherent superposition of (n = 2, 3) sub-levels 


    t_s = t_EWS + t_CLC + t_DLC .

Calculating t_EWS and t_CLC with attosecond accuracy, the authors were able to isolate a temporal shift due to electronic correlation during n=2 and n=3 shake-up ionization of about six and several ten attoseconds, respectively.

The Eisenbud-Wigner-Smith (EWS) contribution, t_EWS, represents the PE phase accumulation during the ionization of helium atoms in the XUV-pump pulse and includes the short-range part of the electron-electron interaction. It was calculated by direct numerical solution of the bound-continuum transition-matrix element of the two-electron system [cf., Eq. amp below] <cit.> using the exterior complex scaling (ECS) method <cit.>.
While this numerical approach is state-of-the-art for helium, its value for transparently revealing physical details of the correlated dynamics during single ionization of helium is limited.

The delay contributions t_CLC and t_DLC are due the simultaneous interaction of the PE with the residual target charge distribution and IR-laser electric field and termed "Coulomb-laser coupling (CLC)" and "dipole-laser coupling (DLC)". 
The CLC term includes the combined effect of the monopole term of the residual charge distribution and IR-laser <cit.>. 
The DLC term represents the effect of the PE interaction with the dipole-component of the charge distribution of the residual ion.
This term tends to be very small or negligible for direct ionization, but provides an important contribution for shake-up ionization, especially for large PE-detection angles relative to the XUV- and IR-pulse-polarization direction <cit.>. 

Applying this scheme to more complex atoms, 
such as experimentally convenient neon and argon gases <cit.>, is computationally challenging and requires concessions (most noticeably mean-field approximations), the accuracy of which (at the one attosecond time scale) still needs to be confirmed.

Understanding the physics behind the individual delay contributions in Eq. t_s for prototypical helium targets is crucial for the development and scrutiny of dependable attosecond streaking models for many-electron targets. It requires a highly differential study of streaked photoemission, going beyond the detection of PEs within a small solid angle about the IR- and XUV-pulse-polarization direction, in addition to the reliable modeling of the target electronic response to the incident XUV and IR pulses. 
For interferometric photoemission from helium, PE-emission-angle-dependent spectra were measured by Heuser et al. <cit.>. These RABITT experiments detected a striking phase shift of the sideband yields for electron-detection angles that are larger than ≈ 50 degrees relative to the polarization direction, even in the direct emission channel (leaving in the ground state). 
A related RABITT study for helium targets was recently carried out by Fuchs et al. <cit.>, who measured the three-dimensional PE-momentum distribution in coincidence with the residual ions, in order to examine experimentally the dependence of RABITT time delays (i.e., relative phase shifts between sidebands in RABITT spectra) on the PE angular momentum and coupling of continuum states in the assisting IR-laser pulse. 
Their measured delays of up to 12 as between outgoing s– and d–wave PEs, further emphasizes the need for highly differential (angle-resolved) data for the comprehensive analysis of (streaked or interferometric) time-resolved photoemission. 

Revisiting the sensitivity of relative streaking delays on the correlated electric photoemission dynamics reported in Ref. <cit.> for direct and shake-up streaked photoemission from helium along the polarization direction of the XUV and IR pulse, we carried out ab initio calculations for angle-resolved direct and (n=2,3) shake-up photoemission.
In Secs. <ref> and  <ref> we review our approach for 
numerically propagating the time-dependent Schrödinger equation (TDSE) for helium exposed to external electromagnetic fields and outline the calculation of two–electron states satisfying single–ionization boundary conditions. 
Sections  <ref> and  <ref> describe our numerical calculation of different contributions (EWS, CLC, and DLC) to the experimentally observable relative streaking delay, while Sec. <ref> briefly discusses noteworthy elements of our computational method, for which further details are given in the Appendices. 
Section <ref> contains our results for and discussion of
streaked photoemission spectra and streaking delays for direct and shake-up ionization (Sec. <ref>),
the EWS delay contribution for XUV single ionization of helium (Sec. <ref>),
the CLC delay contribution (Sec. <ref>), 
the DLC delay contribution in relation to the dipole moment of the residual charge distribution (Sec. <ref>), and a multipole analysis of the evolving residual charge distribution
and its back-action on the PE. Our summary and conclusions follow in Sec. <ref>.
Throughout this work we use atomic units unless indicated otherwise. 




§ THEORY
 



 §.§ Coupled radial equations
 

We describe the single ionization of ground-state para-helium atoms by ultrashort XUV pulses with and without assisting phase-coherent delayed (or advanced) IR-laser pulses by solving the TDSE in full dimensionality,

    tΨ = HΨ = [∑_i=1,2 (H_i + H_Fi ) + V_12] Ψ .

The Hamiltonian H consists of the single-electron Hamilton operators for the interaction of each electron with the atomic nucleus, H_i, the electronic correlation, V_12, and the electron—external-field-interaction terms, H_Fi, given by, respectively,

    H_i     = -_i/2 - 2/r_i ,  
    
    V_12   = 1/ r_2 -  r_1 ,

and 

    H_Fi = ℰ(t)    · z_i        (length gauge)
    
    - A(t)   ·z_i       (velocity gauge) .


The electron position vectors,
r_i  ,  i=1, 2, are chosen so that the atomic nucleus coincides with the origin of our coordinate system. 
The external electromagnetic-field pulses are assumed to be linearly polarized along the quantization (z) axis with cosine-squared electric-field envelopes of pulse length T,

    ℰ(t)    = ℰ_XUV(t-τ) + ℰ_IR(t) 
    ℰ_j(t)    = ℰ_j,0sin(ω_j t) cos(π t/T_j)^2 ,  j = XUV, IR,

and corresponding vector potential A(t) = -∫_-∞^t t'ℰ(t'). The delay (advance) of the IR-pulse center with respect to the center of the XUV pulse is designated by 
τ < 0 (τ >0).


In our applications to photoionization by linearly polarized light,
the magnetic quantum number of the spinless two-electron system, M=0, is preserved and will be omitted below. The relevant set of quantum numbers needed to specify a partial wave, λ = {l_1,l_2,L}, thus comprises the individual angular-momentum quantum numbers l_i according to H_i , i=1, 2, and the total angular-momentum quantum number L.
Expanding the two-electron wave function

    Ψ( r_1,  r_2, t) = 
    𝒮1/r_1r_2∑_λψ_λ(r_1, r_2, t) 𝒴_λ( r_1,  r_2)

in terms of real-valued generalized spherical harmonics

    𝒴_λ( r_1,  r_2) =
           ∑_m_1l_1     l_2     L
     m_1     -m_1     0 Y_l_1, m_1( r_1) Y_l_2, -m_1 ( r_2) ,

we express the TDSE TDSE as the set of coupled equations for the radial partial waves ψ_λ <cit.>

    ∑_λ' H_λ,λ'ψ_λ' = tψ_λ ,

with (r_1,r_2)-dependent angular-coupling-matrix elements

    H_λ,λ' = 𝒴_λH𝒴_λ' .

The square brackets in Eq. genharmonics denote  Clebsch–Gordan coefficients and 𝒮 the symmetrization operator for the two-electron spatial wave function of para-helium. Details about the matrix elements H_λ,λ' are given in Appendix <ref>.






 §.§ Partial-wave analysis of entangled PE states
 

The exact PE states for single ionization of He, Ψ, are solutions of the time-independent Schrödinger equation,

    H_0 Ψ = E Ψ ,
       H_0 = H_1 + H_2 + V_12 ,

with energy eigenvalue E = -Z^2/n_1^2 + k_2^2/2. They are subject to incoming-wave boundary conditions, since the time-reversed ionization process starts with asymptotic states of specific kinetic energy, momentum, or angular momentum. Depending on the chosen boundary condition, Ψ can be expressed in different ways, e.g., as states with a given total angular momentum L, 𝒮|n_1,λ,k_2⟩, or as states, 𝒮|n_1,l_1,m_1, k_2⟩, that asymptotically merge into product states. We expand the total-angular-momentum states in the basis of generalized spherical harmonics [Eq. (<ref>)],

    ⟨ r_1,  r_2|n_1,λ,k_2⟩
       = 1/r_1 r_2∑_l'_1,l'_2ψ_λ',n_1,k_2(r_1, r_2)𝒴_λ'( r_1,  r_2) ,

where λ' = {l_1',l_2',L}. The radial partial waves then satisfy the boundary condition <cit.>

    ψ_λ',n_1,k_2(r_1, r_2) r_2→∞⟶ δ_l_1,l'_1δ_l_2,l'_2 r_1 R_n_1,l_1^(Z=2)(r_1) √(2/π)
       sin[k_2 r_2 - π l_2/2 +1/k_2ln(2k_2 r_2) + σ_l_2 + δ_n_1^λ(k_2)] .

They depend on the radial wave function, R_n,l^(Z=2),  of the residual ion in the bound state ψ_n,l,m^(Z=2)( r) and the energy-dependent Coulomb phase shift,

    σ_l(η) = [Γ(l+1+η)] ,

with Sommerfeld parameter η = -Z/k. The sin factor in Eq. bc derives from the asymptotic form of Coulomb function of the first kind, F_l(kr), and includes the additional energy-dependent phase δ_n^λ(k) due to the short-range electronic interaction during the emission process.

To provide the |l_2,m_2⟩ partial-wave contributions of emitted electrons with asymptotic momentum magnitude k_2, we decouple the outgoing spherical-wave part of |n_1,λ,k_2⟩ in Eq. EL with the unitary transformation

    |n_1,l_1,m_1,l_2,k_2⟩ = |n_1,l_1,m_1,l_2,-m_1,k_2⟩
       = ∑_Ll_1    l_2    L
     m_1    -m_1     0^-δ_n_1^λ|n_1,λ,k_2⟩ ,

where the first equation indicates that we will drop the redundant quantum number m_2=-m_1 from now on.
This allows us to compose entangled two-electron states with asymptotic momentum vector k_2 by combining these partial waves into asymptotic plane-wave states, 

    |n_1, l_1, m_1,  k_2⟩
       = ∑_l_2^l_2/k_2^-σ_l_2 Y_l_2,-m_1^* ( k_2)|n_1,l_1,m_1,l_2,k_2⟩ ,

where l_1 relates to the bound electron and l_2 to the PE. To simplify the notation, 
we drop subscripts and symmetrization operator from now on and refer to the symmetrized scattering states as

    |n, l, m,  k⟩ = 𝒮|n_1, l_1, m_1,  k_2⟩ .






 §.§ XUV only delays
 

Exposure of helium to just the XUV-pump pulse releases PE wave packets whose center is displaced relative to a fictitious free electron that starts moving at the atomic nucleus
with the constant final PE velocity at the instant when the center of the XUV pulse coincides with the atomic center. This "XUV-pulse-only" photoemission delay,

    t_XUV( k) = t_EWS( k) + t_C(k,r) ,

consists of the EWS delay t_EWS and the Coulomb delay 

    t_C(k, r) = 1/k^3[1-ln(2kr)] .

The EWS delay corresponds to the phase accumulation during the absorption of XUV radiation and release of the PE <cit.>. 
The Coulomb delay accounts for ongoing phase accumulations of the PE in the long-range -1/r Coulomb potential of the residual ion. It is identical for all atoms and does not converge as the PE propagates to the detector.




  §.§.§ Dipole transitions


In first-order time-dependent perturbation theory, the photoemission yield and PE-phase information are calculated based on dipole-transition-matrix elements <cit.>. For photoemission by an XUV pulse ℰ_XUV(t) from the initial helium ground state, |i⟩, into asymptotic spherical waves sph_basis and plane waves with momentum k plane, 
the dipole-matrix elements are, respectively,

    M_s   = n, l_1, m, l_2, k(z_1+z_2)i ,
     
    M_p    = n,l,m, k(z_1 + z_2)i

and depend on the energy-dependent phase δ_n^λ. 

For single-photon ionization from the helium ground state, dipole selection rules restrict accessible final-state quantum numbers to L=1, M=m_1+m_2=0 and odd l_1+l_2. Thus, for direct ionization  the only accessible final configuration is (n,l_1,m,l_2)=(1,0,0,1). In contrast, for n=2 shake-up ionization, the five channels 
(2,0,0,1), (2,1,0,0), 
(2,1,0,2), (2,1,± 1,2) 
can be reached. 
Similarly, n=3 shake-up ionization allows the 13 final configurations
(3,0,0,1), (3,1,0,0), (3,1,0,2), (3,1,± 1,2),
(3,2,0,1), (3,2,± 1,1),
(3,2,0,3), (3,2,± 1,3), (3,2,± 2,3).
In general, the number of dipole-accessible final channels is 2n(n-1)+1.





  §.§.§ EWS delay


The EWS delay carries element-specific information on the electronic correlation dynamics during the photo-release. It can be interpreted as the group delay experienced by the PE wave packet relative to a hypothetical reference wave packet <cit.>. 
It excludes the long-range Coulomb delay t_C and can be calculated numerically as the energy derivative of the phase of the dipole-matrix element for entangled spherical-partial-wave and plane-wave final PE states, respectively,

    t_EWS^n,l_1,m,l_2(k)    = E M_s , 
    
    t_EWS^n,l,m( k)    = E M_p .

As a measure for the net delay for direct or shake-up ionization to a specific sub-shell, we define the weighted average of Eq. tEWS:

    t_EWS^n,l( k) = ∑_mM_p^2 t_EWS^n,l,m( k)/∑_mM_p^2 .

Similarly, we average over {l,m} to define t_EWS^n( k).

Instead of calculating matrix elements that yield the delays in first-order perturbation theory [cf., Eq. tEWS], we can retrieve absolute delays for ionization by the XUV-pulse directly from the (non-perturbative) two-electron wave function twoelwf, based on the numerically calculated radial partial-wave functions ψ_λ(r_1, r_2,t). This approach is particularly well suited for our applications in Sec. <ref>, which are based on the numerical propagation of partial-wave functions (cf., Sec. <ref>). 
Projection of the numerically propagated wave function twoelwf on residual hydrogenic states ϕ_n,l,m^(Z=2)( r),
provides the PE probability density for direct (n=1) and shake-up (n>1) ionization 

    P_n,l,m( r,t)
       = 2∫ϕ_n,l,m^(Z=2)( r') Ψ( r',  r, t) (r')^2Ω'r'^2 ,

where the factor of 2 is due to symmetrization of Ψ( r_1,  r_2, t).


Assuming PEs are detected within the detector solid angle Ω_det centered around the direction k, the mean PE distances from the atomic nucleus are now given by integrating Ω over Ω_det as

    r_c^n,l( k, t) = ∑_m∫_Ω_det r P_n,l,m( r,t) r^2Ωr/∑_m ∫_Ω_det P_n,l( r,t)r^2 Ωr .

Similarly, we average over {l,m} to define 

    r_c^n( k, t) = ∑_l,m∫_Ω_det r P_n,l,m( r,t) r^2Ωr/∑_l,m∫_Ω_det P_n,l( r,t)r^2 Ωr .

Omitting the integration over Ω corresponds to an infinitely high PE-detection resolution. 
Fitting the averaged radial PE displacements at times after which


  * the PE wave packet has separated from the charge distribution of the residual ion and 

  * the PE velocity ṙ_c^α( k, t) has converged to the asymptotic velocity ṙ_c^α( k, t_∞) within a preset tolerance 

according to 

    r_c^α( k,t) 
       = ṙ_c^α( k,t_∞) 
    [(t-τ) - (t_EWS^α( k) + t_C(k_f, r_c))] , 
    α   = {n,l,m}, {n,l}, n  ,

we find the EWS delays t_EWS^α relative to the center of the XUV-pulse electric field.






 §.§ Streaked photoemission delay
 



  §.§.§ Absolute and relative streaking time delays


The delay t_XUV is not directly measurable with current experimental techniques, but contributes to the observable photoemission delay in streaked PE spectroscopy, where a probe pulse with an adjustable delay relative to the XUV-pump pulse is introduced <cit.>. The probe pulse, usually in the IR spectral range, maps temporal shifts in the photoemission process to observable shifts in asymptotic PE momenta. Recorded as a function of the pump-probe delay τ, streaked PE spectra allow the extraction of relative streaking time delays Δ t_s by comparing the temporal offset between energetically separable oscillatory photoemission yields (streaking traces). This method for deriving relative streaked photoemission time delays has been applied to electron emission from different atoms <cit.>, different initial electronic levels in the same target <cit.>, and to identify temporal shifts between direct and shake-up photoemission <cit.>. Differences of absolute streaking delays are experimental observables.

Absolute streaking time delays t_s for streaked photoemission from a specific electronic state relate temporal photoemission shifts to the phase of the probe-pulse electric field ℰ_IR(t). Numerical modeling allows the calculation of absolute delays in two distinct schemes, as either "streaked spectral delays (SSDs)" or "PE wave-packet delays (PWDs)". We refer to SSDs as delays that are extracted from streaked PE spectra 
by fitting the center of momentum of a given streaking trace 

    p_f(τ) = -aA_IR(τ + t_s) + b      (a,b > 0)

to the IR vector potential A_IR=-∫_-∞^t t'ℰ_IR(t'). In contrast, we name delays that are extracted from the wave-packet displacement on a numerical grid using Eq. EWS_xuv as PWDs. The accuracy of the SSD t_s relies only on the accuracy of the XUV–IR-delay-dependent PE-momentum spectrum and does not require the calculation of phase derivatives. For peak intensities between 10^8-10^12 W/cm^2, streaking delays and their numerical fitting accuracy are largely independent of the IR-probe-pulse intensity <cit.>.


According to Eq. t_s, the absolute SSD t_s is affected by τ-dependent momentum shifts induced by the probe pulse, in particular the Coulomb-laser and dipole-laser shifts, t_CLC and t_DLC, in addition to the probe-independent shift t_EWS.
For direct photoemission the short-range PWD contribution t_EWS to t_XUV in Eq. t_xuv merges into t_s when a sufficiently weak delayed IR probe pulse is added. In contrast, the long-range Coulomb delay t_C, being sensitive to the phase accumulation along the entire PE trajectory, is affected by the probe-laser field and transformed into the Coulomb-laser-coupling delay t_CLC <cit.>. 
We have verified that t_CLC depends on the PE-energy and -detection direction. For emission along the XUV- and IR-pulse-polarization direction it is identical for direct and shake-up ionization.





  §.§.§ Entangled residual Stark states


For shake-up ionization the probe-laser pulse causes an additional  dipole-laser-coupling temporal shift t_DLC <cit.>, even if the probe-pulse intensity is kept sufficiently low to make the quasi-static Stark effect in the ground-state-helium  entrance channel negligible (as is the case in typical streaking experiments). This delay contribution has no counterpart in t_XUV. In the exit channel the probe-laser mixes degenerate and nearly degenerate states of the residual ion, causing the formation of quasi-static Stark hybrid states with dipole moment d^n( k) as linear combinations of excited states with opposite parity. 
The formation of the residual dipole affects the energy of the two-electron system in the probe-laser electric field and, in turn, the final PE-momentum and streaking delay <cit.>.   

To describe the formation of quasi-static residual Stark states during shake-up ionization, we unitarily transform the scattering states |n,l,m, k⟩ plane to  Stark scattering states with the same (n, m) quantum numbers and total energy -Z^2/n^2 + k^2/2,

    |n,α,m, k⟩ 
    = ∑_l C_α,l^n,m|n,l,m, k⟩ ,

where α denotes the Stark quantum numbers <cit.>. The expansion coefficients {C_α,l^n,m} for the n=2,3 shells of are listed in Tabs. <ref> to <ref> of Appendix <ref>. 


After shake-up ionization by the XUV pulse, once the PE wave packet and residual charge distribution effectively no longer overlap, we represent 
the residual ion in a given 
n channel as a superposition of degenerate Stark states, 

    |*⟩ψ_^(n)( k) = 𝒩∑_α,m c_n,α,m( k)|n,α,m⟩ ,

with normalization constant 𝒩 = (∑_α,mc_n,α,m^2)^-1/2. To obtain the expansion coefficients from the numerically  provided accurate solution Ψ of the two-electron TDSE TDSE for XUV-only ionization, we project at a time t^⋆ = 5.5 fs after the center of XUV pulse

onto the two-electron entangled scattering states ⟨n,α,m, k| in Eq. (<ref>),

    c_n,α,m( k) = 
    ⟨n,α,m, k|Ψ(t^⋆)⟩ .

For photoemission along the XUV-pulse-polarization direction (k =  z), only 
m=0 terms contribute in Eq. bound_Psi. Thus,   
|*⟩ψ_^(n)( k) is rotationally symmetrical about the z axis and its dipole moment,

    d^(n)( k) = -*ψ_^(n)( k) rψ_^(n)( k) ,

aligned with k. In contrast, for PE-detection directions k that are not aligned with z, the combination of terms with m≠0 results in electric dipole vectors vec_dipole that, generally, do not align with z.

A sufficiently weak, quasi-static IR-streaking electric field does not noticeably distort the electronic structure of the initial He atom and residual ion. Its interaction with the XUV-pulse-generated dipole yields the energy shift

    Δ E_dip^(n)( k)     = ℰ_IR(τ) 
    *ψ_^(n)( k) z ψ_^(n)( k)
       = - A_IRτ    d^(n)_z( k)

and corresponding PE momentum shift

    Δ k^(n)   ≈-Δ E_dip^(n)( k)/k
    = -A_IRτd^(n)_z( k)/k .

This equation builds on the assumption that for shake-up emission to a specific residual shell n, PE wave packets corresponding to different asymptotic states
{|n,α,m⟩} in Eq. bound_Psi have the same center of momentum 
k = √(2(ω_XUV - I^(n)_p)).  

Our numerical results in Sec. <ref> below validate this assumption by revealing that, for PE momenta k within the spectral range of the XUV pulse, the relative magnitudes of the coefficients c_n,α,m( k)  keep their proportions (cf. sis2n).

At large PE – distances, the PE has accumulated the EWS phase shift, 
and the PE-momentum shift in the oscillating IR field can be approximated as <cit.>,

    k_f(τ) ≈(k + Δ k^(n))  k
    - A_IR(τ + t_EWS) .

Denoting the polar angle of k_f and k as θ_f and θ, respectively, and defining γ = θ_f-θ, projection of Eq. k_f onto k_f  yields


    k_f(τ)    = kcosγ
       - cosθ_f[A_IR(τ + t_EWS) + t_DLCA_IR(τ + t_EWS)τ]
       ≈ kcosγ - cosθ_f · A_IR(τ + t_EWS + t_DLC) ,

where

    t_DLC( k) = cosγ/cosθ_fd_z^(n)( k)/k≈d_z^(n)( k)/k cosθ .

For the IR-pulse peak intensity considered in this work, γ is less than 1.6^∘, making it negligible in Eq. tDLC. The approximation in Eq. k_f_appr is supported by our typical delays t_DLC being less than 5% of the IR period. 
We note that an alternative approximation is to treat A_IR(τ) as sinusoidal <cit.>, which results in

    t_DLC( k) ≈1/ω_IRtan^-1[ω_IRd_z^(n)( k)/kcosθ] .







 §.§ Propagation algorithm
 
 
We solve the TDSE TDSE using a finite-element discrete-variable representation (FE-DVR) of the partial-wave radial wave functions ψ_λ(r_1, r_2, t) on a 2D adaptive rectangular numerical grid for the radial coordinates (r_1, r_2). For details of our implementation of the FE-DVR method we refer to Ref. <cit.>, references therein, and Appendices <ref> and <ref> of the present work. 





  §.§.§ Length gauge


Expressing the electron interactions with the XUV and IR pulse HFI in the length gauge, we time-propagate the radial TDSE H_couple using the split-operator Crank-Nicolson technique <cit.> with equidistant time steps Δ t, 

    ψ_λ(t+Δ t) = ^- HΔ tψ_λ(t)
          = ^- H_0Δ t/2^- H_intΔ t^- H_0Δ t/2ψ_λ(t) + Δ t^3 .

This scheme assumes the time-dependence of H to be negligible during each time step. In all applications in Sec. <ref> below, we ascertain numerical convergence by selecting sufficiently small time steps for the accumulation of numerical errors in each time step being irrelevant.

The advantageous choice of H_0 = H_1 + H_2 allows the factorization of the propagation along the radial coordinates,

    ^- H_0Δ t/2 = 
    ^- H_1Δ t/2^- H_2Δ t/2,

since H_1 and H_2 commute. 
In the length gauge, the interaction operator H_int = H_F1+H_F2+V_12 in position representation does not contain derivatives. It thus couples all partial waves separately at any given point of the 2D numerical grid (without the need for error-susceptible finite differencing). The split-operator method is therefore particularly suitable for length-gauge calculation.
For all exponentiations we use the Cayley form of the time-evolution operator <cit.>,

    ^- HΔ t = (1+/2 H Δ t)^-1(1-/2  H Δ t) .





  §.§.§ Velocity gauge


For applications in Sec. <ref> below we employ both,  the length and velocity form of H_Fi [cf., Eq. HFI], in order to assess the numerical convergence of our PE-momentum distribution and temporal shifts. While formally quantum mechanics is gauge invariant, approximations introduced in the modeling and the progression of numerical errors tend to affect the accuracy of length- and velocity-gauge calculations unequally. Their comparison is therefore a common tool for assessing the accuracy of numerical models. In particular, as pointed out by Cormier and Lambropoulos <cit.> for the numerical computation of photoemission spectra for strong-field ionization based on partial-wave expansions, the rate of convergence in length and velocity-gauge calculations with the number of included partial waves tends to be different. For above-threshold ionization of hydrogen atoms, the authors obtained converged spectra including 10 partial waves in velocity-gauge calculations, while requiring in excess of 200 partial waves in the length gauge. The physical reason for the very different rate of convergence is a more effective inclusion of the external-field-driven electronic quiver-motion in velocity-gauge final PE states in term of the canonical PE momentum, k -  A(t). In contrast, the same anisotropic dynamics requires a large number of partial waves in the length gauge, where the final-state description hinges on the (time-independent) kinematic PE momentum k.

While this convergence comparison appears to clearly favor the velocity gauge, the split-operator scheme outlined above in part offsets its advantageous inclusion of the canonical PE momentum due to the appearance of the differential operator *z_i in Eq. HFI. In contrast to the length gauge and represented in spherical coordinates, the electron–external-field interaction H_Fi now involves the coupling of two-electron partial waves at adjacent grid points that are addressed in the finite-difference implementation of the *z_i operator on the (r_1,r_2) grid. For this reason, we abandon the split-operator propagation scheme when working in the velocity gauge and, instead, use a Lanczos algorithm <cit.> with an adaptive number of Krylov basis states to keep the numerical error below a set threshold. 

In a Krylov basis {Ψ, HΨ, …, H^K-1Ψ} with K states, where Ψ is the wave function at the current propagation time, H in Eq. TDSE is represented by a tridiagonal matrix

    H =
    [   α_0   β_1      ;   β_1   α_1   β_2      ;         β_2   α_2     ⋱      ;                 ⋱     ⋱ β_K-1;       β_K-1 α_K-1 ] .

The truncation error in the Lanczos algorithm we use to diagonalize H can be estimated as <cit.>

    ϵ_lanc^K≈Δ t^K/K!∏_i=1^K β_i .

We dynamically adjust K at each propagation time step Δ t to satisfy ϵ_lanc^K-1 < 10^-12, rather than using a fixed number of basis states as, e.g., Ref. <cit.>.







§ RESULTS AND ANALYSIS
 



 §.§ Streaking delays
 

Following <cit.>, we assume an IR streaking pulse with a Full-Width Half-Intensity Maximum (FWHIM) of 3 fs and a peak intensity of 4×10^11 W/cm^2. The XUV pulse has a photon energy of 90 eV, FWHIM of 200 as, and peak intensity of 10^12 W/cm^2. 
We use 40 equally spaced XUV–IR-pulse delays for a delay range of 2 IR periods, calculate the center of PE momenta within an infinitely narrow detection cone, and fit to the IR-pulse vector potential to extract the streaking delays according to Eq. fit. Figure <ref> shows our corresponding calculated streaked photoemission spectrum. The streaking traces for direct and n=2 shake-up emission are clearly separated.



Absolute streaking delays t_s for direct and n=2,3 shake-up emission extracted from the spectrum in streaking are show in ang_delay, together with delays for XUV photon energies of 100 and 110 eV. The dependence of t_s on the XUV-pulse photon energy and PE-detection angle for emission is complex, owing to the three distinct terms in Eq. t_s and different (n,l,m) sub-channels. We therefore separately examine contributions to t_s in the following subsections. The most prominent feature, a rapid decrease of t_s with PE-detection angle, is primarily due to the 1/cosθ factor in the DLC term [cf. Eqs. t_s and tDLC]. 







 §.§ XUV-only delays
 

Figure <ref> shows the EWS time delay, t_EWS, extracted from the PE wave packet according to Eqs. r_c_nl and EWS_xuv, for direct ionization (n=1) and all (n,l) shake-up sub-channels of the ionic L and M shells. These results are calculated for emission along the XUV-pulse polarization and XUV photon energies of 90, 100, 110, and 120 eV. Lines in between the small markers are added to guide the eye. For direct (black line) and 2s shake-up (red dashed line) our results in tEWS (a) agree with the phase-derivative calculation of Pazourek et al. <cit.> according to Eqs. tEWS and tEWS_n1l1. For shake-up emission to (2p) our results still agree with Ref. <cit.> with attosecond accuracy and show small discrepancies only at the sub-attosecond level. Generally, in the absence of resonances, EWS delays tend to decrease in absolute value for increasing XUV photon energy, due to the decreasing time left for strong PE interactions. This trend is followed for direct, 2p-, and 3p-shake-up emission. At the sub-attosecond time scale,  the onset for this expected decrease is not reached for 2s-shake-up emission for XUV photon energy up to 120 eV.  



Figure <ref> extends the sub-shell-resolved EWS delays, shown in tEWS  for emission along the XUV-polarization axis, to arbitrary PE-detection directions for XUV photon energies of 90, 100, and 110 eV. 
For 3p and 3d shake-up emission this graph clearly shows the expected decrease in magnitude of the EWS delay for increasing XUV photon energy and a more pronounced dependence on the PE-detection direction in comparison to 2p shape-up emission. Noticeable are also the general increase of the delay magnitude with the principle and angular momentum quantum numbers of the residual shake-up shake-up state, in compliance with the larger extent of the residual charge distributions. 
Delays for 1s, 2s, 3s shake-up emission are angle independent, shown in tEWS, and not duplicated in this graph.




Averaging the EWS photoemission delays for n=2 and 3 according to Eqs. r_c_n and EWS_xuv results in the emission-angle-dependent absolute sub-shell averaged delays displayed in tEWS_ang, in addition to the angle-independent delays for direct XUV photoemission.
The angular dependence is most pronounced for shake-up emission at the lowest XUV photon energy (90 eV). In comparison with the sub-shell-resolved delays in tEWS_sub, the angular dependence of the sub-shell-averaged delays in tEWS_ang
tends to follow more closely the 3p than the 3d delays. 
The less uniform angle dependence of the sub-shell-averaged delays for different photon energies, in particular, the crossing of n=3 delays in tEWS_ang, 
is a result of the sub-shell and XUV-photon-energy- and angle-dependent emission probabilities P_n,l,m( r_2,t) employed as weights in Eq. r_c_n.



To keep the contribution of the EWS delay to the absolute XUV delay t_xuv in perspective, we note that for direct emission by 90 eV XUV pulses the non-convergent Coulomb delay amounts to t_C = -53.8,-57.0,-59.1 (as) at detector distances r_2=0.5,2,5 (m), respectively. The negative long-range Coulomb delay thus constitutes an advance that dominates the EWS delay. 





 §.§ The Coulomb-laser-coupling delay
 



  §.§.§ Definition of CLC delays for hydrogen atoms


The CLC delay accounts for the long-range PE-phase accumulation. For emission along the XUV- and IR-pulse-polarization directions it does not noticeably depend on the target electronic structure and is approximately identical for H and He <cit.>. It is largely a classical quantity, which was found to agree within ∼ 1 as for quantum and classical-trajectory Monte Carlo calculations <cit.>. We assume its long-range classical behavior to apply for any final PE-momentum direction k (i.e, for any PE detection angle θ) and calculate it ab initio,  following Refs. <cit.>, according to

    t_CLC( k) = t^H_100_s( k) - t^H_100_EWS(k) ,

based on the  streaking [t^H_100_s( k)] and EWS delay 

    t^H_100_EWS(k)    = EC( k)z1,0,0
       = σ_l=1E_E=k^2/2

for ground-state hydrogen with outgoing Coulomb waves 

    ⟨ r|C( k)⟩ =
           1/r√(2/π)∑_l,m[^l/k^-σ_lY_l,m^* ( k)] F_l(η, kr) Y_l,m( r) .

We note that t^H_100_EWS is angle independent due to the symmetry of the hydrogen ground state. 

Figure <ref> shows the CLC delay for hydrogen as a function of the PE-detection angle for final PE energies between 17.03 and 85.41 eV, as given by Eq. HCLC. As expected, t_CLC decreases with increasing XUV photon energy (or asymptotic PE kinetic energy). 








  §.§.§ Classical-trajectory CLC calculation


To investigate the angular dependence for PE emission to large angles in CLC, we extended the classical calculation for emission along the XUV-pulse-polarization direction in Ref. <cit.> and performed an angle-resolved classical-trajectory simulation. We assume instant absorption of a single XUV photon with energy E_XUV to provide PEs at aligned initial positions and momenta, r_0 ∝z and k_0 ∝z, respectively, with total energy

    E_PE = E_XUV + E_100^(Z) =  k_0^2/2 - Z/r_0 ,

where E_100^(Z) is the ground-state energy of a hydrogen-like atom. This equation gives k_0 as a function of r_0. We determine r_0 = 0.55 a.u., such that the classically calculated CLC delay, t_CLC^cl, matches the results in CLC at zero degrees, and apply this value for all energies considered in this figure. 

Propagating classical PE trajectories starting at ( r_0,  k_0),  subject to the same IR field, we use for He streaking calculation for a range of XUV–IR pulse delays τ, we obtain τ-dependent asymptotic PE momenta, i.e., classical streaking traces, that allow us to extract the classical streaking delay t_s^H_100,cl. Substitution into Eq. HCLC now yields delays t_CLC^cl that are independent of the PE detection direction, suggesting that the angle dependence of t_CLC in CLC is a quantum-mechanical effect.






 §.§ The dipole-laser-coupling delay
 

Exposure of the residual dipole vec_dipole to the electric field of the IR streaking pulse results in the energy shift
Δ E_dip = - d_z^(n)( k)  ℰ_IR(τ)  
in Eq. E-dip.
This shift affects the entangled PE–residual-ion system and results in the DLC contribution t_DLC to the streaking delay in Eq. t_s.
The effect of t_DLC for PE emission along the XUV-polarization direction is discussed in Ref. <cit.>. In this subsection we examine t_DLC( k) for arbitrary PE detection directions k.

Applying Eq. tDLC, we calculate the coefficients c_n,α,m( k) according to Eq. c2s.

For the XUV-pulse parameters considered in the present work, we find that the ratios of the magnitudes of the complex-valued  coefficients c_n,α,m( k) remain approximately constant over spectral profile of the PE wave packet. This allows us to evaluate these coefficients at the central PE momentum k = √(2(ω_XUV - I^(n)_p)). Figure <ref> illustrates this proportionate change for the example of n=2 shake-up emission. 



The approximate DLC delays t_DLC for n=2 and n=3 shake-up ionization of para-helium and XUV pulses with photon energies of 90, 100, and 110 eV show a striking PE-detection-direction dependence at large angles that is easily traced to the 1/cosθ factor in the dipole model given by Eq. tDLC (tDLC_model). 
Within this model, direct photoemission has no DLC delay, since d_z^(0)=0.
The sign change for shake-up emission at θ≈ 40^∘ indicates a change in the orientation of the residual dipole at larger angles. 






We validate the DLC model results in tDLC_model by computing streaking, EWS, and CLC delays ab initio. The ab initio DLC delays shown in tDLC are obtained from Eqs. t_s and HCLC as 

    t_DLC( k) = t_s( k)  - t^H_100_s( k) 
        -[t_EWS( k) - t^H_100_EWS( k)]  ,

in good overall agreement with the model DLC delays for shake-up emission in tDLC_model. For most angles, this agreement is better than 10 as. It becomes worse at large angles near θ = 90^∘, possibly due to a combination of the model assumptions leading to Eq. tDLC and 
small numerical signal-to-noise ratios (very small yields) at large detection angles. As seen in Tab. <ref>, for the special case of PE emission along the XUV-polarization direction (0^∘), our DLC-model delays are in good quantitative agreement with the correlation delays 
derived from measured streaking delays in Fig. S10 of the supplementary material of Ref. et al. <cit.>.



The DLC model does not explain the angle dependence of the ab initio DLC delay for direct emission shown in the upper graph of tDLC. In addition, at the as timescale, it is at variance with our ab initio calculated delays (Tab. <ref>). 
A possible cause for this variation is that the approximate calculation of t_CLC for hydrogen according to Eq. HCLC does not accurately represent the CLC delay for helium. 
On the other hand, accurately calculating the CLC delay for direct photoemission from helium as 

    t_CLC^He( k) = t^n=1_s( k) - t_EWS^n=1(k) .

for the n=1 curves in tDLC would merely illustrate the difference of the two definitions, t_CLC^He - t_CLC, instead of a deviation from the  physically intuitive DLC model tDLC.






 §.§ Residual multipole analysis
 

During photoemission, the electronic probability density of the residual ion keeps evolving while interacting with the PE. After a sufficiently long time after the photo-release by the XUV pulse, the PE wave-function overlap with the charge distribution of the residual electron is small enough for the conditional electronic probability density,

    P^He^+( r, t) = Ψ( r,  r_PE ,t)^2/∫Ψ( r',  r_PE, t)^2 [3]r' ,

to be associated with the residual ion's electronic charge density, where the PE is assumed as a classical particle located at the maximum of the PE wave-packet probability, r_PE.

The multipole expansion of the electronic potential energy V^He^+ corresponding to the charge density P^He^+, assumed to be confirmed within a volume of radius a < r_PE, results in

    V^He^+( r,t) = ∑_l=0^∞∑_m=-l^l 1/r^l+1 C_l,m(t) Y_l,m( r)    (r > a) ,


    C_l,m(t) = -4π/(2l+1)∫_r⩽ a r^l Y_l,m^*( r) P^He^+( r, t) [3]r .

Using individual terms of this expansion as effective potentials in separate SAE calculations allows the quantitative comparison of the effects of given multipole orders on PE spectra and streaking time delays. Performing such simplified one-electron calculations for the parameters of the present work, we found dominant monopole and dipole interactions, while quadrupole and higher multipole orders are negligible.  

A few snapshots of the residual charge-density evolution are shown in multi_pole
for the ionization of para-helium in 90 eV, 10^12 W/cm^2 XUV pulses with a FWHIM of 200 as (no IR pulse is present). 
We assumed emission along the XUV-polarization direction (more specifically, the +z axis), such that only multipole terms with m=0 are present. The color-coded graphs show residual charge distributions P^He^+( r, t) for increasing classical PE distances r_PE from the atomic nucleus. Residual probability densities for direct and shake-up emission are displayed in Figs. <ref> (a) and (c), respectively. The corresponding coefficients C_1,0(t) of the dipole term in the multipole expansion multipole for direct and shake-up emission are shown in Figs. <ref> (c) and (d), respectively. 



The oscillation period of the contribution for direct emission in multi_pole (b) is about 88 as, corresponding to an energy of 47 eV, close to the energy difference between the undistorted K and L shells of , 40.8 eV. We therefore attribute these small-amplitude oscillations to weak quantum beats of the ground state against transiently populated n=2 levels. The rapid decay of these beats as the PE moves away from the residual ion illustrates the noticeable presence of electronic entanglement, even for direct emission and PE distances way outside the residual charge distribution.

For shake-up ionization the composition and evolution of the residual charge cloud are more complex than for direct emission [multi_pole (d)]. Now the dipole term starts oscillating with a period of about 120 as (35 eV). 0.7 fs after the ionization the oscillation period changes to 530 as (7.8 eV). These two quantum-beat frequencies are close to undistorted K- to L-shell and L- to M-shell excitation energies of , respectively. As expected, the oscillation amplitude for shake-up emission in multi_pole (d) is significantly larger than for direct emission in multi_pole (b). In addition, the oscillations for shake-up ionization persist long after those for direct emission have decayed. These findings are in agreement with an intuitively expected stronger correlation of the PE with electronically excited . 





 §.§ The polarization effect by IR field


In addition to the XUV photoionization dynamics generating polarized residual ions, the electric field of the IR streaking pulse   contributes to the polarization of the residual ion and may thereby affect angle-dependent shake-up photoemission yields and streaking delays. In this section we investigate this effect and show that, for the assumed IR-pulse intensity in the present work, the effect of the residual-ion polarization by the streaking pulse can be neglected.   

For hydrogen atoms this effect was investigated before by Baggesen and Madsen <cit.>,  who pointed out that the polarization of a target atom by the streaking IR-pulse may influence the angle-dependence of streaked photoemission yields and delays. They showed that, for sufficiently strong IR streaking fields, the streaked PE asymptotic energy acquires an additional phase shift relative to 
the IR vector potential A_IR(t), due to the fact that the dipole-laser coupling is proportional to the IR electric field ℰ_IR(t) and thus phase-shifted by 90^∘.

To quantify the distortion of the  residual ion by the streaking field, we calculate its static polarizability <cit.>, 


    β_n, l, 0    = β_n, l, 0^b + β_n, l, 0^c
        =2∑_n',l'^n' nn',l',0zn,l,0^2/E_n'^(0)-E_n^(0)
        +2∑_l'∫_0^∞kk,l',0zn,l,0^2/k^2/2-E_n^(0) ,


based on virtual excitations to hydrogenic bound states, {|n,l,m⟩}, (β_n, l, 0^b) and spherical partial-wave continuum states, |k,l',m⟩, (β_n, l, 0^c) for the nuclear charge Z=2. Our numerical values for the lowest four states of are in close agreement with the scaled theoretical values in Refs. <cit.> (Tab. <ref>). The literature values in Tab. <ref> are calculated for hydrogen atoms and multiplied by 1/16 to account for the scaling β_n,l,0∝ 1/Z^4 <cit.>. 



For the IR-pulse peak intensity assumed in this work of 4× 10^11W/cm^2, corresponding to a peak-electric-field strength of ℰ_IR,0 = 3.4-3 a.u., the distortion of the electric dipole moment of the ground state is 

    Δ d_1,0,0,max = 2 ℰ_IR,0 β_1,0,0 = 1.9-3 ,

where the factor of two derives from the maximal variation of the IR electric field over one IR oscillation period. 
For the n=2 shell of , the largest dipole polarizability is β_2,1,0 = 13.49, resulting in  

    Δ d_2,1,0,max = 2ℰ_IR,0β_2,1,0≈ 9.2 -2 .

The maximal IR-pulse-induced dipole in the n=2 shell is thus about 50 times larger than for the ground state, and the distortion of the ground state may be neglected in comparison with the distortion of excited states by the streaking IR field.

At the same time, the maximal IR-pulse-induced dipole in the n=2 shell in Eq. (<ref>) is about 16 times smaller than the dipole moment of 1.5 of the XUV-pulse-excited n=2 Stark states according to Tab. <ref> in Appendix <ref>. This is consistent with the numerical example 
in Fig. <ref>, where the maximal dipole coefficient for n=2 shake-up emission exceeds the maximal dipole coefficient for direct emission by an order of magnitude. The distortion by the IR streaking pulse of the XUV-photoemission-produced residual dipole is thus a small correction.  

This difference is less pronounced for n=3 shake-up excitation. Here, the  XUV-photoemission-produced residual-dipole moment is 4.5 (cf., Tabs. <ref> and <ref> in Appendix <ref>), exceeding 6.5 times the IR-pulse-induced dipole distortion,

    Δ d_3,1,0,max = 2ℰ_IR,0β_3,1,0≈ 0.69 .

These comparisons indicate that, for shake-up ionization into the K and L shells of , the distortion of the XUV-produced dipole by the IR pulse is negligible, while for n=3 shake-up ionization  further scrutiny is needed in this regard.

We note that once the IR distortion of becomes relevant, the distinction of three separate delay contributions to the streaking delay in Eq. (<ref>)  breaks down, and the EWS term needs to be modified to account for the influence of the IR-pulse during photoemission for small XUV – IR pulse delays. In the present work,  the fact that the three contributions in Eq. (<ref>)  correctly add up to the streaking delay, provides further evidence for the IR-pulse distortion of being negligible.  

Furthermore, the classical  trajectory analysis in Sec. <ref>, excluding the IR-pulse distortion of the residual ion, agrees for emission angles below ≈60^∘ with the quantum result and further supports the very minor role of the IR distortion of the residual ion in the present work.  

This analysis of the distortion of the residual ion by the IR pulse may provide a reference for futures streaking-pulse parameters that break the non-distorting dipole assumption. For a discussion of the IR-pulse distortion of helium in the entrance channel, we refer to Ref. <cit.>. For the description of relative streaking delays, the IR-pulse distortion of the neutral helium target affects direct and shake-up channels equally and thus tends to cancel in the calculation of relative streaking delays. 






§ SUMMARY AND CONCLUSIONS
 



In this work we confirmed the results of previous theoretical and experimental investigations for streaked direct and shake-up photoemission from helium atoms along the IR- and XUV-pulse linear polarization. We performed ab initio and model calculations and extended these earlier investigation to  arbitrary PE detection directions, scrutinizing EWS, CLC, and DLC contributions for direct and (n=2,3) shake-up ionization.

For large PE detection angles we found that CLC delays become angle dependent and tracked this dependence as a purely quantum-mechanical effect that cannot be explained with classical-trajectory simulations. For shake-up ionization we find dominant contributions to the streaking delay due to DLC. Confirming a physically intuitive model for the DLC effect, we validated the model predictions for the DLC delay against ab initio calculations over a large range of PE detection angles. 

Distinguishing the production of polarized residual ions due to the XUV-photoemission dynamics and possible contributions from its time-dependent polarization by the oscillating IR pulse, we estimated the IR-pulse-induced polarization contribution. Our model DLC calculations without including the IR distortion of the XUV-pulse-generated dipole agreeing with the DLC-delay contribution derived from ab initio calculations provide evidence for the IR-pulse distortion of the residual dipole being negligible. 








§ COUPLING-MATRIX ELEMENTS AND POISSON INTEGRAL
 

Each term in the Hamiltonian of Eq. TDSE can be represented as a band-diagonal matrix in an FE-DVR basis, and the computational cost for the propagation by one time step scales linearly with the total number of grid points. 
As basis function, we use polynomials defined by Gauss-Lobatto quadrature, scaled to fit each FE interval for the FE-DVR expansion of the radial partial-wave functions,

    ψ_λ(r_1, r_2, t) = ∑_i,j c_i,jϕ_i(r_1)ϕ_j(r_2) .

These basis functions are localized on the grid points of the 2D numerical (r_1 and r_2) grid and approximately orthonormal,

    ϕ_i(x_j) ∼δ_ij,     ∫ϕ_i(x)ϕ_j(x)x≈δ_ij ,

where {x_j} designates the set of grid points for either one of the radial coordinates, r_1 and r_2.

The numerical efficiency of the FE-DVR method derives from these polynomials being adapted to and defined on the employed numerical grid. 
The flexibility of the FE-DVR method stems from its applicability to any given (possibly unequally spaced) set of FE boundary points.




 §.§ Electron-electron interaction


The electronic interaction term, V_12, in the Hamiltonian is expanded in generalized spherical harmonics,

    V_12 = 4π∑_l=0^∞(-1)^l/√(2l + 1)r_ < ^l/r_ > ^l+1𝒴_l,l^0 ( r_1,  r_2) .

For the angular part of the electronic coupling-matrix elements we have

    𝒴_λ𝒴_l,l^0𝒴_λ'
    = 2l+1/4πδ_L,L'√((2l_1'+1)(2l_2'+1)(2L+1))
           ×[    l l_1'  l_1;    0    0    0 ][    l l_2'  l_2;    0    0    0 ]l     l_1'     l_1
     l     l_2'     l_2
     0     L     L

and calculate the radial matrix elements accurately using Poisson integrals <cit.>,

    ϕ_i ϕ_jr_ < ^l/r_ > ^l+1ϕ_i'ϕ_j' =δ_i,i'δ_j,j'
           ×[2l+1/r_i r_j √(ω_i ω_j) (2  T_l)^-1_i,j + r_i^l r_j^l/r_max^2l+1] ,

where

    (T_l)_ij   = ϕ_i[-1/2md^2/dr^2 + l(l+1)/2mr ] ϕ_j
       = -1/2mϕ_id^2/dr^2ϕ_j + δ_i,jl(l+1)/2mr_i^2 .




 §.§ External-field interaction


The length-gauge external-field-interaction term H_Fi in the TDSE [Eq. TDSE] contributes to the matrix element in Eq. H_couple and is evaluated according to

    𝒴_λℰ(t)· z_i𝒴_λ' = ℰ(t) r_i F_λ,λ'^(i)   (i=1,2) ,

where the angular integrals are

    F_λ,λ'^(i) = 𝒴_λcosθ_i𝒴_λ'
       =√(9(2l'_1+1)(2l'_2+1)(2L'+1)/(4π)^2)[ δ_i,1  l'_1   l_1;     0     0     0 ]
       ×[ δ_i,2  l'_2   l_2;     0     0     0 ][  1 L'  L;  0  0  0 ]δ_i,1    l'_1     l_1
    δ_i,2    l'_2     l_2
     1     L'     L   (i=1,2) ,

and the curly brackets denote the Wigner 9j symbols.

In the velocity gauge, H_Fi couples the angular function as 

    𝒴_λ- A(t) z_i𝒴_λ' = - A(t) G_λ,λ'^(i)   (i=1,2) ,

where

    G^(i)_λ,λ' = F^(i)_λ,λ'r_i + g^(i)_λ,λ'/r_i

and

    g_λ,λ'^(i) = -𝒴_λcosθ_i +sinθ_iθ_i𝒴_λ' .

In particular, the hydrogenic field-interaction matrices are <cit.>)

    g_λ,λ'^(1) =
       δ_l_2, l_2'∑_m_1l_1     l_2     L
     m_1     -m_1     0l_1'     l_2'     L'
     m_1     -m_1     0×
        (δ_l_1',l_1+1 l'_1𝒞_l_1,m_1 - δ_l_1,l'_1+1 l_1𝒞_l'_1,m_1)

and 

    g^(2)_λ,λ' =
       δ_l_1, l_1'∑_m_2l_1     l_2     L
     -m_2     m_2     0l_1'     l_2'     L'
     -m_2     m_2     0×
        (δ_l_2',l_2+1 l_2'𝒞_l_2,m_2 - δ_l_2,l'_2+1 l_2 𝒞_l'_2,m_2) ,

with

    𝒞_l,m
    = √((l+1)^2-m^2/(2l+1)(2l+3)) .






§ STARK STATES
 

The following tables give the matrix element of the unitary transformation between the asymptotic Stark states |n,α,m⟩ in Eq. Starkstate and Z=2 hydrogenic bound states |n,l,m⟩ within a given shell n. 
The first column  shows the eigenvalues of the electric dipole moment d_z^(n,α,m) for each Stark state.





These matrices are calculated by diagonalizing z in each degenerate subspace spanned by all states {|n,l,m⟩} for fixed n,m <cit.>. The subspaces with n=2, m=± 1 and n=3, m=± 2 are non-degenerate and have zero polarization with Stark quantum number α=0. The z component of the dipole in Eqs. (<ref>) and (<ref>) is calculated as

    d^(n)_z( k)    = -*ψ_^(n)( k)zψ_^(n)( k) ,
       = -𝒩^2 ∑_α,mc_n,α,m( k)^2 d_z^(n,α,m) ,

where the coefficients c_n,α,m( k) are defined in Eq. c2s.





§ OPTIMIZATIONS
 




 §.§ Selection rules
 

For para-helium Ψ( r_1,  r_2) is symmetric with respect to electron exchange. Spherical harmonics satisfying the relation

    𝒴_λ( r_1,  r_2) = 
    (-1)^l_1+l_2-L𝒴_λ( r_2,  r_1) ,

implies 

    ψ_l_1, l_2^L(r_1, r_2,t) = (-1)^l_1 + l_2 - Lψ_l_2, l_1^L(r_2, r_1,t)

and eliminates about half of the partial waves in Eq. twoelwf, allowing us to keep only those for l_1 ⩽ l_2. This simplification is compromised for the rectangular (non-quadratic) grids we used, where r_1 and r_2 have different spacial ranges, and l_1 and l_2 require different maximal values for numerical convergence.  

The symmetry properties of the Clebsch-Gordan  and Wigner 9j coefficients imply that  Eq. TDSEcouple couples partial waves with the same odd values of l_1+l_2-L. Since the initial ground state of para-helium (L=0, l_1=l_2) has even l_1+l_2-L, partial waves  with odd l_1+l_2-L do not contribute.





 §.§ Adaptive grid with optional sliding window


We propagate the radial wave function on an adaptive (r_1, r_2) grid, which we determine during the numerical calculation by adding a "detector region" outside the complex absorber potentials at the outer grid boundaries of r_1 and r_2 <cit.> (grid_layout). When the wave-function probability density in the detector region on either side reaches a threshold (typically  5-9 ), the length of that side is automatically extended, keeping the existing grid points and basis polynomials unchanged. For single ionization, one radial coordinate has a fixed length and no detector region is needed along its grid direction. In tests we found that this technique can accelerate the numerical throughput for single ionization 2 to 3 times. 








§ ACKNOWLEDGMENTS


We thank Keegan Finger, who participated in the early stages of this work within the "Research Experiences for Undergraduates (REU)" program of the NSF, as well as Van-Hung Hoang and Aihua Liu for stimulating discussions. 
This work was supported by NSF Grant No. 2110633 (Method development for the assessment of transient response in matter to and electronic states in arbitrary time-dependent fields)
and the Chemical Sciences, Geosciences, and Biosciences Division, Office of Basic Energy Sciences, Office of Science, US Department of Energy under Award DEFG02-86ER13491 (Attosecond electron dynamics). 



