Advanced momentum sampling and Maslov phases for a precise semiclassical model of strong-field ionization

Recollision processes are fundamental to strong-field physics and attoscience, thus models connecting recolliding trajectories to quantum amplitudes are a crucial part in furthering understanding of these processes. We report developments in the semiclassical path-integral-based Coulomb quantum-orbit strong-field approximation model for strong-field ionization by including an additional phase known as Maslov’s phase and implementing a new solution strategy via Monte-Carlo-style sampling of the initial momenta. In doing so, we obtain exceptional agreement with solutions to the time-dependent Schrödinger equation for hydrogen, helium, and argon. We provide an in-depth analysis of the resulting photoelectron momentum distributions for these targets, facilitated by the quantum-orbits arising from the solutions to the saddle-point equations. The analysis yields a new class of rescattered trajectories that includes the well-known laser-driven long and short trajectories, along with novel Coulomb-driven rescattered trajectories. By virtue of the precision of the model, it opens the door to detailed investigations of a plethora of strong-field phenomena such as photoelectron holography, laser-induced electron diffraction and high-order above threshold ionization.


I. INTRODUCTION
Strong-field physics is the study of intense and short laser fields interacting with matter, which has enabled the measurement and manipulation of electrons on an attosecond timescale (10 −18 s).The strong-field process high-harmonic generation (HHG) [1,2] has enabled the production of ultra-short attosecond laser pulses [3,4] and given rise to the field of attosecond physics [5][6][7][8].The process of above-threshold ionization (ATI) [9] (specifically strong-field ionization [10]) has been used in ultrafast imaging procedures, e.g., laser-induced electron diffraction (LIED) [11][12][13][14][15], and photoelectron holography [16][17][18].All of these processes and applications are underpinned by recollision or recombination of the photoelectron.As such, all the above advancements and applications were enabled by breakthroughs in understanding and modelling the photoelectron's return.
The strong-field approximation (SFA) [19][20][21] (see Ref. [22] for a historical overview) and the three-step model [1], enabled significant progress in the description of HHG [2] and ATI [23,24], through the understanding of laser-driven returns of the photoelectron.Recollision models have enabled simple quantification of the electron's maximum return energy [1] (3.17U p ) in terms of the ponderomotive energy (U p ), the quiver energy of a free electron in a laser field.This in turn led to a simple equation for the HHG photon energy cut-off [2] (I p +3.17U p ), I p being the ionization potential of the target, as well as the cut-off for the direct and rescattered photoelectron energy [23] (2U p and 10U p , respectively).It was found that HHG and high-energy ATI energy spectra could be well-described by two types of trajectories undergoing a laser-driven recollision with shorter or longer time before returning [25], known as the short and long trajectories.Such insight led directly to considerable experimental developments in HHG, see, e.g., Refs.[26][27][28].In the case of ATI, recollision also plays a crucial role in LIED [11,12], where laser-driven photoelectron recollision is used to image molecules, or photoelectron holography [16,17], where recolliding wavepackets interfere with those that do not recollide.Considerable success of the such recollision models can be attributed to their simplicity, e.g., the Coulomb-field was neglected in the continuum, often allowing an analytic description.
Experimental breakthroughs [29,30] measuring accurate two-or even three-dimensional momentum (see, e.g., Ref. [31]) distributions for the photoelectron in strongfield ionization (instead of photoelectron energy spectra), revealed that models that neglect the Coulomb field during continuum propagation could not reproduce many key features for linearly polarized fields.Following this, a number of new models were developed that could properly account for the effect of the Coulomb potential for the photoelectron in the continuum, see, e.g., [32][33][34][35][36][37][38][39], and Ref. [18] for a review.The majority of these models took a semiclassical approach, applying different approximations and expansion upon Feynman path integrals.With these approaches, a qualitative description of strong-field ionization was possible, while retaining some descriptive power of recolliding trajectories.
Semiclassical models describing strong-field ionization fall into two categories [18,40], 'forward' methods that propagate many trajectories, binning the final momenta, or 'inverse' methods that find the specific trajectories that lead to each final momentum point.The benefit of forward methods is that it is easy to implement and generalize, however, it requires a huge number of trajectories to converge, which makes analysis of trajectories particularly difficult, as well as being numerically challenging.The inverse method provides much easier and clean trajectory-based analysis, while also being much faster to compute.However, the downside is that the inverse method is more difficult to implement and has been less generalizable.The Coulomb quantum-orbit strong-field approximation (CQSFA) [37,39,41,42] is the primary inverse method that has provided a powerful trajectory-based description of strong-field ionization [43][44][45][46].However, until recently, this method was only capable of qualitative agreement with solutions of the time-dependent Shrödinger equation (TDSE).The primary reasons were that important recolliding trajectories were missing, it lacked crucial phases known as the Maslov phase, and key features like pulse envelopes and solving in full three-dimensions instead of two-dimensions were missing.The Maslov phases had been neglected in all semiclassical treatments of strong-field ionization, until its recent inclusion in Ref. [47].The Maslov phase is a geometric winding number [48] that has an analogy with the Gouy phase in optics [49].
In this work, we address all of these issues, presenting a range of developments in the CQSFA that enables a high-level of quantitative agreement with a TDSE solution.Namely, using a similar methodology to [47], we implement the Maslov phases, using an adaptive Monte-Carlo-style momentum sampling methods we capture all valid trajectories, we implement a sin-squared pulse envelope, and provide solutions in three dimensions.To benchmark these changes, we compare the CQSFA to TDSE solutions for three different atomic targets and obtain exceptional agreement.Then, using a trajectorybased analysis, we identify a new class of rescattered solutions that capture all important recollision dynamics.This includes the well-known long and short laser-driven recolliding trajectories, but also trajectories whose return is dominated by the Coulomb potential, and cases that fit in-between these two extremes.We directly link these trajectories to the final momentum distributions and initial momentum, and use this to devise further improvements to the sampling algorithm.
This paper is organized as follows.In Sec.II, we present the theory for the CQSFA, giving a brief overview in Sec.II A, and proving more detail for the newest developments, the calculation of Maslov phases (Sec.II B), the momentum sampling algorithm in Sec.II C, while in Sec.II D we define the CQSFA orbits.In Sec.III, we present our main results, providing a comparison of the CQSFA with the TDSE for three different atomic targets, and providing an in-depth analysis of the ionization times, trajectories, and their contribution to the photoelectron momentum distributions.Finally, in Sec.IV, we present our conclusions.

