Population transfer to high angular momentum states in infrared-assisted XUV photoionization of helium

An extreme-ultraviolet (XUV) laser pulse consisting of harmonics of a fundamental near-infrared (NIR) laser frequency is combined with the NIR pulse to systematically study two-color photoionization of helium atoms. A time-resolved photoelectron spectroscopy experiment is carried out where energy- and angle-resolved photoelectron distributions are obtained as a function of the NIR intensity and wavelength. Time-dependent Schrödinger equation calculations are performed for the conditions corresponding to the experiment and used to extract residual populations of Rydberg states resulting from excitation by the XUV + NIR pulse pair. The residual populations are studied as a function of the NIR intensity (3.5 × 1010 − 8 × 1012 W cm−2) and wavelength (760–820 nm). The evolution of the photoelectron distribution and the residual populations are interpreted using an effective restricted basis model, which includes the minimum set of states relevant to the features observed in the experiments. As a result, a comprehensive and intuitive picture of the laser-induced dynamics in helium atoms exposed to a two-color XUV–NIR light field is obtained.


Introduction
Modern table-top high-harmonic generation (HHG) sources routinely provide extreme-ultraviolet (XUV) attosecond 6 Authors to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. pulses that can be used in laboratory-scale pump-probe experiments, enabling the study of quantum phenomena with attosecond time resolution [1]. In these experiments, the XUV pulse is usually combined with a near-infrared (NIR) pulse. In an XUV-NIR configuration, attosecond resolution is obtained by using the well-characterized electric field of the NIR pulse as a time reference for the attosecond dynamics [1]. In order to unravel the dynamics, it is necessary to understand the interplay between the interactions of the atom with the NIR and XUV pulses. In this work, we systematically explore the NIR-assisted XUV ionization of helium from low to medium NIR intensities as a function of the NIR wavelength and the relative delay between the XUV pulse, containing frequency components both above and below the helium ionization threshold, and the NIR pulse.
As a model for more complex systems in XUV-NIR experiments, helium has been studied extensively both theoretically and experimentally. In a seminal paper from 2007, Johnsson et al demonstrated that a two-color XUV-NIR field can be used to control the ionization yield of helium on a subcycle timescale, i.e. by varying the time delay between the attosecond pulse train (APT) and the NIR pulse [2]. Observed fast oscillations of the He + yield as a function of the time delay between the APT and the NIR pulse were understood in a time domain picture in terms of an interference between wavepackets created in the presence of the NIR field by consecutive attosecond pulses within the APT. In a complementary frequency-domain picture, these oscillations may be understood as the result of an interference between ionization pathways via different intermediate states dressed by the NIR field [3,4]. In reference [2], it was further shown that the amplitude of the oscillations strongly depends on the intensity of the NIR laser pulse. Later works [5][6][7][8] investigated the intensity-and wavelength-dependence of the He + yield by means of cold target recoil ion momentum spectroscopy (COLTRIMS) and showed that resonant absorption of the XUV and NIR light by Rydberg states has a major influence on the observed dynamics, as it affects the phase associated with each ionization pathway. The energies of these Rydberg states are strongly influenced by the moderately intense NIR field, which induces an AC Stark shift, which can shift the states into or out of resonance with the available XUV frequencies. Several experiments have focused on this effect using the attosecond transient absorption spectroscopy (ATAS) technique, providing time-and frequency-resolved information on the time-dependent dynamics of laser-assisted ionization of helium [9][10][11][12][13][14][15][16][17]. For instance in reference [13] an isolated attosecond pulse was combined with an NIR field of moderate intensity (3 TW cm −2 ) in order to observe the subcycle dependence of the AC Stark shift on the XUV-NIR time delay.
Compared to photoion detection and ATAS, photoelectron spectropscopy (PES) can directly provide information on the population of Rydberg states [18]. PES experiments on helium exposed to a two-color XUV-NIR light field at moderate NIR intensities (few TW cm −2 ) revealed intensity-dependent modulations in the photoelectron angular distributions (PADs) [5][6][7][8], which were attributed to contributions of G-waves (photoelectrons with the angular momentum = 4), hence suggesting the excitation of |1snf ( = 3) Rydberg states. In reference [5], it was proposed that population transfer to these bound states occurs through an NIR-induced two-photon transition from the lower-lying |1s2p state, which is itself excited by the XUV radiation, but no quantitative analysis was provided.
The experimental results mentioned above indicate that signals observed in experiments on NIR-assisted XUV ionization can change significantly when the NIR intensity is tuned from a weak field regime where interactions may be described within first-order perturbation theory (i.e. involving processes where a single NIR photon is absorbed or emitted) to one where a moderately intense NIR field can induce multi-photon transitions. As many attosecond experiments are carried out at these intermediate NIR intensities, it is crucial to understand and study this transition in detail.
In this work we report time-, energy-and angle-resolved photoelectron spectroscopy measurements on NIR-assisted XUV ionization of helium using low (ca. 10 10 W cm −2 ) to intermediate NIR intensities (ca. 8 TW cm −2 ), which are accompanied by detailed time-dependent Schrödinger equation (TDSE) simulations in the single-active electron (SAE) approximation. Based on experimental and computational results we develop an effective restricted basis model that provides an intuitive picture of the photoionization dynamics.
The paper is organized as follows: in sections 2 and 3, we introduce the experimental setup and the numerical methods for solving the TDSE, respectively; in section 4 we discuss the experimental and numerical results and identify the Rydberg states that play a role in the experiment; in section 5, we introduce the effective restricted basis model, which emerges from the conclusions from the previous sections; in section 6 we compare the results of the model with the experiment and the TDSE simulations and provide a unified picture of the dominant pathways in the time-dependent two-color XUV-NIR dynamics underlying the experimental and numerical observations; in section 7, we conclude and provide an outlook.

