Disentangling Long Trajectory Contributions in Two-Colour High Harmonic Generation

This work investigates High Harmonic Generation (HHG) in gas targets, induced by intense two-colour driving fields. We compared classical trajectory Monte Carlo simulations based on a semiclassical model of strong field tunnel ionisation of helium to experimental HHG spectra, and investigated the relative contribution of long trajectories to low harmonic orders. This phenomenon can be found even when the experimental setup is chosen to suppress long trajectories and favour phase matching for short trajectories.


Introduction
The fundamental process of High Harmonic Generation (HHG) was first discovered in the late 1980s [1,2].Starting from ultrashort few-femtosecond (fs = 10 −15 s) driving pulses (typically in the visible or near-infrared wavelength domain), HHG is an extremely nonlinear frequency up-conversion process, creating Extreme Ultra-Violet (XUV) pulses with a pulse duration shorter than a femtosecond (the current world record is 43 attoseconds (as = 10 −18 s) [3]).Since then, HHG-created XUV pulses have found diverse applications in ultrafast science [4][5][6][7].First, the intuitive three-step-model of HHG was developed within a single-atom response [8][9][10], which introduced the concept of short and long trajectories in HHG.Then, the macroscopic medium was also taken into account including propagation and phase matching [11][12][13], which typically favours the short trajectories.Finally, the temporal coherence was explored [14], which showed two distinct spatial regions with different coherence times in the far-field pattern of the high harmonics depending on the two different trajectories.As the general understanding became better, HHG spectroscopy was also proposed [15,16], which combines sub-Ångström spatial with attosecond temporal resolution.Dudovich then proposed a new approach for HHG with a two-colour field in an orthogonal configuration that enables controlling and decoupling the angle of ionisation and the angle of recollision and therefore increasing the dimensionality of HHG spectroscopy [17,18].Here, we explore the contribution of the long trajectories and the parent ion interactions in two-colour HHG using Classical Trajectory Monte Carlo (CTMC) simulations.The work presented in this article focuses on HHG in helium gas.
The intuitive three-step-model of HHG was developed by [8][9][10].In the first, quantum mechanical step, a strong driving field tunnel ionises an electron, which is subsequently accelerated classically by the strong laser field.In the case of small ellipticity or linear polarisation, photoelectrons can be driven back into their parent ion by the laser field.These photoelectrons can then either rescatter, or recombine again, representing the third step.In the case of recombination, a high energy photon is emitted, where the frequency of the high harmonic is determined by energy conservation with ω n = n • ω, I p the ionisation potential and E kin the kinetic energy of the photoelectron when recombining.Atomic units are used throughout this manuscript, unless otherwise specified.The linearly polarised driving field is given by with envelope f (t), maximum field strength F 0 , optical frequency ω, and polarisation axis x.
The standard Strong Field Approximation (SFA) [19] neglecting the ion Coulomb field after ionisation, and dipole approximation [20] neglecting the spatial dependence of the laser field are used, as well as a single active electron assumption [21].
Depending on time t i when a photoelectron is liberated, the trajectories will return at a different time t r (t i ) and carry a different kinetic energy E t i , indicated by the colour code in Figure 1.For most energies, there are two trajectories, a long trajectory that starts soon after the peak of the field (greyed out in Figure 1), and a short trajectory that starts later in the cycle (full colours).The energy cutoff of the returning photoelectrons is provided by the single trajectory separating the domains of short and long trajectories.The cutoff energy is calculated to be which corresponds to the maximum possible electron kinetic energy at the time of recombination [22].Depending on the focus position relative to the gas jet for HHG, the contribution of the long trajectories in the HHG signal can be changed [14] and, with the right balance, quantum path interference between the short and long trajectory was observed for the first time in HHG [23] in excellent agreement with theory [24], and then extended to many other gas targets [25].During the derivation of this quantum model by Lewenstein et al. [22], several assumptions of the classical trajectory three-step-model are confirmed.The contributing trajectories must return to their parent ion again, which also imparts conditions on the initial momentum of the trajectory right after tunnel ionisation.Furthermore, in a recent experimental study, it was found that the spatial excursions of HHG trajectories are also well explained by the semiclassical model [26].Therefore, it is valid and helpful to investigate the HHG process in greater detail using Classical Trajectory Monte Carlo (CTMC) simulations.
Since the strong field tunnel ionisation rate declines exponentially with decreasing field strength [19], long trajectories are preferentially created.However, due to the spreading of the wave packet during the propagation, the probability of recombination decreases with increasing excursion time.In addition, long trajectories show more efficient phase matching in the converging beam, when the gas interaction region is placed before the focus of the driving laser beam [27].For short trajectories, on the other hand, the phase matching is more efficient in the diverging beam after the laser focus, and it is easier to achieve than for short trajectories.Optimizing the gas jet parameters after the beam focus, combined with the quantum spreading of the wave packet, one can work against the preference towards long trajectories from the ionisation rate.In practice, therefore, the optimised HHG spectrum consists almost exclusively of short trajectory contributions [18,27].

