Pulse length effects in long wavelength driven non-sequential double ionization

The classical trajectory Monte Carlo model [43, 44] has been widely used to study double ionization (DI) of atoms in strong laser fields [15, 45, 46]. The laser field in our classical model is expressed as $\mathbf{E}(t)=E_{0}f(t)\sin(\omega t)\,\hat{\mathbf{x}}$, where the pulse envelope $f(t)=\sin^{2}(\frac{t\pi}{\tau})$, $\tau$ is the pulse duration, $E_{0}$ is the laser peak amplitude, and $\omega$ is the laser central frequency. The Keldysh parameter [39] in this work is smaller than 1, and thus we can use the ADK theory [47] to describe the tunneling step of the first laser-ionized electron. The tunneling rate can then be expressed as

$$W_{0}\left(t\right)=4\left(\frac{4I_{p}^{(1)}}{E\left(t\right)}\right)^{\frac{2}{\sqrt{2I_{p}^{(1)}}}-1}\exp\left[\frac{-2\left(2I_{p}^{(1)}\right)^{\frac{3}{2}}}{3E\left(t\right)}\right],$$

(1)

where the first ionization potential is fixed to $I_{p}^{(1)}=0.579$ a.u. The initial conditions are set as follows: the position of the first tunneling electron is approximated by $\mathbf{x}_{0}=-\frac{I_{p}^{(1)}}{|E(t)|}\hat{\mathbf{x}}$, its transverse momentum is approximated by a Gaussian distribution with width $\sqrt{\frac{|E(t)|}{2\sqrt{2I_{p}^{(1)}}}}$ and the longitudinal momentum is assumed to be zero, while the initial condition for the second bound electron is determined by a canonical ensemble [48], in which its energy equals the second ionization potential of Ar ($I_{p}^{(2)}=1.016$ a.u.). The sampling number at different tunneling times $t_{i}$ is assigned according to the weight $W_{0}(t_{i})$. Having set the initial ensemble, the classical equations of motion are then solved:

$$\frac{d\mathbf{r}_{i}}{dt}=\frac{\partial H}{\partial\mathbf{p}_{i}},\,\frac{d\mathbf{p}_{i}}{dt}=-\frac{\partial H}{\partial\mathbf{r}_{i}}.$$

(2)

Here, the Hamiltonian of the two-electron system can be written as

$$H=\sum_{i=1}^{2}\left[\frac{1}{2}p_{i}^{2}-\frac{2}{\sqrt{r_{i}^{2}+s}}\right]+\frac{1}{\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|}+\sum_{i=1}^{2}\mathbf{r}_{i}\cdot\mathbf{E}(t),$$

(3)

where the parameter $s$ is set to $0.1$ to avoid the Coulomb singularity at the origin. In our classical model, the bound electron in the initial ensemble is marked as $e_{1}$ and the tunneling electron is marked as $e_{2}$. After $e_{2}$ tunnels out at $t_{i}$, we define the instant at which $e_{2}$ returns back and the two electrons are closest as the recollision time $t_{re}$. Subsequently, we define the instant at which the energy of the bound electron becomes larger than 0 and the electron is not trapped afterwards at the ionization time $t_{ion}$ of $e_{1}$. In total, more than $10^{5}$ double ionization events are collected to guarantee the convergence of the results.

The quantum mechanical solution for two-electron systems interacting with a long-wavelength strong laser field is a formidable numerical problem, even when the electron spin is neglected. To overcome this difficulty, reduced-dimensional models, which have long proved their usefulness in the single active electron (SAE) case [9, 49, 50], have made their way toward multi-electron studies [51, 52, 53, 54, 55, 56, 57]. In the case of two-electron systems, the simplification involves a reduction from 3 spatial dimensions per electron (3+3) to 1 spatial dimension per electron (1+1).