II. THEORY
In this section, a short deviation of the CQSFA is presented in order to highlight some small deviations from the original formulation [37,39].The following section provides some details on the calculation of the Maslov phase [47].Throughout, we will work in atomic units (a.u.) unless otherwise stated.
A. The CQSFA Our starting point is the S-matrix amplitude [37] describing the transition from initial system state |ψ 0 ⟩, to a final scattering state |ψ p f ⟩ with asymptotic momentum p f .The Hamiltonian has been separated into a field-free part H 0 and time-dependent interaction part H I (t) and U (t, t 0 ) is the time-evolution operator for the entire system.We consider atomic systems within the singleactive electron approximation (SAE), such that the fieldfree Hamiltonian takes the form where p is the canonical momentum of the electron and V (r) is the effective SAE potential.For all targets considered in this text, the potential is of the form where Z is the charge of the parent ion and with coefficients from Ref. [50].We work exclusively within the electric dipole approximation, for which the interaction Hamiltonian in the length gauge is given by where the electric field is related to the vector potential To proceed, we approximate the final continuum eigenstate as a field-dressed plane wave |ψ p f (t)⟩ ≃ |p f + A(t)⟩ ≡ | pf (t)⟩, and insert the identity in Eq. ( 1) with I p the ionization potential of the target, and the ionization matrix element where |ψ 0 ⟩ is the time-independent initial state.In this way, we capture the time-evolution following ionization using the propagator ⟨ pf (t)|U (t, t ′ )| p0 (t ′ )⟩ from an initial dressed momentum state | p0 (t ′ )⟩ to a final dressed momentum state | pf (t)⟩.The semiclassical nature of the CQSFA stems from approximating this propagator by a semiclassical expression.
In the CQSFA, the propagator in Eq. ( 8) is written in terms of path integrals, to which semiclassical approximations can be applied.This is done by time slicing, where U (t, t ′ ) is split into N + 1 operators using its composition property, each propagating a time ∆t.Between these operators the identity is resolved in dressed momentum states and position states, allowing one to write [51] where (10) The limit N → ∞, is now taken, in order to turn the product of integrals into path integrals.
The above formulation used the transition matrix element in a momentum representation, propagating from an initial momentum p0 to a final pf .It is also possible to work in a 'mixed' representation, where propagation instead is from an initial position r 0 to a final momentum [38].The two forms of the matrix elements can be related as (11) which, together with Eq. (10), shows that the 'sliced' action in a mixed representation includes an additional boundary term Returning to Eq. ( 8) we see that if we pull out an exponential dependence of d, Eq. ( 9), we can write the transition amplitude in the mixed representation now including the p0 integral in the product.In the limit N → ∞, we obtain the path integrals with the mixed representation action and Hamiltonian The fact that we have an action corresponding to a mixed representation is the key point of this section.The effect of this is that when we apply the saddle-point approximation (SPA) to the path integrals, in order to obtain a semiclassical description, we obtain a stability matrix of the form ∂p(t)/∂p 0 , instead of the ∂p(t)/∂r 0 that is obtained in a momentum representation [38,52].Next, we perform the SPA on the time and path integrals [37,51], thereby obtaining the desired transition amplitude where ν s is the Maslov index, to be properly introduced in the following section.The sum over s is over all paths connecting the final and initial momentum and fulfilling the classical equations of motion, which follow from stationary action in time, position and momentum [37] respectively.The subscript s indicates that a quantity is connected to the classical trajectory Eqs. ( 20)-( 21) given by initial momentum and time, p 0 and t s , related through Eq. (19).In this way s functions as a collective index s = (p 0 , t s ) of the parameters uniquely identifying a given trajectory (Eq.( 19) can have several solutions for a given p 0 ).Equations ( 19)-( 21) are referred to as the saddle-point equations (SPEs).The t s is usually related to the time of ionization, with the first SPE, Eq. ( 19), expressing energy conservation at ionization.For a vanishing potential p s (τ ) becomes constant, and Eq. ( 19) reduces to the SPE from the lowest order SFA [37,39], providing a link between the ionization time and final momentum of the electron.This is in agreement with the fact that the potential is neglected in the lowest order SFA.When considering SFA trajectories in the following, we are thus considering solutions to Eq. ( 19) and (20) for a constant momentum value, p 0 = p f , with p f the final momentum value of the photoelectron.
J s in Eq. ( 18) is the Jacobian determinant in the mixed representation in general evaluated at some time τ along the classical trajectory, Eqs. ( 20) and (21).Note that in Eq. ( 18) J s (t) is to be evaluated at asymptotic final time t, where the field is over and the electron far from the parent ion.The Jacobian matrix is also referred to as the stability matrix, since it describes the trajectory's sensitivity to the initial momentum value.At points along a trajectory where J s (τ ) = 0, so-called focal points, there is an insensitivity to changes in initial conditions of the trajectory in some direction.Effectively, this means that several trajectories (with different initial conditions) will pass through the same point, and we thus have a bunching of trajectories at the focal points.If the final momentum of a trajectory happens to be in a focal point, this bunching effect will be apparent in the final transition amplitude Eq. ( 18) and cause a breakdown of the semiclassical approximation in those points [53].Surfaces in the final momentum distribution containing these focal points are known as 'caustics' and are clearly visible in calculated photoelectron momentum distributions (PMDs) as a line of high transition amplitude.The ionization times in Eq. ( 18) are generally complex t s = t Re s + it Im s , so the integral of the action, Eq. ( 16), becomes a complex contour integral.We pick a contour composed of two parts; first we integrate from t s to the real time axis, parallel to the imaginary time axis.This part is interpreted to correspond to the tunneling of the electron out of the potential.During the integration of this first contour, the momentum is assumed to be constant.Thereafter, the action is integrated from t Re s to infinity, corresponding to propagation in real time.In total, the action, Eq. ( 16), is thus split into tunneling and propagation parts [54].In the tunneling part of the action, an integral over the potential appears.In practice this will be replaced by a Coulomb correction factor [47,55,56]: where κ = 2I p .The reasoning for this substitution comes from Perelomov-Popov-Terent'ev (PPT) theory [55,56], which is widely used to explain the ionization rates of atoms.Specifically, by regularizing the integral and expanding in the path to second order, one obtains a finite integral, which resembles the expression that has been substituted.A similar result was also recently derived in a CQSFA setting in Ref. [40].In order to propagate t → ∞, after the pulse is over and the electron is sufficiently far from the residual ion, we analytically propagate the equations of motion and compute the action, using the same procedure as outlined in [38].
In order to calculate the ionization matrix element, Eq. ( 9), we use the asymptotic expansion of a Coulomb wave function where Y m ℓ (θ, ϕ) is a spherical harmonic, κ = 2I p and C ℓm an expansion coefficient.This expansion has previously been shown to work well in the length gauge [57], also when comparing with experimental data for molecules [58].Following the regularization procedure outlined in Ref. [57], we find the matrix element evaluated in the SFA time saddle points, t s , to be with the abbreviations The constant of proportionality is neglected, since we renormalize the simulation results.For the simulations performed here, we only employ one spherical harmonic in our initial state expansion, choosing C 0,0 = 1 for an s-state (hydrogen and helium) and C 1,0 = 1 for a p-state (argon), aligned along the laser polarization direction.
Finally, we assume that the electron starts from r s (t s ) = 0 at the start of the tunneling, in which case we can simply set the phase factor r 0 • p0 in Eq. ( 13) to zero, and neglect the boundary term in the action, i.e., neglect the last term in Eq. ( 16).

