Attoclock with bicircular laser fields as a probe of velocity-dependent tunnel-exit positions

We consider ionization of helium under the influence of a strong laser pulse with an electric field E(t) by using a single-active-electron (SAE) description in dipole approximation. The SAE description approximates the relevant electron dynamics in attoclock-like experiments with helium very well [54], compare also the reviews [37, 38]. To this end, the TDSE,

$$\begin{equation}\mathrm{i}\frac{\partial }{\partial t}\psi (\mathbf{r},t)=\left(\frac{1}{2}{\left[\mathbf{p}+\mathbf{A}(t)\right]}^{2}+V(\mathbf{r})\right)\psi (\mathbf{r},t),\end{equation} \tag{ 1 }$$

is solved in velocity gauge for a given vector potential $\mathbf{A}(t)=-{\int }_{-\infty }^{t}\mathbf{E}({t}^{\prime })\mathrm{d}{t}^{\prime }$ and a binding potential V(r). In 3D, the effective potential V for the helium atom is chosen as by Tong and Lin [55], but with the singularity removed using a pseudopotential for the 1s state with a cutoff radius r_cl = 1.5 a.u. [56]. In order to reduce the numerical workload, we additionally study the dynamics in a 2D helium model. To this end, the potential is further softened by replacing $r\to \sqrt{{r}^{2}+0.34}$ such that the correct ionization potential I_p = 0.90 a.u. is obtained. The ground state (serving as initial state) is calculated as eigenstate of the short-time field-free time-evolution operator for a time step Δt = 0.02 a.u. as described in [57]. This method of finding the initial state is beneficial for obtaining reliable PMDs even in regions of low signal in momentum space.

We solve numerically the TDSE using the split-operator method on a Cartesian grid [58]. In 3D, we follow the scheme presented in [59] and divide the three-dimensional configuration space in an inner region and an asymptotic outer region. This approach is related to the mask method introduced in [60, 61]. The binding potential is fully included on the inner grid that spans 409.6 a.u. in each dimension with a spacing of Δx = 0.4 a.u. At the edge of the inner grid, the binding potential is turned off smoothly. In the asymptotic outer region the interaction of the electron with the ionic potential is neglected. Hence, this wave function is represented in momentum space at all times and it is propagated until a final time t_f by means of Volkov states. The inner region has an absorbing boundary covering a distance of 50 a.u. from the edge of the grid. The absorbed parts of the inner wave function are not discarded, but added coherently to the wave function of the outer region. Thus, the momentum-space wave function $~{\psi }(\mathbf{p},{t}_{\text{f}})$ of the ionized system is collected in the outer region.

In 2D, we use the same method as in 3D (with a grid spanning 682.7 a.u. and using Δx = 0.17 a.u.) when the laser parameters are such that the vector potential A(0) is large. For sufficiently low vector potentials A(0), however, we choose a large inner grid spanning 4096 a.u. in each direction, which is enough to contain the wave function till the end of the laser pulse (at time t_p). To obtain the momentum representation $~{\psi }(\mathbf{p},{t}_{\text{f}})$ of the outgoing photoelectron wave packet at a sufficiently large time t_f ≫ t_p, we remove the localized bound states with a mask function in the range r < 20 a.u. Afterwards the remaining part of the wave function ψ(r, t_p) is further propagated using the eikonal approximation

$$\begin{equation}~{\psi }(\mathbf{p},{t}_{\text{f}})\approx \frac{{\text{e}}^{-\mathrm{i}\frac{{p}^{2}}{2}({t}_{\text{f}}-{t}_{\text{p}})}}{2\pi }\int {\mathrm{d}}^{2}\mathbf{r}\enspace {\text{e}}^{-\mathrm{i}\varphi (\mathbf{r},\mathbf{p})}\psi (\mathbf{r},{t}_{\text{p}})\end{equation} \tag{ 2 }$$