Experiment
We use the same experimental set-up for generating highharmonics and measuring photoelectron momentum distributions of helium as described in reference [18]. Figure 1 shows a schematic diagram of the experimental set-up from the top view and the side view, parts (a) and (b), respectively. In short, 35 fs, 4.2 mJ NIR laser pulses at a repetition rate of 1 KHz are generated by a dual-stage femtosecond amplifier system (Komodo-Dragon, KMLabs). We split the laser beam by a 50% beam splitter into two arms, beam A and beam B. Beam A is further split into two parts by a half-moonshaped fused silica plate (FS-1) with a thickness of 1 mm in the proportion 70% (beam C) to 30% (beam D). For the NIR intensity dependence measurements beam D was used. Its intensity was controlled by closing an iris (not shown in the figure). For the XUV-NIR delay dependence measurement, for practical reasons, beam D was blocked and beam B was introduced using an insertable half-moon mirror HM (the lower 30% of the beam profile equivalent in shape to . The side views α and β correspond to the regions α and β in the top view. FS1 and FS2 are fused silica plates, HM is a half-moon shaped mirror, Ge is a germanium mirror, TM is a toroidal mirror, VMI is a velocity map imaging apparatus, Jet A is a pulsed gas jet to generate high-harmonics, Jet B is a pulsed gas jet to inject helium gas to the interaction region of the VMI. The HM mirror is inserted to perform XUV-NIR delay dependence measurements (see text). (c): the measured high-harmonic spectra as a function of the relative grating distance in the compressor of the laser system. The NIR wavelength is given in the top axis. The horizontal line appearing around 23 eV is due to a background signal caused by the CCD camera.
beam D). Both beams (beam C and beam B/D) are focused into an argon gas jet ( jet A) to produce high-harmonic emission in the XUV region. Beam C is reflected from germanium plates at a Brewster angle to suppress the co-propagating NIR beam, while beam B/D is reflected from a silver-coated surface and provides the NIR dressing field for the experiment. Beam B/D is transmitted through a second half-moonshaped, fused silica plate (FS-2), which blocks the XUV light generated by beam B/D in the jet. The XUV-NIR delay is adjusted by a translational stage placed in the path of beam B. XUV and NIR pulses have parallel polarizations and are both focused by a gold-coated toroidal mirror (ARW optical corporation) into a helium gas jet. The photoelectron momentum distribution resulting from the ionization of helium is measured by velocity map imaging (VMI) apparatus. The image is captured by a CCD camera. The XUV pulse is dispersed by a flat-field grating placed after the VMI apparatus and the spectrum is measured by a microchannel plate. We scan the photon energy of the high-harmonics in the XUV pulse by changing the grating distance in the compressor of the Ti: Sapphire amplifier. In figure 1(c), we present the measured high-harmonic spectra as a function of the grating distance in the compressor of the laser. As the grating distance is changed, the photon energy of high harmonics shifts monotonically. We calibrate the NIR intensity in the interaction region of the VMI photoelectron spectrometer by assuming that non-resonant two-color photoionization by the XUV-NIR pulses causes a ponderomotive energy shift, U p = |F NIR | 2 /4ω 2 NIR , in the photoelectron spectrum, where F NIR is the electric field strength and ω NIR is the angular frequency of the NIR pulse in atomic units. We also assume that the NIR pulse has no chirp in the range where the photon energy is scanned in this experiment.

Numerical simulations
We solve numerically the exact TDSE for a helium atom in the SAE approximation in the presence of the electric field of the combined XUV and NIR laser pulses. The code used here is described in reference [19]. In the calculations, the XUV and NIR pulses have a sine-squared envelope that is 4 and 40 NIR cycles long, respectively. Note that a 4-NIR-cycle long XUV pulse is used in order to obtain harmonics with a bandwidth corresponding to the experimentally observed one. The XUV pulse consists of four equally intense (3.5 ×10 10 W cm −2 ) odd harmonics of the fundamental NIR frequency, ranging from the 13th (H13) to the 19th (H19) harmonic. Moreover, calculations were performed where only one of these harmonics was included. The TDSE was solved by expanding the wavefunction on a 2D grid containing 150 000 radial points equally spaced by Δr = 0.15 au and including orbital angular momenta up to max = 16. Angle-resolved photoelectron spectra were evaluated up to E max = 0.4 au (ca. 11 eV) with a step of ΔE = 10 −4 au. An advantage of the code developed in [19] is that it allows us to calculate the photoelectron spectrum both at positive and 'negative' energies. While positive energies correspond to the production of photoelectrons that can be detected experimentally, the 'photoelectron spectrum' at negative energy permits an extraction of the residual occupation of Rydberg states of the atom at the end of both pulses. This allows us to track the residual population in Rydberg states as a function of the experimental parameters, i.e. the NIR intensity, the NIR and XUV wavelength and the relative time delay between the two pulses. Figure 2 shows a selection of experimental photoelectron momentum distributions for ionization of helium by the twocolor XUV-NIR field for different NIR intensities (top row) and wavelengths (bottom row). All experimental PADs shown in figure 2 and in the rest of the manuscript are slices through the 3D momentum distributions obtained from the experimental VMI images by performing an Abel inversion with the BASEX algorithm [20]. The rings in the first image of the first row of figure 2 correspond to ionization of helium atoms by harmonics H17 and H19. With the increasing NIR intensity new features appear at lower photoelectron momenta corresponding to ionization of Rydberg states prepared by harmonics H13 and H15. For higher NIR field intensities the angular distributions change revealing an eight-lobe structure. The observed features are labeled in figures 3(a) and (b), where distributions for λ NIR = 788 nm and I NIR = 3.5 ± 0.17 TW cm −2 are shown for the relative XUV-NIR delay of 0 fs (figure 3(a)) and 360 fs (figure 3(b)), where positive time delays indicate that the XUV pulse arrives before the NIR pulse. Based on the detected momentum of the photoelectrons, we can assign a specific ionization process to each feature. The outer feature at |p| 0.40 au (corresponding to a kinetic energy of 2.16 eV) is the result of direct ionization of the helium ground state |1s 2 by H17 (26.74 eV). We label it as H17. The feature at 0.63 au (5.40 eV) corresponds to ionization of the helium ground state by H19, which we label as H19. The angular distributions of these two features do not change appreciably with time delay because the NIR pulse is not directly involved in these ionization channels (i.e. they only change to the extent that ionization by H17 or H19 is implicated in other ionization mechanisms). Their angular distributions correspond to the emission of a P-wave, i.e. the formation of a continuum state with orbital angular momentum = 1. The sideband feature (labeled as S18) visible in figure 3(a) between these two ionization channels is the result of absorption (emission) of an NIR photon accompanying the absorption of an H17 (H19) photon. In figure 3(b) the sideband is not visible since there is no temporal overlap between the XUV and NIR pulses. Features at lower momenta in figures 3(a) and (b) correspond to ionization of the helium atom by H13 (20.45 eV) or H15 (23.60 eV) in combination with one or more NIR photons. Based on the measured photoelectron momenta of 0.073, 0.232 and 0.275 au (corresponding to photoelectron kinetic energies of 0.073, 0.73 and 1.03 eV, respectively), we assign these features to NIR ionization of n = 3, 4, 5 Rydberg states populated by the XUV pulse, where n is the principal quantum number of the helium Rydberg state |1sn . The PADs reflect the angular momentum of the outgoing photoelectrons and, as such, depend on the number of NIR photons implicated in the ionization processes. Four-lobed patterns at τ = 360 fs ( figure 3(b)), which include a local maximum in the angular distribution at 90 degrees with respect to the XUV + IR polarization axis, indicate the participation of Dwaves (i.e. photoelectrons with orbital angular momentum = 2). D-wave photoemission can arise in two-photon XUV-NIR ionization of helium: owing to the selection rules for dipole transitions, XUV-only ionization from the ground state |1s 2 can only create P-waves, while absorption of an additional NIR photon creates a superposition of S-and D-waves.