Apart from the technical advantages of such models, the dimensional reduction comes with additional freedom to select the axes (geometry). The most straightforward choice of the reduced-dimensional axes in the case of linear polarization laser fields is the so-called Eberly or Rochester (EB) model [58, 59, 60]. In this approach, each electron is left on one axis exactly parallel to the polarization axis. This choice makes it easier to compare with experimentally available information. However, the downside is that the e-e Coulomb interaction is overestimated. This problem can be remedied in the context of reduced-dimensional models with a different choice of the simulation axes. One such choice is given by the classical Eckhardt-Sacha (ES) model [61], where each electron travels on a one-dimensional track inclined with respect to the polarization axis at a constant angle. The choice of the inclination angle in this model is derived from the existence of a saddle in the electrons potential energy during simultaneous escape, and as such is better for modeling the (e,2e) NSDI ionization pathway (see Fig. 1(a)).

The laser parameters in this work, i.e., moderate laser intensity, long wavelength, and long laser pulse duration, all result in a vast expansion of the wavefunction encoded in the position-space representation, which poses a difficulty for traditional grid-based methods. The common method to deal with such problems is the t-SURFF technique, which has recently been successfully extended to two-electron systems [62, 63], although at shorter wavelengths. Here, we employ a different method, which involves coherent transfer of the (multi-)electron wavefunction [64, 65, 66]. In this approach, the electron wavepacket, once sufficiently far from the core, is coherently transferred into separate, non- or semi-interacting regions in which the evolution is performed in momentum or mixed representations, preventing its absorption and hence the loss of information. In the following, we describe the details of the method.

The physical space is divided into four identically sized grids (regions) - the main region ($N$) hosting the initial wavefunction and the auxiliary, two single-ionized ($S^{1}$, $S^{2}$) and one double-ionized ($D$) regions (see Fig. 1(b)). The $S^{1}$, $S^{2}$, and $D$ regions store wavefunctions that evolve in parallel and similarly to the $N$ region using the split operator method, except that the evolution in the coordinate corresponding to the freed electron is performed solely in the momentum space without the Coulombic interaction (Volkov evolution). The free (ionized) electron direction is marked by the upper script in the $S$ region symbol. In the case of the $D$ region, the wavefunction evolves without any Coulombic interaction in the momentum space. The coherent wavepacket transfer process is cascading, i.e., occurs in two steps: In the first step, parts of the wavefunction from the $N$ region, which occupy the space past half of the grid size in the $i$ direction ($i={1,2}$), are subtracted from it and transferred to the $S^{i}$ region, stored in mixed representations: $\psi(r_{2},p_{1})$ for $S^{1}$ and $\psi(p_{1},r_{2})$ for $S^{2}$ (light green and red parts of the $N$ region in Fig. 1(b)). In the second step, parts of the wavefunction from the $S^{i}$ region extending past half of the grid size in the $i$ direction are transferred to the $D$ (doubly-ionized) region stored in the $\psi(p_{1},p_{2})$ representation (light blue parts of the $S^{1}$ and $S^{2}$ regions in Fig. 1(b)).

The description above can be summarized by Eqs. (4)-(6) in which the Hamiltonians $H_{N}$, $H_{S}^{i}$, $H_{D}$ correspond to the regions $N$, $S^{i}$, $D$, respectively:

	$\displaystyle\hat{H}_{D}$	$\displaystyle=$	$\displaystyle\sum_{i}\frac{1}{2}\hat{p}_{i}^{2}+A(t)\hat{p}_{i}$		(4)
	$\displaystyle\hat{H}_{S}^{i}$	$\displaystyle=$	$\displaystyle\hat{H}_{D}-\sum_{j\neq i}\frac{2}{\sqrt{\hat{r}_{j}^{2}+\epsilon}}$		(5)
	$\displaystyle\hat{H}_{N}$	$\displaystyle=$	$\displaystyle\hat{H}_{D}+H_{ee}-\sum_{i}\frac{2}{\sqrt{\hat{r}_{i}^{2}+\epsilon}}.$		(6)

The electron-electron interaction part $H_{ee}$ depends on the applied model described above. For the EB model, the interaction is described by $H_{ee}^{EB}$ (Eq. (7)), while for the ES model, it takes the form $H_{ee}^{ES}$ (Eq. (8)), i.e.

	$\displaystyle H_{ee}^{EB}$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{(\hat{r}_{1}-\hat{r}_{2})^{2}+\epsilon}}$		(7)
	$\displaystyle H_{ee}^{ES}$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{(\hat{r}_{1}-\hat{r}_{2})^{2}+\hat{r}_{1}\hat{r}_{2}+\epsilon}}.$		(8)