with $\varphi (\mathbf{r},\mathbf{p})=\mathbf{p}\cdot \mathbf{r}+{\int }_{{t}_{\text{p}}}^{{t}_{\text{f}}}\mathrm{d}{t}^{\prime }V(\mathbf{r}+\mathbf{p}({t}^{\prime }-{t}_{\text{p}}))$. From the momentum-space wave function $~{\psi }(\mathbf{p},{t}_{\text{f}})$ of the ionized wave packet the PMD is obtained as modulus square: $w(\mathbf{p})=\vert ~{\psi }(\mathbf{p}){\vert }^{2}$. While the results from the 2D TDSE simulations are fully converged with respect to variation of the numerical parameters, the convergence is more difficult to assess for the time-consuming 3D results. In 3D, the numerical capabilities prevent us from further enlarging the used grids in position space. However, we have alternatively solved the TDSE with the generalized pseudospectral method [62, 63] and have calculated the PMD by projection on Coulomb waves. For the lowest intensity used, we have verified that the results from the two numerical methods are in very good agreement. We estimate the numerical error for the presented 3D results to be smaller than 10 percent.

The essence of the classical backpropagation method is first to propagate fully quantum-mechanically the initial state forward in time by solving the TDSE till a final time t_f [31–33]. Afterwards the liberated wave packet is transcribed to a classical phase-space distribution, which is then used as initial distribution for a swarm of classical trajectories evolving backwards in time. In previous implementations, e.g. [29, 31–33], the position-space quantum mechanical wave packet has been used to initiate the classical trajectories and the local-momentum method has been applied to assign a local momentum p(r) to each position r. Here, we follow a slightly different path and we use the momentum-space wave function $~{\psi }(\mathbf{p},{t}_{\text{f}})$ to initiate the classical backpropagation. To this end, we express the wave packet $~{\psi }(\mathbf{p},{t}_{\text{f}})=R(\mathbf{p},{t}_{\text{f}})\mathrm{exp}[\mathrm{i}S(\mathbf{p},{t}_{\text{f}})]$ by means of a real-valued amplitude and phase. For a given momentum p, the corresponding position is obtained by the 'local-position method'

$$\begin{equation}\mathbf{r}(\mathbf{p},{t}_{\text{f}})=-{\nabla }_{\mathbf{p}}S(\mathbf{p},{t}_{\text{f}}).\end{equation} \tag{ 3 }$$

Naturally, this method only works well when interference has a negligible effect on the final electron momentum distribution. For binding potentials with finite support and sufficiently long times t_f, it can be shown that the obtained phase-space distribution is equivalent to the one from previous implementations based on the local-momentum method (see section 3.4 for further details).

To better understand the tunneling dynamics, we identify the exit point of the under-the-barrier motion by propagating the swarm of classical trajectories backward in time following Newton's equation. For each trajectory, the backpropagation is stopped at a time t₀, when the velocity in the direction of the instantaneous electric field E(t₀) vanishes (velocity criterion, see references [31, 32]). This results in an initial distribution w_ini(t₀, v_⊥, v_z ) that can be parameterized by an ionization time t₀, a velocity component in the polarization plane v_⊥ and a velocity component in light propagation direction v_z . The position of the electron at the ionization time t₀ can be interpreted as tunnel-exit position. By construction, the quantum-mechanical momentum distribution at time t_f is exactly reproduced by classical trajectories that are launched according to the distribution w_ini(t₀, v_⊥, v_z ) at the velocity- and time-dependent tunnel-exit positions.

An attoclock-like configuration with quasi-linear fields provides a clean setup to study the momentum dependence of the attoclock shift. To this end, we use a field consisting of a superposition of two counter-rotating circularly polarized fields with the following vector potential

$$\begin{equation}\mathbf{A}(t)=-\frac{2}{3}\frac{{E}_{\text{peak}}}{\omega }\left[\left(\begin{matrix}\hfill \mathrm{cos}(\omega t)\hfill \\ \hfill \mathrm{sin}(\omega t)\hfill \end{matrix}\right)+\frac{1}{4}\left(\begin{matrix}\hfill -\mathrm{cos}(2\omega t)\hfill \\ \hfill \mathrm{sin}(2\omega t)\hfill \end{matrix}\right)\right].\end{equation} \tag{ 4 }$$