Experimental results
At zero relative delay τ = 0 (figure 3(a)), the PADs associated with ionization of n = 4, n = 5 are substantially different compared to the case of τ = 360 fs: an eight-fold pattern appears, with local maxima along the polarization axis as well as at angles of 45, 90 and 135 degrees with respect to the polarization axis, suggesting the involvement of G-waves ( = 4). Continuum states with this angular momentum can be reached by absorption of one XUV photon and three NIR photons. Since more complex ionization processes involving up to three NIR photons are observed at the overlap of the XUV and NIR pulses, in the following we will mainly focus on this case (τ = 0 fs).
To investigate the NIR wavelength dependence of the PADs associated with features n = 4 and n = 5 in the second row of figure 2, the NIR wavelength is varied between λ NIR = 785 nm and 796 nm while keeping the NIR field strength fixed at F NIR = 0.009 au (3.0 TW cm −2 ). Note that the NIR wavelength also determines the XUV photon energies via ω XUV = qω NIR , where q is the harmonic order. As the NIR wavelength is red shifted, the ionization features of the Rydberg states move toward lower momenta as expected. It is also clear that the NIR wavelength has a strong influence on the relative intensities of different features: at 785 nm the n = 5 feature is the most intense, while at 795 nm the n = 3 feature becomes the most prominent one. This is due to the red shift of the H15, which results in the increase of the excitation probability of the |1s3p state in comparison to the |1s4p and |1s5p states. To gain a further insight the angular distributions are analyzed in term of anisotropy parameters β K [21,22]. Each feature is integrated in momentum space in the narrow range around the n = 4 and n = 5 photoelectron peaks and the resulting angledependent distributions are fitted to the expression [21,22] where θ is the angle between the photoelectron emission direction and the laser polarization axis, σ 0 is the total cross-section, P K are Legendre polynomials and β K are the anisotropy parameters. The sum in equation (1) is restricted to a maximum of K = 8, because in the PADs we do not observe photoelectrons with angular momentum higher than = 4. Figure 4 reports the experimentally obtained anisotropy parameters β K associated with the features n = 4 (figure 4(a)) and n = 5 (figure 4(b)) as a function of the NIR wavelength. Note that K = 8 terms in equation (1) result purely from Gwaves (i.e. photoelectron continuum states with = 4). The other terms contain contributions from all -waves. Consequently, in order to investigate the relative strength between the G-waves and other contribution to the PADs, it is sufficient to focus on the NIR wavelength dependence of the anisotropy parameters β 8 . While for the Rydberg feature n = 4 (figure 4(a)) a slight increase of G-wave contributions is apparent for shorter wavelengths, the corresponding PADs are mainly dominated by D-wave contributions over the whole NIR wavelength region probed in the experiment, as indicated by the small values β 8 parameter. On the other hand, the PADs associated with the feature n = 5 show a clear wavelength dependence as the β 8 anisotropy parameter increases from 0.6 (at 796 nm) to 1.0 (at 785 nm) toward shorter NIR wavelengths. The observed difference in the NIR wavelength dependence of the PADs for n = 4 and n = 5 features indicate that the relative strength of the D-and G-wave contributions in the angle-resolved photoelectron spectrum results from intermediate Rydberg states resonant with the XUV and NIR frequencies used in the experiment.
Since the NIR field induces AC Stark shifts on all states, which results in changes of the resonance positions, we investigate the intensity dependence of the PADs of features n = 4 and n = 5. In figures 4(c) and (d) we report the NIR field strength F NIR dependence of the anisotropy parameters β K at a fixed NIR wavelength of λ NIR = 788 nm. For both features the β 8 parameter grows as the NIR field strength is increased. This indicates that at higher NIR field strengths the multi-photon transitions leading to the excitation of these high-angular momentum states = 4 become dominant over the one NIR photon transitions leading to excitation of = 2 states.

