Numerical investigation of the sequential-double-ionization dynamics of helium in different few-cycle-laser-field shapes

Philipp Wustelt,1,2,* Max Möller,1,2 Markus S. Schöffler,3 Xinhua Xie ( ),3 Vaclav Hanus,3 A. Max Sayler,1,2 Andrius Baltuska,3 Gerhard G. Paulus,1,2 and Markus Kitzler3,† 1Institute of Optics and Quantum Electronics, Friedrich-Schiller-University Jena, D-07743 Jena, Germany 2Helmholtz Institute Jena, D-07743 Jena, Germany 3Photonics Institute, Technische Universität Wien, Vienna, Austria (Received 1 November 2016; published 14 February 2017)


I. INTRODUCTION
Studying the quantum dynamics of multielectron systems is a fundamental issue.One of the most elementary processes is the removal of both electrons from a helium atom.Double ionization of helium by strong laser fields serves as a benchmark process for the understanding of light-matter interactions in complex systems and has been investigated extensively both theoretically and experimentally over the past few decades [1].
There are two primary processes leading to laser induced double ionization.The first is nonsequential double ionization, which can be observed for close-to-linearly polarized laser fields.In this process, a field-driven recollision of a tunnel ionizing electron is responsible for the knockout of the second electron [2][3][4].This process can be efficiently controlled by manipulating the shape of a few-cycle laser field, using, e.g., the duration or carrier-envelope phase (CEP) of the pulse [5].The second process, sequential double ionization (SDI), becomes dominant as the ellipticity, , of the laser field increases.In this case, the second electron is emitted via an additional tunnel ionization step.As the removal of the first electron results in a more tightly bound second electron, the second ionization requires a higher laser intensity.For elliptically polarized laser pulses the ion momentum distribution (IMD) resulting from SDI exhibits a characteristic four-peak structure [6,7].This is well understood within the framework of the semiclassical model of strong-field ionization [2,8] under the assumption of two independent, sequential ionization events.
This assumption is questionable, however, if the SDI process is squeezed into an extremely short time interval, which can be done by using very short pulses.The resulting SDI dynamics has been probed in a recent benchmark measurement of the CEP-dependent IMD from helium produced by SDI in intense close-to-circularly polarized few-cycle pulses [9].It was demonstrated that the four-peak structure of the IMD, typically observed for longer pulses, may turn into a six-peak structure for few-cycle fields.This behavior has been reproduced by a purely sequential-ionization model, where the relative timings of the two-electron emission events are linked only by the interdependence of the field-governed ionization rate [10].However, certain differences between measured and simulated results remained [9], and only chirp-free pulses were considered.
Here we investigate how different laser parameters, such as different orders of chirp, and the resulting pulse shape influence the dynamics of SDI and the measured IMD.We focus on determining the parameter ranges for which the aforementioned four-peak structure turns into a six-peak structure.Specifically, we investigate the influence of the pulse parameters: CEP (φ 0 ), pulse duration (τ 0 ), peak intensity (I ), and second-(φ (2) ) and third-order (φ (3) ) spectral phase on the IMDs.
Chirped pulses are a well-established tool in molecular and atomic quantum control with femtosecond pulses [11][12][13][14][15].In contrast, the influence of chirp is rarely considered for extremely fast processes such as strong-field tunnel ionization on the attosecond time scale.However, a good understanding of the influence of chirp and the other aforementioned parameters is particularly important for evaluating differences between experiments and theoretical models on SDI in fewcycle fields [9,16].One reason for the lack of investigations of the influence of these parameters is the difficulty in precisely measuring and controlling them.
Here, our approach is to separately observe the effects of each parameter in the resulting IMD and underlying SDI dynamics.The paper is organized as follows.First we discuss our model and how it differs from existing theories (Sec.II, with details of the simulations in the Appendix).Next, in Sec.III, we demonstrate the validity of our model by reviewing and extending the quantitative comparison with recently published data on SDI of helium [9].Finally, in Sec.IV, the model is used to explore the effects of the aforementioned pulse parameters.