For the particular field-strength ratio 2:1 of the fundamental to the second harmonic, the laser field is approximately linearly polarized near the peaks of the electric field [29]. In particular, around t = 0 we can write

$$\begin{align}\mathbf{A}(t)& ={A}_{x}(0){\mathbf{e}}_{x}-\frac{{E}_{\text{peak}}}{{\omega }_{\text{eff}}}\enspace \mathrm{sin}({\omega }_{\text{eff}}t){\mathbf{e}}_{y}+\mathcal{O}({t}^{4}),\\ \mathbf{E}(t)& ={E}_{\text{peak}}\enspace \mathrm{cos}({\omega }_{\text{eff}}t){\mathbf{e}}_{y}+\mathcal{O}({t}^{3}).\end{align} \tag{ 5 }$$

Hence, the electric field resembles a linearly polarized field along the y-axis with an effective frequency ${\omega }_{\text{eff}}=\sqrt{2}\omega$ and a peak field strength E_peak related to the time-averaged intensity of the field by $I=\frac{5}{9}c{{\epsilon}}_{0}{E}_{\text{peak}}^{2}$. The adiabaticity of the ionization process can be quantified by the Keldysh parameter $\gamma =\sqrt{2{I}_{\text{p}}}{\omega }_{\text{eff}}/{E}_{\text{peak}}$. For numerical calculations, the vector potential of equation (4) is multiplied with an envelope f(t) = cos⁴(ωt/(2n_p)) of n_p cycles duration. To avoid the appearance of ATI rings in the relevant part of the momentum distribution, we use three- or five-cycle pulses. The characteristic shapes of vector potential and electric field are depicted in figure 1(a).

Zoom In Zoom Out Reset image size

The PMD obtained from the solution of the 3D TDSE for helium is shown in figure 1(a). It consists of three separated cigar-like regions of high probability that belong to the three maxima of the electric field per optical cycle. For this calculation, an effective wavelength of 800 nm is used, i.e. the fundamental wavelength corresponding to ω is 1131 nm. The global maximum corresponds to the region of almost linear polarization around the peak of the laser pulse at t = 0. In a short-range potential, the global maximum of the distribution would be located at the corresponding negative vector potential −A(0) (indicated as red dot in figure 1). In a long-range potential, the probability distribution is centered in p_x -direction at p_x = −A_x (0) (compare also figure 2(b)), but the whole cigar is shifted towards positive momenta in p_y -direction. This global attoclock shift Δp_max ≈ 0.27 a.u., represented by the maximum of the total distribution, has been investigated in [29].

Zoom In Zoom Out Reset image size

In addition to the global attoclock shift, we can analyse the shift along p_y for each momentum in the p_x –p_z -plane individually. We refer to this quantity as the momentum-dependent attoclock shift. As an example, figure 1(b) shows the signal in the p_z –p_y -plane (at constant p_x = −A_x (0)) with each column being normalized independently. We determine the most probable momentum Δp_y using a Gaussian fit for each momentum p_z . The result is show as black line in figure 1(b). At p_z = 0, this shift is equivalent to the global attoclock shift Δp_max. However, for non-zero p_z we find smaller shifts reaching values as small as ≈0.16 a.u. at p_z = 0.6 a.u. corresponding to only ≈60% of the global attoclock shift.

The momentum-dependent attoclock shift as a function of p_x and p_z is depicted in figure 1(c). Again, at each p_x , p_z , the shift is obtained from a Gaussian fit to the p_y distribution. For the three-cycle laser pulse used here, the maximal attoclock shift of Δp_y ≈ 0.285 a.u. is not found exactly at p_x + A_x (0) = 0 but at a slightly smaller value p_x + A_x (0) ≈ −0.2 a.u. The momentum-dependent attoclock shift decreases as a function of p_x and p_z when going away from this point. We attribute the slight asymmetry of the momentum-dependent attoclock shift as a function of p_x + A_x (0) to a tiny rotation of the probability distribution (by only ≈2° in the polarization plane) that is mostly caused by the deviation from a purely linearly-polarized field near t = 0. Hence, the asymmetry diminishes for longer laser pulses that would be usually used in an experiment, see figure 2(a).

