Ultrafast laser-induced strong field phenomena play a pivotal role across a broad spectrum of fundamental areas within modern physics. For instance, the first step in the three-step model of high harmonic generation (HHG)^1,2, a cornerstone of attosecond science, relies on strong field-driven tunnel ionization. Similarly, powerful techniques like ‘attoclock’ employ this phenomenon to probe intricate details of quantum tunnelling^3,4. Strong field induced Coulomb explosion, fragmentation, tunneling and multiphoton ionization in the velocity map imaging (VMI)⁵ and cold target recoil ion momentum spectroscopy (COLTRIMS)⁶ experiments are used to study the variety of most fundamental physical and chemical processes in atoms and molecules^7,8. In the wake-field acceleration experiments^9,10, the strong laser field allows to accelerate particles with a record gradients, enabling the construction of high performance accelerators of much smaller size than conventional ones or even at the photonic chip scale^11,12,13.

Accurate measurement, optimization and control of laser intensity distribution and pulse duration in the interaction region are critically important for most of the strong field-driven experiments.

Upon initial examination, the problem may appear straightforward: to directly measure or estimate fundamental laser pulse characteristics, such as energy, width, and duration, and subsequently deduce the peak intensity. However, it has long been recognized and experimentally validated^2,14,15 that such a seemingly simple approach usually results in unacceptably large errors, frequently exceeding 50%, which can significantly impact the interpretation of experimental results. Several factors can contribute to such errors, including imperfect laser beam quality, aberrations in the focusing optics, distortion of ultrafast pulses due to spatial chirp (e.g., caused by misalignment of the compressor grating), pulse dispersion in the generation medium, and more.

Furthermore, high-power laser systems often are prone to temporal fluctuations in pulse intensity and variations in their spatial and temporal profiles, compounding the issue due to the highly nonlinear nature of strong-field phenomena. Therefore, a reliable in situ approach to measure and control these parameters in the focal volume becomes critically important for all areas involving strong-field laser physics.

To address this challenge, various methodologies have been reported. Typically, intricate, specialized, and costly techniques have been employed to achieve this objective, including mass spectrometry¹⁶, velocity map imaging (VMI)¹⁷, and cold target recoil-ion momentum spectroscopy (COLTRIMS)¹⁵. Furthermore, the utilization of these approaches necessitates very low gas densities, generally incompatible with HHG experiments, and often extremely low transverse temperatures, which further complicate the experimental setup, requiring larger system dimensions.

For high gas density measurements, particularly critical in strong-field physics, such as attosecond science, the estimation of laser peak intensity has often relied on the high harmonic generation (HHG) cutoff^18,19. However, in practice, determining the cutoff is nontrivial, as contrary to an idealized and clearly defined cutoff, a series of low-intensity harmonics frequently appear in the experimentally measured HHG spectrum^2,20. As a result, determining the cutoff harmonic order is somewhat subjective. Furthermore, it has been demonstrated that macroscopic effects, arising from averaging and phase matching, can distort the single-atom HHG cutoff response²¹, leading to inaccurate intensity estimates and a lack of access to pulse duration estimation within the interaction region.

In this study, we introduce a novel in situ measurement and control approach that solves this long-standing critical problem, opening the door to unprecedented accuracy in strong-field, in particular, attosecond science experiments. By conducting independent measurements in helium (He), the gas with a highest known ionization potential, I_p = 24.6 eV, and argon (Ar), I_p = 15.76 eV²², with exact the same beam parameters in the focal volume, we have identified indelible and meaningful signatures of double ionization, and barrier-suppression ionization (BSI).

Furthermore, we realized a dynamic theoretical model that incorporates single ionization, sequential double ionization (SDI), non-sequential double ionization (NSDI), non-sequential excitation (NSE), recollision-induced excitation with subsequent ionization (RESI), and BSI phenomena for strong field ionization. This model demonstrates excellent agreement and consistency with our experimental results. Importantly, while demonstrated here using argon as a test case, this approach seamlessly extends to the experiments with other atomic gases, including neon (Ne), krypton (Kr), and xenon (Xe).

Ion detection scheme

The schematic of our experimental setup is shown in Fig. 1. A linearly polarized laser pulse from a Ti:sapphire amplifier, centered around λ = 800 nm with a pulse duration of 30 fs (measured at the laser output), was focused by a spherical f = 500 mm focusing mirror M2 through the vacuum chamber window (W) onto a gas jet (GJ). To maintain precise control over the pulse energy, a motorized λ/2 waveplate (RWP) and a thin-film polarizer (TFP) were employed. The focal plane on the gas jet was monitored via an image-conjugated microscope, which comprised a beam-picking mirror (FM), an objective lens (OL), and a CCD imager.

The ion detector mesh ({rm{I}}{{rm{D}}}^{-}) is connected to the negative voltage of the high voltage power supply (({{rm{V}}}_{{rm{PS}}})) through the shunt resistor (({R}_{1}=1,MOmega)). The output is monitored by an oscilloscope (OSC) connected to the circuit via a DC decoupling capacitor (({C}_{1}=100,{pF})) and a resistor (({R}_{2}=10,MOmega)). The input capacitance of the scope (15 pf) and input impedance ((1,MOmega)) are not shown in the figure but are accounted for in the calculations

Full size image

The ion detector^19,23,24 is composed of a stainless steel grounded plate (ID^G) connected mechanically and galvanically to the body of the pulsed valve, and the metal mesh (ID^–) attached to the (ID^G) with dielectric rods and connected to the negative high-voltage, V_PS via a simple passive electronic circuit shown in Fig. 1. In this experiment, we employed a pulsed valve (Parker series 9 solenoid valve) that is widely recognized for its versatility and widespread use in numerous strong-field experiments. It is important to note that our approach is not limited to this specific valve and can seamlessly adapt to other gas jet sources as well. We applied a voltage of V_PS = −500 V and set the distance between the grounded plate and the mesh to d = 80 mm. Both the grounded plate and the mesh are circular, each with an 80 mm diameter.