In addition to the efficient multi-region evolution, we also collect information about the quantum current fluxes streaming through boundaries located in the $N$ region (i.e. before the coherent transfer takes place; see Fig. 1(c)). Such boundaries are located at distances proportional to the classical quiver amplitude $E_{0}/\omega^{2}$ and inherit their naming from the previously defined grid regions, except for an addition of the $+/-$ sign signifying whether the flux was along the field axis (correlated) or perpendicular to it (anti-correlated) towards the $D$ region.

We record the 3D PMD taken in coincidence with the $\rm{Ar}^{2+}$ ions. Since the dead time particle detector is about 10 to 20 ns [13], we detect only one electron in the first quadrant [67, 68]. The other electron momentum can be reconstructed following the momentum sum conservation of the three involved particles coincidence, i.e.

$$P_{x}^{Ar^{2+}}=-(P_{x}^{1}+P_{x}^{2}),\;\;\;\;P_{y}^{Ar^{2+}}=-(P_{y}^{1}+P_{y}^{2}),\;\;\;\;P_{z}^{Ar^{2+}}=-(P_{z}^{1}+P_{z}^{2}).$$

(9)

By extracting the doubly charged ion longitudinal momentum $P_{x}^{Ar^{2+}}$ and one of the electrons $P_{x}^{1,2}$, we can deduce the other electron longitudinal momentum $P_{x}^{1,2}$. In order to ensure momentum conservation and realize a safe reconstruction of the missing particle, we filter out the number of coincidence events per laser pulse to be exactly one. In Figs. 2(a)-(c), we show the ($e_{1},e_{2}$) correlated PMD along the laser polarization axis and under different laser pulse durations. The PMD from the shortest pulse duration (see Fig. 2(a)) shows correlated momenta with similar amplitude along the diagonal $P_{x}^{1}=P_{x}^{2}$. When increasing the pulse length, the two outgoing electron momenta start displaying an asymmetry as shown in the PMD from Fig. 2(b). Finally, at longer pulse durations (see Fig. 2(c)), the ($e_{1},e_{2}$) correlated momentum map clearly exhibits a cross-like shape [25]. As a consequence, the experimental PMDs show a strong dependence on the pulse duration at a fixed laser intensity.

The mechanisms of NSDI under long-wavelength laser pulses have already been studied with CTMC [28, 30]. The rescattering electron returns with a large momentum because of the high ponderomotive energy, however, emitting the bound electron just requires a small part of such rescattering energy. Thus, there is a clear asymmetry between the two electrons. With the help of CTMC, we can clearly understand the underlying dynamics. However, previous works [28, 29] did not consider the influence of the pulse duration. In Figs. 3(a)-(c), we show that the PMDs have clear differences under different pulse durations and the results are consistent with the experimental results (see Fig. 2).

To understand the influence of the pulse duration on the NSDI mechanisms, we show the correlated PMD without symmetrization in Figs. 3(d)-(f). In these PMD spectra, $e_{2}$ is the tunneling electron. The tunneling electron tends to have smaller momenta, and this can be understood by a classical approximation. The final momentum of the tunneling electron can be approximated by $p_{re}-A(t_{re})$, where $p_{re}$ and $A(t_{re})$ are the momentum and the vector potential at the recollision time, respectively. Right after the recollision, the tunneling electron still carries a large remaining momentum, which has an opposite direction compared to the vector potential, and thereafter it is slowed down by the laser field, leading to a small momentum in the final state. On the contrary, the bound electron acquires a small part of the rescattering energy, and the total energy of it after recollision is around zero. Thus, the bound electron is ionized with a small $p_{0}$ and the final momentum can be approximated by $-A(t_{0})$, where $p_{0}$ and $A(t_{0})$ are the momentum and the vector potential at the ionization time of $e_{1}$, respectively. For a high rescattering energy, the direct (e,2e) channel dominates, thus, $t_{re}$ and $t_{0}$ are close to each other and are around at the maximum of the vector potential. Hence, the final momentum of $e_{1}$ ($-A(t_{0})$) is larger than $e_{2}$ ($p_{re}-A(t_{re})$). The differences in Figs. 3(a)-(c) can be summarized by pointing out two main factors: the momentum of $e_{1}$ becomes larger as the pulse duration increases, while the momentum distributions of $e_{2}$ (the tunneling electron) only have slight changes. We will analyze in detail why the pulse duration has such an influence on the PMD of the two electrons next.