Numerical results
PADs obtained via numerical solution of the TDSE under the same conditions as the experimental PADs are shown in figures 3(c) and (d). Note that in figure 3(d) for numerical stability we choose the relative delay of 58 fs, which is smaller than the experimental delay of 360 fs. The delay in the experiment was chosen to ensure there is no overlap between the XUV and NIR pulses and avoid possible contamination by experimental NIR pre-pulses. In the calculations the relative time-delay between the XUV and NIR pulses, although smaller than in the experiment, is such that the XUV pulse fully arrives before the NIR pulse as in figure 3(b). All features observed in the experiment can be clearly identified in the numerical results. We note that the relative intensities of the n = 3 and n = 4 Rydberg features are not reproduced well in the TDSE simulation. This can be attributed to the model shape of the XUV spectrum used in the simulation. The spectrum used in the simulation does not perfectly reproduce the tail of the experimental spectrum, which strongly enhances the signal associated to the ionization of the n = 3 feature (|p| < 0.2 au) compared to the n = 4 and n = 5 transitions. This is due to the oscillator strength of the transition from the ground state to the |1s3p state, which is larger than the ones corresponding to transitions to the |1s4p and |1s5p states. As described above, the TDSE calculations allow us to track the residual population in the bound states after the end of both pulses as a function of the relative time delay τ , the carrier wavelength λ NIR (or, equivalently, the carrier frequency ω NIR ) and the NIR field strength F NIR . Figure 5(a) shows the residual population in selected Rydberg states calculated for a range of carrier NIR wavelengths λ NIR (760-820 nm) at a fixed NIR field strength of F NIR = 0.004 au (0.56 TW cm −2 ) and for a relative time delay of τ = 0 fs. Note that the NIR carrier wavelength determines the photon energies of the high-harmonics peaks in both the experiment and the TDSE simulations. For each Rydberg state the residual populations peaks at a certain wavelength, indicating that a resonant transition is involved. For the |1s3p , |1s4p and |1s5p Rydberg states the residual populations peak respectively at NIR wavelengths of 805, 783 and 773 nm: this is due to the energy of H15 becoming resonant with the energy of the respective Rydberg state 15ω NIR = E np . The wavelength dependence of the residual populations in the |1snf states requires a more detailed analysis. Since they are not directly coupled to the ground state, transitions to the |1snf states require the absorption of more than one photon. The resonance conditions for these high angular momentum states should be evaluated in terms of a combined XUV and NIR multi-photon transition. In reference [5], it was suggested that the |1snf Rydberg states are populated via the |1s2p state by the absorption of H13 and two NIR photons. To test this assumption, we performed TDSE calculations for XUV pulses containing only H13 or only H15. In figures 5(b) and (c) we report the residual populations in the Rydberg states |1s4f and |1s5f for XUV pulses composed of only H15 (dark shaded areas) or only H13 (lighter shaded areas). The residual population obtained for XUV pulses containing all four harmonic orders are also shown for comparison as dashed lines. According to figures 5(b) and (c), the residual populations in the |1snf states are reproduced when the XUV pulse is only composed of H13 and are lower by almost two orders of magnitude when the XUV pulse is only composed of H15. Therefore, figure 5 confirms that the population transfer from the ground state to the |1snf states occurs predominantly by means of a multi-photon absorption process where the absorption of H13 is accompanied by the absorption of two NIR photons.
The NIR field induces AC Stark shifts on all states, resulting in changes of the resonance positions and the strengths of the transitions. To investigate the role of AC Stark shifts in the two-color dynamics, the residual populations of the Rydberg states are studied as a function of the NIR intensity. This is done by fixing the NIR wavelength at λ NIR = 788 nm-where both the |1snp and |1snf states are efficiently populated in the presence of the NIR field as shown in figure 5(a)-and varying the intensity of the NIR field from zero up to 4 TW cm −2 , for a relative delay of τ = 0 fs. Figure 6(a) shows the residual population in the same states as in figure 5(a) as a function of the NIR field strength F NIR . At low NIR field strengths below 0.002 au (0.14 TW cm −2 ), the residual population at the end of the pulses is limited to the states that are optically coupled to the ground state via the absorption of one XUV photon, i.e. the |1snp states. Since in the one-photon regime the ionization rate of the bound states scales quadratically with the NIR field strength, the residual populations of the |1snp states in figure 6(a) fall off quadratically with the field strength F NIR up to approximately F NIR = 0.005 au (8.8 TW cm −2 ). As the NIR field strength is increased, a qualitative change of the residual populations in the Rydberg states occurs. The residual populations in the |1s4f and |1s5f states becomes larger than in the |1s5p and |1s4p states for field strengths higher than F NIR = 0.005 au (0.88 TW cm −2 ) and F NIR = 0.007 au (1.72 TW cm −2 ), respectively. This change in the dominant Rydberg character is reflected in the experimental data of figures 4(c) and (d), where we observe that G-waves dominate over D-waves only above an NIR field strength of 0.01 au ( 3.51 TW cm −2 ), which is in agreement with the previous literature [5][6][7][8]. Above F NIR = 0.01 au, the residual populations of all Rydberg states drop significantly as the NIR pulse increasingly depopulates the Rydberg states by ionization.

