Polarization in Strong-Field Ionization of Excited Helium

We analyze how bound-state excitation, electron exchange and the residual binding potential influence above-threshold ionization (ATI) in Helium prepared in an excited $p$ state, oriented parallel and perpendicular to a linearly polarized mid-IR field. Using ab initio B-spline Algebraic Diagrammatic Construction (ADC), and several one-electron methods with effective potentials, including the Schr\"odinger solver Qprop, modified versions of the Strong-Field Approximation and the Coulomb-Quantum Orbit Strong-Field Approximation (CQSFA), we find that these specific physical mechanisms leave significant imprints in ATI spectra and photoelectron momentum distributions. Examples are changes of up to two orders of magnitude in the high-energy photoelectron region, and ramp-like structures that can be traced back to Coulomb-distorted trajectories. The present work also shows that electron exchange renders rescattering less effective, causing suppressions in the ATI plateau. Due to the long-range potential, the electron continuum dynamics are no longer confined to the polarization axis, in contrast to the predictions of traditional approaches. Thus, one may in principle probe excited-state configurations perpendicular to the driving-field polarization without the need for orthogonally polarized fields.


I. INTRODUCTION
Above-threshold ionization (ATI) is a strong-field phenomenon in which an atom absorbs more photons than are energetically required for it to ionize. Since its first observation [1], ATI has been employed to provide unprecedented insight into the interaction of intense laser radiation with atoms or molecules. The high-energy photoelectrons released into the continuum, together with subfemtosecond resolution, make ATI a powerful tool for probing real-time electron dynamics. Applications of ATI directly associated with imaging and target reconstruction are, for instance, laser-induced electron diffraction (LIED) [2] and ultrafast photoelectron holography [3]. Such imaging methods combine high photoelectron currents and the possibility of retrieving phase differences without the need for elaborate interferometric schemes. See, e.g., [4] for a review.
For that reason, orbit-based models have been used for over two decades due to their huge predictive power. By relating quantum transition amplitudes to electron paths, which may interfere for a given photoelectron energy, they provide a physical picture of ATI as a laser-induced rescattering process ( [5][6][7]; for a review see [8]). Attosecond resolution is guaranteed in the rescattering process and ensuing phenomena because it takes place within a fraction of a field cycle, which, typically, is a few femtoseconds. Established approaches such as the strong-field approximation (SFA) [9][10][11] consider the target as a source term and employ a structureless continuum, i.e., field-dressed plane waves. This has allowed the transition amplitudes associated to strong-field phenomena to be written as a laser-dressed Born-type series [5][6][7], which considerably reduces the numerical effort involved. Besides the easy implementation, an obvious advantage is that the interaction with the core is well-defined. This makes rescattering processes straightforward to identify. Furthermore, the SFA can be easily associated with interfering electron orbits, if formulated in conjunction with saddle-point methods [12][13][14].
For that reason, the insight provided by the SFA was vital for identifying the overall structure and interpreting key features in ATI spectra. The ATI spectrum consists of a low-energy region, up to 2U p , where U p is the ponderomotive energy, and of a high energy region extending to approximately 10U p . The 2U p cutoff corresponds to the 'classical' kinetic energy limit for 'direct' ATI electrons, which are freed into the continuum via tunnel or multiphoton ionization, and reach the detector without further interaction with the target. The ionization mechanism that predominates depends on the Keldysh parameter for the system in question, γ = I p /2U p , where I p is the ionization potential. A Keldysh parameter of greater than one implies multiphoton ionization, while tunnelling is dominant if γ ≤ 1 [9].
The high-energy region is mainly dominated by high-order ATI (HATI), which is characterized by a plateau that does not significantly change across the energy distribution, falling off at the cutoff point of approximately 10U p . The widely accepted explanation of this feature is that the high-energy photoelectrons result from the ionized electron being driven back by the laser field to scatter elastically off its parent ion. The rescattered electron is then accelerated further in the field to energies up to 10U p , which can be derived from a classical model [15]. Alternatively, the freed electron may recombine with a bound state of the parent ion, which leads to high-order harmonic generation (HHG), or recollide inelastically with the core, releasing other electrons. The latter processes give rise to laser-induced nonsequential double and multiple ionization (NSDI, NSMI) [16,17].
Nonetheless, the aforementioned picture is an oversimplification. For example, ATI peaks have a substructuretermed Freeman resonances -caused by the ponderomotive shifts of states that produces resonant enhancements [18]. This effect can also cause broadening of the peaks and even large energy shifts. Furthermore, recent discoveries have cast doubt on the distinction between tunneling and multiphoton regime based on the Keldysh parameter -for instance, photoelectron spectra obtained from solving the time-dependent Schrödinger equation (TDSE) show features consistent with multiphoton absorption even for laser intensities that correspond to the tunnelling ionization regime [19]. In addition to that, specific features of the plateau can vary from system to system. Some atomic targets, such as krypton, have a dropping plateau while others, including xenon, have a flatter plateau. Specific laser intensities can also cause resonance-like enhancements to the spectra. A variety of theoretical models have been proposed in recent years to account for these surprising observations, including Floquet theory [20], channel-closing theory [21,22] and models based on the analysis of TDSE solutions [23].
Moreover, dynamic effects such as charge migration [24][25][26][27][28], bond breaking [29], nuclear degrees of freedom [30] and multielectron resonances [31,32] are expected to be important for larger systems. Extended targets also imply a structured continuum, or a blurred bound-continuum distinction. Modeling such dynamics is not an easy task even for the simplest possible case, i.e., a single-electron continuum under the joint influence of the laser field and a long-range potential. The residual potential blurs the distinction between "direct" and "rescattered" as the electron now may be lightly deflected, or undergo a soft scattering by the core [33]. This goes beyond the clear-cut distinction imposed by Born-type approaches such as the SFA [34,35]. All the above suggests the driving-field polarization as a tool to assess what the SFA leaves out. In fact, experimental and theoretical studies have been carried out that demonstrate the dependence of recollision outcomes -such as HATI or HHG -on the polarization of the incident radiation [36,37]. These studies show that the yield is increasingly suppressed as the ellipticity of the laser pulse is increased. The generally accepted explanation for this observation is that an elliptical laser will imbue the ionized electron with a non-zero drift velocity that reduces the chances of recollision with the parent ion by the electron [38]. Subsequently, it was demonstrated that substituting laser polarization for atomic polarization yielded the same order of magnitude suppression in the production of HHG [39]. The present work will extend this line of inquiry to the dependence of the production of HATI on the atomic polarization. We will focus for simplicity on the helium atom prepared in an excited state of ns 1 (n + 1)p 1 interacting with a linearly polarized field. The polarization of the p orbital can then be oriented along the laser polarization axis or in the plane perpendicular to it.
In case the atomic polarization is oriented parallel to the laser polarization, we expect there to be no suppression to the process outlined above, wherein a continuum electron is driven by the laser field to scatter off the parent ion. However, if the atomic polarization is oriented perpendicular to the laser polarization, ionization in the direction of the laser polarization is precluded due to the atomic orbital symmetry. Therefore, ionization is only possible if the electron acquires a nonzero drift velocity in the plane perpendicular to the laser polarization. This should greatly reduce the probability of a successful recollision with the core, as the electron's drift velocity, unimpeded by interaction with the laser, will displace it away from the parent ion. Thus, in the high intensity, low frequency regime for which the recollision model is applicable, we expect the atomic polarization to significantly affect the HATI spectrum. On the other hand, residual central potentials may favor ionization or rescattering off the polarization axis, which deviates from the standard predictions. For instance, in [40], rescattering in circularly polarized fields may occur due to the interplay of the laser field and the long-range potential.
In this work, we consider excited helium parallel and perpendicular to the driving-field polarization, for which we calculate photoelectron spectra and momentum distributions. We use approaches that incorporate the residual binding potential and/or the core dynamics: the B spline algebraic diagrammatic construction (B-spline TD-ADC), the one-electron Schrödinger solver Qprop [41,42], the SFA and the Coulomb quantum-orbit strong-field approximation (CQSFA). They will allow for an assessment of what the simple rescattering picture leaves out. The B-Spline TD-ADC is an ab-initio method that solves the 3D atomic many electron time-dependent Schrödinger equation for a neutral system [43]. Specifically for the two-electron model studied in this work, it accounts for electron exchange and recollision-induced excitations. In the strong-field context, the B-Spline ADC has been used in high-order harmonic generation [44] and pump-probe spectroscopy [45]. The CQSFA is an orbit-based method that incorporates the Coulomb potential and the external laser field on equal footing. It has a high predictive power as it allows for quantum mechanical pathways to be switched on and off at will. In the context of ultrafast photoelectron holography, it led to possibly the most thorough investigations of how holographic patterns form [33,[46][47][48]. This work is organized as follows. In Sec. II, we discuss the theoretical methods used in this work. Subsequently, in Sec. III, we present our results for initial perpendicular and parallel polarized states. Specifically, we investigate angle-integrated spectra (Sec. III), and photoelectron momentum distributions (Sec. III B). Finally, our main conclusions are presented in Sec. IV.

II. THEORY
In order to calculate ATI spectra or photoelectron-momentum distributions, we must compute the time-dependent wave function and project it onto an asymptotic continuum state with a well-defined momentum. Below, we will discuss several ways to do so. We will use atomic units throughout, unless otherwise stated, and the dipole approximation.

A. B-Spline algebraic diagrammatic construction
We calculated the ATI spectra of excited, polarised, helium atoms using the time-dependent (TD) B-spline algebraic diagrammatic construction (ADC) ab initio method [43-45, 49, 50]. Within the TD B-spline ADC approach, the 3D many-electron time-dependent Schrödinger equations (TDSE) for the neutral helium atom interacting with the intense mid-IR laser field, given by is solved by making the ansatz for the time-dependent many-electron state of neutral helium. In where Ψ N GS represents the ground state of neutral helium, while the basis states Ψ N j refer to the correlated many-electron configuration states of the ADC theory for N -electron neutral systems [51].
The total time-dependent Hamiltonian in Eq. (1) for the time-evolution of the system interacting with the pulse readsĤ HereĤ N 0 is the field-free many-electron ADC Hamiltonian describing the neutral system. The laser-atom interaction is described in the length gauge asD z E(t), where the z component of the dipole operator isD z = N j=1ẑ j and the summation over the j index runs over all the N electrons of the atom.
The single-particle basis set used in this approach consists of spherical harmonics Y l,m (θ, ϕ) for the angular part and B-spline functions B i (r ) for the radial coordinate. The single particle basis functions used in this calculation are therefore expressed as The initial excited states used in the present simulations are |Ψ N (t = 0) ⊥ = |1s2p x and |Ψ N (t = 0) = |1s2p z in the perpendicular and parallel set-ups, respectively. In this work, we have used the lowest level of the ADC-hierarchy compatible with a correct description of the strong-field ionization of excited helium by the mid-IR laser pulses, i.e. ADC(1). Within ADC(1), the configuration manifold used to describe the ionisation of helium by the laser pulses, via solving the TDSE, consists of the singly excited one-hole-one-particle (1h − 1p) configurations. The 1h − 1p manifold allows one to describe ionization as well as excitation of helium atom in any of the singly-excited bound states. The time propagation of the unknown coefficients α j (t) of the B-spline ADC(1) many-electron wave-function is performed by means of the Arnoldi-Lanczos, algorithm [44,45,52]. The B Spline ADC within ADC(1) has been successfully used in the strong-field regime to model the intensity-dependent interference minimum that is present in the high-order harmonic spectra of CO 2 [44].

B. One-electron models
The methods below approximate the system by a one-electron system, whose evolution is given, in atomic units, by the single-electron time-dependent Schrödinger equation (TDSE) where the Hamiltonian describes the joint influence of the binding potential and the external field. Thereby, gives the field-free one-electron atomic Hamiltonian andr andp denote the position and momentum operators, respectively. The coupling with the field is given by the interaction Hamiltonian H I (t) =r · E(t) and H I (t) = p · A(t) + A 2 /2 in the length and velocity gauges, respectively, where E(t) = −dA(t)/dt is the electric field of the external laser field and A(t) the corresponding vector potential. Eq. (5) equation is either solved numerically using the freely available software Qprop [41,42] or semi-analytically using the SFA or the CQSFA. Throughout, we employ the effective potential where f (r) = a 1 e −a2r + a 3 re −a4r + a 5 e −a6r (8) and r(τ ) = r(τ ) · r(τ ). The coefficients are chosen as a 1 = 1.231 a.u., a 2 = 0.662 a.u., a 3 = −1.325 a.u., a 4 = −1.236 a.u., a 5 = −0.231 a.u. and a 6 = 0.480 a.u. [53]. These parameters were obtained by fitting to a numerical potential computed by self-interaction free density functional theory [54,55].
For Qprop, the TDSE is solved in the velocity gauge, while for the semi-analytic approaches the length-gauge Hamiltonian is used. This is mainly for practical reasons: the numerical solution of the TDSE converges faster for the velocity gauge, while the length gauge gives better results for ATI in the SFA [56] and CQSFA [4]. For Qprop we refer to [41,42], while the semi-analytic models will be summarized below. In all one-electron models, for the perpendicular initial p state, the valued real p x orbital is used via the corresponding coherent superposition of states with m = ±1.
In Qprop it is possible to compute strong field ionization for the initial states |ψ 2p±1 , where the subscript ±1 refers to the quantum number m. The 'x' oriented state may be computed by the following superposition, One can simply take the amplitude components (both real and imaginary) for the initial states |ψ 2p±1 and compute the coherent superposition above (9) to receive |ψ 2px amplitude, which is needed to provide the yield for the photoelectron momentum distributions.
A convenient starting point for the semi-analytic expressions is the Schrödinger equation (5) in integral form, namely where U a (t, t 0 ) = exp[iH a (t − t 0 )] is the time-evolution operator associated with the field-free Hamiltonian, and the time evolution operator where T denotes time-ordering, relates to the full Hamiltonian H(t) evolving from an initial time t 0 to a final time t. Alternatively, one may construct the integral equation as where U (V ) (t, t 0 ) is the Volkov time-evolution operator, associated with the Hamiltonian We wish to calculate the transition amplitude ψ p (t)|U (t, t 0 )|ψ 0 from a bound state |ψ 0 to a final continuum state |ψ p (t) with momentum p. Using Eq. (10) and the orthogonality between continuum and bound states, it can be written as with |ψ 0 (t ) = exp[iI p t ] |ψ 0 , where I p is the ionization potential and we have set t 0 = 0. One should note that, apart from considering a single active electron, no approximation has been made in the time propagation described by Eq. (14). Some approximations and/or further assumptions will be employed next.
Physically, the SFA relies upon the assumption that all bound states apart from the initial state of the system do not need to be considered. That is, the state of the system at a time t can be written as where |ψ p are the continuum states, and |ψ 0 (t) is the population in the ground state, often written |ψ 0 (t) = a(t) |ψ 0 . If the laser field is not intense enough for the ionization to reach the saturation level, a widespread approximation is a(t) = e iIpt , which means that the ground-state depletion can be neglected.

Transition amplitudes
For ATI, the lowest nonvanishing term of the Born-type series is obtained if the full time evolution operator is replaced by the Volkov time evolution operator U (V ) (t, t ) in Eq. (14). The subsequent term will include a single interaction in Eq. (12). This will lead to two terms that give the amplitude for producing a photoelectron with a particular well-defined momentum, The term M Dir (p) can be physically interpreted as a direct contribution from electrons that reach the detector with no further interaction with the core, while M Resc (p) is associated with electrons that are driven back to the parent ion by the laser field and are subsequently elastically rescattered. The expression for the direct term is where θ is the polar angle, H I (t) is the length-gauge atom field interaction and the action reads The rescattered term is given by The three steps in the rescattering process delineated in Eq. (19) are: (i) ionization at time t to a Volkov state with canonical momentum k; (ii) propagation in the laser field from t to t, where t > t ; (iii) scattering off the ion to a final Volkov state with canonical momentum p at time t, where V can be a model potential for the ionic system in question. It is very common to calculate the SFA transition amplitudes employing saddle-point methods, as they simplify the numerical effort and provide a very intuitive, clear picture in terms of rescattering electron trajectories (for details see, e.g., [58]). However, the standard saddle-point method selects the rescattering processes that are exactly along the driving-field polarization. This will be problematic for initial bound states aligned perpendicularly to the field polarization. A discussion of how to compute these matrix elements explicitly and how to implement saddle-point methods in the present context is provided in the appendix.
For early discussions of how to incorporate rescattering in ATI see Refs. [5] and [6]. Further detail about the SFA and strong-field ionization is provided in the tutorial [59] and in the review articles [8,60]. Introducing the term M Resc (p) is also known as the "improved strong-field approximation" [61].

Bound-state transitions
The standard formulation of the SFA method assumes that the active electron occupies a specific ground state and that it can only transition into continuum states through interaction with the laser. The influence of other bound states on the dynamics of the system is typically not considered.
However, if the initial state of the system is excited to a point where it is within energetic reach of other excited states, this assumption may become problematic. In the case of the helium atom excited into the 2 1 P state, particular attention should be paid to the transition dipole between the 2 1 P state and the 2 1 S state. Here, the energy gap is relatively small, approximately 0.0221 a.u (0.6 eV), and the transition dipole moment is large, approximately 2.9 a.u. In principle, other possible bound-state transitions, such as 2 1 P -1 1 S or 2 1 P -3 1 S can be more safely ignored as the energy difference that characterizes these transitions is significantly larger. Excitation has also been incorporated in the SFA in the context of NSDI [62][63][64][65][66][67]. For coherent superpositions of states in the SFA which, however, are only coupled via the continuum see, e.g., [68].
To incorporate this into the SFA model above, we modify the system's state space vector (see Eq. (15)) as follows: This would modify the first term in Eq. (16) to where the index j runs over the excited states to be included in the calculation. The same modification is made to the expression for the term M Resc (p) in Eq. (16). The model now describes ionization from a series of strongly coupled bound states. We determine the time-dependent behaviour of the bound states in a laser field with the system of coupled equations: where I pj is the ionization potential for the excited state j and d j−i is the matrix element of the dipole transition from one excited state to another with the orientation of this element determined by the spatial polarization of the excited orbitals. If the number of excited states under consideration is two and the quantity 2E 0 · d j−i (where E 0 is the maximum field strength) is much greater than the bound-state energy difference between the states I Pj − I Pi , then the system of equations in (22) can be solved perturbatively [69]. The solution for our helium system is whereĨ p is the average energy of the two excited states and the 0 in the superscript indicates that this solution is to zeroth order. This means that the energy difference is neglected and the states are considered as degenerate. The condition for this solution 2E 0 · d 2 1 S−2 1 P I P 2 1 P − I P 2 1 S is not valid in the case of a finite pulse at the edges of the envelope, but remains applicable provided an alternative condition is met: that ω I P 2 1 P − I P 2 1 S where ω is the field frequency. It can be readily seen from Eqs. (23) and (24) that the solution for the perpendicular configuration is This is because the 2 1 P -2 1 S transition is not possible due to the symmetry constraints of the dipole transition in the case of the perpendicular configuration.

The Coulomb quantum-orbit strong-field approximation
Instead of constructing iterative schemes around Eq. (10), one of which is the SFA, one may also employ pathintegral methods and construct approaches that incorporate the field and the binding potential on equal footing. One of such approaches is the Coulomb quantum-orbit strong-field approximation (CQSFA) [35,46]. It is constructed from the transition amplitude which is obtained by inserting a closure relation in the initial field-dressed momentump 0 = p 0 + A(t ) in Eq. (10) and is exact for a one-electron system. Eq. (25) takes the system from the initial photoelectron bound state |ψ 0 (t ) to its final momentum state |ψ p (t) = |p f (t) , where the tilde indicates field-dressing.
t is the vector potential, denote the electron's initial and final field-dressed momenta, respectively. The time-evolution operator U (t, t ) is associated to the full Hamiltonian, and thus includes the laser field and the binding potential. We take the interaction Hamiltonian H I (t ) to be in the length, gauge.
Applying time-slicing techniques [70,71] gives where D p and Dr are the integration measures for the path integrals [35,70], and the prime denotes a restriction.
The action and the Hamiltonian in Eq. (26) reads as and respectively. In Eqs. (27) and (28), the intermediate momentum p and coordinate r have been parameterized with regard to the time τ andp = p(τ ) + A(τ ). One should note that no Born-type expansion was used in the derivation of Eq. (26). This implies that the CQSFA cannot be viewed as a field-dressed perturbative series with regard to the binding potential, and that the distinction between direct and rescattered electrons is blurred. This is even clearer when Eq. (26) is solved using saddle-point methods, i.e., we search for values of t , r and p that render the action (27) stationary. This gives Eq. (29) yields the energy conservation at tunnelling, and Eqs. (30) and (31) give the subsequent electron propagation. We compute Eq. (25) along a two-pronged contour that starts in at the complex time t = t r + it i , goes parallel to the imaginary time axis, i.e., t to t r , and, subsequently, along the real time axis from the real ionization time t r to the final time t → ∞ [72][73][74][75]. The first and second arms of the contour are associated with under-the-barrier dynamics and continuum propagation, respectively.
In the under-the-barrier part of the contour, we kept the momentum constant and equal to p 0 [4,48]. This reduces Eq. (29) to In the continuum propagation, the action can be simplified according to the procedure discussed in [35,46,53,76]. The CQSFA transition amplitude obtained using the saddle-point approximation reads where p s and r s and t s are determined by Eqs. (30), (31) and (32), respectively. The term in brackets gives the stability of a specific trajectory, and In practice, we use ∂p s (t)/∂p s (t s ) instead of the stability factor in Eq. (33). This may be obtained employing a Legendre transformation and will not affect the action if the electron starts from the origin. The continuum propagation is performed by solving the inverse problem, i.e., given a final momentum p f (t), we find the initial momentum p 0 at the tunnel exit. This is defined as r 0 (t r ), where A further approximation used here is to take a real tunnel exit where r 0z indicates the component of the tunnel trajectory along the field polarization direction. This will simplify the problem and lead to real trajectories in the continuum. For discussions of how to solve the full complex problem see [48,77]. The term p + A(t s )|H I (t s )|ψ 0 contains the influence of the initial bound state, i.e., its geometry, and will be important for the present work. Explicitly, where and φ α (r) is given by a Gaussian basis set For details see [68].
In the CQSFA, one may identify four main types of orbits: 1. Type 1 orbits behave like short, direct SFA trajectories that are slowed down by the Coulomb potential. They leave in the direction of the detector and their transverse momentum does not change direction.
2. Type 2 orbits are field-dressed Kepler hyperbolae whose start times are displaced in half a cycle from those in type 1 orbits. They start on the "wrong" side and are turned towards the detector. They are lightly deflected by the potential and can be loosely associated with the long SFA orbits. Their transverse momentum does not change its sign during continuum propagation.
3. Type 3 orbits are also field-dressed Kepler hyperbolae starting in the same half cycle as type 2 orbits, but interact more strongly with the core. Their transverse momentum changes direction during continuum propagation. They have no counterpart in the SFA as they depend on the interplay between the driving field and the binding potential and are not embedded in a Born-type method. In previous work, we have approximated orbit 3 within an analytic piecewise model in which rescattering is incorporated. However, mimicking the effect of orbit 3 and accurately reproducing the spider was only possible by including the acceleration caused by the Coulomb potential during the continuum propagation and the Coulomb phase [47]. This is a substantially different scenario than that encountered with the SFA, be it direct or rescattered, or in approaches focused on short-range potentials [78]. In the limit of vanishing residual potential for the continuum propagation, orbits 2 and 3 become degenerate and tend to the long direct SFA orbit [35]. 4. Type 4 orbits start in the same half cycle as type 1 orbits, but go around the core before reaching the continuum. Their behaviour is close to that of SFA rescattered trajectories, and they may reach energies of around 10U p . However, one should notice that the process they undergo is different and their propagation in the continuum is not restricted to the field-polarisation axis.
One should note that the energies of orbits 2 and 3 may go beyond the direct ATI cutoff of 2U p , due to the influence of the residual potential. This contributes to blur the distinction between direct and rescattered, as mentioned above. A limitation of the CQSFA is that it does not properly account for processes involving excited bound states. The basis chosen in Eq. (25) by using the closure relation inp 0 works well for continuum states but does not reproduce excitation processes accurately. Still, it does take into account highly excited states which are strongly coupled with the field and behave similarly to the continuum. The four types of orbits stated above were first identified in [34].

III. RESULTS
In the following, we will analyse how different polarisations in the initial bound states influence angle-integrated photoelectron spectra and photoelectron momentum distributions, with emphasis on several counterintuitive features. We used a bandwidth-limited mid-IR pulse of intensity I = 3.2 × 10 13 W/cm 2 , frequency ω = 0.775 eV, that corresponds to a wavelength λ = 1600 nm, linear polarisation along the z direction, with a ponderomotive energy (U p ) of 7.65 eV. Unless otherwise stated, the binding energy used for this initial state is 3.286 eV (I p = 0.1238 a.u.).
For this parameter range, we have verified that the electron in the ground state orbital 1s is not appreciably affected by the IR pulse in the TDSE solvers.
For the B-Spline ADC method and the SFA, the pulse used is described by a cosine squared envelope for the vector potential of amplitude A 0 and polarization vectorê z , and with a pulse duration T set to 10 cycles (ca. 50 fs). The calculations have been performed using a linear B-spline knot sequence [43] with a radial box radius of R max = 500 a.u. and N b = 625 radial B-splines. The maximum angular momentum employed in the monocentric expansion of Eq. (3) was l max = 30. Convergence of the results with respect to the basis set parameters has been checked. A complex absorbing potential (CAP) with starting radius R CAP = 400 a.u. has been used to absorb the wavefunction and avoid its reflections from the grid boundary. Each and every B-spline ADC(1) energy-dependent photoelectron spectra presented here were calculated by coupling the TD B-spline ADC method to the t-surff technique.
For the CQSFA, for simplicity we consider a linearly polarized monochromatic wave, whose vector potential is described by This is a good approximation for long enough pulses. In practice, we take the field to be over 20 cycles long, and include ionization events within up to four laser cycles, that is, the range of ionization times t lies within 0 ≤ t ≤ 8π/ω. A restricted temporal unit cell is necessary for practical purposes, but may introduce some boundary effects. For a complete discussion see [79]. For Qprop, we have employed a four-cycle flat-top pulse with additional linear half-cycle turn on and off for the comparison to the CQSFA and SFA models. In the comparison of Qprop with the B-spline ADC model, we have matched the appropriate pulse of sin 2 , whose vector potential is given by Eq. (40).
A. ATI Spectra

One and two-electron Schrödinger solvers
We will start by discussing the angle-integrated ATI spectra obtained from the numerical solutions of the TDSE employed in this work, calculated with the diagrammatic B-Spline ADC method and with the freely available oneelectron Schrödinger solver Qprop [41,42]. These results are presented in Fig. 1, for the parallel and perpendicular initial-state configurations (upper and lower panel, respectively). For the parallel configuration, the spectra obtained with the B-Spline ADC method exhibits a long ramp-like structure extending from near the ionization threshold up to the high photoelectron energy of 10U p , while the Qprop computations exhibit a low-energy region followed by a ramplike decrease from 2U p to 4U p (see pink shaded region), with a flat plateau up to 10U p . The high-energy suppression, in the case of TD B-Spline ADC, is caused by the "potential" that the excited/ionised electron experiences, which has contributions from both the Coulomb and exchange terms calculated with the Hartree-Fock orbitals of the helium atom. These exchange terms effectively repel the rescattered electron if it gets too close to the core, and thus render rescattering less effective. This influence increases with the photoelectron final energy, as higher-energy electrons must approach the core more closely. These ramp-like features are not altered significantly if the bound-state transition 2 1 P − 2 1 S is removed from the B-Spline ADC computation (see blue and green curves).
On the other hand, if the initial 2p state is aligned perpendicular to the driving-field polarization (lower panel in Fig. 1), the computations lead to strikingly different results. For clarity, we have included the parallel-polarized spectrum without the 2 1 P − 2 1 S transition. One should note that, due to selection rules, this transition is forbidden in the perpendicular configuration. Up to the photoelectron energy of roughly 2.5U p , there are at most subtle differences, with the Qprop spectrum being slightly larger than the B-Spline ADC spectra. Subsequently, one sees a steep ramp-like suppression for the (one-electron) Qprop computation up to 4U p succeeded by a second ramp of much gentler slope. The B-Spline ADC spectra are significantly less suppressed for energies higher than 2U p , regardless of configuration. Thus, in contrast to the parallel-aligned case, electron exchange appears to enhance rescattering. The results suggest that exchange effects, which are included in the B-spline ADC simulation, lead to a decrease of the dependence of the yield on the direction of the driving field with respect to the initial atomic polarization, especially in the (approximately) 2U p − 6U p energy region. Interestingly, for the perpendicularly polarized initial state the B-Spline ADC computation exhibits multiple ramps instead of the monotonic behavior identified for parallel alignment. Thereby, one may identify two ramps, for photoelectron energy ranges 2U p ≤ E k ≤ 4U p and 4U p ≤ E k ≤ 6U p , respectively. Furthermore, in comparison to the parallel-aligned case (blue curve), the B-Spline ADC spectrum  [41,42], while the green (solid) and the blue (dotted) curves have been computed with the diagrammatic B-Spline ADC method including and excluding the bound-state transition 2 1 P − 2 1 S, respectively. For the Qprop calculation a 4-cycle sin 2 pulse was used, while for the B-Spline ADC a 10-cycle sin 2 pulse was taken and pulse length of 53 fs. The upper and lower panels have been calculated for z (parallel) and x (perpendicular) orientation of the initial 2p bound states of helium with regard to the driving-field polarization. The blue curve in the lower panel has been computed for parallel polarization, but removing the bound-state transitions, the most important among them being 2 1 P − 2 1 S. Each curve has been normalized to its largest value to facilitate a qualitative comparison. The pink and violet shaded areas highlight the photoelectron energy ranges 2Up ≤ E k ≤ 4Up and 4Up ≤ E k ≤ 6Up, respectively. computed for the perpendicular configuration (green curve) exhibits a decrease of one order of magnitude for energies higher than 4U p even in the absence of the 2 1 P − 2 1 S transition. Still, this behaviour occurs for much larger energies than those in the Qprop spectrum. The remarkably different behaviour from Qprop suggests a multielectron character for the second ramp as well.

Comparison with orbit-based methods
For that reason, we will now compare the outcome of Qprop with that of semi-analytical, orbit-based methods. We will start by discussing the spectra obtained with the strong-field approximation (SFA), which is shown in Fig. 2. In all cases, the energy regions in the SFA spectra can be clearly associated with direct and rescattered ATI, with a distinct cutoff at photoelectron energy 2U p followed by long plateaus up to the rescattered ATI cutoff of 10U p . The direct contributions dominate up to approximately 4U p , with the plateau prevailing for higher energies. For perpendicular alignment, the plateau is much more suppressed than in the parallel-aligned case. This is consistent with the picture put across by the SFA, that rescattering occurs for a small angular range around the polarization axis. Because there is not much probability density with which the returning electronic wave packet can overlap, the resulting signal will be suppressed. This adds up to an overall suppression that occurs for similar reasons at the instant of ionization. Since tunneling is highly directional and occurs predominantly along and in the vicinity of the polarization axis, the overall signal will be around two orders of magnitude weaker in the perpendicular case. We have, however, normalized the yield in the lower panels of the figure in order to perform a qualitative comparison.
In comparison with the Qprop calculation, the SFA plateau is strongly suppressed for parallel orientation, although they follow the same trend (see upper panel in Fig. 2). In fact, it is noteworthy that there are only two orders of magnitude difference between the direct and rescattered signal in the Qprop case, while for the SFA this difference amounts to four orders of magnitude. Moreover, the behavior in the direct region is much flatter than for the SFA, with a ramp for photoelectron energies between 2U p and 4U p (see pink shaded region). Interestingly, including the 2 1 P − 2 1 S bound-state transition in the SFA as discussed in Sec. II B 1 leads to subtle changes in the plateau, but does not alter its overall shape or intensity.
For perpendicular-aligned initial states (lower panel in Fig. 2), the Qprop and SFA computations are much more similar, with a steep decay in several orders of magnitude for 2U p ≤ E k ≤ 4U p (dubbed 'mid-energy ramp') followed by a plateau. Still, in that energy region the TDSE result always follows that of the SFA from above, and exhibit a ramp-like structure. One should note that this mid-energy ramp appears in a single-active electron setting. This is a key difference with regard to the behavior seen for the second ramp in the B-Spline ADC computations; see discussion of Fig. 1. The mid-energy ramp is associated with hybrid orbits, which are missed by the SFA, and will be analyzed next.
In Fig. 3, we plot the spectra obtained with Qprop against those from the CQSFA. The spectra obtained with the CQSFA exhibit an almost constant, much higher plateau, whose onset occurs for energies immediately higher than 2U p . This is in stark contrast with the ramp-like behavior obtained in the previous cases. Furthermore, a direct comparison of the parallel and perpendicular configurations shows a suppression of at least one order of magnitude for the plateau in the latter case. This is much less than that observed in all other methods, and it is also curious that the parallel and perpendicular spectra exhibit very similar shapes. This suggests that, for the parameter range employed, and in particular the ionization potential I p = 0.1238 a.u. of the excited states, the CQSFA orbits are behaving in an unusual way, which is much more dictated by the binding potential than by the laser field.
A further investigation is presented in Fig. 4, in which, in addition to the previous curves, calculated for excited helium, we have included a CQSFA computation for a slightly larger ionization potential (I p2 = 0.1750 a.u.) and 2p initial states parallel and perpendicular to the laser-field polarization. This is quite different from the behavior we have observed if the correct ionization potential (I p1 = 0.1238 a.u.) is taken. Thereby, one can see that, for the perpendicular and parallel-oriented cases, the spectra practically overlap up to energies up to 2U p . However, for the perpendicular case there is a long ramp in the photoelectron energy region 2U p ≤ E k ≤ 4U p followed by a strongly suppressed ATI plateau. This behavior is in better qualitative agreement with the Qprop outcome than that obtained for the correct ionization potential. Next, we will analyze these behaviors in more detail using photoelectron momentum distributions.

B. Photoelectron momentum distributions
In the following, we will have a closer look at the electron dynamics by inspecting photoelectron momentum distributions (PMDs). For simplicity, we restrict ourselves to a one-electron scenario employing the CQSFA and using Qprop [41,42]

units)
FIG. 2: Angle-integrated photoelectron spectra computed for excited helium (Ip = 0.1238 a.u.) using Qprop (orange and brown dotted curves for parallel and perpendicular cases, respectively) and the strong-field approximation (SFA) (remaining curves), for the same peak-field intensities and wavelengths as in the previous figures. For the Qprop calculation a flat-top four cycle pulse with an additional half a cycle turn on and off was used, while for the SFA we have taken only a single cycle. Each curve has been normalized to its largest value to facilitate a qualitative comparison. For the perpendicular configuration, the SFA spectrum is one order of magnitude smaller than its parallel counterpart in the direct region (energies smaller than 2Up and the rescattered plateau is two orders of magnitude smaller. The top panel (parallel alignment) also shows all scattering (Direct and Rescattered) contributions in the SFA for the bound coupled state, 2s2p (purple wide dashed curve). The pink and violet shaded areas highlight the photoelectron energy ranges 2Up ≤ E k ≤ 4Up and 4Up ≤ E k ≤ 6Up, respectively.
helium (upper and lower rows, respectively) in a four-cycle trapezoidal pulse. There is an overall suppression of the signal in the perpendicular configuration, but in the present discussion we will focus on qualitative features and the physics behind them. All distributions exhibit well-defined ATI rings and several holographic structures, such as a spider-like pattern around the polarization axis, fan-shaped fringes close to the ionization threshold and spiral-like patterns perpendicular to the driving-field polarization. These structures have been discussed in previous publications [33,[46][47][48] and have been observed in numerical computations and experiments (for details see the recent review article FIG. 3: Angle-integrated photoelectron spectra computed for excited helium (Ip = 0.1238 a.u.) in initial 2p states oriented parallel (z) and perpendicular (x) to the laser-field polarization, using Qprop (dotted orange and brown lines) and the coulomb quantum-orbit strong-field approximation (CQSFA) (purple and red lines, dashed and solid respectively), for the same peakfield intensities and wavelengths as in Fig. 1. For the Qprop calculation, a flat-top four cycle pulse with an additional half a cycle turn on and off was used, while for the CQSFA, we have considered ionization events within one cycle of a monochromatic field. Each curve has been normalized to its largest value to facilitate a qualitative comparison. The pink and violet shaded areas highlight the photoelectron energy ranges 2Up ≤ E k ≤ 4Up and 4Up ≤ E k ≤ 6Up, respectively. [4]). There are, however, some qualitative differences between the CQSFA and the TDSE computations. First, the holographic patterns in the TDSE outcome are asymmetric with regard to the p ⊥ axis and the ATI rings are less sharp. This asymmetry is very likely caused by a strong bound-state depletion, which is present in the TDSE case. Depletion will affect the probabilities associated with different ionization events such that the later events will be suppressed, and this will cause distortions in the patterns, making the distributions more similar to those from a few-cycle pulse. As such, we have verified that the initial excited-state population is quickly reduced for the parameter range employed (not shown). This effect has not been incorporated in the CQSFA. The CQSFA also has some asymmetry for a different reason. Despite using a monochromatic field, we must restrict the ionization events to a finite temporal range, which introduce some artificial asymmetry depending on how this range is chosen. This can be eliminated by the incoherent averaging procedure described in [79]. the Coulomb quantum-orbit strong-field approximation (CQSFA) (purple and red curves, wide dashed and solid respectively) computed for excited helium (Ip1 = 0.1238 a.u.) using the same parameters as in Fig. 3, but with two additional (blue and green, dot-dahsed and small dashes respectively) curves computed using the CQSFA for a slightly larger ionization potential (Ip2 = 0.1750 a.u.). Each curve has been normalized to its largest value to facilitate a qualitative comparison, except for those computed for the larger ionization potential Ip2. The pink and violet shaded areas highlight the photoelectron energy ranges 2Up ≤ E k ≤ 4Up and 4Up ≤ E k ≤ 6Up, respectively. In panels (c) and (f), we have considered the CQSFA and a slightly larger ionization potential (Ip = 0.175 a.u.) as a test case. The upper and lower rows have been computed for initial bound states aligned parallel and perpendicular to the driving-field polarization, respectively. For Qprop, we have used a four-cycle flat-top pulse with an additional half cycle turn on and off. For the CQSFA, we have approximated the driving field by a monochromatic wave. In all cases, we considered a driving-field peak intensity I0 = 3.2 × 10 13 W/cm 2 (Up = 7.65 eV), and wavelength λ = 1600 nm. The thick circular dashed line indicates the direct ATI cutoff 2Up and the yield in each panel has been normalized to its maximum value, but there are roughly between one and two orders of magnitude difference between the parallel and perpendicular configuration. Furthermore, the yield for the perpendicular case is suppressed by a factor of 10 2 for CQSFA and 10 for Qprop, respectively. The green dashed rectangle on panels (d) and (f) indicate the suppressed rescattering ridge for a particular perpendicular momentum, whereas the dashed triangle shows off axis suppression in panels (b), (c), (e) and (f). Finally, the red squares seen on all panels show the spilling beyond the 2Up cutoff as mentioned above.
Another noteworthy aspect in the PMDs are superimposed circular regions of radius 10U p centered at (p , p ⊥ ) = (±2 U p , 0), which overlap at (p , p ⊥ ) = (0, 2 U p ) along the perpendicular momentum axis. The boundaries of such regions determine a ridge, which is a well known indicator of rescattering. The rescattering ridge is present in the CQSFA computation, regardless of the orientation of the initial orbital [see middle panels in the figure], but is strongly suppressed for the TDSE in the case of perpendicular orientation [see lower left panel]. This suppression is not isotropic, and remnants of the ridge can be seen around the axis p ⊥ perpendicular to the driving-field polarization. If, on the other hand, we consider a slightly larger ionization potential for the CQSFA, we also observe a very faint rescattering ridge for the perpendicular case (see right column in the figure). In addition to that, there are suppressions of the photoelectron signal around the p ⊥ axis, which are stronger for perpendicular alignment. Subtler features are a spilling of the spider-like structure beyond the direct 2U p cutoff and a wedge-like suppression (shown on Fig. 5 by the dashed square and triangle, respectively.) along the p axis for the perpendicular bound-state configuration in all approaches.
We will next analyze the above features in terms of CQSFA orbits. With that aim in mind, in Fig. 6, we plot PMDs obtained using individual CQSFA orbits. The plots illustrate the momentum region they occupy, and how they are influenced by the parallel and perpendicular configurations. We will commence by looking at the parameters in Figs. 5(b) and (e), which, apart from an overall suppression of roughly two orders of magnitude in the perpendicular case, look qualitatively similar.
Orbits 1 and 2 lead to approximately elliptical PMDs, whose major axis is along p and which are mostly constrained by the 2U p direct ATI cutoff. This is expected as both orbits behave primarily as direct SFA orbits. Spilling beyond this cutoff occurs mostly along the driving-field polarization and is caused by the presence of the Coulomb potential. For perpendicular orientation, there is a wedge-like suppression along the p axis for the contributions of orbit 1 [see Fig. 6(b) in comparison with Fig. 6(a)], whose angular range broadens as |p | increases. This effect can be understood as follows: if an electron propagates along orbit 1, its final and initial momentum will be similar, and, for high enough final momenta, close to the standard SFA [46]. Thus, if the electron is freed along the field-polarization axis, it will continue to propagate along this direction. For small angles, its propagation direction will not change significantly. The 2p orbital is strongly directional, so that if it is oriented perpendicular to the laser-field polarization, tunneling along the p axis will be hindered. The wedge-like shape comes from the fact that this is not a perfect mapping, and the final and initial momenta will only coincide for large values of p . Still, the above discussion will approximately hold. In contrast, as orbit 2 is essentially a field-dressed hyperbola, the binding potential plays a stronger role and the transverse momentum mapping is non-trivial. In this case, the wedge-like feature is absent.
The PMDs associated with orbits 3 and 4, shown in the bottom rows in Fig. 6, occupy momentum regions within a broader angular range than the previous single-orbit distributions. This is due to them undergoing a stronger deflection and interaction with the core. Previous studies have shown that orbit 3 has a hybrid character [33,47], with no counterpart in the SFA, while orbit 4 behaves essentially like a rescattered orbit. The PMD associated with orbit 3 mainly occupies the momentum region defined by the overlapping rescattering ridges, with a slight spilling near the polarization axes. In both parallel and perpendicular configurations, there is a suppression along p = 0, which, however, is more pronounced in the latter case. Orbit 4 leads to a PMD with clear rescattering ridges for both parallel and perpendicular configurations and a caustic structure at the boundary. This indicates that, for the present parameter range, the dynamics associated with orbit 4 are not predominantly located along the polarization axis, as  Fig. 5 and a single cycle, but considering a slightly higher ionization potential (Ip = 0.175 a.u.). The first columns panels (a) and (c) have been computed for parallel case. Whereas perpendicular alignment with regard to the driving-field polarization, is the second column, panels (b) and (d) and the circular lines indicate the direct ATI cutoff at energy 2Up across all panels. The two upper panels have been normalized to their highest value, but there is a difference of one to two orders of magnitude between the parallel and perpendicular cases. In the lower panels the same scale was used.
would be expected in an SFA treatment of rescattered trajectories. On the other hand, if we consider a slightly larger ionization potential, there is a much stronger suppression of the ridge for orbit 4 and of the signal around the p ⊥ axis for orbit 3. This is illustrated in Fig. 7.
These are great examples of how rescattering may no longer be restricted to the polarization axis due to the Coulomb potential. For the correct value of the ionization potential associated with the excited state, the CQSFA tunnel exit is located in a region for which the Coulomb potential is dominant. This suggests that orbit 4 will behave in a much less directional way than what one expects within the SFA framework. It will be able to probe orbitals oriented perpendicular to the driving-field polarization, and the rescattering ridge will be present. On the other hand, for higher ionization potentials, the tunnel exit will be located in the region for which the field is dominant. This will reduce the angular range for which rescattering may occur, and will cause a strong suppression for photoelectron energies higher than 2U p . This issue is illustrated in Fig. 8, which shows that this transition depends quite critically on the interplay between the ionization potential and the ponderomotive energy. For the ionization potential I p = 0.15 a.u., which is slightly larger than that of excited helium (orange curve), the binding potential dominates. This causes orbit 4 to leave in a perpendicular direction to the field, until approximately 2 a.u. Subsequently, the orbit leaves and returns to the core region propagating in the expected direction, that is, parallel to the driving-field polarization. Still, the large initial transverse momentum acquired by orbit 4 makes it much less directional than what is predicted by the SFA. This means that, for that parameter range, this orbit will also probe orbitals with perpendicular alignment. This explains the similar CQSFA spectra in the previous section, and the prominent rescattering ridges in Fig. 6(h) and Fig. 5(e). In contrast, if the tunnel exit is located in a region for which the field is dominant (blue curve), orbit 4 will behave in a much more directional way, practically along the polarization axis. This will cause an overall suppression for initial bound states oriented perpendicular to the driving-field polarization. Nonetheless, there are still momentum components in the perpendicular direction and the orbit does not leave parallel to the field. This, together with the presence of intermediate orbits, leads to remnants of the rescattering ridge in the PMDs.
Therefore, even in a single-electron setting, the residual binding potential influences rescattering. It gives rise to orbits of hybrid character, and, in some instances, makes ionization and/or recollision deviate from the driving-field polarization axis. Hybrid orbits are absent in Born-type approaches such as the SFA, and will contribute for the anisotropic suppression in the rescattering ridge reported in this work. Ionization, and in some cases propagation, perpendicular to the laser-polarization axis may occur and will partially hinder the suppression of the rescattered ATI signal. One should note, however, that the CQSFA results for excited helium strongly deviate from the TDSE computations. Interestingly, the TDSE agrees more with the CQSFA if a higher (artificial) ionization potential is taken in the CQSFA.

IV. CONCLUSIONS
In this work, we have probed excited helium in initial 2p-state configurations parallel and perpendicular to a mid-IR, linearly polarized field. We use different numerical methods, namely the B-spline algebraic diagrammatic construction (ADC) [43] and the one-electron Schrödinger solver Qprop [41,42], whose outcomes were compared to those of two semi-analytic methods: the Strong-Field Approximation (SFA) and the Coulomb quantum-orbit strongfield approximation (CQSFA). Our main goal was to understand the limitations of the rescattering model in its simplest form, i.e., that dictated by the SFA, and what subtleties must be taken into consideration. Overall, it is noteworthy that different computations may lead to spectra whose shapes and intensities are quite distinct.
The SFA predicts a sharp decrease in the ATI signal for energies above the direct ATI cutoff 2U p for all cases, and a strong suppression in the high-energy ATI plateau, which extends up to the energy of 10U p , for the perpendicular initial state. This is due to tunnel ionization and recollision being located around the field polarization axis. This implies that both processes are suppressed due to the initial-state geometry. First, there will not be a substantial probability density to undergo tunneling and, upon the act of rescattering, there will not be a significant overlap of the returning wavepacket with the orthogonally oriented bound state. Furthermore, in the SFA no changes that may occur in the core and alter its geometry are incorporated. This sharp decrease is similar to that resulting from changing the driving-field polarization, instead of that of the initial bound state. Any deviations from those patterns mean that the dynamics are more intricate, either because ionization or recollision may occur off axis, or because processes involving the core must be taken into consideration.
The results presented in this paper reveal a subtler picture, with the angle-integrated ATI spectra exhibiting ramplike structures in some computations. This is most striking for the B-spline ADC computation, which goes beyond a single-active electron and has the least degree of physical approximations. The B-spline ADC spectra exhibit multiple ramps and a smooth transition from direct to rescattered ATI, both for parallel and perpendicular orientation. A ramp-like structure is expected if there are many superimposed processes and/or substantial depletion. This will tend to weaken effects caused by the initial bound-state geometry or suppress high-energy features by removing parts of the bound-state population via several channels. It could also be caused by mechanisms rendering rescattering less effective.
It is noteworthy that the B-Spline ADC spectra behave in strikingly different ways from those obtained with the one-electron computations. The spectra from Qprop and the analytical models exhibit a much more obvious transition from the direct to the rescattered regimes, while the B Spline ADC spectra present a continuous ramp. Removal of a specific bound-state excitation channel (2 1 P − 2 1 S) in the B-Spline ADC influences the spectra. However, the changes are too subtle to justify all the discrepancies. The current results suggest that in the parallel-aligned case the presence of electron exchange renders rescattering less effective by introducing repulsive effects, while for the perpendicular-aligned configuration multielectron effects enhance tunneling and low-energy rescattering. Because exchange is non-local, it renders both effects less directional, which decreases the suppression in the perpendicular case. Interestingly, a suppression of one order of magnitude between parallel and perpendicular polarised initial p states has also been observed for ATI spectra computed for argon using the time-dependent density functional theory (TDDFT) [80]. At the ADC(1) level at which our simulations were done, we are close in accuracy to the TDDFT, and both computations incorporates electron exchange. This suggests that rescattering in ATI is a delicate effect that probes the fine details of the core structure and core-photoelectron interaction including its non-local (exchange) component. The influence of electron exchange will also be energy dependent. The closer to the core the rescattering electron gets, the higher the photoelectron energy will be. This also means that exchange effects will become increasingly important, leading to a ramp.
Still, even for single-electron computations, the spectra exhibit quite different features, depending on what is included or left out. For instance, in the parallel-aligned case, the SFA spectrum shows a much stronger drop in signal for the plateau than that of Qprop. In the perpendicular-aligned scenario, the Qprop outcome follows the SFA spectrum from above in the direct region, and decays sharply after the direct ATI cutoff. In contrast, in the CQSFA the plateau is overestimated unless an artificially large ionization potential is taken. Moreover, it varies much less dramatically with the initial states' orientation than all other methods considered in this work.
The features mentioned above can be understood in greater depth by looking at angle-resolved photoelectron momentum distributions (PMDs) and by analyzing the electron orbits involved. This analysis has been performed for the CQSFA and Qprop. The agreement is good for the main features, such as the rescattering ridge and a suppression near the polarization axis that occurs for the perpendicular configuration. However, the holographic patterns are much more irregular for Qprop than for the CQSFA. Furthermore, in the perpendicular configuration the rescattered ATI signal is much more suppressed for Qprop. This includes most of the signal beyond the energy of 2U p and the rescattering ridge. A better agreement is obtained if the ionization potential is artificially increased in the CQSFA computation.
The irregularities in the Qprop holographic patterns are explained by the huge amount of bound-state depletion, which favours certain events over others and compromises the contrast in the quantum interference patterns. Furthermore, there may be other mechanisms through which the electron is freed in the continuum, such as over-thebarrier ionization and coupling with highly excited states. These features are not incorporated in the CQSFA. The over-enhancement of the rescattered contributions in the perpendicular case is due to the Coulomb potential being dominant upon ionization, for the CQSFA orbit type that leads to the rescattering ridge. This introduces a large perpendicular momentum component in the orbits, so that they initially tunnel in a perpendicular direction and may return farther away from the polarization axis.
Further studies indicate that a tunnel exit in a region for which the field dominates, which can be achieved with a slightly larger ionization potential, will alter this behavior. In this case, the orbit will be more localized around the polarization axis, so that tunneling and rescattering will be suppressed for perpendicular configurations. Interestingly, this suppression is closer to the behavior observed for the TDSE computation. The TDSE outcome will be influenced by ac Stark shifts, additional ionization channels, bound-state broadening and over-the-barrier ionization, which are not included in the CQSFA. Support for this assumption is provided by further TDSE computations using weaker binding energies, which exhibit a rescattering ridge for the perpendicular-aligned initial states (not shown). Some combination of these effects likely leads to the rescattering behaviour observed in the CQSFA when a larger ionization potential is used. Discussions of how to include Stark shifts in the SFA are provided in [81,82]. For a proposal of how to measure Stark shifts in excited helium, see [83,84].
If, on the one hand, this suggests that the CQSFA may need modifications for highly excited, loosely bound states, on the other hand this is a nice illustration of how the interplay of the binding potential and the external laser field may lead to counterintuitive behaviors and move rescattering or ionization away from the driving-field direction. This, together with the core dynamics effects observed for the B-spline ADC spectra, shows that there may be other possibilities for imaging than those based exclusively on the field exploring bound-state geometry or driving-field shapes. Finally, one should note that preparing targets in excited and oriented states is within experimental reach. For instance, one may use the Magneto-optical trap recoil ion momentum spectroscopy (MOTRIMS) technique (see, e.g. [85,86] in the context of electron collisions and [87] for multiply ionized Rubidium in a strong elliptically polarized field). Moreover, one may prepare the 2 1 P state of helium either by resonant excitation from the metastable 2 1 S state [88,89], or by a direct extraction from the ground state [90,91]. The polarization of the excited state can be controlled as described in [92], and the lifetime of the 2 1 P state is several orders of magnitude larger than the timescales involved in the present work. Usual excitation results in a superposition of the ground and excited state helium. However, for the parameters employed in this work, the ground state is too strongly bound to result in any appreciable ATI. Therefore, we are probing the excited-state population and whatever population stays in the ground state is simply lost for us.

V. ACKNOWLEDGEMENTS
We would like to thank Dr Bridgette Cooper for her insights and tutorials in using GAMESS and Dr Emilio Pisanty for his knowledge on colour schemes to enhance plots without the loss of key features.

Appendix. SFA Matrix elements
In this Appendix, we provide the main steps for calculating the integrals I 1 = dre ik·r V (r) (42) and I 2 = dre −ik·r r cos θψ 0 (r) (43) in Eqs. (17) and (19), for the effective binding potential (7). To avoid singularities in the Fourier transform of this potential, it is multiplied with a damping factor: V (r) → V (r) × e − r . Numerical simulations show that this factor does not distort the shape of the ATI spectrum, although it can shift its overall magnitude; for this study we have chosen = 1.0. For other methods to treat this singularity see, e.g., Ref. [93]. Using the spatial representation Ψ 0 (r) = R n0l0 (r)Y l0m0 (Ω) of the initial state, where Ω is the solid angle, Eqs. (42) and (43) can be simplified as I 1 = −4π 1 2 + k 2 + a 1 2 + a 2 2 + k 2 + 2a 3 a 4 ( 2 + a 2 4 + k 2 ) 2 + a 5 2 + a 2 6 + k 2 I 2 = (4π) 3/2 √ 3 lm i −l Y lm (Ω k ) drr 3 R n0l0 (r)j l (kr) dΩY * lm (Ω)Y 10 (Ω)Y l0m0 (Ω, ) where we have inserted the damping factor into Eq. (7), substituted cos θ by a spherical harmonic and used the multipole expansion for the exponential function in three dimensions. In Eq. (45), the solid angle Ω is defined in coordinate space and the solid angle Ω k with regard to the intermediate momentum k. The spatial radial integral can be solved analytically for an initial state described by a hydrogenic wave function. For the ground state (1s) we use Similarly, for 2s we obtain For the excited state 2p, where l 0 = 1, the sum over l runs from 0 to 2, giving the following three expressions in ascending order of l The angular integral of the product of three spherical harmonics is given by the product of Clebsch-Gordan coefficients (CG), whose general form reads as dΩY * l1m1 (Ω)Y l2m2 (Ω)Y l3m3 (Ω) = (2l 2 + 1)(2l 3 + 1) 4π(2l 1 + 1) CG(l 2 , l 3 , l 1 ; m 2 , m 3 , m 1 )CG(l 2 , l 3 , l 1 ; 0, 0, 0).
If the initial state is perpendicularly polarized, then the angular part Y l0m0 (Ω) in Eq. (45) will have some combination of m 0 = 1, −1. Therefore, the only nonzero contribution to the rescattering term will come from the element in the summation for which m = 1, −1. One should note that, so far, the rescattered ATI transition amplitude (19) incorporates all possible paths for this three-step process are integrated over, all possible ionization times for the ionization (when the laser is active) and rescattering processes, and all the possible intermediate Volkov states with canonical momentum k. The integrals in Eq. (19) can be calculated using saddle-point methods, which simplify the numerical effort involved and are easily relatable to an intuitive, orbit-based picture [14,58,94]. However, depending on the problem at hand care must be taken. If the saddle approximation is used in all three dimensions, then k x = 0 and k y = 0, this constrains the vector k to be fixed in the z-direction, leaving no contribution at all and not enabling meaningful comparisons to the numerical simulation results to be made. Physically, this implies that only those rescattering processes that occur exactly along the laser polarization axis are included in this approximation, and the remaining processes are selected out. For discussions in the context of molecules see [95].
We therefore use the saddle point approximation for just two dimensions in the integral and numerically integrate over the dimension of the atomic polarization (chosen in the x direction). The action can be expanded in three dimensions, so that S(k; t, t ) = 1 2 t t dt (k z + A z (t )) 2 + k 2 x + k 2 y Each variable in the integrand that we approximate is substituted by its stationary value obtained with the saddlepoint approximation, obtained setting the derivative of the action is equal to zero. This gives k st y = 0.
for ∂S(k; t, t )/∂k z = 0 and ∂S(k; t, t )/∂k y = 0, respectively. For each dimension in which we use the approximation, we make the following substitution where α is a small parameter to avoid a singularity when t → t . If the saddle point approximation is used in three dimensions, then this parameter is not necessary unless the initial state has p symmetry. This is because the exponent in the last integral of Eq. (19) goes to zero in the limit when t → t and the spherical Bessel function of such an argument tends to zero for all orders apart from for j 0 -this function only appears in the summation for an initial p state. Thus, modifying (45) with (46), (47) and (48), we obtain E(t )e iI P t e −iS(k sp z ;t,t ) e iS(p,t) dk x e −iS(kx;t,t ) × 1 (2π) 3 dr e i(k−p).r V (r ) × 1 (2π) 3/2 dre −i(k+A(t )).r cos θψ 0 (r) where k = k sp zẑ + k sp yŷ + k xx . This specific approach is being used both for parallel and perpendicular bound-state orientation. This integral was also performed in [96].