Disentangling Long Trajectory Contributions in Two-Colour High Harmonic Generation

Abstract

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.

1. 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

$${\omega }_n=E_{\mathrm{kin}}+I_{\mathrm{p}},$$

(1)

with ${\omega }_n=n·\omega$, $I_{\mathrm{p}}$ the ionisation potential and $E_{\mathrm{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

$$\begin{array}{cc}\hfill {\mathit{F}}_x\left(t\right)& =f\left(t\right)F_{0}cos\left(\omega t\right)\widehat{x},\hfill \end{array}$$

(2)

with envelope $f\left(t\right)$, maximum field strength $F_0$, optical frequency $\omega$, 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\left(t_i\right)$ 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

$$E_{\mathrm{cutoff}}=n_{\max }\omega =3.17U_{\mathrm{p}}+I_{\mathrm{p}},$$

(3)

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].

2. 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\omega }/I_{\omega }=\xi$ acts along the orthogonal direction y

$${\mathit{F}}_y(t,\phi )=f\left(t\right)F_0\sqrt{\xi }cos(2\omega t+\phi )\widehat{y},$$

(4)

slightly displacing the trajectories, depending on the two-colour delay $\phi$.

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 $\phi (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

$${\phi }_{\max }\left(n\right):=\phi (t_i,v_{i,y}=0).$$

(5)

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 ${\phi }_{\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].

3. 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 $\phi$ 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\left(t\right)={cos}^2(\Omega t)$ envelope for the field, $\Omega$ was chosen to yield a 6 fs FWHM in intensity, and the carrier wavelength was $\lambda =800$ nm.

For the subsequent propagation, classical Newton’s equations were solved in two versions, using the MATLAB (R2016a, MathWorks^® 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 Figure 3B and Figure 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.

4. 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,

$${\mathit{k}}_n=n·{\mathit{k}}_1+\mathit{K},$$

(6)

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. Quantities in phase matching.

z	longitudinal position in direction of beam propagation
r	radial position away from the beam axis
$t_i$	exit time
$t_r$	return time
$\tau =t_r-t_i$	excursion time
${\mathit{k}}_n=n\frac{\omega }{c}{\widehat{k}}_1$	radiated harmonic
${\mathit{k}}_1=\frac{\omega }{c}\widehat{z}+\nabla arg\left(\frac{1}{b+2iz}exp\left\{-\frac{\frac{\omega }{c}{r}^2}{b+2iz}\right\}\right)$	fundamental Gaussian beam
$S_{\mathit{p}}(t_r,t_i)$	classical action for electron wave packet propagation
$\mathit{K}(z,r)=\nabla \left(n\omega t_r-S_{\mathit{p}}(t_r,t_i)\right)$	intrinsic atomic phase
N	total number of interacting targets along the axis
j	gas jet diameter

. 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.

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

$$\delta k=k_n-|n·{\mathit{k}}_1+\mathit{K}(z,r=0)|.$$

(7)

For an interaction region on axis ($r=0$) in forward propagation along the axis, this mismatch can be estimated as

$$\delta k(z,n,\tau )=n\frac{\omega }{c}-\left|n\left(\frac{\omega }{c}-\frac{2b}{b^2+4z^2}\right)+\frac{F_0^2{z}_0^{2}z}{2{\omega }^2{(z^2+z_0^2)}^2}·\tau \right|,$$

(8)

where we used approximations described in Ref. [27].

Figure 5A shows the phase mismatch $\delta k$ of the wave vectors for two sample excursion times ${\tau }_{\mathrm{short}}=25$ au and ${\tau }_{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 $\left|\delta k\right|$).

Starting from the phase mismatch $\delta k(z,n,\tau )$ Equation (8), one can estimate the total amplitude efficiency with the help of a geometric series

$$GS(z,n,\tau ,N)=\sum _{t=0}^{N-1}exp\left(i\delta k(z,n,\tau )\frac{j}{N}·t\right)\frac{1}{N}.$$

(9)

This geometric series oscillates around a value for low N, and converges for larger $N\gtrsim 40$. The limit of $N\to \infty$ is given by

$$GS(z,n,\tau )=\frac{i(1-exp(ij\delta k(z,n,\tau )\left)\right)}{j\delta k(z,n,\tau )}.$$

(10)

A typical XUV pulse contains several orders of magnitude more photons than $N=40$ per harmonic order, which justifies taking the limit of $N\to \infty$. Then, the absolute square

$$PM(z,n,\tau )=GS(z,n,\tau )G{S}^{\ast }(z,n,\tau )=\frac{2}{j^2\delta k^2(z,n,\tau )}-2\frac{cos\left(j\delta k\right(z,n,\tau \left)\right)}{j^2\delta k^2(z,n,\tau )}$$

(11)

yields the efficiency for the intensity of the high-order harmonic radiation. For perfect phase matching

$$\underset{\delta k\to 0}{lim}PM(z,n,\tau )=1,$$

(12)

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\mathrm{mm}$ [43].

5. CTMC Simulations

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

5.1. Displacement Gate

Scanning the relative phase $\phi$ between the second harmonic and the fundamental allows for finding the ${\phi }_{\max }\left(n\right)$ 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

$$n=\frac{1}{\hslash \omega }\left(E_{\mathrm{kin}}\left(t_r\right)+I_{\mathrm{p}}\right).$$

(13)

Then, the summed up probability weights of all returned trajectories contributing to harmonic order n from the CTMC simulation with two-colour phase $\phi$ 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

$$I(n,\phi )=C\left(n\right)cos\left[2(\phi +{\phi }_{\max }\left(n\right))\right]+B\left(n\right)$$

(14)

function, where $C\left(n\right)$ describes the contrast of displacement gate and $B\left(n\right)$ 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.

The retrieved ${\phi }_{\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 ${\phi }_{\max }\left(n\right)$ 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 ${\phi }_{\max }\left(n\right)$ curves extracted from the different CTMC simulations show remarkable agreement with the experimental spectrum for harmonics $n\ge 25$. A constant shift between the experimental and CTMC values can be attributed to the unknown absolute phase difference in the experiment [18].

5.2. Long Trajectories Contribution

In Section 5.1, it was found that the maximizing phase ${\phi }_{\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 ${\phi }_{\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 ${\Gamma }_{\mathrm{PM}}$ of a returning trajectory is given by

$${\Gamma }_{\mathrm{PM}}=\Gamma (t_i,p_i)·PM(n,\tau ),$$

(15)

where $\Gamma (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\mathrm{mm}$, with a jet diameter of $j=0.95\mathrm{mm}$ [43].

Figure 6 contrasts two displacement gate plots. In Figure 6A, only short trajectories are considered, with their probability weights given by $\Gamma (t_i,p_i)$. In Figure 6B, the HHG spectra were calculated by summing up the corrected weights ${\Gamma }_{\mathrm{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 ${\phi }_{\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,\tau )$ 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 ${\phi }_{\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 ${\phi }_{\max }\left(n\right)$ 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\left(n\right)$ 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.

6. 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 ${\phi }_{\max }\left(n\right)$ 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\ge 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.