Effective restricted basis model
In order to obtain more qualitative insight into the dynamics of two-color XUV-NIR ionization of helium, an effec- tive restricted basis model is constructed. The model is built by expanding the wavefunction |Ψ(t) in terms of a restricted basis which includes the bound states and the continua of the helium atom that can be populated by one XUV photon and up to three NIR photons, i.e. including states with orbital angular momenta up to max = 4. The experimental and numerical PADs show features associated with Rydberg states with angular momentum = 1 and = 3 and a principal quantum number of n 5. Therefore, we restrict the expansion of the wavefunction |Ψ(t) in terms of the |1snp and |1snf states to n 5. The PADs do not show any features associated with |1sns and |1snd states, but these states serve as off-resonant intermediates. Therefore we expand the wavefunction on these series without any restriction on the principal quantum number. Hence, the following ansatz for the wavefunction is introduced: (2) where the system is initially in the ground state, i.e. 1s 2 |Ψ(−∞) = 1. c n (t) and c (t) are the time-dependent amplitudes of the bound states and the continua, respectively. For shortness of notation we indicate the amplitudes of the bound states only by the quantum numbers of the excited electron c n (t) and continuum states with energy by c (t), assuming that the second electron always occupies the 1s shell. This wavefunction is inserted in the TDSE i∂ t |Ψ(t) = H(t)|Ψ(t) , where the time-dependent Hamiltonian in the length gauge is given byĤ(t) =Ĥ 0 −d · E(t, τ ). HereĤ 0 describes the field-free atom,d is the dipole operator and E(t, τ ) is the two-color field, with relative delay τ between the XUV and NIR pulses. The XUV and NIR pulses, are assumed to have Gaussian envelopes with full-width at halfmaximum (FWHM) of 5 fs and 50 fs and, like in the numerical TDSE simulations of section 3, the XUV pulse is composed of a sum of equally intense (3.5 × 10 10 W cm −2 ) harmonics qω NIR with q = 13, 15, 17, 19. Both pulses are linearlypolarized and their polarizations are taken to be parallel. The transitions considered in the model are restricted by the dipole selection rules for linearly polarized light Δ = ±1, Δm = 0. For the XUV field we consider couplings of the ground state |1s 2 to the |1snp series by H13 and H15 as well as to the |1s p continuum by H17 and H19. In figure 7 we depict the relevant states and transitions considered in the effective model.
It is possible to use simple considerations in order to exclude some of the NIR-induced transitions from the model. Due to the relatively weak NIR intensity used in the experiment, multi-photon NIR-only transitions from the ground state and the AC Stark effect of the ground state are ignored. Hence the NIR field influences only the Rydberg states and the continua. Moreover, we exclude NIR-induced continuumcontinuum transitions in our treatment due to the small probability of an electron absorbing an NIR photon when far away from the ion core.
At the intermediate intensities considered in this work the NIR field induces substantial AC Stark shifts on both the Rydberg and continuum states. AC Stark shifts of Rydberg states in helium were investigated most recently in reference [23]. Following reference [23], we model the AC Stark shift as a ponderomotive shift U p (t) = |F NIR (t)| 2 /4ω 2 NIR , except for the |1s2p states. The AC Stark shift of the |1s2p state is known to have a complex NIR intensity dependence [23]. We extract its intensity-dependent AC Stark shift from the numerical results of reference [23]. The insertion of the ansatz (2) in the TDSE leads to a system of coupled differential equations for the bound c n (t) and continuum c (t) amplitudes. Note that in the ansatz (2) the expansion of the wavefunction |Ψ(t) on the |1sns and |1snd states involves two infinite sums (since there is no restriction on the principal quantum number n); this in principle would lead to an infinite number of coupled equations, which would make the model numerically untreatable. In order to avoid treating the two Rydberg series explicitly, it is possible to use an adiabatic elimination procedure [24][25][26]. The procedure is based on the fact that the NIR-induced transitions leading to excitation of |1sns and |1snd Rydberg states are off-resonant (i.e. the energy separation between the Rydberg states of these and the neighboring series are much larger/smaller than one NIR photon). Consequently, their amplitudes (c ns (t) and c nd (t)) oscillate quickly compared to the ones of the neighboring Rydberg series (c np (t) and c nf (t)); in this procedure the approximationċ ns/d (t) 0 is then made and the resulting adiabatic expression for the discrete amplitudes is substituted in the remaining coupled differential equations. The explicit treatment of the amplitudes for the continuum states c (t) would also make the model computationally expensive, due to the appearance of integrals over the photoelectron energy in the system of coupled equations. By taking the approximation of a flat continuum, i.e. a slow dependence of the dipole matrix elements between bound and continuum states 1sn |d|1s on the photoelectron energy , it is possible to remove the continuum amplitudes from the coupled equations as well [24][25][26]. We refer the reader to appendix A for a more detailed description of the adiabatic elimination procedure for the off-resonant and continua states.
By applying the adiabatic elimination of the off-resonant and continuum states and using the rotating wave approximation (RWA), a system of coupled differential equations for the amplitudes of the ground state |1s 2 and the |1snf and |1snp Rydberg states is obtained. Denoting the vector of amplitudes as c(t) = {c 1s 2 (t), c 2p (t), c 3p (t), c 4p (t), c 5p (t), c 4f (t), c 5f (t)}, the system of differential equations may be cast in a matrix form: here H 0 is a diagonal matrix of the field-free energies of the discrete states, which are taken from the NIST database of atomic energy levels [27]. U p (t) is a Stark shift matrix, which we model as explained previously; Ω(t) accounts for the one-photon dipole coupling of the ground state to the |1snp series; Γ(t) is a complex-valued matrix which appears in the adiabatic elimination procedure (see appendix A), which accounts for one-photon coupling of the Rydberg states to the continua, as well as two-photon couplings of the Rydberg states to each other via V-type and Λ-type transitions; finally, V(t) takes into account the coupling of the |1s2p state to the |1snp and |1snf states via absorption of two NIR photons. The matrix elements of the matrices mentioned above are given by here we used the short notation d ik = i|d|k and d i = i|d|1s for the dipole matrix elements. The formal expression for the matrix element V 2p,n (t) is obtained by secondorder perturbation theory. As explained in more detail in appendix B, we manually adjust the strength of the matrix elements of V(t) in order to reproduce the intensity-dependent behavior of the residual populations obtained by numerical solution of the TDSE in figure 6(a). The strength of the matrix elements is reported in table B1 of appendix B. The diagonal elements of the matrix Γ(t) are the decay rates of the bound states toward the continua while the offdiagonal elements correspond to the two NIR-photon coupling of two Rydberg states with Δ = ±2 via the continuum and via off-resonant bound states reachable by the selection rules for linearly polarized light. The parameter q ij is the Fano parameter [26,28], which governs the ratio between the transition strength of the two-photon coupling via the offresonant bound states (V-type transition) and the continuum states (Λ-type transition). The Fano parameters are adjusted according to the matrix elements of the matrix V(t) (see appendix B for more detail) and are reported in table B1 of appendix B. All remaining dipole matrix elements are calculated using the central potential model for helium as in [28,29]. In this approach, the excited state wavefunction is taken as hydrogen-like with effective quantum numbers n, aimed at reproducing the field-free energy of the Rydberg state. This allows for analytical evaluation of all dipole matrix elements. As shown in appendix C, the central potential method is in excellent agreement with the existing literature [30,31] for the range of bound and continua states considered in this paper.
After the adiabatic elimination of continua, the amplitudes c (t) are no longer directly determined. After numerical integration of the differential equations (3), they can be calculated by inserting the resulting discrete amplitudes c n (t) in the formal solution for the continuum equations in the RWA where ω and ω n are the kinetic energy of the electron in the continuum and the energy of the Rydberg state,s respectively. This expression takes into account the timedependent ponderomotive shift of the continua due to the NIR field via φ(t )= t t 0 dτ U p (τ ), where, as before, U p (t) = |F NIR (t)| 2 /4ω 2 NIR . The energy-resolved photoelectron distributions are given by the sum over all angular momenta of the squared amplitudes of the continua at large detection times P = |c (t → +∞)| 2 . The PADs can be found by expressing the outgoing photoelectron wavefunction in partial waves corresponding to the different continua [32]: where θ is the angle between the momentum of the photoelectron and the laser polarization axis, δ = arg[Γ(1 + − i/k)] (where Γ is the Gamma function) is the Coulomb phase shift of the th partial wave with momentum k = √ 2 and P is the Legendre polynomial of order . In order to facilitate the comparison with the PADs obtained in the experiment and the TDSE calculations, we produce PADs resolved in the momentum of the photoelectron f(k, θ) = (d /dk)f( , θ) in figures 3(e) and (f ).