II. SEMICLASSICAL MODEL OF SEQUENTIAL DOUBLE IONIZATION
The SDI model used here and in Ref. [10] is detailed in the Appendix.In short, we assume a laser pulse with Gaussian spectrum defined by the parameters: peak intensity, I ; transform-limited full-width-at-half-maximum (FWHM) pulse duration of the intensity envelope, τ 0 ; center wavelength, λ 0 ; carrier-envelope phase, φ 0 , and second-order, φ (2) , and third-order, φ (3) , spectral phase.This defines the field, E(t), and the corresponding vector potential, A(t).We use the ionization rate proposed by Tong and Lin [17], which is based on the static field ionization rate [18] and adjusts this rate with a prefactor to correct overestimation in the high-intensity overbarrier ionization regime.This rate is used to solve numerically the coupled rate equations for the time-dependent population of the two charge states [19].
A classical Monte Carlo simulation is used to generate an ensemble of combinations of ionization times, t 1 and t 2 , which are subsequently used to calculate the asymptotic electron momenta, p 1 = −A(t 1 ) and p 2 = −A(t 2 ).The resulting ion momentum of the doubly charged ion, p ion , is given by the sum, p ion = −(p 1 + p 2 ).Three-dimensional focal volume intensity averaging is used to account for the extended laser-target overlap in the measurement [20].
Our model assumes strictly sequential double-ionization dynamics.Forces between the electrons as well as between the electrons and the ion are neglected.However, performing the Monte Carlo simulation according to the numerical solution of the coupled rate equations, which includes the exponential dependence of the tunneling rate on the field strength, accounts for the fact that the probability for ionizing the second electron is sensitive to the time-dependent population of the first charge state.This interdependence is strongly affected by the subcycle shape of the field and couples the relative timing of the two ionization events.This becomes more important the shorter the pulse is, since the time intervals during which the first and second ionization events take place start to overlap due to the rapidly increasing intensities of such pulses.
For SDI in an elliptically polarized field, the ab initio approach would be direct integration of the two-electron time-dependent Schrödinger equation.However, this approach is computationally extremely demanding and thus has not been realized to our knowledge.Moreover, this approach yields only limited insight into the underlying dynamics.
For these reasons, semiclassical models have been the preferred approach (see the recent review Ref. [1]).In our model we neglect several aspects to simplify the theory while covering the relevant effects.Specifically, we do not include the time-dependent modulation of the ionization rates induced by the oscillation of the remaining bound wave function caused by the removal of the first electron [19,21,22].Additionally, we neglect forces between the two electrons and ions as well as between electrons [23][24][25].However, the model reproduces the main features observed experimentally and yields a great deal of insight into the underlying dynamics.Note that recently a fully classical calculation which considers similar laser parameters [26] yielded comparable results for the ion momentum distributions corresponding to SDI of argon with elliptically polarized laser fields ( = 0.75).