B. The Maslov Phase
The full CQSFA transition amplitude, Eq. ( 18), contains the phase-factor e −i πνs 2 known as the Maslov phase with ν s the Maslov index.This phase factor has not been included in previous versions of the CQSFA, but as we will show it is paramount for achieving excellent agreement with the TDSE with a semiclassical model for ATI.Recently, the Maslov phase was included in a semiclassical strong-field model [47], similar to the CQSFA.In this section, we elaborate on this procedure, based on Refs.[53,59].
First, a brief note on what the Maslov index represents.While its consequence is readily tractable from Eq. ( 18), changing the phase of the transition amplitude of each trajectory, the understanding of when it arises and how to calculate it is not so simply understood.There, however, exists an analogy from optics [47]: When focusing a Gaussian beam, it asymptotically obtains a π phaseshift when passing through its focus, known as the Gouy phase shift.Similarly, the ionic potential may act as a focusing lens and focus a manifold of trajectories with infinitesimally small deviations in initial conditions to a single point.When and where these focal points are in space and time depend upon both the ionic potential, the properties of the laser-electron interaction potential, and consequently also on each trajectory's initial conditions.At these focal points, for which J s (τ ) = 0, the Maslov index may instantaneously change its value in integer steps.The Maslov phase is hence a time-dependent function, ν s = ν s (τ ), and the contributions of each focal point along the trajectory must be calculated.
To calculate J s (τ ), we require the equations of motion of the stability matrices, which are in the form of a Jacobi initial-value problem.The derivation of the equations of motion and their initial conditions are detailed in Appendix A. We find that the time-evolution of the stability matrices in the mixed representation are governed by the following coupled matrix-differential equations d dτ where H[V (r s (τ ))] denotes the Hessian of the potential, which, as noted, must be evaluated along the classical trajectories.The initial conditions of the stability matrices read For Eq. ( 4) the Hessian reads where 1 denotes the 3 × 3 identity matrix and ⊗ the tensor product.Propagation of Eqs. ( 28) and ( 29) along with the SPEs, Eqs. ( 19)-( 21), then allows determination of J s (τ ) through Eq. ( 22).
With J s (τ ) obtained for any time τ along the trajectory, we can determine the times τ i for which J s (τ i ) = 0, that is, the focal points.Since ν s only changes in these points, its value at the final time t can be written as the sum of these changes at the focal times where ν s,0 is the initial value at the start of the trajectory, ∆ν s (τ i ) the change at a given focal point (to be determined below), and the sum runs over all focal points τ i ∈ [Re(t s ), t].The initial value ν s,0 is only nonzero if the trajectory happens to start upon a focal point [59].When ν s (t) is to be used in the final transition amplitude, Eq. ( 18), t is again to be seen as an asymptotic time, t → ∞, but practically large enough that the laser field is over, and the photoelectron is far from its parent ion.At the focal points, the determinant of the stability matrix will be zero, J s (τ i ) = 0, and consequently ∂p s (τ i )/∂p 0 will have some number, m, of linearly independent eigenvectors with corresponding eigenvalues of 0, denoted d (k) , that is, m solutions to the eigenvalue problem with λ k = 0 for k = 1, . . ., m.Given the d (k) eigenvectors at these times, the change in Maslov index at a focal point ∆ν s (τ i ) is given by the following formula where Π(τ i ) is a matrix with the following components where we have defined From our experience, setting ν s,0 = 0 for all cases will not notably change the result due to the occurrence of focal points at the start of the propagation being extremely rare.