The inclusion of the mesh, instead of solid component, is instrumental in minimizing the gas flow resistance, thereby enabling efficient evacuation of the residual gas by a vacuum turbomolecular pump located beneath the vacuum chamber. The pulsed gas valve was installed on a micron-precision (XYZ) vacuum positioner, ensuring that the laser beam intersected the jet at the center, approximately 100 µm below the injection flange. This intentional positioning was chosen for optimal conditions for attosecond science experiments (HHG), facilitating in situ measurements of HHG^2,19. This straightforward construction not only enhances the versatility of the ion detector but also simplifies its deployment across a wide range of experimental configurations. Consequently, our approach can be easily applied in a diverse array of strong-field experiments.

When the laser pulse interacts with the gas jet, the gas molecules and atoms undergo ionization within the pulse’s relatively short duration, typically on the order of tens of femtoseconds, through the strong-field ionization mechanisms. As a result, the gas medium’s resistance drops, enabling current to flow between the ground plate and the mesh, recorded as the voltage evolution on an oscilloscope as illustrated in Fig. 1. Utilizing a professional analog electronic circuit simulator LTspice, we performed the calculations to determine the transient response of the circuit, and fit it to the experimentally measured voltage trace, with a typical result shown in Fig. 2a.

a Example of measured (dotted blue curve) vs. calculated (solid red curve) ion signal temporal evolution for argon jet. The laser pulse energy was 0.6 mJ, and the gas backpressure was 7.5 Bar. The best fitted parameters were found ({R}_{{ID}}left({t}_{0}right)=13.27,MOmega) and ({tau }_{{ID}}=30.74,mu s) with goodness of fit ({R}^{2}=0.98). b Total charge (integrated current, modeled in LTspice and shown by dotted blue curve) vs. measured ion signal. The goodness of fit (solid red line) ({R}^{2}=0.999)

Full size image

The laser ionized gas was modeled as a time-dependent resistor R_ID(t) = R_ID(t₀) exp[(t − t₀)/τ_ID] in parallel with a capacitor, representing the ion detector geometry. Transmission lines with defined impedances and lengths were included to account for cable effects. The initial resistance R_ID(t₀), ionizing pulse arrival time t₀, and the decay constant τ_ID were fit to experimentally measured data.

The minor “bump” observed 1 millisecond after the primary pulse results from the ionization of residual gas induced by a successive laser pulse from our 1 kHz amplified laser system. We systematically varied the laser intensity by adjusting the Rotating Waveplate (RWP) in Fig. 1, and recorded the voltage trace for each intensity setting. For each measured voltage trace, we used the LTspice to compute the total charge, which corresponds to the integral of the current between the ground plate and the mesh, and confirm a linear dependence between the ion signal and the total charge as shown in Fig. 2b. The term “ion signal” can refer to either the maximum voltage measured on the oscilloscope trace or the integral of the voltage trace (as used in our analysis), both result in the linear dependence between the ion signal and the total charge.

The measured spatial intensity distribution of the beam in the focus is shown in Fig. 3. The input beam has been apertured down to 9.25 mm before the focusing mirror (M2) with an iris to improve the beam quality and to extend effective Raleigh range in the focus. By scanning the image-conjugated microscope (depicted as a dashed box in Fig. 1), we validated that the beam width and shape do not change within the measured 0.4 mm extension (FWHM) of the gas jet.

The inner black and the outer red line represent the (gamma =1) contours for helium and argon gases respectively, emblematically illustrating the regions of tunneling (inner part) and multiphoton (outer part) regimes. The intensity was calibrated using the energy per pulse of 0.4 mJ and ({tau }_{{FWHM}}=43) fs. The spatial intensity distribution in the focal plane was imaged using a microscope as illustrated in Fig. 1

Full size image

Strong Field Ionization Dynamics

Before delving into intricacies of strong-field ionization processes, it is essential to emphasize the clear qualitative distinctions between the perturbative and strong-field regimes of nonlinear responses of matter to a laser field.

In the perturbative regime, at relatively low intensities (typically ≤ 10¹³ W/cm²), the interaction between the laser’s electric field and the electrons within an atom or molecule is significantly weaker than the Coulomb interaction between the electrons and the nucleus. As a result, the material’s polarizability can be expanded as a converging power series in the laser electric field strength^25,26. Many hallmark nonlinear optical phenomena arise in this regime²⁷, such as second- and third-harmonic generation^25,28, parametric down-conversion²⁹, self-phase modulation^30,31, and supercontinuum generation³², among others.

In contrast, in the strong-field regime, where the optical field strength becomes comparable to the Coulomb interaction between the electron and nucleus (typically ≥ 10¹⁴ W/cm²), the polarizability power series expansion no longer converges. This regime is dominated by ionization processes, which play a pivotal role in our measurements. Single-electron ionization of neutral atoms or ions represents the most fundamental aspect of strong-field ionization and has been extensively studied for several decades^3,4,33,34,35. Keldysh’s seminal contribution³³ laid the foundation for our understanding of strong-field ionization in the so-called ’transparency region’, where the ionization potential, I_p, exceeds the energy of the incident laser field, (hslash {omega }_{L}).