III. COMPARISON BETWEEN SIMULATION AND MEASUREMENT
Good agreement between the measured [9] and the simulated SDI distributions (see Fig. 1) is obtained for a pulse Here and throughout the manuscript the measured IMDs were rotated such that the major polarization axis is parallel to the x axis as it is assumed in the simulations.(e), (g) Same as panel (c) but for CEP = 0 (e) and CEP = π/2 (g), where results with a CEP ±5 • around the nominal value were summed.(d), (f), (h) Corresponding momentum distributions simulated by the CTMC method with τ 0 = 4.38 fs and an intensity of I = 9.6 × 10 15 W/cm 2 at the position of the gas target.CEP dependence of the measured (i) and simulated (j) momentum distributions along p z (with all other directions integrated over).Note that the arbitrary offset of the CEP is adapted in the simulation to fit best to the count rate maxima with small momenta in the experiment.duration of τ 0 = 4.38 fs (φ (2) = 0, φ (3) = 0, λ 0 = 750 nm, = 0.95) in the simulation.This pulse duration is somewhat shorter than the one measured in the experiment, but is still within the uncertainty of the measured values (τ 0 = 4.5 fs, φ (2) ≈ 0, φ (3) ≈ 0, λ 0 = 750 nm, = 0.95).
The peak intensity in the simulation required to reproduce the measured data was identified by matching the cutoffs of the IMDs in experiment and simulation based on the relation between the electron momenta and the laser field vector potential, p 1,2 = −A(t 1,2 ).Varying the peak intensity in the laser focus, best agreement was found for I 0 = 7.4 × 10 16 W/cm 2 .Note that in the experiment, the laser focal position was outside the gas jet (for details see the Appendix).Therefore the intensity on target was considerably lower than in the focus.Using the notation of Ref. [9], the intensity on target within the limited spatial extension of the gas target that selects a small part of the focal volume in the measurement is on average equivalent to an intensity of I = 9.6 × 10 15 W/cm 2 .This result is in very good agreement with the experimental value of I = 1 × 10 16 W/cm 2 , obtained by an independent calibration with the method described in Ref. [27].
Particularly striking is the agreement for the CEP-integrated IMD in the polarization plane of the field [see Figs.1(a)-1(d)].The main features observed in the measurement, including the splitting of the four-peak structure into a six-peak structure, are reproduced by the simulations.Even the CEP-resolved results, that are compared in Figs.1(e)-1(h), agree quite well.However, subtle differences are noticeable as indicated by the arrows in Figs.1(g)-1(h).These deviations transfer into a corresponding feature in the CEP-resolved momentum distribution along the z axis, p z [compare Fig. 1(i) with Fig. 1(g) and note the arrows].In the following, we will investigate the dependence of these deviations on several parameters and thus obtain insight into their origin.

IV. INFLUENCE OF THE PULSE PARAMETERS ON MOMENTUM DISTRIBUTIONS
Now that we have established the validity of the model by reproducing measured IMDs both integrated over all CEPs and for certain values of the CEP, we investigate the influence of the pulse parameters (I , τ 0 , φ (2) , φ (3) ) on the simulated momentum distributions.Specifically, to obtain further insight into the SDI dynamics in few-cycle laser pulses, we examine how the pulse parameters influence the main features of the IMD and for which parameter values the best match between the measured and simulated data is achieved.
CEP-resolved SDI is particularly sensitive to the pulse parameters as the first ionization typically occurs at significantly, often more than an order of magnitude, lower field strength.Thus, the precision with which one needs to know the pulse shape for an accurate modeling is increased compared to that required if only single ionization is modeled, in particular, when the initial ionization step occurs early on the rising edge of the pulse, the shape of which is very sensitive to the laser parameters.In contrast, the second ionization step predominantly happens within the smooth and flat region around the pulse peak.In addition, the temporal shift of spectral components introduced by chirp of different orders can strongly affect the ionization dynamics.
These considerations are particularly important in view of the experimental facts that (i) nonlinear spectral broadening and subsequent recompression may easily result in a complicated spectral phase with pre-and postpulses [28] and (ii) in difficult to precisely characterize complicated pulse structures.The typical method to generate few-cycle pulses starts with spectral broadening of the incoming pulses by the ionization-induced plasma nonlinearity as well as Kerr nonlinearity in a gasfilled hollow-core fiber [29,30].After spectral broadening, the linear and nonlinear chirps introduced by the process are compensated for by dielectric multilayer mirrors in order to produce a few-cycle field.
Note, however, that this compensation is typically not perfect and optimized for minimizing second-and third-order chirp, which leaves some level of nonlinear chirp.Moreover, the resulting spectra are typically far from Gaussian and exhibit strong spectral modulations in intensity [31].This often results in pronounced pre-and/or postpulses, which can easily reach amplitudes of ∼ 10% of the main pulse.Therefore, knowledge of the spectral phase can be expected to become important for processes like SDI, which is sensitive to such low-intensity pulse structures due to the early emission of the first electron.