C. The inversion problem
In order to calculate the transition amplitude, Eq. ( 18), we sum over all saddle-point solutions fulfilling Eq. ( 19)- (21).Given that this transition amplitude is formulated via a mixed representation path integral, this means we have an initial condition r s (t s ) = 0 and a boundary condition p s (τ → ∞) → p f .Thus, to find solutions, one must find all the initial momenta p 0 that lead to a particular final momentum p f .This is known as the inversion problem.Due to the difficulty of finding all of these solutions, many models take an alternative approach, where many starting trajectories are forward propagated and the final momentum binned, so that trajectories that end in the same final bin are summed coherently.To resolve interference, this requires very small bins and typically in excess of a billion trajectories.Additionally, it has also been shown in Ref. [60] that these 'forward' approaches do not yield the correct sampling weight in terms of the Jacobian J s , leading to 1/|J s | instead of the correct 1/ |J s | computed by inverse approaches (see Eq. ( 18)).
In earlier iterations of the CQSFA [39,41,42], not all solutions to the SPEs were included.This method assumed that there were only 4 solution per laser cycle and used an iterative procedure to find them.Later, additional solutions were identified [42] but were not incorporated into the CQSFA.The approach taken here, is to use a Monte-Carlo-style sampling method, to find all the initial momenta that lead to any particular final momentum, see Fig. 1.This works by generating a partially random guess for an initial momentum p 0,g , which (propagating via the SPEs) leads to a particular final momentum p f,g (p 0,g ).Then, in analogy to the 'forward' method, we bin this solution, associating it with the closest final momentum grid point p f,j to p f,g (p 0,g ), where j denotes the index of the grid.However, now this method differs from the 'forward' method, as we change the initial momentum to make the resulting final momentum approach that of the grid point.More precisely, the initial guess momentum p 0,g is used as a starting point to find a slightly different initial momentum p ′ 0,g so that the final guess momentum and the grid point momentum are essentially equal, by solving |p f,g (p ′ 0,g ) − p f,j | < ϵ, where ϵ is taken to be a very small value of the order ϵ ∼ 10 −12 .When sampling several initial guesses, care must be taken that any new initial guesses p 0,g do not lead to a solution that has already been binned.
If enough initial momenta are provided by a random number generator over a large enough range of p 0 , all unique solutions can be found in this way.However, this algorithm can be very inefficient if solutions are clustered very densely in small areas of p 0 momentum space, which is the case for rescattered trajectories.An example of the set of all initial momentum obtained for a Helium simulation is given in Fig. 1.The rescattered trajectories have a negative initial momentum coordinate perpendicular to the laser field polarization p 0⊥ < 0, which form densely clustered points that are more difficult to sample.
To address this, we implemented an algorithm that can detect narrow regions of densely clustered solutions, and then sampling can be more concentrated in these regions.These regions can be detected after performing an initial sampling over the full region depicted in Fig. 1.Then a grid is placed over the region p 0⊥ < 0, which is coarse enough so that if a dense cluster of points intersects a box in the grid, the initial sampling will have found at least one solution in this box.The clusters can then be intensively searched by restricting guesses only to boxes in the grid with at least one solution already present.
A similar approach was recently implemented in Ref. [40], however, in this case the sampling was done with a Gaussian distribution.This is efficient for the direct trajectories that do not undergo rescattering, but less so for the rescattered trajectories.Ref. [40] also did not consider a pulse envelope or the Maslov phases that we consider here.

D. Orbit types
In previous versions of the CQSFA the division of the quantum orbits into 4 types [61] was explicitly used in the simulations [37,42].For a laser linearly polarized along the z direction, the classification compares the tunnel exit and the final momentum along the polarization direction, p f ∥ , to the final and initial perpendicular momentum p f ⊥ and p 0⊥ , through the signs of the products r 0∥ p f ∥ and p f ⊥ p 0⊥ , as summarized in the legend of Fig. 2. We briefly describe the 'standard' trajectories of each orbit type, also shown on the same figure.Orbit type 1 (r 0∥ p f ∥ > 0, p f ⊥ p 0⊥ > 0) corresponds to direct trajectories, not revisiting the ion core.Type 2 (r 0∥ p f ∥ < 0, p f ⊥ p 0⊥ > 0) describes trajectories driven past the core by the field, but without much interaction.Type 3 (r 0∥ p f ∥ < 0, p f ⊥ p 0⊥ < 0) corresponds to trajectories performing a field-driven forward scattering of the ion core.Last, type 4 (r 0∥ p f ∥ > 0, p f ⊥ p 0⊥ < 0) describes trajectories performing a field-driven backwards scattering of the core.
Orbit types 1 and 2 only interact weakly with the ion potential, and can thus be related to the trajectories found in the lowest order SFA (see the text following Eq.( 21)).
On the other hand, orbit type 3 and 4 require the presence of the Coulomb potential.Since the semiclassical trajectories here are restricted to a plane we use a 2D definition of the orbit types, but a simple generalization to 3D trajectories is possible [62].We note that other trajectories than these 'standard' versions exist, such as multi-pass trajectories, having multiple close encounters with the ion core, and directly re-colliding trajectories, where the ionization happens at low fields and the scattering is primarily driven by Coulomb forces.While these were neglected in previous CQSFA versions [42], the current version of the CQSFA will include all trajectory dynamics, assuming sufficient sampling.