The nature of the ionization process is characterized by the adiabaticity or Keldysh parameter, denoted as (gamma). This parameter is defined as (gamma =2pi frac{{text{T}}_{text{tunnel}}}{{text{T}}_{text{L}}}), where ({text{T}}_{text{L}}=frac{2pi }{{omega }_{text{L}}}) represents the laser field’s period, and ({T}_{rm{tunnel}}) is referred to as the “tunneling time”. The tunneling time³³ is defined as a ratio of twice the length of the potential barrier, ({L}_{rm{barrier}}) divided by a classical velocity, ({{rm{v}}}_{{rm{e}}}), of an electron moving under the barrier. The electron’s velocity is given by the virial theorem: ({{rm{v}}}_{{rm{e}}}=sqrt{2{I}_{p}/{m}_{e}}), where ({m}_{e}) is the electron’s mass. Assuming a point-like potential for electron-nuclear interaction and approximating the total potential with a triangular barrier solely created by the laser field, one gets ({L}_{rm{barrier}}={I}_{p}/e{E}_{L}).

Consequently, for the tunneling time, we get ({T}_{rm{tunnel}}=frac{2{{rm{L}}}_{{rm{barrier}}}}{{{rm{v}}}_{{rm{e}}}}=frac{sqrt{2{I}_{p}{m}_{e}}}{e{E}_{L}}), and for the Keldysh parameter, we obtain:

For a given laser frequency, a weaker laser field, ({E}_{L}), results in a longer potential barrier length, ({{{L}}}_{{rm{barrier}}}), leading to an increase in the Keldysh parameter, γ. In the limit of (gamma gg 1) or ({E}_{L}ll frac{sqrt{2{I}_{p}{m}_{e}}{omega }_{L}}{e}), the laser field undergoes multiple oscillations during the tunneling time, defining the multiphoton regime. In this regime, the ionization rate, denoted as ({W}_{i}), exhibits a power-law dependence on the electric field, a well-known feature in perturbative non-linear optics²⁷: ({W}_{i}propto {E}_{L}^{2K}), with (K=leftlfloor {I}_{p}/{omega }_{L}+1rightrfloor) representing the minimum number of absorbed photons required to conserve energy, known as multiquantumness. Conversely, in the limit of (gamma ll 1) or ({E}_{L}gg frac{sqrt{2{I}_{p}{m}_{e}}{omega }_{L}}{e}), the driving laser field experiences minimal change during the tunneling time, defining the well-known tunnelling regime with ionization rate for hydrogen atom given in Gaussian units³⁶ as ({W}_{i}=frac{4{m}_{e}^{3}{e}^{9}}{{E}_{L}{hslash }^{7}}exp left[-frac{2{E}_{a}}{3{E}_{L}}right]), where ({E}_{a}=frac{{m}_{e}^{2}{e}^{5}}{{hslash }^{4}}) is the electrical field experience by an electron situated at the Bohr radius ({a}_{0}=frac{{hslash }^{2}}{{m}_{e}{e}^{2}}) in the hydrogen atom.

It is essential to emphasize that in our experiment, as illustrated by the measured profile in Fig. 3, the intensity varies across the beam so that at the central region of the beam, ionization may predominantly occurs in the tunneling regime, while as we move away from the center, the field strength gradually decreases, transitioning towards the multiphoton ionization regime. This intermediate regime is particularly rich and interesting as it contains signatures of both tunnelling and multiphoton ionization mechanisms. Helium’s higher ionization potential restricts the dominant tunneling ionization to a smaller area compared to argon.

Shortly after Keldysh’s seminal work, Perelomov, Popov and Terentev in their work^34,37,38 addressed the effect of Coulomb interaction on ionization rate. They obtained Eq. (4) for the optical cycle-averaged ionization rate ({W}_{rm{PPT}}left({E}_{L},{omega }_{L}right)) for the linear polarized driving laser field (from now on and throughout the rest of the paper, unless otherwise explicitly stated, we will use atomic units, (hbar=e={m}_{e}=1)).

The refinement of strong field ionization models along with the experimental effort for its validation is an ongoing active field of research^{16,35,39,40,41,42,43,44,45,46,47,48}. It has been experimentally shown that ionization rate predicted by the PPT model is excellent for atoms up to the Keldysh parameter γ ~ 3–4^16,45 and even valid for small molecules⁴⁶. Another widely used ionization model is the Ammosov, Delone, Krainov (ADK) that can be obtained from PPT model in the limit of γ→0. Main advantage of the ADK is its relative simplicity. However, it has been shown experimentally that the validity of the ADK formula is limited to the tunneling regime with Keldish parameter γ ≤ 0.5⁴⁹. In this work, we employ the PPT model for calculating the cycle-averaged ionization rate ({W}_{01}) Eq. (4) including the BSI correction Eq. (18).

After a neutral atom is ionized by a strong laser, with a probability rate ({W}_{01}), resulting in the formation of a singly charged ion, there exists a probability that the remaining ion will undergo further ionization in the laser field, ultimately yielding a doubly charged ion. This phenomenon is referred to as sequential double ionization. The probability rate, denoted as ({W}_{12}), can be calculated using the same PPT ionization Eq. (4), but with the consideration of the ionization potential, ({I}_{p}), for the singly charged ion, which is higher than that of the neutral atom. Additionally, the residual charge is equal to Z = 2. Consequently, these factors lead to a reduced ionization rate when compared to single-electron ionization.

When an electron gets ejected from the parent ion, there is a probability that the electron wavepacket, driven by the oscillating laser field in the continuum, recollides with the parent ion. There are four possible outcomes of this process. Firstly, the recolliding electron may experience elastic scattering, or it may contribute to the HHG and generation of attosecond pulses¹. This is what makes our approach particularly useful in attosecond science, as it enables accurate in situ characterization of the spatially resolved intensity and duration of strong-field femtosecond pulses within the interaction volume. However, these processes do not contribute to ionization in the medium after the laser pulse is gone.