Two-Colour HHG
Here, we build on concepts first introduced in Ref. [18], and related to Refs.[28,29].Shafir et al. [18] employ a two-colour HHG method to extract the instant of time when electrons exit from their bound state and enter the continuum.Figure 2 illustrates the underlying concept.The strong fundamental driving field acts along the longitudinal coordinate x, as described in Equation (2).It dominates the ionisation rate and the major component of the photoelectron trajectories.A weak co-propagating second harmonic field with relative intensity I 2ω /I ω = ξ acts along the orthogonal direction y slightly displacing the trajectories, depending on the two-colour delay ϕ.
In Figure 2, examples of short trajectories are drawn, always with an initial transverse momentum v i,y = 0.For each set of initial conditions (t i , v i,y ), a specific two-colour delay ϕ(t i , v i,y ) can be found for which the transverse displacement cancels out.The trajectory then returns to the parent ion at return time t r and can contribute to the HHG spectrum.The transverse velocity v i,y = 0 is the most probable initial condition both in the adiabatic case [30,31] and in the non-adiabatic (almost) linear case [32,33].Therefore, the HHG signal of harmonic order n is maximised exactly for This displacement gate method as developed by the Dudovich group gives access to the initial and return times for different harmonic orders in the HHG spectrum [18,28,34].To perform the analysis and extract the times, Shafir et al. [18] neglected the effect of the ion Coulomb potential on the trajectories of the photoelectrons, according to SFA, allowing them to analytically describe the photoelectron trajectories.Then, combining the requirement that a trajectory must return to the origin after time t r with the extracted ϕ max from the measurement and setting v i,y = 0, one can solve for the initial and return times t i and t r [18].