The RESI mechanism [33, 69] is always important when considering NSDI. Also, multi-return dynamics, in which the electron returns without a recollision due to the transverse displacement and the recollision may happen after several return times, is relevant when considering NSDI under long-wavelength laser pulses [28]. Therefore, it is necessary to examine whether the two mechanisms are the reasons behind the experimental PMD spectra. For a longer-duration laser field, the tunnel-recollision process may contain more return times than a shorter pulse. To exclude the influence of these events, in which the number of returns is larger than the maximum value under the shorter pulse, in Figs. 4(a) and (b) we select the DI events in which the time difference between tunneling and recollision ($\Delta t_{1}$ is directly related to the multi-return trajectories [28]) is smaller than 3T. However, we still can see a clear difference in the PMDs under different pulse durations. Thus, those multi-return events, in which the number of returns is larger in the longer-duration laser field, do not account for the experimental results. In Figs. 4(c) and (d), we select DI events in which the time difference between recollision and emission of the bound electron is smaller than $0.1T$ (excluding the influence of RESI). Here, the correlated PMD spectra still show a clear difference. Therefore, there is no direct relation between the experimental results and the RESI mechanism.

In Fig. 5(a), we show the energy distribution of $e_{1}$ right after the rescattering of $e_{2}$ for two different pulse durations. Both distributions are located around 0 a.u., and therefore we conclude that the initial momentum of $e_{1}$ is small at the ionization time $t_{ion}$, and the final momentum of $e_{1}$ can be approximated by $-A(t_{ion})$ in both cases. In Fig. 5(b), when the laser FWHM changes from 4T to 14T, the distribution of $|A(t_{ion})|$ changes from 2 a.u. to 2.5 a.u. (the most probable value). This change can help us to directly understand why the momenta of $e_{1}$ become larger for longer pulse durations in Figs. 3(d)-(f). Then we just need to consider why $A(t_{ion})$ tends to become larger as the pulse duration increases. There are two possible factors: the phase of the vector potential at $t_{ion}$, or the envelope function $f(t_{ion})=\sin^{2}(\frac{t_{ion}\pi}{\tau})$ (e.g. the instantaneous laser intensity). In Fig. 5(c), we show the phase distribution of $\omega\,t_{ion}$. If the phase distribution accounts for the dynamics behind Fig. 5(b), the phase should be closer to 0 or $\pi$ to obtain a larger instantaneous vector potential when we increase the pulse duration. However, such a tendency is not observed in Fig. 5(c), thus, we can not claim the phase distribution as the reason behind the PMD of Fig. 5(b). In Fig. 5(d), the distributions of the value of $f(t_{ion})$ show a large difference under different laser pulse durations. When the laser FWHM is 14T, the distribution of $f(t_{ion})$ is located around 1.0. However, when the pulse duration is shorter, the distribution of $f(t_{ion})$ lies in the interval [0.7,1.0]. Thus, for the case of a shorter pulse, the instantaneous vector potential tends to be smaller than that of a longer pulse because of the envelope-induced intensity difference. Therefore, we can understand the PMD of $e_{1}$ under different pulse durations through Fig. 5(d).