The recolliding electron wavepacket may also cause either the ionization of the parent ion, which is referenced as non-sequential double ionization (NSDI), with a rate of ({W}_{rm{NSDI}}left(Eleft(tright)right)) or ion’s excitation to a higher energy state, denoted as 1*, with a rate ({W}_{rm{NSE}}left(Eleft(tright)right)), and its subsequent ionization (RESI)^50,51,52, with a rate ({W}_{{1}^{* }2}left(Eleft(tright)right)). This ({W}_{{1}^{* }2}) ionization rate is calculated again using the PPT ionization model with the ionization potential of singly charged ion in the excited state.

To account for NSDI and the RESI processes in ionization dynamics, we implement a semiclassical model for the calculations of the ({W}_{rm{NSDI}}left(Eleft(tright)right)) and ({W}_{rm{NSE}}left(Eleft(tright)right)) rates⁵². In the nutshell, we calculate the probability of ionization (Pleft({t}_{i}right)) in the time window (left[{t}_{i},{t}_{i}+{dt}right]) within an optical cycle of period T, (Pleft({t}_{i}right)=Wleft(Eleft({t}_{i}right)right)bar{P}left(0le tle {t}_{i}right){dt}=Wleft(Eleft({t}_{i}right)right)exp left[-{int }_{0}^{{t}_{i}}Wleft(Eleft({t}^{{prime} }right)right)d{t}^{{prime} }right]{dt}), where (Wleft(Eleft({t}_{i}right)right)) is an instantaneous ionization rate Eq. (19) and (bar{P}left(0le tle {t}_{i}right)=exp left[-{int }_{0}^{{t}_{i}}Wleft(Eleft({t}^{{prime} }right)right)d{t}^{{prime} }right]) is the probability that no ionization occurs in the (left[0,{t}_{i}right]) time window. For each ionization time, we calculate the classical trajectory of electron within an oscillating field (Eleft(tright)=Aleft(tright)cos left(omega tright)), the return time t_r and return energy ({E}_{k}left({t}_{r}right)). A(t) is a slowly varying envelope of the pulse that is assumed to be constant within a time scale of one optical period (the so-called slowly varying envelope approximation (SVEA)⁵³). In semiclassical picture, upon ionization the electron wavepacket emerges from the tunnel barrier with a momentum coordinates, ({p}_{perp }) perpendicular to classical trajectory direction (in the atomic units)^54,55:

This wavepacket spreads as it propagates in the continuum during the ({t}_{r}-{t}_{i}) time window. Upon recollision, we calculate the overlap of the recolliding wavepacket with the respective energy ({E}_{k}left({t}_{r}right))-dependent cross-sections: ({sigma }_{rm{NSDI}}left({E}_{k}left({t}_{r}right)right)), the cross-section of electron impact ionization of ion, and ({sigma }_{rm{NSE}}left({E}_{k}left({t}_{r}right)right)), the cross-section of electron impact excitation of ion^56,57. We determine the probability of the process of interest by examining the overlap, as derived in the section “Non-sequential double ionization” of the Methods. The rates ({W}_{rm{NSDI}}) and ({W}_{rm{NSE}}) are found by normalizing the respective probabilities by half the period of the driving laser field.

When laser field approaches or exceeds the critical value of ({E}_{rm{BSI}}={I}^2_{p}/left(4Zright)), the ionization potential ({I}_{p}) reaches or surpasses the potential barrier created by the combination of Coulomb attraction to the atomic core with charge Z and the instantaneous laser field. In this so-called barrier-suppression ionization (BSI) region, it is well-known that strong-field ionization models^35,38 overestimate the ionization rate^49,58. To account for this discrepancy, we use a correction factor⁵⁹ for ionization rate calculations Eq. (18), as detailed in the section “Barrier-suppression ionization correction” of the Methods.

Taking into account all the aforementioned processes, we derive the following set of coupled equations that describe the ionization dynamics:

where ({N}_{0}left(vec{r},tright)), ({N}_{1}left(vec{r},tright)), ({N}_{{1}^{* }}left(vec{r},tright)), and ({N}_{2}left(vec{r},tright)) are the densities of neutral, singly ionized, singly ionized in the excited state and doubly ionized atoms/molecules, respectively.

To comprehensively elucidate the ionization dynamics within the focal volume in both space and time, we solve the system of Eq. (3) taking into account the measured intensity cross-section of the laser beam in the focus as shown in Fig. 3, as well as the change of gas jet density along the laser propagation direction within the jet. The gas density distribution, ({N}_{0}left(vec{r},t=0right)) within the jet has been measured using the laser induced strong-field ionization gas jet tomography²⁴ by scanning the laser beam across the slice of interest and measuring the ion signal, as depicted in Fig. 4. The spatial intensity distribution of the laser beam couples into the ionization dynamics equations via spatial dependence of respective ionization rates, (Wleft(Eleft(vec{r},tright)right)).

The argon gas jet was measured and tomographically reconstructed with a back pressure of 7.5 bar and a laser intensity of 0.5 mJ. The maximum ions charge density is ({rho }_{{ions}}^{{MAX}}=3times {10}^{15},{rm{cm}}^{-3}), and the maximum gas density is ({N}_{0}^{{MAX}}=8.6times {10}^{15},{rm{cm}}^{-3}). Gas jet is visualized as cut-outs within isosurfaces of constant gas density

Full size image

Figure 5 illustrates the typical ionization dynamics for helium and argon at the focus of a femtosecond laser pulse. For this analysis, we assume a typical for a femtosecond laser (f^{2}(t)={mathrm{sech}}^{2}(t/tau)) temporal intensity pulse profile⁶⁰, with a FWHM pulse duration of 43 fs and pulse energy of 2 mJ. The rationale for choosing these pulse parameters becomes evident from the discussion below. For pulse durations in this range, we can safely assume space-time invariance^61,62,63, i.e. (Eleft({vec{r}}_{perp },tright)={E}_{rm{peak}}varepsilon left({vec{r}}_{perp }right)fleft(tright)), where ({E}_{rm{peak}}) is the electric field peak amplitude, (varepsilon left({vec{r}}_{perp }right)=sqrt{{I}_{N}left({vec{r}}_{perp }right)}), with ({I}_{N}left({vec{r}}_{perp }right)) being the measured (as depicted in Fig. 3) normalized focused beam intensity distribution.