CTMC Method
CTMC simulations offer a more direct access to the dynamics and trajectory parameters of photoelectrons that generate high harmonics.The absolute phase between the two colours ϕ is known, whereas, in the experiment, only a phase relative to an unknown constant shift can be directly determined [18].Through a comparison of the resulting HHG spectra with experimental data, the CTMC also offers the possibility to investigate the relative contributions of short and long trajectories to the high harmonic spectrum.
The initial conditions of the CTMC are based on the adiabatic Ammosov Delone Krainov (ADK) theory [30,31].A cloud of 150000 trajectories is launched during a quarter cycle of the electric field oscillation, starting at the peak of the field at the centre of the ultrashort pulse, until the zero-crossing of the field.The probability weight is defined by the ADK ionisation rate, and each trajectory receives a random initial transverse momentum (in the plane orthogonal to the instantaneous laser field direction) following the ADK momentum distribution [19,35,36].The longitudinal initial momentum is assumed to be zero, following a standard approach [19,32].The exit point for each trajectory was given by the radius calculated from a parabolic coordinates approach [35,[37][38][39], and taken exactly in the direction of the instantaneous laser field.The laser parameters are chosen to imitate the experiment on helium [18].As was done for all presented theories in Ref. [18], we calibrated the pulse amplitude such that the retrieved HHG cutoff energy matched with the experiment, resulting in F 0 = 0.1139 au.The fundamental IR pulse had a f (t) = cos 2 (Ωt) envelope for the field, Ω was chosen to yield a 6 fs FWHM in intensity, and the carrier wavelength was λ = 800 nm.
For the subsequent propagation, classical Newton's equations were solved in two versions, using the MATLAB (R2016a, MathWorks R Inc., Natick, MA, USA) native ode45 function.This 4th order Runge-Kutta scheme employs 5th order error estimates, and adapts the time-steps to fulfill a relative error tolerance of 10 −7 .The first version of the equations of motion are according to SFA, where the force term only includes the electric field.The second version propagates in the superposed forces of the electric field pulse, the ion Coulomb force and the field induced dipole on the ion.The latter corresponds to the TIPIS model (Tunnel ionisation in Parabolic coordinates with Induced dipole and Stark Shift) [39,40].More details of the CTMC implementations can be found in [35,41].
Purely classical trajectories only experience scattering off the Coulomb potential, but of course do not recombine during the numerical calculation.Thus, possible recombination events must be filtered out of the recorded trajectories a posteriori.From our CTMC simulations, we find the best qualitative agreement between the extracted kinetic energies and experimental HHG spectra if the recombination condition "returning to the individual exit radius" is used (see Figure 3A).These findings agree with Ref. [42], where Xiong et al. also conclude that a recombination distance equal to the exit point was necessary to yield the correct description for short trajectories in the near-threshold part of the HHG spectrum.When checking the kinetic energy at a closer specific radius, then short trajectories gain additional kinetic energy from being accelerated by the laser field while travelling closer to the ion (see Figures 3B and 4).Long trajectories do not experience this acceleration because they return to the ion when the field is already driving them away again, decelerating them slightly when they have to approach closer than their individual exit point.To illustrate this key difference between short and long trajectories, Figure 4 plots the radius over time for two exemplary trajectories.The red dotted curve is a long trajectory, which was ionised right when the electric field (blue dashed) reaches its maximum value, corresponding to the shortest possible exit radius (grey dot-dashed line) in the simulation.Compare also the sketch in Figure 1 for short and long trajectories within the SFA and three-step-model.The green dotted curve represents a short trajectory that was launched just shortly before the electric field changes sign.The comparatively small field strength this late in the quarter cycle results in a larger exit radius.Evidently, the short trajectory must first be accelerated towards the ion in order to reach a smaller recombination distance condition.If this radius is chosen as recombination condition R, then short trajectories travel a relatively large distance from their initially larger exit radius towards the ion nucleus until they cross this radius, while being accelerated by the laser field (blue dashed) that has already turned around.Long trajectories, however, can be decelerated by the field during their approach to the ion.

Macroscopic Propagation
Both the semiclassical three-step-model [8][9][10] as well as the quantum theory [22] describe the single atom response of HHG.In order to have efficient production of high-order harmonic radiation, phase matching conditions for constructive interference at harmonic order n, must be satisfied over the complete propagation distance of the driving laser field through the target interaction region.All new parameters and quantities are defined in Table 1.Equation ( 6) takes account of the Gouy-phase of a Gaussian beam, the phase of the electron wave packet (described by the classical action) and the phase due to the time passed until the recombination.On the other hand, it neglects dispersion effects induced by the neutral and ionised gas.
fundamental Gaussian beam S p (t r , t i ) classical action for electron wave packet propagation K(z, r) = ∇ nωt r − S p (t r , t i ) intrinsic atomic phase N total number of interacting targets along the axis j gas jet diameter Experimental setups are typically optimised for short trajectory phase matching on axis (see also Ref. [18] and their supplementary material).In that case, the phase matching condition for long trajectories is best fulfilled for emitted wave vectors diverging from the optical axis [27].These contributions to the final HHG spectrum can easily be filtered out spatially, as was done in Ref. [18].However, there still might be a small contribution of long trajectories to the low harmonics even when looking along the axis.The phase mismatch is then given by For an interaction region on axis (r = 0) in forward propagation along the axis, this mismatch can be estimated as where we used approximations described in Ref. [27].
Figure 5A shows the phase mismatch δk of the wave vectors for two sample excursion times τ short = 25 au and τ long = 80 au at harmonic order n = 22, and the gas jet positioned after the focus of the beam (z > 0).In an optimal position of the jet after the focus, long trajectories are very efficiently, but not completely, suppressed along the axis (corresponding to large |δk|).
Starting from the phase mismatch δk(z, n, τ) Equation ( 8), one can estimate the total amplitude efficiency with the help of a geometric series This geometric series oscillates around a value for low N, and converges for larger N 40.The limit of N → ∞ is given by A typical XUV pulse contains several orders of magnitude more photons than N = 40 per harmonic order, which justifies taking the limit of N → ∞.Then, the absolute square yields the efficiency for the intensity of the high-order harmonic radiation.For perfect phase matching the estimate yields perfect efficiency of constructive interference, as expected.Figure 5B shows the phase matching efficiency for varying positions of the gas jet after the focus, with a jet diameter of j = 0.95 mm [43].

CTMC Simulations
In this section, we compare our CTMC calculations with the experimental results published by Shafir et al. in Ref. [18].

Displacement Gate
Scanning the relative phase ϕ between the second harmonic and the fundamental allows for finding the ϕ max (n) for each harmonic order n in the spectrum, for which the recorded intensity generated at that harmonic order is maximised.Out of the complete cloud of trajectories started in the CTMC, every returning trajectory is assigned their harmonic order n with Then, the summed up probability weights of all returned trajectories contributing to harmonic order n from the CTMC simulation with two-colour phase ϕ is shown in Figure 6A.These displacement gate plots are normalised to the maximum intensity within each harmonic order individually.The normalisation is to show the maximum intensity within each harmonic order more clearly in the plot.
The horizontal line cut of each harmonic order can be fitted with a function, where C(n) describes the contrast of displacement gate and B(n) is the constant background intensity.As an example, Figure 7 shows the line out for harmonic order 31 of the dataset displayed in Figure 6A, together with the fitted function.
(A) (B)  The retrieved ϕ max for the two different simulation cases are shown on top of the experimental displacement gate plot in Figure 8. Neglecting the Coulomb potential (cyan) does not have any significant effect on the HHG spectra for most of the harmonic spectrum.Comparable ϕ max (n) are found both in the CTMC with Coulomb and the CTMC without Coulomb cases, and only around the cutoff energy do the two curves diverge.In addition, the ϕ max (n) curves extracted from the different CTMC simulations show remarkable agreement with the experimental spectrum for harmonics n ≥ 25.A constant shift between the experimental and CTMC values can be attributed to the unknown absolute phase difference in the experiment [18].

Long Trajectories Contribution
In Section 5.1, it was found that the maximizing phase ϕ max in the CTMC agrees reasonably with the experimental results, except for low harmonics slightly above the HHG threshold n < 25.For these low harmonic orders, the ϕ max extracted from CTMC shows an opposite slope from the experimental displacement gate plot.This is the range of energies where most long trajectories are initiated since the field is the strongest.
In order to estimate the efficiency of the HHG build up, the phase matching efficiency Equation ( 11) is calculated for all numerically computed returning trajectories.The final phase matching corrected probability weight Γ PM of a returning trajectory is given by where Γ(t i , p i ) is the ionisation probability at time t i with initial transverse momentum p i .For the following calculations, the jet centre is assumed to be located at z = 5 mm, with a jet diameter of j = 0.95 mm [43].Figure 6 contrasts two displacement gate plots.In Figure 6A, only short trajectories are considered, with their probability weights given by Γ(t i , p i ).In Figure 6B, the HHG spectra were calculated by summing up the corrected weights Γ PM from Equation (15), including also long trajectories.The effect for higher harmonic orders is relatively small.However, including long trajectories does yield a change in the slope of the ϕ max curve, similar to that found in the experiment (see Figure 9).Figure 6B is much noisier than the plot focusing on only short trajectories Figure 6A.While there are more trajectories included in Figure 6B, the increased number is not enough to completely compensate for the additional statistical fluctuation due to the introduced complex phase for all trajectories.
The shape and the value of PM(n, τ) strongly depends on the detailed parameters of the beam geometry, such as the spot size on the focusing mirror, the exact geometry of the gas jet in the interaction region with the laser beam, the minimal beam waist, the confocal parameter of the focused beam, and the jet position relative to the focus.Performing the CTMC analysis with just slightly different beam geometry leads to differences in the extracted ϕ max larger than their respective 95% confidence intervals for the fits in the extraction process (shown as the band covering the data points in displacement gate plots).Considering this, it is difficult to predict accurately the ratio of short and long trajectories corresponding to the experimental conditions.However, including long trajectories based on phase matching efficiency always shifted the ϕ max (n) curve to larger phases, with the strongest influence on the lower harmonics.Taking into account long trajectories therefore leads to qualitatively better agreement with the experimental displacement gate, compared to the version where all long trajectories are excluded a priori.Additionally, the oscillation contrast C(n) of the displacement gate plot especially for the lower harmonic orders is reduced by adding long trajectories.Compare again the two plots in Figure 6. Figure 10 highlights that the displacement gate contrast in the CTMC with only short trajectories is much stronger (blue) than the contrast in the experiment of [18] (red).In the case when long trajectories are considered as well (green), the contrast drops to comparable values.

Conclusions
We have shown that CTMC simulations, whether they take the ion Coulomb force into account or not, essentially reproduce the three-step-model of HHG [10].The extracted ϕ max (n) is virtually independent of the ion Coulomb force.Consequently, the reconstructed initial and return times in Ref. [18] are expected to be robust, and the approximations used for the calculations are valid.If exclusively short trajectories are used for the analysis, then the experimental displacement gate can be reproduced for harmonic orders n ≥ 25.Long trajectories, however, were needed for qualitative agreement between the CTMC simulations and the experiment for low harmonic orders.Their contribution can be estimated by looking at the phase matching efficiency.The displacement gate contrast strongly depends on whether long trajectories are considered or not.Agreement is only achieved again if phase matching is considered for all trajectories.

Figure 1 .
Figure 1.Photoelectron trajectories contributing to High Harmonic Generation (HHG).The colour code represents the kinetic energy of the trajectory at the return, where red is low and dark blue is high energy.Long trajectories are greyed out because typically the short trajectories dominate the spectrum.The black solid line represents the driving laser field −F x (t).

Figure 2 .
Figure 2. Short trajectories in two-colour HHG displacement gate.Some trajectories are displaced and miss the parent ion on the return, and others can come back to where they originated.The colour code corresponds to the kinetic energy carried by the electrons at the return, blue = high energy, red = low energy.

Figure 3 .Figure 4 .
Figure 3. Kinetic energy of single classical trajectory calculations with Coulomb.(A) recorded at the exit radius for each trajectory; (B) recorded at a radius of 3 au from the parent ion.A fixed radius yields spurious energy for trajectories which are ionised late in the cycle, and low harmonics in the spectrum can not be recovered.Dark blue stars (light blue circles) represent short (long) trajectories, respectively.Red dots show the experimental data from Shafir et al. black the values from quantum stationary phase calculation[32] and grey the three-step-model[10]. Values for the last three provided from[18].The values extracted from the Classical Trajectory Monte Carlo (CTMC) simulation in (A) reproduce the three-step-model so closely that the grey line is barely visible.The CTMC is of course based on this analytical model, except that it takes the Coulomb potential into account.

Figure 5 .
Figure 5. (A) phase mismatch for an example long and short trajectory, depending on the position of the gas jet after the laser beam focus; (B) phase matching efficiency for two example trajectories depending on the position of the gas jet z.

Figure 6 .
Figure 6.CTMC with Coulomb simulation, the two displacement gate plots show the influence of long trajectories.(A) only short trajectories with extracted ϕ max (n) shown in white; (B) both long and short trajectories are considered.

Figure 7 .
Figure 7. Intensity oscillation of harmonic order 31 while scanning the two-colour phase ϕ.Blue dots are extracted from Figure 6A, the red line is the fit of Equation (14).

Figure 9 .
Figure9.The extracted ϕ max (n) curves for the case of only short trajectories (white), and all trajectories with phase matching (green) are compared to the experimental displacement gate plot (colour bar)[18].

Figure 10 .
Figure 10.The displacement gate contrasts C(n) extracted from fitting Equation (14) to the experimental data (red [18]), the CTMC TIPIS model (Tunnel ionisation in Parabolic coordinates with Induced dipole and Stark Shift) with only short trajectories (blue) and the CTMC TIPIS model considering all trajectories (green).The error bars show the 95% confidence interval of the Equation (14) fit.

Table 1 .
Quantities in phase matching.