The observed momentum-dependent attoclock shifts can be explained by means of a classical adiabatic model consisting of two steps [20, 64, 65]: (i) laser-induced ionization and (ii) acceleration of the electron in the laser field as well as the ionic potential. The ionization step launches an electron at time t₀ with an initial velocity v₀. In the adiabatic limit γ → 0 the initial velocity v₀ is perpendicular to the electric field E(t₀) at the instant of tunneling, v₀ ⋅ E(t₀) = 0 [64, 66]. Since the electric field is approximately linearly polarized along y near t = 0, the v_y -component vanishes for ionization times t₀ around t = 0 and the setting is nearly isotropic in the v_⊥–v_z -plane which is perpendicular to the electric field E(t₀). Hence, in this case we can write the initial velocity as v₀ ≈ v_⊥ e_x + v_z e_z .

After ionization, the acceleration of the electron by the electric field and the Coulomb force of the parent ion can be described by Newton's equation. It maps the initial conditions (t₀, v₀) to the final momenta p(t₀, v₀). If we first consider the potential-free case, the motion in the laser field only results in p = −A(t₀) + v₀. For our particular bicircular field (4), we find that in the vicinity of t = 0 the release time t₀ is mapped to the p_y -component of the final electron momentum and the initial velocity v₀ ≈ v_⊥ e_x + v_z e_z is mapped to the p_x - and p_z -components [29]. Hence, a change of the initial time t₀ leaves the p_x - and p_z -components approximately unchanged and a change in the initial velocity leaves the p_y -component approximately unchanged. Thus, our bicircular field offers a clean setup to study effects of an initial velocity v₀.

The influence of the Coulomb force leads to an additional momentum change Δp_C that is essential to explain the attoclock shift, because it modifies the mapping to the asymptotic momentum of the electron,

$$\begin{equation}\mathbf{p}({t}_{0},{\mathbf{v}}_{0})=-\mathbf{A}({t}_{0})+{\mathbf{v}}_{0}+{\Delta}{\mathbf{p}}_{\text{C}}({t}_{0},{\mathbf{v}}_{0}).\end{equation} \tag{ 6 }$$

We assume that the electron trajectory starts at an initial position $\mathbf{r}({t}_{0})={\mathbf{r}}_{0}=-{r}_{0}\hat{\mathbf{E}}({t}_{0})$, i.e. at the exit of the tunnel, where $\hat{\mathbf{E}}({t}_{0})=\mathbf{E}({t}_{0})/E({t}_{0})$ with the electric field strength E(t₀) = |E(t₀)|. During its motion in our quasi-linear field the electron is driven away from the parent ion and does not come back close to it (in contrast to pure linearly polarized fields). Hence, we assume that, in our case, it is a good approximation to treat the Coulomb field as a perturbation and evaluate its contribution to the final momentum by integration of the Coulomb force along a trajectory r_L(t) governed solely by the laser field [64]:

$$\begin{equation}{\Delta}{\mathbf{p}}_{\text{C}}=-{\int }_{{t}_{0}}^{\infty }\mathrm{d}t\enspace \nabla V({\mathbf{r}}_{\text{L}}(t)).\end{equation} \tag{ 7 }$$

The Coulomb force decreases rapidly with increasing distance from the ion. Hence, in the adiabatic limit of vanishing Keldysh parameter γ → 0, the change in momentum is acquired in a time much shorter than an optical cycle of the driving light field. In the relevant time window around t₀, the driving field can be assumed to be constant and the light-field-driven trajectory can be approximated as