A. Dependence on pulse duration
Before we discuss the dependence of SDI on fine variations of the pulse duration in the few-cycle regime, we compare SDI by few-and multicycle pulses.For multicycle pulses, the ionization probability of each ionization step is distributed over many laser cycles and widely spreads in time [see Fig. 2(c)].Thus, the number of electrons emitted in the positive or negative direction is almost identical and the momentum distributions for both single and double ionization are symmetric even for fixed CEP [see Fig. 2(e)].This results in a symmetric double ionization distribution [see Fig. 2(a)].Further, as shown in Fig. 2(c), the long and slowly varying pulse envelope removes any correlation between the phase of the laser field at the times of the first and second sequentialionization steps, i.e., the classical times at which the first and second electrons are born in the continuum are unrelated.
In contrast, for few-cycle pulses the ionization probabilities of both ionization steps are squeezed into a single cycle or even into a single half cycle of the field [see Fig. 2(d)].As the circularly polarized field rapidly changes in amplitude and direction, i.e., during the outwards spiraling of the electric-field vector, the observed photoelectron momentum distribution takes the shape of a spiral and becomes asymmetric [see Fig. 2(f)].The shorter the pulse, the stronger the interdependence of the ionization times, t 1 and t 2 , and of the corresponding momenta, p 1 and p 2 , which leads to an asymmetric IMD [see Fig. 2(b)].For the limit of a single contributing half cycle and a fixed CEP, the first electron emission would be confined to a very small angular sector and the range of instants of the second ionization event and the resulting emission directions are extremely sensitive to the shape of the instantaneous electric field.
Simulating the pulse duration dependence from the fewcycle to the multicycle regime in Fig. 3, demonstrates how the main features of the IMD depend on the pulse duration, τ 0 , still assuming φ (2) = 0, φ (3) = 0.For example, as shown in Fig. 3(a), a cut through the IMD around p x = 0 starts at a four-peak-structure for very short pulses, changes into a six-peak structure around the pulse duration used in the experiment (τ 0 ≈ 4.5 fs), and changes back to a four-peak structure for longer pulse.As shown in Fig. 3(b), this behavior manifests itself in a pulse-length-dependent IMD width, which we characterize using the variance, μ 2 of p z .To further illustrate the pulse-length dependence, examples of the CEPintegrated two-dimensional IMDs for selected pulse durations are shown in Fig. 3(c).
Extracting both electron emission times from the simulation [see Fig. 3(e)] demonstrates that for durations around 4.38 fs the ionization probability for the second electron extends over two laser cycles, causing the six-peak structure in the IMD.For longer pulses, the second electron emission is distributed over more cycles and the six-peak structure turns back to a four-peak structure.
Analyzing the individual electron momenta, p z1 and p z2 , yields the τ 0 dependence of the CEP-integrated momentum coincidence map of the first and second emitted electron [see Fig. 3(d)].In these maps, parallel emission of both electrons into the same hemisphere [outer ring in the IMD in Fig. 3(c)] is located in the first and third quadrants, while the antiparallel electron emission case (inner ring) is located in the second and fourth quadrants.While only a very narrow region of antiparallel momentum combinations of p z1 and p z2 (second and fourth quadrants) is possible for 3.5 fs, a splitting for the antiparallel emission starts to appear around 4 fs and vanishes for τ 0 longer than 5 fs, which coincides with the observations of the six-peak structure.The structures in the first and third quadrants corresponding to parallel two electron   emission are independent of the pulse length [see Fig. 3(d) in the coincidence maps], which leads to a shape of the outer ring in the IMD that is independent of the pulse duration [see Fig. 3(c)].This result might be a starting point for explaining measurements of electron coincidence data for the double ionization of argon by 7-fs pulses [16] without the need to invoke any electron-electron correlation effects.An additional observation in Fig. 3 is that t 1 , and thus p 1 , is remarkably insensitive to the pulse duration as this ionization step is fully saturated, while t 2 is more susceptible to changes of τ 0 and the CEP enabling subcycle control of SDI [see Fig. 3(f)].