Discussion
The photoelectron angular distributions obtained from the effective restricted basis model for an NIR field intensity of 3.51 TW cm −2 and a carrier wavelength λ NIR = 788 nm are compared to the experimental and TDSE results in figures 3(e) and (f ) for a relative delay of τ = 0 fs and τ = 58 fs, respectively, where, like in the TDSE calculations, the positive time-delay is chosen such that the XUV pulse fully arrives before the NIR pulse. The PADs clearly reproduce the features observed in the experimental and numerical results. The residual populations in the Rydberg states |c n (∞)| 2 for a relative XUV-NIR time delay τ = 0 fs are shown in figure 5(d) as a function of the carrier wavelength λ NIR for a fixed field strength F NIR = 0.004 au. They reproduce the TDSE results in figure 5(a) qualitatively and quantitatively, which confirms that the model is able to capture the essential transitions that are responsible for the combined XUV-NIR ionization. Small deviations are observed at longer NIR wavelengths above λ NIR = 800 nm. TDSE calculations (not reported here) show that in this region other Rydberg states such as the |1s5s state are also populated, which were removed in the model by the adiabatic elimination procedure. This indicates that the assumption for the adiabatic elimination, i.e. the off-resonant nature of transitions toward the |1s5s state, breaks down at these longer wavelengths. In figures 5(e) and (f ), we report the residual population in the |1snf states for an XUV pulse composed of harmonic H13 (light shaded areas) or H15 (darker shaded area) only. In the case of an XUV pulse consisting of H13, the residual populations for the |1snf states are much larger, confirming that these states are populated via the |1s2p state. While in both cases the restricted basis model overestimates the contributions from V-type and Λ−type transitions to the residual populations in the |1snf states for H15, the deviations are larger for the |1s4f state. This might be a result of the adjustment of the strength of the twophoton couplings in the model, which aimed at reproducing the intensity dependence of the residual populations at the fixed wavelength λ NIR = 788 nm. Figure 6(b) shows the residual populations resulting from the model as a function of the field strength F NIR for a carrier wavelength of λ NIR = 788 nm and for a relative delay of τ = 0 fs. The main features observed in the TDSE results are reproduced, although small deviations persist especially at higher field strengths. Similar to the TDSE results, the residual population in the |1snf states becomes larger than the population in the |1s5p state at a field strength of ca. F NIR = 0.005 au (0.9 TW cm −2 ); this indicates that, above this NIR field strength, the pathway leading to the |1snf states via the |1s2p state dominates over the one-photon excitation of the |1snp states by H15. At even larger field strengths F NIR > 0.01 au (3.51 TW cm −2 ) the model predicts larger residual populations in the Rydberg states compared to the TDSE results. This is expected, as at these larger field strengths the Keldysh parameter γ = ω √ 2E b /F NIR (where E b is the binding energy) approaches unity for the Rydberg states, indicating that the NIR field cannot be treated perturbatively. As a consequence, the model underestimates the ionization rate of the Rydberg states. By using the model, we can explain why the appearance of G-waves occurs predominantly during the overlap between the XUV and NIR pulses. In figure 8, we report the timedependent populations |c n (t)| 2 of the Rydberg states for the case of a relative delay between the XUV and NIR pulses of τ = 0 fs ( figure 8(a)) and τ = 58 fs ( figure 8(b)) for an NIR field with a field strength of F NIR = 0.01 au and a wavelength λ NIR = 788 nm. The populations in the |1snp states display fast oscillations during the XUV pulse, where each of the oscillations correspond to one of the attosecond pulses composing the XUV attosecond pulse train. The amplitude of the oscillations depends on two factors: the absorption crosssection of the |1snp Rydberg states (which scales with respect to the principal quantum number of the bound state as n −3 ) and the detuning between the harmonic frequencies and the energy of the Rydberg states in the presence of the NIR field. The n −3 dependence of the absorption cross-section of the |1snp Rydberg states results in a smaller population transfer to the higher lying states with respect to the lower lying ones. The effect of the relative detuning of the harmonic fre- quencies to the Rydberg state energies may be instead appreciated by looking at the relative amplitude of the oscillations of each |1snp state in figure 8(a): the population in the |1s4p state displays an almost smooth time-dependence due to its energy being almost resonant with H15, while the population of the other |1snp states varies sharply at each oscillation around t = 0 fs.
In figure 8(a) at the end of the NIR pulse the population is largely in the |1snf states while the |1snp states are almost completely depleted. For the n = 3, 4, 5 members of the |1snp Rydberg series the depletion is largely due to one-photon ionization toward the |1s s and |1s d continua. The time-dependent behavior of the |1s2p state at τ = 0 fs requires a more detailed analysis. Since the H13 is off-resonant with respect to the |1s 2 − |1s2p transition, the population in the |1s2p state is present only transiently during the XUV pulse. Therefore, a limited time window is available for the NIR field to transfer the transient |1s2p population to the high-lying |1snf states via the absorption of two NIR photons. In figure 8(a), two-photon Rabi oscillations with a period of 14 fs between the |1s2p (blue line) and |1snf (magenta and orange lines) populations can be observed on the falling edge of the NIR field. Similar Rabi oscillations between these states were also observed in the literature [23,33]. When the two pulses are not overlapping in time, the excitation of |1snf states and the subsequent one-photon ionization toward the |1s g continuum is much less efficient. In figure 8(b) we show the time-dependent populations during the NIR pulse for τ = 58 fs. Prior to the onset of the NIR pulse the XUV pulse induces substantial population in all |1snp states (notice that the population of the |1s3p and |1s4p states even goes offscale in figure 8(b)). On the rising edge of the NIR pulse, one-photon processes dominate and a large drop of population in the |1snp states is observed due to ionization toward the S-and D-continua. Around the peak of the NIR pulse the two-photon coupling between |1snp and |1snf states becomes appreciable. In the insets of figure 8, we show the photoelectron angular distributions associated with the ionization of n = 5 Rydberg states, obtained by integrating the PADs within an interval of 0.30 eV around the photoelectron peak corresponding to one-photon ionization of n = 5 Rydberg states; the switching between G-and D-wave is evident as the NIR-XUV delay is varied between τ = 0 fs and τ = 58 fs.
To complete the analysis we investigated the effect of the XUV pulse chirp on the resulting PADs. We introduced positive and negative chirp leading to approximately twice longer XUV pulses as compared to the unchirped pulses used in the main calculations. The chirp of the XUV pulse does not change the photoelectron distributions appreciably, apart from small variations of the relative intensities of Rydberg-mediated ionization channels due to the redistribution of the XUV pulse intensity under the envelope of the NIR probe pulse.
Let us now turn to the question of why the off-resonant pathway to the |1snf states via the |1s2p state (accompanied by the absorption of an H13 photon) dominates over other pathways coupling the |1s4p and |1s5p states (produced by the absorption of an H15 photon) to the |1snf states by Λ−type transitions and V-type transitions involving the same number of NIR and XUV photons (see figure 7). For the explanation, it is sufficient to invoke the Fano propensity rules [34]. The propensity rules state that absorption (emission) of a photon tends to increase (decrease) the orbital angular momentum . That is, in absorption (emission) the dipole matrix element 1sn |d|1sn is larger (smaller) when = + 1 ( = − 1). These propensity rules in atomic transitions are valid for bound-bound transitions as well as for bound-continuum and continuum-continuum transitions [34,35]. The Λ− and V-type transitions violate the propensity rules and, therefore, the associated transition probabilities are strongly suppressed with respect to the H13 + 2ω pathway. To show this quantitatively, we report in table 1 the radial parts of the dipole matrix elements calculated using the effective hydrogen-like helium potential mentioned above and described in detail in appendix C for the Λ− and Vtype pathways involving the |1s4p and |1s5p states and the 2p + 2ω IR pathway from the |1s2p to the |1snf states via the |1s3d state. The second, third and fourth columns of table 1 show the absolute value of the radial dipole matrix elements between the initial state and the |1s3d state, between the |1s3d state and the |1s4f state and between the |1s3d state and the |1s5f state, respectively. The last two columns of table 1 correspond to the products of the matrix ele- Table 1. Absolute values of the radial part of the dipole matrix element (in atomic units) of the pathways coupling the |1snp and |1snf states, as obtained in the effective restricted basis model. The Λ and V symbols indicate, respectively, the pathways via the continuum |1s d and via the |1s3d state (see figure 7). For the calculation of the bound-continuum matrix elements, we have used an NIR wavelength λ NIR = 788 nm.
Pathways ments of each dipole transition for the V-, Λand 2p + 2ω NIR transitions and are a good indicator for the relative strength of the different pathways. The matrix element for the 2p + 2ω NIR pathway is at least four times stronger than that of all other pathways. Moreover, recent work [23] has clearly shown that the AC Stark shift of the |1s2p state due to the NIR field-for intensities similar to the ones used here-tends to lower the energy of the Rydberg state, shifting it toward resonance with H13. The transition probability is therefore enhanced also by the AC Stark shift of the |1s2p state, which is included in the present model.