$$\begin{equation}{\mathbf{r}}_{\text{L}}(t)={\mathbf{r}}_{0}+{\mathbf{v}}_{0}(t-{t}_{0})-\frac{1}{2}\mathbf{E}({t}_{0}){(t-{t}_{0})}^{2}.\end{equation} \tag{ 8 }$$

At the peak of the electric field (corresponding to t = 0) the derivative of the electric field $\dot {\mathbf{E}}(0)$ vanishes and, hence, the leading-order correction to equation (8) also vanishes. Thus, for times near t = 0, the light-field-driven trajectory of equation (8) is a good approximation. Inserting equation (8) in equation (7), the momentum change due to the Coulomb force in a bare −Z/r potential is given by (to second order in the initial velocity v₀)

$$\begin{align}{\Delta}{\mathbf{p}}_{\text{C}}=& \left[\frac{\pi }{2}\frac{Z}{\sqrt{2{r}_{0}^{3}E({t}_{0})}}-\frac{3\pi }{16}\frac{Z{\mathbf{v}}_{0}^{2}}{\sqrt{2{r}_{0}^{5}{E}^{3}({t}_{0})}}\right]\hat{\mathbf{E}}({t}_{0})\\ & -\frac{1}{2}\frac{Z}{{r}_{0}^{2}E({t}_{0})}{\mathbf{v}}_{0}.\end{align} \tag{ 9 }$$

The momentum correction of equation (9) can be viewed as an expansion of the exact result in terms of the potential strength, represented by the charge Z, and the Keldysh parameter γ. The second term of equation (9) leads to some focusing in the velocity coordinate v₀ in the sense that the momentum distribution gets narrower in p_x - and p_z -directions compared to the distribution of initial conditions. In contrast, the first term of equation (9) points in the direction of the electric field E(t₀). In the quasi-linear field, electrons that are liberated at the peak of the electric field (corresponding to t₀ = 0) are thus deflected in p_y -direction. Within this simple model the momentum-dependent attoclock shift introduced in section 3.1 is given by Δp_C,y of equation (9). Obviously, this shift strongly depends on the tunnel-exit position r₀ of the electron and, hence, the attoclock can be used as a 'nano-ruler' to investigate the tunnel-exit position.

In previous classical models it has usually been assumed that the tunnel-barrier width r₀ depends only on the time t₀ via the electric field strength E(t₀) (compare, e.g. the reviews [37, 38] on the attoclock). In the simplest approach it is assumed that the electron tunnels adiabatically in a 1D cut along the direction of the laser field [67]. This is referred to as field-direction model (FDM). After a separation of coordinates a one-dimensional potential barrier is formed by an atomic potential −Z/r and the electric field of the laser. The other directions lead only to an increased effective ionization potential ${I}_{\text{p}}^{\text{eff}}={I}_{\text{p}}+\frac{1}{2}{\mathbf{v}}_{0}^{2}$ such that the tunnel-exit position is given by [20]

$$\begin{equation}{\mathbf{r}}_{0}=-\frac{{I}_{\text{p}}^{\text{eff}}+\sqrt{{({I}_{\text{p}}^{\text{eff}})}^{2}-4ZE({t}_{0})}}{2E({t}_{0})}\hat{\mathbf{E}}({t}_{0}).\end{equation} \tag{ 10 }$$

The barrier width as a function of the v_z -component of the initial velocity is shown in figure 3(a). Here, a non-zero initial velocity leads to an increased width of the tunnel barrier. If the influence of the central potential at the tunnel exit is neglected (Z = 0), the tunnel-exit position simplifies to the result for a triangular barrier (referred to as ${I}_{\text{p}}^{\mathrm{e}\mathrm{ff}}/E$ model)

$$\begin{equation}{\mathbf{r}}_{0}=-\frac{{I}_{\text{p}}+\frac{1}{2}{\mathbf{v}}_{0}^{2}}{E({t}_{0})}\hat{\mathbf{E}}({t}_{0}).\end{equation} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