III. RESULTS AND DISCUSSION
In this section, we present and discuss results obtained with the model and compare to calculations performed using the freely available time-dependent Schrödinger equation (TDSE) solver, Qprop [63].Qprop works within the electric dipole and SAE approximation, and can calculate photoelectron spectra using the i-SURFV method.The radial step size and time step used in simulations were adjusted to ensure the correct ground state energy in imaginary time propagation and convergence of the PMDs, down to ∆r = 0.01 a.u. and ∆t = 0.0025 a.u.for argon.The size of the angular grid was chosen as r max = 100 a.u. with the t-SURFF boundary at 300 a.u.Angular momenta up to ℓ = 100 were considered for argon and hydrogen, and up to ℓ = 150 for helium.
For all calculations, both the TDSE and CQSFA, the coefficients for the effective SAE potentials for Eq. ( 5) are taken from Ref. [50].Throughout we take the vector potential to be a linear sin 2 -field of the form where ẑ is a unit vector in the z polarization direction, U p is the ponderomotive potential, ω is the carrier frequency of the field, N c is the number of cycles and φ is the carrier-envelope phase (CEP).All the results in the following sections will be for a 800 nm (ω = 0.057 a.u.) N c = 2 or 4 cycle pulse, with φ = π/2.The intensity was varied for each target, using 1.3 × 10 14 W/cm 2 (U p = 0.29 a.u., γ K = 0.94) for hydrogen, 4×10 14 W/cm 2 (U p = 0.88 a.u., γ K = 0.72) for helium and 2 × 10 14 W/cm 2 (U p = 0.44 a.u., γ K = 0.81) for argon, to ensure low Keldysh parameters, γ K = I p /2U p , without having over the barrier dynamics.
The systems here all have cylindrical symmetry.Therefore, it is sufficient to consider a 2D slice of the results, so we consider only results in the zy-plane for x = 0. Furthermore, the symmetry causes all the semiclassical trajectories to be confined within a plane, meaning sampling over initial momenta also only has to be performed in the zy-plane.We note, however, that the calculations are performed with a 3D model which has implications for determining the correct value of J s (τ ), Eq. ( 22), [47] and also allows for calculations for systems without cylindrical symmetry.In the following sections we will refer to the y and z coordinates as r ⊥ and r ∥ , the coordinate perpendicular and parallel to the laser polarization respectively.Similarly, we denote the parallel (p z ) and perpendicular (p y ) momentum components p ∥ and p ⊥ .Unless otherwise stated, the momenta used in the following will all refer to the final momenta of the semiclassical trajectories.

A. Validation of the model
In Fig. 3, we compare CQSFA and TDSE simulations performed for hydrogen, helium, and argon.The simulation parameters for helium are the same as used in Ref. [47], in order to enable direct comparison.For all targets, we see a very good agreement between the TDSE and the CQSFA with the Maslov phase included, Fig. 3(a)-(c), in the lower energy region before the presence of caustics.As mentioned in the section following Eq.( 22), the caustics are surfaces containing focal points for the classical trajectories at the final asymptotic time.Because of the breakdown of the semiclassical model in these points, the caustics are generally visible in the PMDs as sharp lines of discontinuity in the signal, e.g. in Fig. 3(a) at p ∼ (0, −0.8) a.u.
When we include the Maslov phase into the CQSFA, it is evident that especially the location of the interference fringes along p ∥ at p ⊥ ∼ 0 are very well-reproduced.The fact that the Maslov phase is indeed necessary to obtain this agreement is illustrated in Fig. 3(d)-(f).Here, the CQSFA with and without the Maslov phase is shown on FIG. 3. Photoelectron momentum distributions for H, He and Ar subject to a two-cycle sin 2 -pulse with a carrier-envelope phase set to φ = π/2 [see Eq. ( 38)].The wavelength is set to 800 nm, corresponding to ω = 0.057, for all calculations.Due to the cylindrical symmetry, [H(t), L ∥ ] = 0, only half of the PMD is required.The horizontal dashed white lines in each subplot signify the boundary between the upper and lower parts.In the first row, we compare TDSE results (upper half) for H (a), He (b) and Ar (c) to the improved CQSFA, including the Maslov phase (lower half).In the second row, we then compare the CQSFA calculations with (without) the Maslov phase in the bottom half (top half) for H, He and Ar in (d), (e) and (f), respectively.The intensity of the laser is for each target given by 1.3 × 10 14 W/cm 2 for H, 4 × 10 14 W/cm 2 for He and 2 × 10 14 W/cm 2 for Ar.These parameters ensure that we are within the tunneling regime and the Keldysh parameter γK satisfying γK < 1 for all targets.The single colorbar is representative for all PMDs.
each half, and a shift of interference fringes is clearly seen along the p ⊥ = 0 axis, see panel (e) box 1 for He.
Another important structure affected by the inclusion of the Maslov phase is the fringes at p ∥ > 0 almost parallel to the p ∥ axis, commonly known as the 'spider-legs' [18].This interference pattern is found for a wide range of targets, and is visible in all the PMDs in Fig. 3, most clearly for He.As seen on panels (d)-(f) the Maslov phase results in a 'lift' of the central fringe of the spider (panel (e) box 2 for He), increasing its width to a size in much better agreement with the TDSE results.This effect is further illustrated in Figs.4(b), (d), and (f), showing the signal along a slice of the PMDs of Fig. 3 for p ⊥ = 0.1 a.u.While the CQSFA results both seem to agree with the TDSE simulations for low p ∥ , the inclusion of the Maslov phase clearly creates a difference for p ∥ > 0.5 a.u., increasing agreement with the TDSE results significantly.In Figs.4(a), (c), and (e), a slice of the signal along p ∥ = 0.4 a.u.also illustrates how the location of the other spider legs are affected by the Maslov phase.For all targets the agreement is better with the inclusion of the Maslov phase, perhaps most clearly for hydrogen, where quite large shifts in the interference minimum can be seen (compare e.g. the gray dashed line at p ⊥ = 0.25 a.u. to the blue at p ⊥ = 0.30 a.u. in Fig. 4(a)).This illustrates the necessity to take the Maslov phase into account, if one wants to perform analysis based on the spider structure.
Still, the CQSFA and TDSE results do not agree completely.One particular region where the results seem to disagree, is for high values along the p ∥ axis, where the CQSFA results drop off too fast.This is clearly illustrated for both hydrogen and helium in panels (b) and (d) in Fig. 4, and can also be seen in the PMDs of Fig. 3.The disagreement seems to be worst for hydrogen.The reason for this could be the fact that the hydrogen simulation is performed for the lowest intensity of all the simulations, and in connection also the simulation with the highest Keldysh parameter (γ K = 0.94).To test the range of validity of the CQSFA, we performed several simulations with values of Keldysh parameters in the range [0.7, 1.1].In general, testing suggests that the CQSFA becomes better for longer wavelengths [62], higher intensities and lower Keldysh parameters, a trend in conformity with the well performing helium simulation which has the highest intensity and lowest Keldysh parameter (γ K = 0.72).The CQSFA also performs well for intensities corresponding to over-the-barrier dynamics.
Another source of discrepancy between the CQSFA and TDSE results is, as mentioned, the caustics.The effect of those are worse for the hydrogen and argon simulations, as seen in panel (a) and (c) of Fig. 3, where they cut off most of the structure at high energies, an effect also visible as a sharp drop-off on panels (a) and (f) of Fig. 4. Since the caustics represent classical turning points, their position can be pushed outwards from the center by increasing, for example, the intensity or the number of cycles of the pulse.
The results also show that the argon simulation is less converged than the hydrogen and helium simulations, see Figs. 4(c) and (f), where the blue and gray lines showing the Ar CQSFA results are clearly less smooth than the corresponding curves for He and H, Figs.4(a)-(d).The reason for this is that the effective potential used for argon gives rise to many additional (and seemingly irrelevant) solutions to the SPEs, meaning significantly more sampling of orbits is required to reach convergence.It is therefore worth noting, that the inclusion of f (r) in Eq. ( 4) is not critical for the current simulations.Just using the plain Coulomb potential (but the correct I p for the targets) yields results almost identical to those in Fig. 3, with some small difference in the signal value for the high energy rings, while convergence is reached with less sampling.It could be, however, that the effective potential plays a larger role for other systems or field parameters, such as fields with more cycles where higher energies are reached.
Overall, the large degree of agreement between the CQSFA results and the TDSE simulations shown in Fig. 3 illustrates the accuracy and validity of the semiclassical model.