B. Dependence on peak intensity
As shown in the previous section, the existence of the sixpeak structure requires that the second electron emission is distributed over more than one optical cycle, while the first electron emission is saturated within essentially a single half cycle.If the intensity is fixed, this condition can be met by adjusting the pulse duration.Alternatively, with a fixed pulse duration, a similar behavior can be achieved by altering the intensity, as is shown in Figs.4(a) and 4(b) in comparison with Fig. 3(a).These figures display a cut through the IMDs around p x = 0.In the measured IMDs, the transition from four to six and back to four peaks for increasing I is similar to that seen for increasing τ 0 .Like the pulse duration dependence, these observations for the intensity dependence can be understood by acknowledging that increasing intensity enables more half cycles to contribute to the ionization of the second electron.In contrast to the measured results, in the simulated IMDs the six-peak structure also survives for very high intensities (see Fig. 4).The reason for the discrepancy between measured and simulated IMDs for high intensities remains unclear, when only chirp-free pulses are considered.
Another observation is that the intensity-dependent variance, μ 2 , increases monotonically with I , because the second  (2) ) scale.(e) Vector potential of a laser pulse with negative and positive φ (2) .(f) Selected ion momentum distributions for equal positive and negative φ (2) .emitted electron experiences a larger vector potential, since the second ionization step takes place at higher intensities.

C. Dependence on second-order spectral phase
A more subtle knob to control subcycle emission in SDI is the second-order chirp, φ (2) , of the few-cycle fields.Figure 5 compares the φ (2) dependence of the measured and simulated IMDs at fixed peak intensity along p z and the corresponding variance, μ 2 of p z .In the measurement φ (2) was varied by adding glass to the beam path and adjusting the laser-pulse energy to keep the peak intensity constant [9].The added glass leads to increased pulse durations, but leaves the spectrum unchanged.Because of the increasing φ (2) , the laser frequency becomes time dependent.
The measured IMDs feature a transition between four and six peaks, which is well reproduced by the simulation [see Fig. 5(d)].The transition can be explained using similar arguments as those used for explaining the pulse duration dependence.That is, the first ionization step is saturated and the second step is distributed over an increasing number of optical cycles as φ (2) is increased and the pulse is stretched.
In the simulation, we also considered negative φ (2) .The resulting IMDs and μ 2 evolution reveal an interesting effect.For negative φ (2) [see Fig. 5(c)], the variance, μ 2 , is considerably smaller than for positive φ (2) of the same absolute value.This goes along with qualitatively different features in the chirp dependence of the cut through the IMD as well as the IMD itself [see Fig. 5(d)].
The effect can be explained by recalling that the chirp induces a time-dependent frequency sweep, which translates into time-dependent electron momenta, because of the mapping of ionization time into the electron momenta via p = −A.If the leading edge of the laser pulse is blueshifted (φ (2) < 0), the first emitted electron experiences a smaller vector potential and thus acquires smaller electron momenta.In contrast, the second electron emitted around the center of the pulse, where the carrier frequency remains unchanged as compared to unchirped pulses, does not experience such a shift in momentum.Conversely, for a redshifted (φ (2) > 0) leading edge, the first electron obtains larger momentum due to the larger vector potential.This effect of the chirp φ (2) on A is visualized in Fig. 5(e), which shows A for a chirped pulse with Gaussian, and thus symmetric, envelope of the electric field, E(t).Due to the relation, A(t) = − t −∞ E(t )dt , the envelope for the vector potential is asymmetric.As the ion momentum, p ion , is the sum of both electron momenta, it follows that if p 2 stays more or less the same, while p 1 is controlled by the frequencies present at the leading edge, the distribution of the resulting p ion is wider or narrower for pulses with positive or negative φ (2) values, respectively, in agreement with Fig. 5(c).