The attoclock shift as a function of the p_z -component of the momentum extracted from the numerical solution of the TDSE for an intensity of 4 × 10¹⁴ W cm⁻² is shown in figure 3(b) together with the classical momentum shift Δp_C,y that is estimated by equation (9) for different choices of the tunnel exit. In the classical model, larger tunnel-exit positions lead to a weaker influence of the Coulomb force on the outgoing electron and hence smaller attoclock shifts. Even though the classical momentum change Δp_C,y for the FDM tunnel exit is overall slightly too small, the variation of the attoclock shift with momentum p_z is well reproduced. However, if the velocity dependence of the tunnel-exit position is artificially neglected by setting ${I}_{\text{p}}^{\text{eff}}={I}_{\text{p}}$, the variation of the attoclock shift as a function of p_z is drastically reduced. The difference between the shifts from the classical model including the velocity-dependent tunnel exit and the attoclock shifts from numerical TDSE calculations are at least partially caused by the adiabatic choice of the tunnel-exit position, the perturbative evaluation of the Coulomb correction and the approximated field-driven trajectory.

A simple formula for the momentum change of the classical model as a function of the initial velocity can be obtained (to second order in v₀) by using the ${I}_{\text{p}}^{\text{eff}}/E$-tunnel-exit position of equation (11), leading to

$$\begin{equation}{\Delta}{p}_{\mathrm{C},y}=\frac{\pi ZE({t}_{0})}{{(2{I}_{\text{p}})}^{3/2}}-\frac{9\pi }{4}\frac{ZE({t}_{0})}{{(2{I}_{\text{p}})}^{5/2}}{\mathbf{v}}_{0}^{2}.\end{equation} \tag{ 12 }$$

The first term describes the global shift at v₀ = 0 or equivalently at the final momenta p_x + A_x (0) = 0 and p_z = 0, whereas the second term leads to the decrease of the shift as a function of p_x + A_x (0) and p_z . Within this model, two thirds of the velocity-dependent term of equation (12) result from the dependence of the tunnel-barrier width r₀ of equation (11) on the velocity and one third results from the linear velocity dependence of the trajectory (second term in equation (8)).

The classical model introduced in this section rests on two basic assumptions. First, both tunnel ionization and the subsequent Coulomb-induced momentum change of the outgoing electron are treated adiabatically, i.e. the electron velocity at the tunnel exit has no component along the ionizing field and the field is assumed to be constant during the short phase of Coulomb-induced acceleration. Second, the idea of the model is that the attoclock shifts in the final momentum distribution, which are obtained by peak search, correspond to ionization at the field maximum. We have attempted to refine this model by sampling the entire distribution of initial times and velocities in a classical trajectory Monte Carlo (CTMC) simulation as it has been successfully used to model the global attoclock shifts, e.g. in references [25, 39]. Here, the trajectories are weighted according to the Ammosov–Delone–Krainov rate [68] and for all trajectories, Newton's equation of motion are solved numerically with full inclusion of the external electric field and the Coulomb force. Afterwards, the attoclock shifts are extracted from the obtained final momentum distributions. In general, the CTMC method is expected to work well in the setup considered in the present work because interference between trajectories does not play a role. However, as far as the momentum-dependent attoclock shifts are concerned, we have found that the agreement with the TDSE results is not fully convincing as the momentum dependence is overestimated. We believe that this deficiency may be due to slight inconsistencies in the model such as the ad-hoc combination of initial distributions from quantum mechanical theories with the classical Newton's equation. Therefore, we do not show results from the CTMC method. In contrast, the simple classical model, consisting of the equations above, works surprisingly well as shown by our results below. Therefore, the discussion in the remainder of this paper focuses on the idea that the attoclock shifts correspond to ionization at the field maximum without taking the full distribution of initial conditions into account.