B. Trajectory dynamics
The high level of accuracy of the CQSFA adds a new level of credibility to orbit-based analyses.Given how well the semiclassical model reproduces the TDSE results, it seems likely that the sampled trajectories and their properties could shed light onto key aspects of the underlying physics.Studying the different types of trajectories found in the new model, along with their impact on the PMDs, could thus lead to an improved understanding of the ATI process.
To investigate the different types of trajectories and their effects on the transition amplitude, we plot in Fig. 5(a1) the ionization time as a function of final parallel momentum for orbits sampled in the helium simulation shown in Fig. 3(b), at p ⊥ = 0.25 a.u.(the CQSFA result in Fig. 3(b) is shown for p ⊥ < 0, but because of the cylindrical symmetry the result for p ⊥ > 0 is identical).Each vertical slice in the plot shows all saddle point times for a particular final momentum point in the PMD.The corresponding SFA solutions (see the text following Eq.( 21) for definition) are shown in dashed black lines.Clearly, the majority of orbit 1 (blue) and 2 (orange) solutions lie along the SFA solutions, indicating their close resemblance to the SFA trajectories.Between each pair of SFA 'lines' one sees additional structure composed of orbit type 3 (green) and 4 (red).These structures contain the rescattering trajectories, which will be studied in detail.
The rescattering structure with the most impact on the transition amplitude is the framed structure in Fig. 5(a1), containing orbits with ionization times near the highest peak of the electric field.A zoom-in of this structure is seen in panel (b1) and a collection of corresponding trajectories displayed in panel (b2), to show the evolution in trajectory dynamics that the structure represents.The lowest orbit 4 line, labeled by box A, contains field-driven scattering trajectories; the trajectory has a turning point at the maximum of the electric field, then it returns to the ionic core and scatters after spending around one laser cycle in the field.These orbits correspond to the 'standard' orbit 4 in Fig. 2, and were used in previous CQSFA models (see, e.g., Fig. 16 in Ref. [18]).Box B marks a turning point in the orbit 4 structure.Here the field maximum occurs after the trajectories turning points and the return to the ionic core now happens within the same cycle as ionization.In this way, the orbits from boxes A and B are similar to the 'long' and 'short' orbit from the rescattered SFA [64] and HHG [25].As we move to ionization times further from the electric field maximum, these short orbits behave less like the SFA orbits, as the Coulomb interactions become more dominant, making the trajectory turning point happen sooner compared to the field maximum, see the trajectory from box C. At box D the scattering happens before the next field maximum, while at box E we have rescattering before the laser amplitude changes sign.Interestingly, and quite unlike normal SFA orbits, these orbits are ionized at times approaching 0 electric field.For the same reason, however, their contribution to the transition amplitude is small.It is possible that these trajectories play a more prominent role for more complex targets or laser fields.In such cases, they could be an interesting probe of the system potential with which they interact strongly.
It is worth noting that other trajectory types than the ones just discussed are present in Figs.5(a1) and (b1).For example, the orbit 1 and 4 structure found between box C and E on Fig. 5(b1) represents multipass trajectories, rescattering with the parent ion several times.Even 'chaotic' trajectories, where the electron moves in a complicated and unpredictable manner, can be found in the less connected parts of Fig. 5(a1), primarily for lower times (Re(t s ) ∈ [0, 90] a.u.), since these trajectories needs time interacting with the field to develop.Due to their extreme sensitivity to initial conditions, however, these 'chaotic' trajectories should have a large value of J s , Eq. ( 22) and in turn a small contribution to the final transition amplitude.
We also note that the orbit structure shown in Fig. 5 can become more complicated for other potentials.One example is the argon potential used to obtain the results in Fig. 3(c) and (f), where a lot of additional solutions become possible, complicating the sampling as previously mentioned.This however also shows the possibility for additional types of trajectory dynamics to be involved in the description of more complicated targets.