D. Influence of third-order spectral phase
Variations of φ (2) lead to a time-dependent frequency within the pulse, but the pulse envelope remains symmetric.This is no longer the case when introducing third-order dispersion (TOD), φ (3) .A positive TOD produces an envelope that exhibits a sequence of postpulses with decreasing intensity in addition to an asymmetric main pulse with a steep falling edge.Changing the sign of the TOD reverses the temporal order of this envelope.Despite having only a fraction of the amplitude of the main pulse, as previously discussed, these satellite pulses influence the SDI dynamics by changing the ionization times significantly.
Figure 6(d) shows simulated IMDs for a scan of the TOD with all the other parameters being kept constant (I = 9.6 × 10 15 W/cm 2 , τ 0 = 4.38 fs, φ (2) = 0).Varying TOD in the range of ±180 fs 3 leads to FWHM durations of the main pulse up to 8 fs.For increasing magnitude of TOD, the six-peak structure quickly vanishes and only a four-peak structure is observed.This shows that the measured six-peak structure cannot be due to prepulses or postpulses.
Similar as for the variations of φ (2) , the IMDs for different signs of the spectral phase, φ (3) , are different.This is not surprising as φ (3) causes an asymmetric field envelope and thus also an asymmetric vector potential.Analyzing the φ (3)  dependence of the IMDs' widths in Fig. 6(c) shows that the maximum of the width is shifted towards negative values.This might be explained by a steeper slope of the leading edge of the main pulse's envelope that occurs for a small amount of FIG. 6. Simulated third-order spectral phase (φ (3) ) dependence of He 2+ momentum distributions for a laser pulse with τ 0 = 4.38 fs (φ (2) = 0): (a) Same as Fig. 1(d) but for a Gaussian pulse with φ (3) = 70 fs 3 , (b) CEP dependence of momentum distributions along p z for φ (3) = 70 fs 3 (with all other directions integrated over), (c) the variance (second-order moment) in p z as a function of φ (3) , (d) momentum integrated over |p x | < 2 a.u. as a function of φ (3) [compare to Fig. 5(d)], (e) same as panel (a) for φ (3) = −70 fs 3 , (f) same as panel (b) for φ (3) = −70 fs 3 , (g) envelope of the intensity for φ (3) = −70 fs 3 , (h) ionization times of the first and second emitted electron (integrated over all values of the CEP), and (i) momentum distribution of the radial momentum, p r = p 2 z + p 2 x , of the first and second electron (integrated over all values of the CEP).
However, for stronger negative TOD, starting around values of −50 fs 3 , the prepulses play a significant role.As Fig. 6(h) shows for φ (3) = −70 fs 3 , the first ionization step is separated between the main pulse and the prepulse, which results in a bimodal momentum distribution of the first emitted electron [see Fig. 6(i)].Interestingly, this splitting of the electron momentum distribution is not translated into an additional splitting of the IMD, since only the typical four-peak structure is observed in the corresponding CEP-integrated IMD [see Fig. 6(e)].This can be explained by the fact that p ion is the convolution of the distributions of p 1 and p 2 .As a consequence of the convolution in the final distribution of p ion , the splitting of p 1 is smeared out.
Another interesting effect shown in Figs.6(e)-6(i) is that the increased relevance of prepulses for larger negative TOD coincides with the appearance of additional peaks in the CEP-resolved IMD.These peaks are very similar to the CEP-dependent features observed in the measurement, but not reproduced by the simulations using the Fourier-transformlimited pulses that we described in Sec.III (see Fig. 1).Comparing Figs.1(i) and 6(f) reveals that for φ (3) = −70 fs 3 a similar feature spreads over a wide range of CEPs in the IMD at p z ≈ 5 a.u.(note the arrows).However, for the CEP-integrated IMD obtained with negative φ (3) , the six peaks merge nearly completely into a four-peak structure.Thus, for the range of parameters investigated here, it is not possible to obtain full agreement between simulation and measurement for both the CEP-resolved and the CEP-integrated IMDs by varying single parameters.This indicates that the features observed in the experiment that are marked by arrows in Fig. 1 are likely to be caused by very small deviations of the phase from a perfectly flat spectral dependence, which highlights the delicate dependence of the IMDs on very fine changes of the pulse parameters and their interplay.Obtaining insight into elusive effects caused by electron-electron interaction [9,16] from measured IMDs thus necessitates superb characterization of the pulse parameters.