In the following, we employ the models of the attoclock shift for ionization occurring at the peak of the electric field, i.e. at t₀ = 0 and E(t₀) = E_peak. In this case, the simple formula of equation (12) suggests that the classical momentum change in the Coulomb field is a linear function of the electric field strength E_peak. For a systematic comparison between the models and the TDSE, we perform a series of calculations with varying peak field strength while keeping the vector potential A_x (0) fixed (by adjusting the wavelength). For sufficiently weak fields, we find indeed a linear field-strength dependence of the global attoclock shift Δp_max extracted from TDSE, see figure 4(a). The slightly increased slope observed in TDSE (compared to equation (12)) is partially caused by the approximated tunnel-exit position of equation (11) (compare the result using the FDM tunnel exit) and partially by non-adiabaticity: the Keldysh parameters γ are in the range between 0.28 and 0.57. For very strong fields, depletion effects play an important role so that the peak of the ionization rate shifts to earlier times and hence the attoclock shift from TDSE decreases [28, 31, 32].

Zoom In Zoom Out Reset image size

To quantify the relative variation of the attoclock shift as a function of the momentum by a single figure of merit, we define the relative difference M between the global attoclock shift and an outer attoclock shift

$$\begin{equation}M=\frac{{\Delta}{p}_{\mathrm{max}}-{\Delta}{p}_{\text{out}}}{{\Delta}{p}_{\mathrm{max}}}\end{equation} \tag{ 13 }$$

with Δp_out being the attoclock shift at a given value of v_out = p_x + A_x (0). In general, other measures could also be used, e.g. the curvature or the width of the momentum dependence, that would all allow for similar conclusions. For v_out between 0.3 a.u. and 0.5 a.u., the same characteristic behavior of M is observed in the TDSE calculations. We present the results for v_out = 0.4 a.u. as a function of the field strength E_peak in figure 4(b). The simple formula of equation (12) predicts that the relative difference M of the attoclock shift is independent of the field strength, see the gray horizontal line in figure 4(b). However, we find that the relative difference M from TDSE increases monotonically as a function of the electric field. This observation is in agreement with the classical momentum change of equation (9) when combined with the FDM tunnel exit of equation (10). Nevertheless, we note that for all considered field strengths, the numerical solutions of the TDSE suggest in the adiabatic limit a slightly smaller variation M compared to the classical model of equation (9).

The main ingredient of the classical ad-hoc model presented in section 3.2 is the dependence of the tunnel-exit position r₀ on the initial velocity v₀ of the electron. The classical backpropagation introduced in section 2.2 offers the opportunity to extract such a correlation between the electron's initial conditions from ab initio calculations. Here, we again restrict ourselves to trajectories that depart from the tunnel barrier at the peak of the electric field strength, i.e. we analyse only those backpropagation trajectories for which the backpropagation stops at t₀ = 0. Additionally, due to the nearly isotropic setting in the plane perpendicular to the direction of the instantaneous electric field it is sufficient to study the dependence along one coordinate axis of the initial velocity.

The y-component of the initial position parallel to the electric field E(t₀) is shown in figure 5(a) as a function of initial velocity v_⊥ in the polarization plane for an intensity of 10¹⁵ W cm⁻². The backpropagation result for a long-range potential is in excellent agreement with the FDM of equation (10). As long as the ionization process is sufficiently adiabatic and over-the-barrier ionization is unimportant, this good agreement is also observed for a broad range of intensities, see figures 4(c) and (d). Increasing non-adiabaticity leads to slightly smaller exit positions of the electron compared to the static tunneling limit in agreement with references [31, 69].

Zoom In Zoom Out Reset image size