When it comes to the tunneling electron $e_{2}$, the momentum distributions only have slight changes under different pulse durations. In Fig. 6(a), we can see that the tunneling electron tends to have a larger momentum right after the recollision for longer pulse durations. In Fig. 6(b), the distributions of the vector potential at the time of recollision $t_{re}$ resemble those in Fig. 5(b). Thus, the final momentum of $e_{2}$ can be approximated by $p_{re}-A(t_{re})$. Because both $p_{re}$ and $|A(t_{re})|$ tend to become larger as the pulse duration increases, the ultimate momentum of $e_{2}$ does not necessarily increase. To gain a better quantitative understanding, we calculate the expected energy of $e_{2}$ right after the recollision and show the results in Table 1. Here, $p_{re}$ is approximated by $\sqrt{2E}$ ($E$ is the energy of the tunneling electron $e_{2}$ right after the recollision). Estimates are consistent with the simulation results of Figs. 3(d)-(f). All results suggest that the instant laser field at rescattering, which is different due to the different values of the envelope function, distinctly affects the PMD spectra. Thus, the mechanism behind the evolution of the experimental results for different pulse durations can be attributed to an envelope-induced intensity effect.

The calculated PMDs of Argon exposed to 3100-nm pulses of variable length for the ES and EB models are shown in Fig. 7. Evolving the wavefunction simultaneously on 4 separate grids, using the Hamiltonians from Eqs. (4)-(6) and allowing the system to relax through $0$-field post-propagation, accumulates the PMD information in the $D$ region presented in the plots. As stated before, in the case of the EB model, the $p$ axes are parallel to the polarization axis, while they are oblique in the ES model. Convergence near the $p$ axes is slow, and a cross-shaped gap is visible on all plots. In some of the previous works at shorter wavelengths (e.g., [66]), an additional final coherent transfer (masked extraction and addition) steps from the $S$ and $N$ regions successfully completed the data without polluting the data with bound states. For the wavelength considered in this work, applying such a technique with a fixed mask appears to introduce bound states for smaller FWHM pulse lengths (4T and 6.6T), while improving the dataset only slightly for both models at 13.5T FWHM pulse length (results not shown). This feature is indicative of a substantial population of excited bound states at shorter pulse lengths. The general changes of the ES model observed in PMDs as the pulse length increases are consistent with the experimental data and the previously discussed classical results. For the EB model, the change in PMDs is less visible. The observed difference between the results obtained with the EB and ES models may be attributed to the fact that the ES model is better suited to modeling the (e,2e) NSDI ionization pathway.

The mechanism extracted in the classical case, confirming the PMD changes under different pulse durations due to different distributions of envelope values at the time of ionization (see Fig. 6(d)), also finds its counterpart in the quantum case. To show this effect, we shall analyze the probability fluxes (see Fig. 1(c) for definitions and Fig. 8 for results) that correspond to single ionization and various double ionization scenarios. The flux describing single ionization involves all possible events, that is, electrons that will return and recollide with the parent ion and electrons that will fly away to the detector. Regardless of the pulse length, single ionization peaks are located around the maximum of the pulse and throughout the pulse duration (and for some time after the pulse is gone). To infer the role of rescattering in DI, we calculated two fluxes populating the $D$ region, namely the one that directly crosses from the $N$ region and the other that comes from the $S$ region. Additionally, we are able to calculate fluxes to each of the PMD quadrants separately, i.e., to focus on both electrons escaping in the same direction (subscript $+$) or electrons escaping in opposite directions (subscript $-$). Looking at the envelope-aligned probability fluxes towards the $D$ region for different pulse lengths (see Figs. 8(c)-(h) and (c’)-(h’)) we notice a significant and consistent shift towards later times for shorter pulses. The shift of the maximum flux probability for the DI ($N\rightarrow D_{+}$) towards later times in case of shorter pulses confirms the envelope-induced intensity effect.

In essence, the envelope-induced intensity effect comes from the tunnel-recollision dynamics. In Fig. 9(a), we show the distribution of the time difference between tunneling and recollision. In both pulse durations, the time differences are most likely located around 1.2T, indicating that the multiple-return recollision process dominates in NSDI under 3100-nm laser fields [28] and the traveling time distributions are similar under different pulse durations. In Fig. 9(b), the tunneling time distributions under the scaled time are both located around the same value (center of the pulse). Considering the travel time after tunneling, the recollision most likely occurs 1.2T later. Thus, the recollision happens further back in the laser envelope with a smaller laser intensity for a shorter pulse in Fig. 9(c). Such an envelope-induced intensity effect under long-wavelength linear polarized lasers connects the PMD and the tunnel-recollision dynamics, a point that has not been discussed before.