V. SUMMARY
In conclusion, we have compared measured ion momentum distributions after double ionization of helium by close-tocircularly polarized few-cycle fields to classical trajectory Monte Carlo simulations, which assumed a purely sequentialionization model.The simulations emphasize that both ionization steps in SDI are linked by the interdependence of the field-governed ionization rate, i.e., that the probability for emission of the second electron is sensitive to the first ionization step.
Using this simple model we have demonstrated remarkably good agreement with the observations in the measurements [9] over a very large range of laser-pulse parameters: CEP, pulse duration τ 0 , peak intensity I , and second-order φ (2) and third-order φ (3) spectral phase.The parameter dependence of the unusual six-peak structure observed in the measured ion momentum distribution is reproduced in a very specific and narrow parameter range and it shows that it is very sensitive to the envelope and specific waveform of the ionizing laser field.
Extracting detailed information about the individual electron momenta using simulated coincidence maps yielded detailed insight into the ionization dynamics.Scanning the pulse duration we showed that the six-peak structure can be observed when the first ionization step is quickly saturated and thus localized in time while the second ionization step is distributed over more than one optical cycle.
Investigating the influence of the laser field's second-and third-order spectral phase demonstrated a strong impact on the ion momentum distribution for both types of chirp.This finding sets high requirements on the precision of the pulse characterization when using the described measurements for testing different models of SDI.Further, our results highlight opportunities for controlling two-electron emission on a subcycle time scale using circularly polarized few-cycle fields.that has the Fourier-transform-limited FWHM pulse duration in intensity τ 0 .The prefactor, F 0 , is connected to peak intensity F 0 = √ I in atomic units at the maximum of a Gaussian laser spot in three dimensions (see Ref. [33]). To where a subscript indicates that parameters for respective transitions are used.The factor nlm coincides with the ADK (Ammosov-Delone-Krainov) rate [18], while TBI is an empirical correction factor, which is regularly used to account for the fact that the ADK rate overestimates the ionization rate at the high field strengths used in the measurement [17,34].The quantity κ is related to the ionization potential, κ = √ 2|I P |, and C nl measures the amplitude of the electron wave function in the tunnelling region.Here, we follow the notation used in Ref. [35].Table I summarizes the parameters used in the calculations.
To account for the number of electrons that contribute to the ionization signal, the ionization rate needs to be weighted by the actual population [19].In helium, this means for the first ionization step He⇒He + = 2 TBIHe (t), (A12) where TBIHe (t) is the ionization rate with the parameters given in Table I.The factor of 2 is dropped for the second ionization step as there is only one electron left.
From the solution of the rate equations in the laser field one obtains the time-dependent populations, P He (t), P He + (t), and P He 2+ (t), for a given ellipticity , peak field strength F 0 , pulse duration τ 0 , and carrier-envelope phase φ 0 .
For the ionization times, t He⇒He + 1 and t He + ⇒He 2+ 2 , one wants to create a nonuniform random distribution of times that both conforms to the rate equations [see Eq. (A7)] and maintains causality, i.e., ensures that the first electron is ionized before the second with t He⇒He + 1 < t He + ⇒He 2+ 2 .To do this, we first solve the coupled rate equations to determine the time-dependent ionization probability, dP He + /dt, for the first ionization step.Second, the nonuniform random ensemble of ionization times, t He⇒He  , is chosen to adhere to this ionization probability.Third, for each and every initial ionization time in the ensemble, t He⇒He + 1 , the rate equations are solved again, but with the starting conditions P He (t He⇒He + 1 ) = 0, P He + (t He⇒He + 1 ) = 1, and P He 2+ (t He⇒He + 1 ) = 0.This yields P He 2+ (t) from which dP He 2+ /dt is calculated.Finally, the second ionization time, t He + ⇒He 2+ 2 , is randomly chosen for each and every initial ionization time following the aforementioned ionization rate.This procedure ensures that the ensemble of ionization time pairs adheres to the rate equations and that t He⇒He + 1 < t He + ⇒He 2+ 2 .The final momentum, p ion , of the He 2+ ion after double ionization is calculated as the sum of the momenta of the emitted electrons p ion = −(p 1 + p 1 ).Neglecting the initial velocity after tunneling, any interaction between the two electrons, as well as the Coulomb influence of the parent ion, the electron momenta are given by the vector potential at the instant of ionization, p e = −A(t), where A(t) = − t −∞ E(t )dt .Therefore, the final momentum of the He 2+ ion is given by the sum of the vector potentials at the two ionization times, i.e., p ion = A(t He⇒He   ).Finally, the simulations take into account the three-dimensional distribution of the laser focal geometry in the experiment [9].For this the calculation is repeated for different maximum field strengths.The contribution at each field strength to the ensemble of ion momenta is weighted with the corresponding ionization probability and the relative abundance of this intensity value in the three-dimensional Gaussian focus in accordance with the experimental conditions [20].
The simulated results were obtained with parameters chosen to be the same to the experimental conditions.The laser ellipticity was = 0.95.The center wavelength was 750 nm.For the CEP-independent calculations the CEP is randomly chosen for each ionization event.For comparison of the simulations with the experiment the focal geometry was assumed to be Gaussian with a Rayleigh length of 200 μm.In the CEP-resolved experiment the gas target thickness was 20 μm at a distance of 150 μm from the laser focal point.This was taken into account in the simulations.Note that due to the shift of the target away from the focal position the effective intensity in the target area is strongly reduced compared to the focal peak intensity.The CEP averaging due to the Gouy phase shift across the Rayleigh length of the laser beam was neglected in the calculation.This is justified because the gas jet was much shorter than the Rayleigh length [9].