C. Orbit-resolved PMDs
In this section, we present and analyze PMDs made from only one or two orbit types.This is the first time that such PMDs have been produced with the CQSFA for few cycle pulses.As will be seen, the structures in the orbit-resolved PMDs can be directly related to Fig. 5 panel (a1) and especially (b1).This comparison is aided by the fact that the turning points in the orbit structure (e.g.box B in panel (b1)) correspond to a point on a caustic in the PMDs.This is the case since Re(t s ) is parameterized by the initial momentum p 0 , see Eq. (19).Thus, a turning point in Re(t s ) can be expressed as Hence, a focal point, J s (t) = 0, and a caustic, will lead to a divergent turning point in Re(t s ), as can be observed in Fig. 5(b1) box B. In Fig. 6, the PMDs for orbit types 3 and 4 for a 2cycle and 4-cycle pulse is shown.As we already saw in Fig. 5, the somewhat arbitrary definition of the orbit types, mean that orbit types 4 and 3 changes into the other at p ∥ = 0.In order to get a continuous PMD, we stitch together the contributions of the two orbit types on each half of the momentum distribution about p ∥ , as indicated in the panels.Some of the features in Figs.6(a), (b) can be directly related to the distribution of Re(t s ) and final parallel momentum p ∥ , Fig. 5(a1), with the high signal in Fig. 6(b) originating from the orbits in the central 'branch', Re(t s ) ∈ [90, 140] a.u. as shown in Fig. 5(b1), and the lower signal in (a) originating from the branch below, Re(t s ) ∈ [50, 90] a.u.These are the only two orbit branches expected to correspond to a strong signal, since only the photoelectrons that correspond to these solutions enter the continuum at high field values, while also having adequate time to subsequently propagate in the electric field.
We see that the momentum value where the caustic appears on the left side of Fig. 6 The difference in transition probability in panels (a) and (b) of Fig. 6 is clearly an effect of the two-cycle pulse used.Indeed, in PMDs for a 4-cycle pulse, as can be seen on panels (d) and (e) of the same figure, the value of transition amplitude for a given orbit type is about the same for both p ∥ > 0 and p ∥ < 0. The effect of 2 extra field cycles is largest for orbit 4, where the high-energy interference rings are now present for both positive and negative p ∥ .Furthermore, one notices an additional caustic within the orbit 4 signal in panel (d) (starting at p ∼ (1.9, 0) a.u.) also carrying over to the orbit 3 signal.This caustic is similar to the 2-cycle caustic in panel (b), and is indeed linked to a structure similar to that in Fig. 5(b1), but for the last parts of the 4-cycle pulse.The additional fine interference fringes found within the boundary of the caustic (compared to the 2-cycle case) thus represents interference between trajectories originating from two different cycles of the pulse.Other than these inter-cycle interference fringes, along with a shift of the caustic to higher p ⊥ values, the orbit 3 signal is not much affected by the increase in field cycles.
Figure 6 (c) and (f) shows the total signal from orbits 3 and 4, that is, the transition amplitude, Eq. ( 18), for orbit 3 and 4 coherently summed before taking the norm square.For the two cycle pulse, panel (c), we only have a relatively small amount of interference between the two orbits, since the difference in signal of the orbit types is large (compare Fig. 6 (a), (b)).The main interference effect is the emergence of fine interference fringes within the area of the left orbit 4 signal.For the 4 cycle pulse, Fig. 6(f), interference between orbit type 3 and 4 is more pronounced, since the difference in signal for the two orbit types is less.Especially, we see the emergence of the so-called spiral structure [65] in the central parts of the PMD.The results here show that patterns coming from interference with orbit 4 depends most strongly on the few-cycle nature of the pulse, which was to be expected from the single orbit type PMDs.

D. Connections to initial momenta
The orbits found within the structure in Fig. 5(b1) capture most of the orbit types 3 and 4 signal in the PMD.It is therefore interesting to investigate where in the initial momentum space these orbits originate from, and to optimize sampling of the initial conditions.In Fig. 7, we consider the initial momentum values for the orbits found in the central structure (Fig. 5(b1)) across all values of final momenta.Furthermore, the specific initial conditions for the orbits shown in Fig. 5(b1), with p ⊥ = 0.25 a.u., are marked with black dots.We also include the boxes A-E from Fig. 5(b2), now mapped to initial momentum space.Clearly, judging by the density in the boxes in Fig. 7, the initial conditions in box A, B, C, and D, that are giving rise to the high signal structure in the PMDs Fig. 6, are located within a small area, p 0∥ ∈ [0; 1.5] a.u., p 0⊥ ∈ [0; −0.2] a.u.The other parts of the figure, such as the orbits for p 0∥ > 2.5 a.u., containing box E, were seen to have very little influence on the PMDs, even though they are densely populated in initial conditions.In addition, we have plotted the initial momentum values for the p ⊥ = 1 a.u.orbits, as seen in Fig. 5(b1) in orange.Here, it is seen how the orbits giving the most signal in the PMDs (the ring structure close to box A, B, C) is even more localized in initial momentum space.
Comparing with the total amount of sampled initial conditions, Fig. 1, one sees that it is only a very small amount of the sampled trajectories that matters for the resulting PMDs.In particular, there are many densely populated areas of initial conditions that are low in probability, meaning a clustering approach to sampling quickly leads to unnecessary oversampling.For more complicated potentials, e.g., for molecular targets (and to some degree even for the Ar potential used here) this could make sufficient sampling very difficult.One way to counter this would be to sample in an iterative manner: First a set of random initial conditions are chosen, their trajectories propagated and the parameters important for their contribution to the transition amplitude, primarily Im{S( ps , r s , t s )} and the Jacobian of the stability matrix J s (t), see Eq. (18), are determined.Based on these parameters, one can estimate which regions of initial conditions should be sampled in greater detail, and prevent sampling in regions of low probability density.