For helium (a–c) and argon (d–g). The temporal evolution of the spatial ion distributions is visualized in the Supplementary movie. This movie provides further insights into the dynamics observed in the figure

Full size image

Figure 5a depicts the time evolution of the laser intensity, the density of the neutral helium atoms, and the density of the singly ionized He atoms at the center of the driving laser beam. Due to the very high ionization potential of helium, double ionization is negligible for the considered laser pulse, and ionization of neutral atoms with the ({W}_{01}) rate is the dominant process. Consequently, for helium, the sum of neutrals and singly ionized atoms’ density remains constant. It is evident that the ionization becomes significant before the peak of the laser pulse. In this example, ~20% of the initial helium gas gets ionized, 10 fs before the pulse reaches its peak intensity (at this point in time, the intensity reaches 85% of its peak). After the laser pulse has passed, ~80% of the gas at the beam’s center has been ionized. This creates a distinct “hole-burning” effect in the spatial distribution of neutral atoms, as shown in Fig. 5b, complemented by the generation of spatially “molehill”-like localized singly ionized helium atoms, whose normalized density is visualized in Fig. 5c.

Figure 5d presents the time evolution of the ionization process for argon at the beam’s center, driven by the same laser pulse used in the helium example. Notably, while helium reached 70% ionization only past the peak intensity of the laser pulse on its falling slope, argon achieves this same level by the time the femtosecond pulse reaches just 40% of its peak intensity on the leading edge of the pulse. This significant difference clearly arises from argon’s lower ionization potential compared to helium.

These observations underscores the fact that strong-field ionization exhibits high sensitivity not only to the peak laser intensity but also to the temporal evolution of the driving pulse, particularly within the time window around the peak where the majority of the gas undergoes ionization.

Another major consequence of argon’s lower ionization potential compared to helium is that double ionization plays a significant roles in the ionization dynamics of argon. This is in striking contrast to helium, where the impact of these processes on ionization was negligibly small. Initially, as the laser pulse intensity increases, strong-field single ionization begins to dominate the ionization dynamics. This phenomenon leads to an increase in the density (({N}_{1}left(tright))) of singly ionized (A{r}^{+1}) ions, naturally accompanied by a decrease in the density of neutral atoms (({N}_{0}left(tright))), as illustrated in Fig. 5d. With the increasing density of (A{r}^{+1}) ions and laser intensity, double ionization gain prominence, leading to the double ionization of (A{r}^{+1}) ions (depicted as the dotted black curve labeled (A{r}^{+2}) in the figure).

This phenomenon depletes the transient population of singly ionized (A{r}^{+1}) ions at the beam’s center. Notably, the duration of the transient (A{r}^{+1}) population at the beam’s center is only 19 fs (FWHM), significantly shorter than the 43 fs laser pulse duration. While the high laser intensity at the center drives double ionization, it weakens with increasing distance further away from the center, leading to single ionization there as the dominant process.

This spatial dependence results in the characteristic “molerun”-like distribution of (A{r}^{+1}) ion density across the focal plane once the laser pulse has passed, as illustrated in Fig. 5f. This phenomenon is accompanied by a localized “molehill” peak of (A{r}^{+2}) ions (Fig. 5g) and a significantly more depleted “hole-burning” compared to helium (Fig. 5e).

Experimental results

We now turn to the experimental results of strong field ionization measurements as a function of laser pulse energy. Figure 6a and c show the measured ion signal for helium and argon, respectively, at various back pressures in the pulsed valve, ranging from 1.5 bar to 7.5 bar. As expected, increasing back pressure leads to a higher gas density in the focal volume of the gas jet, consequently enhancing the ion signal for a given gas. Additionally, at the same pulse energy, the argon signal surpasses that of helium due to the difference in their ionization potentials.

a In helium b In helium normalized c In argon d In argon normalized

Full size image

To gain deeper insights, we normalized the same measurements as shown in Fig. 6b and d. Surprisingly, our experiment reveals that the normalized ionization curves for both helium and argon exhibit no discernible dependence on the gas density within the jet. Notably, apart from a single outlier point at 4.5 bar for argon, all ionization curves remain strikingly close within the expected range of random measurement error. We attribute this lone outlier among the ten measurements in both helium and argon to an anomaly in the driving laser stability and exclude it from further consideration.

It’s important to note that after the laser pulse is over, the multitude of recombination processes may occur within the ionized plasma, such as three-body recombination, radiative recombination, and others⁶⁴. Accompanied by space charge effects, such as Debye shielding⁶⁵ and ambipolar diffusion⁶⁴, which disrupts the charge current and separation of ions and electrons within the laser-induced plasma subjected to the electric field of an ion detector, these recombination processes can decrease the fraction of ions and electrons reaching the detector plates and consequently diminish the measured ion signal.

One might anticipate a pronounced, generally nonlinear influence of gas density on the measured ion signal due to the interplay of these processes, potentially altering the ionization curve shapes. However, our experimental results clearly indicate that it is not the case – the normalized ionization curves exhibit remarkable resilience and remain independent of gas density in the jet within the studied range.