FIG. 1 .
FIG. 1. Measured and simulated He 2+ momentum distributions (integrated over all values of the CEP) along the p z (a) and p x (b) directions and (c) measured momentum distribution in the polarization plane, integrated over all values of the CEP, for a pulse duration τ 0 = 4.5 fs and an intensity of I = 1 × 10 16 W/cm 2 in the position of the gas target.Here and throughout the manuscript the measured IMDs were rotated such that the major polarization axis is parallel to the x axis as it is assumed in the simulations.(e), (g) Same as panel (c) but for CEP = 0 (e) and CEP = π/2 (g), where results with a CEP ±5 • around the nominal value were summed.(d), (f), (h) Corresponding momentum distributions simulated by the CTMC method with τ 0 = 4.38 fs and an intensity of I = 9.6 × 10 15 W/cm 2 at the position of the gas target.CEP dependence of the measured (i) and simulated (j) momentum distributions along p z (with all other directions integrated over).Note that the arbitrary offset of the CEP is adapted in the simulation to fit best to the count rate maxima with small momenta in the experiment.

FIG. 2 .
FIG. 2. Comparison between long and short pulses: Simulated results for long (30 fs) (left) and short pulse (3.5 fs) (right) with CEP = [−5 • , 5 • ].(a) and (b) Ion momentum distribution.(c), (d) E z , E x , and envelope of laser field and ionization times for both electrons.(e), (f) Momentum distributions of the first and the second emitted electron.

FIG. 3 .
FIG.3.Simulated pulse duration dependence of the CEP integrated He 2+ momentum distributions (I = 9.6 × 10 15 W/cm 2 , φ(2) = 0, φ(3) = 0).(a) IMD along p z for |p x | < 2 a.u. as a function of pulse duration, (b) pulse duration dependence of the variance μ 2 of p z , (c) ion momentum distribution for different pulse durations, (d) coincidence map along p z of the momentum of the first vs the second emitted electron, and (e) ionization times for all (solid line) CEPs and for CEP of zero (dashed) and π/2 (dotted) of the first and the second emitted electron.

FIG. 4 .
FIG. 4. Intensity dependence of He 2+ momentum distributions for laser pulses with τ 0 = 4.38 fs (φ (2) = 0, φ (3) = 0) integrated over all values of the CEP and |p x | < 2 a.u.for measurement (a) and simulation (b).Note, for each intensity the distribution is normalized to the maximum.The intensity scale displays I at the position of the gas target as given in the measurement [9].
model the ionization dynamics, first the time-dependent population of the different charge states, i.e., P He (t), P He + (t), and P He 2+ (t), are calculated by solving the set of coupled rate equations