IV. CONCLUSION AND OUTLOOK
In this paper we have presented an improved version of a semiclassical model for ATI, the CQSFA, and show-cased its capabilities in reproducing results from the TDSE and ability to disentangle the intricate electronic dynamics on an attosecond timescale.To systematically verify the accuracy of the model we have analyzed the PMDs for three atomic targets, hydrogen, helium and argon, ionized by a two-cycle sin 2 pulse of linearly polarized light, and found excellent agreement.We have fully included the Maslov phase in the calculation of the transition amplitudes, which drastically improves the agreement with TDSE solutions.In earlier versions of the CQSFA, one found deviations in the holographic patterns in the PMDs, which have now been reconciled.Secondly, through a new, Monte-Carlo style sampling of the initial conditions, we solve the inversion problem by randomly sampling the initial momenta and propagating the trajectories forward in time.With this approach, we find 'new' types of quantum trajectories, that were not previously included in the CQSFA.These make a new class of rescattered trajectories that go from laser-driven returns that ionize near the peak of the laser field (including the long and short solutions) to Coulomb-driven return that ionize near zero field.Third, we have also expanded the model to accommodate few-cycle pulses of arbitrary shape, which is crucial for comparisons to stateof-the-art few-cycle ATI experiments, and solve in three dimensions.
The high level of accuracy of the CQSFA yields a new level of credibility to trajectory-based analysis within the field of strong-field physics.We have illustrated that it is possible to extract the trajectories most relevant for atomic ionization and relate them directly to interference structures in the PMD, enabling further optimization of the sampling algorithm.This will enable the use of a trajectory-based analysis for more complicated targets, where the new trajectories highlighted could play a role.Furthermore, we imagine the CQSFA being used for strong-field simulations where direct TDSE solutions becomes unpractical, e.g., ionization of larger molecules, or longer wavelength, such as Ref. [62].Access to accurate PMDs and trajectory analysis could give insight into attosecond dynamics, and even function as a sensitive chiral probe [66].

FIG. 1 .
FIG. 1.Initial momentum values sampled for a simulation of helium ionized by a two cycle sin 2 laser pulse with intensity 4 × 10 14 W/cm 2 and wavelength of 800 nm.The darker color indicates high density of sampled points.The red box marks the region of final momentum values investigated in the simulation, i.e., all the shown initial conditions lead to a final momentum value within this box.Note the difference in scale between the p 0∥ and p 0⊥ axis.Figs.3(b) and (e) shows the corresponding PMD.

FIG. 2 .
FIG. 2. An example of the four different orbit types in the CQSFA for an asymptotic final momentum p f = (0.6, 0.25) a.u., taken from the helium simulation of Fig. 1.The arrowheads indicate the direction of travel of the photoelectrons and are separated by 0.1 laser pulse periods in time.The legend indicates the definition of each orbit type in terms of the sign of the product between initial position and final parallel momentum, and final and initial perpendicular momentum.The trajectories shown here are the ones usually captured in previous CQSFA models.

FIG. 4 .
FIG. 4. Transition amplitude signal for the TDSE simulations (orange lines) and the CQSFA with Maslov phase (blue line) and without (gray striped line), through different slices of the PMDs in Fig. 3. Left column [(a), (c), (e)] shows the signal for p ∥ = 0.4 (a.u.) for hydrogen, helium and argon respectively.The right column [(b), (d), (f)] shows the signal for p ⊥ = 0.1 (a.u.) for the same order of elements.The same laser parameters as in Fig. 3 are employed.All signals are normalized according to their respective maximum values in the plotted region.

FIG. 5 .
FIG. 5. (a1):Real time of the temporal saddle points as a function of final parallel momentum p ∥ for the orbits used in the helium simulation of Fig.3(2-cycle sin 2 pulse with intensity of 4 × 10 14 W/cm 2 and wavelength of 800 nm) at p ⊥ = 0.25 a.u.The different orbit types(1)(2)(3)(4) are plotted in different colors, and the dashed line indicates the ionization times obtained using standard SFA for the same final momentum values.(a2): Electric field shown as a function of the real part of the ionization time.(b1): Zoom-in of orbits marked with a gray box in panel (a1).The light gray orbits in this panel are orbit 3 and 4 for p ⊥ = 1 a.u. and are not related to the other panels.(b2): Real space trajectories of one orbit from each box in (b1).The trajectory color and line style match the corresponding box in the insert.Note that from box D an orbit 3 trajectory is plotted.The crosses mark points along the trajectories with electric field maxima.
(b) (p ∥ ∼ −1.75 a.u.), matches very well with the orbit 4 turning point in Fig. 5(b1) box B, indicating that the circular interference pattern within this caustic originates from the interference of the long and short re-scattered orbit 4 (Box A and C), which is well known.The orbit 3 turning point in Fig. 6(b1), close to Box D, on the other hand, does not result in a visible caustic in the orbit-resolved PMD in Fig. 6(b) on the right-hand side (at p ⊥ = 0.25 a.u.).This indicates that the upper part of this main structure (containing box E), associated with Coulombdriven recolliding orbits 3 and 4 solutions, correspond to a very small transition amplitude.This situation changes at p ⊥ ∼ 0.6 a.u., where a caustic appears on the orbit 3 side of Fig. 6(b) (p ∥ ∼ 1.75 a.u.).This is due to a slight change in the structure from Fig. 5(b1), where at p ⊥ > 0.6 a.u. the orbit 3 line near box D detaches from the orbit 3 solution above and instead connects to the orbit 3 solution below that lies along the SFA solution, forming a new loop structure together with the orbit 4 solutions in Box A, B, and C.This is shown by the light gray orbits in the same panel, indicating the orbit 3 and 4 structure for p ⊥ = 1 a.u.The formation of this ring structure increases the importance of the orbit 3 solutions around the right turning point, both leading to the caustic at p ∼ (1.8, 0.6) a.u. and the appearance of interference fringes above p ⊥ ∼ 0.6 a.u. on the orbit 3 side of Fig. 6(b).

FIG. 6 .
FIG. 6. Orbit resolved PMDs for He of orbit type 3 and 4. Panel [(a), (b), (c)] shows the results for a two-cycle pulse, parameters identical to that of Fig. 3. Panel [(d), (e), (f)] shows results for a 4-cycle pulse, with all other parameters identical.Panels [(a), (b), (d), (e)] shows the PMDs for the single orbits 3 and 4, with each half of a plot, separated by a dashed line, corresponding to the orbit type indicated by the number in the top corner.Panels [(c), (f)] shows the PMDs for the orbit 3 and 4 solutions combined, i.e., the orbit 3 and 4 signal is summed coherently before taking the norm square.