Considering that the ion signal is proportional to the total ionized charge as illustrated in Fig. 2b, we calculate and fit the theoretically calculated ion yield, ({Y}_{rm{ions}}), to the measured normalized ion signals shown in Fig. 6b, d for every back pressure for both helium and argon. The ion yield is determined by the integral of ion charge density, ({rho }_{rm{ions}}left(vec{r}right)) (as depicted in Fig. 4), over the focal volume: ({Y}_{!!rm{ions}}={iiint }_{!!{V}_{rm{focal}}}{rho }_{rm{ions}}left(vec{r}right){dV}). The ion charge density, after the ionizing laser pulse has passed, is given by ({rho }_{rm{ions}}left(vec{r}right)=left({N}_{1}left(vec{r}right)+{N}_{{1}^{* }}left(vec{r}right)right)+2{N}_{2}left(vec{r}right)), where ({N}_{1}left(vec{r}right)), ({N}_{{1}^{* }}left(vec{r}right)) and ({N}_{2}left(vec{r}right)) represent singly and doubly ionized densities calculated at the end of the pulse by solving the set of Eq. (3), as has been described.

Note that we measure both the total energy of the pulse as well as spatial intensity distribution in the focal volume, as shown in Fig. 3. For the measured charge and gas density distribution in the generated plasma, our analysis shows that the spatial distortion of the beam in the interaction region is negligible. Therefore, for a specific pulse shape, the only free parameter in the fitting procedure is the pulse duration, ({tau }_{rm{FWHM}}), or equivalently, the peak intensity ({I}_{rm{peak}}). More details about the fitting procedure are provided in the Methods.

Figure 7a shows a typical result of this fitting procedure in helium at a back-pressure of 7.5 bar. The optimal fit yields a full width at half maximum pulse duration of ({tau }_{rm{FWHM}}=41.7pm 0.5) fs. The upper horizontal axis in the figure displays the pulse peak intensity corresponding to this pulse duration.

a For helium, the extracted pulse duration is ({tau }_{rm{FWHM}}=41.7pm 0.5) fs with a goodness of fit ({R}^{2}=0.997). b For argon the extracted pulse duration is ({tau }_{rm{FWHM}}=43.3pm 1.6) fs with a goodness of fit ({R}^{2}=0.989). Both intensity and (gamma) were calculated using ({tau }_{rm{FWHM}})

Full size image

It’s noteworthy that increasing the ionizing pulse energy generates a plasma with higher density and larger volume within the gas jet by the time the driving laser pulse ends. This, in turn, can potentially amplify recombination and space-charge effects, as previously discussed, and subsequently diminish the measured ion signal at higher intensities. In the presented version, our theoretical model does not incorporate these macroscopic recombination or space-charge phenomena. Therefore, one might anticipate a systematic error that increases at higher intensities, potentially distorting the shape of the measured signal and causing a deviation from the theoretically calculated curve.

However, our results unequivocally demonstrate the opposite—the normalized measured ionization curves exhibit remarkable agreement with our theoretical model, characterized by a parametric fit error of ±0.5 fs (95% confidence interval) and goodness of fit of ({R}^{2}=0.997), even at high pulse energies and gas back-pressure. In principle, this pulse duration uncertainty encompasses both statistical errors and systematic errors. By repeating the measurements, we find a statistical errors to be within 1 fs range for single-pulse energy scan across 150 points spanning the energy range. To isolate systematic errors arising from theoretical model imperfections, we analyzed and presented in Fig. 7 data averaged over 100 ion signal measurements for each pulse energy, effectively reducing statistical error to 0.1 fs.

Only at low pulse energies ((< 0.8) mJ) and correspondingly very low ionization yield do we measure ionization yield in helium higher than predicted by the theory. This is attributed to the ionization of the residual gas impurities in helium, though not visible on the scale of the plot and having negligible effect on the fitting, as these parasitic contributions are very low compared to the dynamic range of the measured ion signal. Given the high ionization potential of helium, at low intensities, residual gases such as nitrogen, oxygen, and water vapor mainly contribute to the ionization signal, surpassing the minor contribution from ionization of helium itself.

Figure 7b depicts a representative result of the analogous fitting procedure in argon at the same back-pressure of 7.5 bar. The best fit results in a pulse duration of ({tau }_{text{FWHM}}=43.3pm 1.6) fs. In contrast to helium, where single-electron ionization dominates across the measured intensity range due to its high ionization potential, double ionization phenomena comes into play in argon. Our calculations show that SDI plays a major role, followed by NSDI, while RESI contribution was negligible in the relevant intensity range. Additionally, BSI has been found to play significant role as well. This complex interplay also explains why, while the fitting results for argon remain good with a pulse duration fit error of ±1.6 fs and ({text{R}}^{2}text{=}0.989), they fall slightly short of the exceptional goodness of fit achieved for helium. Remarkably, despite this intricate interplay of ionization processes in the model, we obtained a pulse duration in argon that closely matches that measured in helium.

It is crucial to include these strong field ionization phenomena in the computations; neglecting them leads to an inconsistent pulse duration in argon compared to the helium result, indicating a deficiency in the comprehensive ionization model. For instance, not accounting for sequential double ionization in argon results in a pulse duration of 74 fs, a stark >70% deviation from the helium (42 fs) or well-modeled argon (43 fs) result. Omitting BSI only in argon leads to 58 fs, still very significant 38% deviation.

At the same time, we validated that neglecting SDI (or any double ionization process) in helium ionization calculations made, as expected, no difference to the pulse duration result. Neglecting NSDI doubles the systematic error in argon compared to helium, but this effect is significantly less important than the SDI contribution. This is because the contribution of ({text{W}}_{12}) dominates over ({text{W}}_{text{NSDI}}) in argon already at the peak intensity of (5times {10}^{14},{rm{W}}/{rm{cm}}^{2}).

Figure 8 summarizes the pulse durations extracted for both helium and argon across various back pressures. Notably, the average pulse duration across all gases and pressures is ({tau }_{text{FWHM}}=43.1pm 1.6) fs.