The initial conditions from backpropagation can be used to start classical trajectories (propagating forward in time following Newton's equation). In analogy to the model in section 3.2, we determine numerically the momentum change Δp_C introduced in equation (6) due to the influence of the Coulomb force for these trajectories. The value of Δp_C,y for v₀ = 0 as well as the corresponding relative variation are also shown in figures 4(a) and (b). Besides the plotted results for t₀ = 0, we have confirmed that the values for the momentum change Δp_C,y depend only weakly on the exact initial time t₀ in the vicinity of the peak of the electric field strength (not shown). For sufficiently low intensities (such that effects like depletion can be neglected) the global attoclock shifts extracted from TDSE follow indeed nearly the classical momentum change Δp_C,y from backpropagation. Here, the TDSE solution is consistent with the assumption of vanishing ionization times. However, for large intensities the classical momentum change Δp_C,y is larger than the global attoclock shift from TDSE. This indicates that the most probable ionization time is negative in this case which is in agreement with reference [29]. Nonetheless, the results from backpropagation reproduce the variation M of the attoclock shift very well, see figure 4(b). These theoretical findings support the interpretation that the momentum-dependent attoclock shifts are strongly affected by the velocity-dependent tunnel-exit positions.

Additionally, we have also studied the tunnel-exit positions for a 2D short-range potential

$$\begin{equation}V(r)=\frac{-1.37\enspace \mathrm{exp}(-{\left(r/3\right)}^{4})}{\sqrt{{r}^{2}+0.34}}\end{equation} \tag{ 14 }$$

with the correct ionization potential of helium. The y-component of the initial position from backpropagation is in very good agreement with the simple ${I}_{\text{p}}^{\text{eff}}/E$ estimate of equation (11) and a quadratic dependence on v_⊥ is clearly visible, see figure 5. This is expected as the strong-field approximation (SFA) (in saddle-point approximation) [66, 70, 71] typically describes the ionization dynamics for short-range potentials very well. Within this theoretical framework, the momentum-space based backpropagation can be carried out analytically (to exponential accuracy) resulting in the well-known exit position

$$\begin{equation}{\mathbf{r}}_{\text{SFA}}=\mathrm{R}\mathrm{e}{\int }_{{t}_{s}}^{\mathrm{R}\mathrm{e}\enspace {t}_{\text{s}}}\enspace \mathrm{d}t\left[\mathbf{p}+\mathbf{A}(t)\right],\end{equation} \tag{ 15 }$$

where t_s is the complex-valued saddle-point time. In the adiabatic limit this expression reduces to the estimate of equation (11) [33].

Both the FDM, equation (10), and the SFA, equation (15), imply that electrons tunneling at the peak field strength (t₀ = 0) appear in the direction of the instantaneous field E(t₀) (here the y-direction). For short-range potentials, this assumption is consistent with the backpropagation results as the obtained x-component of the tunnel-exit position (that is perpendicular to the electric field E(t₀) at the instant of tunneling) vanishes, see figure 5(b). Surprisingly, for long-range potentials, the backpropagation results in non-zero x-components of the tunnel-exit position for v_⊥ ≠ 0. The x-components take values as large as 1 a.u. for v_⊥ = −0.6 a.u. This correlation of initial momentum and initial position is not limited to the initial velocities in the polarization plane but we have observed it in any direction that is perpendicular to the direction of the laser electric field at the instant of tunneling. Moreover, we also have found this correlation for other lights fields (e.g. circularly polarized light) which shows that this finding is not limited to our bicircular field. In additional calculations, we could show that the absolute value of these additional position offsets depends only weakly on the intensity of the field (not shown). Hence, we expect that for sufficiently weak laser fields this component is unimportant compared to the (large) tunnel-barrier width in field direction (y-component).

The backpropagation results presented above have been calculated using the method introduced in section 2.2 which is based on the momentum-representation of the wave function. Compared to the previous implementations [32, 33] based on the position representation, our implementation is beneficial for two reasons: (i) only a relatively small region of the position space needs to be represented explicitly. This is especially beneficial for long-wavelength driving fields and in 3D calculations. (ii) The result for the tunnel-exit position is in good approximation independent of the time t_f when the forward quantum propagation is stopped. A comparison between the two backpropagation versions in figure 5(b) shows that the result using the position-space implementation [32, 33] depends on the final time t_f. For long times t_f → ∞ it converges slowly against the momentum-space result. In the opinion of the authors, the momentum-space implementation appears to be the correct result for the backpropagation scheme.