Conclusions
In this work we have systematically studied the time-resolved NIR-assisted XUV photoionization of atomic helium at moderate NIR intensities up to several TW cm −2 for a range of NIR wavelengths. We have employed velocity map imaging spectroscopy to record the NIR wavelength and intensity dependence of the angle-and energy-resolved momentum distribution of the photoelectrons created in two color XUV-NIR photoionization for different time delays between the XUV and NIR pulses. The results show a dependence of the PADs on the relative delay between the XUV and NIR pulses; at large delays when the XUV comes first, onephoton ionization of the |1s3p , |1s4p and |1s5p Rydberg states mainly contributes to the photoelectron spectra and D-waves are observed in the PADs. At zero delay when the two pulses are overlapping in time, the angular distributions change dramatically and point to the population of the |1s4f and |1s5f states and their one-photon ionization toward the |1s g continuum. The G-wave contributions to the PADs are enhanced for shorter NIR wavelengths and higher NIR intensities. where c n (t), c (t) are the time-dependent amplitudes of the bound and continuum states, ω n are the field-free energies of the Rydberg states (in atomic units) and is the continuum energy. The sum over the bound states principal quantum number is taken over n = + 1, . . . , ∞ and we restrict the orbital angular momentum to the maximum value of = 4, i.e. we consider the bound states up to the |1sng series and continua states up to |1s g . The initial condition is that the system is in the ground state |Ψ(−∞) = |1s 2 . In the above ansatz, we have set the energy of the ground state as the zero energy, i.e. ω 1s 2 = 0. The upper limit of the continuum energy is set at +∞, since for the photon frequencies hereby considered resonances such as doubly excited states lie well above the energies that are reached in the experiment. Inserting the ansatz in the TDSE leads to a system of coupled integro-differential equations for the discrete and continuum amplitudes. In order to simplify the description, we make use of the adiabatic elimination of the continua as described in [24][25][26]; the latter consists in formally integrating the equations for the continua and substituting the resulting expressions into the equations for the discrete states. Making use of the rotating wave approximation together with the assumption that the continuum is flat, that is the bound-free dipole matrix element 1sn |d|1s is a slow function of the continuum energy = ω n + ω NIR in the range spanned by the bandwidth of the laser field with frequency ω NIR , it is possible to define an effective decay rate and AC Stark shift of the bound state due to the coupling to the continuum: where P indicates the Cauchy principal value of the integral and ω i is the field-free energy of the Rydberg state in atomic units. In the numerical simulations, we approximate the Stark shift S i (t) by the ponderomotive shift U p (t) = |F NIR (t)| 2 /4ω 2 NIR for all Rydberg states except for the |1s2p state for which we extract the field dependence of the Stark shift from [23]. Note that the adiabatic elimination of the continua also leads to off-diagonal elements describing the coupling of the Rydberg states via the common continua (equation (9) in the main text). For the off-resonant bound states, we proceed analogously by noticing that when the detuning between two bound states Δ ij = E i − E j + ω NIR is substantially larger than the bound-bound dipole matrix element i|d| j , one can approximateċ i (t) 0. This results in the adiabatic expression for the c i (t) amplitude which may then be substituted in the remaining differential equations. The resulting system of equations is the one in equations (3)-(12) of the main text.