Pulse duration ({tau }_{rm{FWHM}}) extracted by fitting the theoretical model to the measured ionization curves for helium (red crosses) and argon (blue squares). The error bars indicate a parametric fit error (95% confidence interval)

Full size image

As each data point in the plot is shown after averaging reducing statistical error to 0.1 fs, the error bars indicate the pulse duration parametric fit error. The average value for pulse duration for the measurements in He is (42.3pm 1.5) fs and in argon (44pm 1) fs. Thus we find that each of the systematic errors due to imperfections in theoretical model, the variations in the back pressure and usage of different gases does not exceed 2 fs or is less than 5% of the measured pulse duration. The main source of statistical error in our setup was caused by beam pointing stability issues—this can be further improved, for instance, by implementing active beam stabilization techniques⁶⁶. However, for most practical purposes, the demonstrated statistical error (~1 fs), even without averaging, is already satisfactory and comparable to the systematic errors.

The pulse duration at the output of our regenerative amplifier is 30 fs; however, in situ measurements using our scheme reveal a pulse length 40% longer. Since peak intensity is inversely proportional to pulse duration, relying solely on the output pulse duration would result in a significant error in estimating the peak intensity within the interaction region. This once again underscores the importance of in situ measurements for both pulse duration and peak intensity, especially in strong field and attosecond science experiments.

To assess the impact of temporal pulse shape on our results, we considered both the typical sech² and Gaussian profiles used to characterize femtosecond laser outputs. Additionally, we modeled the influence of strong dispersion-induced pulse distortions, by propagating 30 fs pulse through 50 mm fused silica. Notably, the estimated laser peak intensity remained remarkably robust, showing a maximum variation of only 3% across different pulse shapes and distortions.

It is important to stress that the measured ion signal depends not only on the peak intensity, but also on the pulse duration. As the strong field ionization is highly nonlinear process, the shape around the peak of the pulse becomes particularly important. Interestingly, for sech², Gaussian, and dispersion-distorted pulses with the same FWHM, the pulse intensity varies by no more than 2% of its peak within the upper half of the intensity profile, where most ion signal is generated. However, the main variation in intensity between these pulses occurs around their off-peak regions, where intensity is low and impact on the strong field ionization is minor. This further underscores an important advantage of our technique for strong-field ionization experiments: it is inherently sensitive to the pulse characteristic most crucial for the process, namely the intensity profile around the peak, where the majority of strong field ionization phenomena take place.

The robustness of our method to strong dispersion distortions is crucial for practical applications, as the pulse’s spectral phase can be significantly affected by various factors. These include, among others, dispersion from optics and propagation in the air, suboptimal compressor grating settings, and potential distortions from target gases. While the target gas distortions were not significant in the specific experiment being discussed, the method’s resilience to such factors remains important for ensuring accurate and reliable results in diverse real-world scenarios. By validating this robustness, the method can be more effectively applied across different systems and environments, enhancing its overall usefulness.

Utilizing the presented ionization model, our method is valid within the slowly varying envelope approximation (SVEA)⁵³, i.e. (left|{dA}left(tright)/{dt}right|ll {omega }_{L}left|Aleft(tright)right|), indicating that the variation of the electric field envelope within an optical cycle is negligible. However, it is important to stress that satisfying the SVEA condition is crucial only within a time window around the peak of the laser pulse, where the ionization rate, and therefore the electric field, are substantial, thereby relaxing the SVEA requirement. For example, in our experiments using a laser with a central wavelength of 800 nm, our analysis shows that these criteria are well satisfied for pulses with durations ≥ 10 fs. One interesting direction for future research involves refining the ionization model to expand its validity to the near-single-cycle limit of femtosecond laser pulses, with controlled carrier envelope phase (CEP)⁶⁷, and taking into account the possible impact of space-time coupling⁶³.

In conclusion, the accuracy of strong-field experiments, particularly in attosecond science, should no longer be constrained by uncertainties in the femtosecond laser beam characteristics within the interaction region. Spatially resolved, in situ intensity measurements of the driving laser, combined with the measurement of spectrally resolved HHG wavefronts²⁰, pave the way for resolving a long-standing issue in attosecond science experiments: the spatial averaged measurements across the focal volume⁶⁸. This approach holds the potential to significantly enhance experimental precision, for instance, in imaging orbitals^18,69, tracing molecular dynamics⁷⁰, identifying and time resolving tunneling wave packets ⁷¹, following Auger decay⁷², and probing polar molecules with high harmonic spectroscopy² – enabling access to a single-molecule response in a strong laser field⁶⁸.

Helium’s exceptionally high I_p allows the single-electron PPT model to adequately describe its ionization within the relevant intensity range. However, most other molecules and even noble gas atoms possess significantly lower I_p—similar or even lower than argon. To comprehensively elucidate the ionization dynamics of such systems, accounting for double ionization and BSI becomes crucial. Frequently, ion measurements have been interpreted solely through single-electron strong-field ionization models, considering only peak intensity ionization rates; our findings clearly illustrate the inadequacy of this approach.

Furthermore, the implications of our findings extend beyond improving the fidelity of strong-field and attosecond science experiments. The observed robustness to pressure variations and remarkable agreement between our measurements and the single-atom ionization model for helium suggest that ion detection approach can serve as a valuable experimental tool for validating and refining theoretical models of strong-field ionization for more complex atoms and molecules.

Strong-field ionization rate PPT

In this section, we explicitly summarize the formulas for the optical cycled averaged ionization rate, employing the Perelomov, Popov, Terentev (PPT) formulas³⁴, which are used in this work with certain corrections:

Where:

({E}_{L}) is the laser field,

({omega }_{L}) is the laser cycle frequency,

({n}^{* }) is the effective principal quantum number:

({l}^{* }) is the effective orbital quantum number:

({E}_{0}) is the characteristic atomic electric field:

where k is the number of above threshold ionization photons (ATI)²², starting with the minimum number of ATI photons, ({k}_{min }) (this formula for clarity is written in SI units):

Here, ({U}_{p}) is a ponderomotive energy, which is the cycle-averaged kinetic energy of the electron motion in the field. The ({k}_{min }) is rounded up to the nearest integer here, whereas in the original PPT paper³⁴, ({k}_{min }) was not rounded, which leads to unphysical results.

The simplified cycled averaged ADK formula for the ionization rate, ({W}_{rm{ADK}}left({E}_{L},{omega }_{L}right)), is often used³⁵ (can be obtained from ({W}_{rm{PPT}}left({E}_{L},{omega }_{L}right)) in the limit γ → 0 as shown in³⁴ Eq. (59)):

We explicitly emphasize the use of the full ({W}_{rm{PPT}}) cycle averaged ionization rate Eq. (4), omitting the ADK approximation, for single electron ionization in this work.

For argon gas with an orbital angular momentum quantum number (l = 1), we calculated the ionization rate by averaging over all possible values of its magnetic quantum number (m), representing the projection of the angular momentum along the direction of the applied field. This methodology follows the approach outlined by Larochelle et al. in their investigation of the Coulomb effect in multiphoton ionization of rare-gas atoms⁴⁵ and is consistent with the averaging procedure detailed in equation (61) of PPT model³⁴.

Barrier-suppression ionization correction

The BSI correction to the ionization rate, described in the Eq. (18)⁵⁹, has been implemented. The relevant parameter values a₁, a₂, and a₃ are provided in Table I. It has been demonstrated that this modified ionization rate formula aligns well with ({text{E}}_{text{L}}le 4.5{text{E}}_{text{BSI}}), falling within the range the peak laser field used in this work.

(Note, that we corrected a typo in the original paper⁵⁹, where a minus sign appears in the Eq. (18)).

The BSI formula we use is a further development of the correction proposed by Tong et al.⁷³. It has been noted that their original correction for the tunneling ionization fits well for ({text{E}}_{text{L}}le 2.5{text{E}}_{text{BSI}}). However, for larger values of ({E}_{L}), the ionization rate decreases for argon, which is inconsistent with the expected physical behavior.

Non-sequential double ionization

Non-sequential Double Ionization (NSDI)^1,74 is a process in which a singly charged ion undergoes secondary ionization due to recollision with the electron released earlier in the same laser cycle.

The initial step of non-sequential double ionization (NSDI) involves the ionization of an atom within the laser cycle. Consequently, the instantaneous ionization rate³⁴ is used:

where time dependency comes from the field, (E(t)=A(t)cos (omega t)), where (A(t)) is the amplitude of the specific cycle under the slowly varying envelope approximation.

As discussed in the Main section of the manuscript, we calculate the probability of ionization, (Pleft({t}_{i}right)), in the time window (left[{t}_{i}, {t}_{i}+{dt}right]) within an optical cycle of period T.

Electrons which return to their parent ion possess a probability of either ionization or excitation, determined by the electron’s kinetic energy and the respective energy dependent cross-section for ionization/excitation by collision (({sigma }_{{rm{NSDI}}/{rm{NSE}}}))^56,57. This probability is computed as an overlap between the electron’s wavepacket spread in the laser field and the corresponding cross-section. The calculations are outlined as follows:

Where, N is the normalizing factor,

∆r_⊥ is the electron wavepacket spread,

({R}_{c}) is the collision radius,

Finally, the NSDI rate and NSE rate are,

If the electron, following its collision with the parent ion, excites rather than ionizes the ion, the ion itself may undergo ionization from the laser field while in its excited state. This phenomenon is termed recollision-induced excitation with subsequent ionization (RESI)^50,51,52. From the ion’s excited state, its ionization rate is calculated by Eq. (18), with the proper values for the ion’s I_p and Z = 2.

Procedure for pulse duration extraction by fitting to the ionization curve

In our experiment, we measure the laser’s pulse energy (PE), the normalized spatial intensity distribution in the interaction region (({I}_{N}left({vec{r}}_{perp }right))) as illustrated in Fig. 3, and the ion signal as a function of PE (ionization curve).

The PE is related to the laser field amplitude (Eleft(vec{r}_perp,tright)) by:

where ({E}_{rm{peak}}) denotes the laser’s pulse peak field. The pulse peak intensity (({I}_{rm{peak}})) can be expressed as:

We then relate the peak intensity (({I}_{rm{peak}})) to the pulse energy (PE) using equation:

Here, (c) is the speed of light in vacuum, ({varepsilon }_{0}) is the vacuum permittivity. Note the ({int_{-infty }^{infty}};{f}^{2}left(tright){dt}={Cs}times {tau }_{rm{FWHM}}), where ({text{C}}_{text{S}}) is a constant depending on the normalized temporal pulse envelope (f(t)) (e.g., for f(t) = sech (2arcsech((sqrt{0.5}))t/({tau }_{rm{FWHM}})), ({text{C}}_{text{S}}=1/text{arccosh}left(sqrt{2}right))). Therefore, the pulse duration ({tau }_{rm{FWHM}}) is only free fitting parameter.

For each ({tau }_{rm{FWHM}}), we solve the set of ionization dynamics rate equations and calculate the total ion signal ({Y}_{rm{ions}}) (see details in main paper). We use the Matlab Curve Fitting Toolbox to fit the ionization curve to the calculated ({Y}_{rm{ions}}) and determine the optimal ({tau }_{rm{FWHM}}), the fitting error, and the goodness of fit ({R}^{2}). Subsequently, we use Eq. (27) to extract the ({I}_{rm{peak}}).