Appendix B. Matrix elements of V(t) and Fano parameters q ij
The time-dependent matrix H eff (t) describing the system of coupled differential equations (3) was briefly described in the main part of this article. In this appendix we describe the adjustment of the matrix elements of the matrix V(t) mentioned in the main text as well as the Fano parameters q ij appearing in the expression for the matrix elements of the matrix Γ(t). The expression for the matrix elements of V(t) is given in second-order perturbation theory and rotating waveapproximation by where the index i runs over the members of the |1sns and |1snd series for transitions to |1snp final states and |1snd series for transitions to |1snf final states. In the above expressions τ is the relative delay between XUV and NIR. Due to the poor convergence of the sum i resulting from perturbation theory, we take these matrix elements as fit parameters in our calculations, with the aim of reproducing the contribution of G-waves in the experiment and TDSE calculations of section 4. We vary the strength between the V 2p,4f and V 2p,5f matrix elements while keeping the ratio of the matrix elements toward final |1snf and |1snp states fixed at V 2p,nf /V 2p,np = 10. The order of magnitude roughly reproduces the ratio between matrix elements from the |1s3d state toward final |1snp and |1snf states (see the second and third columns of table 1). We find optimal agreement for V 2p5f = 60 au and V 2p4f = 30 au. The matrix elements are reported in table B1. The parameter q ij (equation (10) in the main text) is the Fano parameter [26,28], which gives the ratio between the transition strength of the two-photon coupling via the offresonant bound states (V-type transition) and the continuum states (Λ-type transition). From the above two-photon coupling matrix elements V 2p,n between the |1s2p and upper lying |1snf and |1snp states for n = 4, 5, we estimate the V-type couplings between |1snp and |1snf states by restricting the sum in equation (18) to the |1s3d state, since all other |1snd members of the Rydberg series are far away from resonance with the considered transitions

Appendix C. Effective potential for He
We follow the approach described in [28,29] and assume a central potential model, where the excited electron wavefunction is taken as a hydrogen-like wavefunction. The effective quantum numbers in the model potential are different than the ones obtained in the typical quantum defect theory (QDT) as, in our case, the radial quantum number is forced to have integer values while the orbital angular momentum can assume non-integer values. We refer the reader to the literature for further details on the adopted model [29]. For the relevant states in He considered in this article, the quantum defects quickly become small for increasing principal quantum number n and orbital angular number . We take the field-free energies E n for the Rydberg states and the ionization potential I p from the NIST database and find the effective principal, orbital and radial quantum number from the expressions The hydrogenic wavefunctions describing the Rydberg and continua states are then defined by the effective quantum numbers. When considering the radial part of the dipole matrix element for bound-bound and bound-free transitions, this approach allows one to find an analytical expressions in terms of the Appell hypergeometric function F 2 [28].
To test the quality of the approximation by the effective wavefunctions, we compare them with the literature in figure C1: the agreement is excellent. As mentioned before, since the quantum defects decrease for increasing n and , we expect an even better agreement for the |1snf and |1s states with = 0, 2, 4. Furthermore, we find excellent agreement with the partial photoionization cross-sections for the |1s2p and |1s3p states toward the |1s s and |1s d continua found in literature [31] and the transition strength for bound-bound transitions between the first few Rydberg states of the |1sns and |1snp series when these effective