Electronic nonadiabatic dynamics in enhanced ionization of isotopologues of hydrogen molecular ions from the exact factorization perspective
Received
14th December 2016
, Accepted 14th February 2017
First published on 9th March 2017
It was recently shown that the exact potential driving the electron's dynamics in enhanced ionization of H_{2}^{+} can have large contributions arising from dynamic electron–nuclear correlation, going beyond what any Coulombicbased model can provide. This potential is defined via the exact factorization of the molecular wavefunction that allows the construction of a Schrödinger equation for the electronic system, in which the potential contains exactly the effect of coupling to the nuclear system and any external fields. Here we study enhanced ionization in isotopologues of H_{2}^{+} in order to investigate the nuclearmassdependence of these terms for this process. We decompose the exact potential into components that naturally arise from the conditional wavefunction, and also into components arising from the marginal electronic wavefunction, and compare the performance of propagation on these different components as well as approximate potentials based on the quasistatic or Hartree approximation with the exact propagation. A quasiclassical analysis is presented to help analyse the structure of different nonCoulombic components of the potential driving the ionizing electron.
I. Introduction
The phenomenon of chargeresonance enhanced ionization (CREI), predicted more than twenty years ago,^{2–4} is a prominent example of the complex coupling of electronic motion, ionic motion, and strong laser fields. At a critical range of internuclear separations, the ionization rate of a molecule in a laser field can be orders of magnitude greater than the rate from the constituent atoms. The ionization rate has been explained by a quasistatic argument involving instantaneously frozen nuclei in the pioneering works of ref. 3–9 for which the timedependent Schrödinger equation (TDSE) is solved for various clamped nuclear (cn) configurations , i.e., 
 (1) 
where 
 (2) 
Here and are used to collectively denote the electronic and nuclear coordinates, _{e} is the electronic kinetic energy operator, and Ŵ_{ee} is the electron–electron repulsion. Furthermore, contains the electron–nuclear interaction, which, for a diatomic molecular ion with one electron, as will be considered here, has the form −Z_{1}/r − R_{1} − Z_{2}/r − R_{2}, i.e. a Coulombic doublewell. The critical internuclear separation for CREI can be qualitatively explained by the following argument: at stretched geometries in a given static field, the energy levels in the upfield atom Starkshift upwards while the inner barrier from the internuclear Coulombic potential also grows (Fig. 2 in ref. 3). Provided that the field is turned on fast enough such that any population in the upfield level (LUMO) does not significantly tunnel back to the downfield atomic level (HOMO), the molecule can rapidly ionize over both the inner and outer barriers, which gives rise to the ionization rate observed to be enhanced by orders of magnitude compared to the atomic rate. By requiring that the Starkshifted LUMO level exceeds the top of both the inner and outer fieldmodified Coulombic barriers, one finds R_{c} = 4.07/I_{p} for the critical internuclear separation for CREI. The analysis can be generalized to the case of a laser field, where the field's period is shorter than the tunneling time.^{2,3,10,11}
It has been pointed out that, in actuality, the underlying assumption in this picture of the electron adiabatically following the field is not quite adequate, as the ionization tends to occur in multiple subcycle ionization bursts, not at the peak of the field cycles.^{10,11} Moreover, when applied to an experiment where the molecule's initial geometry is far from where the CREI is expected to occur, the premise of the quasistatic picture can become a little shaky: in particular to observe CREI, the molecule must be stretched to the CREI region rapidly enough such that appreciable electron density still remains largely unionized. In fact, in many experiments, CREI is not observed because too little electron density reaches the CREI region over the course of the applied field.^{12–14} Typically a large fraction of the nuclear density remains near equilibrium, while a part of it dissociates. So, (i) representing its potential on the electron as a Coulombic doublewell is not appropriate, and (ii) the fragment velocities can be comparable to the electronic velocities in some dissociating channels questioning the very notion of the quasistatic description of nuclear–electron interactions. In ref. 1, it was shown that the nuclear dynamics can contribute significant components to the exact potential that arises from correlated electronnuclear motion. Neglecting these contributions leads to severe errors in the prediction of the electronic motion.
In this paper, we analyse the exact potential in more detail by providing a complementary decomposition to the one that was studied in ref. 1, comparing with different approximations, studying the massdependence, and investigating a quasiclassical treatment. The rest of the paper is organized as follows: first, a brief review of the exact factorization approach is presented in Section II, with a focus on the reverse factorization introduced in ref. 15. This provides us with a TDSE for electrons evolving on a single timedependent potential that accounts for the coupling to the nuclear system and any external field in an exact way. We present two different ways of decomposing this potential and some approximate potentials based on conventional approximations. In Section III, we study a onedimensional model for the (a)symmetric isotopologues of H_{2}^{+} subject to a linearly polarized laser field for two different situations, compared with dynamics of different approximate potentials. In Section IV, we take a first step to analyse quasiclassically the structure of various components of the potential driving the ionizing electron. Finally, some concluding remarks are presented in Section V.
II. Exact factorization approach
The exact factorization of the timedependent electron–nuclear wavefunction introduced in ref. 15–19 enables the rigorous definitions of exact potentials acting on the nuclear subsystem and the electronic subsystem in coupled electron–ion dynamics. These potentials follow from writing the full molecular wavefunction as a single product , or , where partial normalization conditions on the conditional wavefunctions and , respectively, render each factorization unique up to a gauge transformation. In the first product the equation for the nuclear wavefunction, , is a TDSE, while in the second product the equation for the electronic factor is a TDSE. The TDSEs in both cases contain a vector potential and a scalar potential that in the former case provides us with the exact manybody nuclear density and current density of the complete electron–nuclear system while in the latter case it leads to the exact manybody electronic density and current density of the whole system. When the electronic dynamics is particularly of interest, as in our present study of ionization, we focus on the second factorization, and investigate the potential that appears in the TDSE for .
If we consider the case of one electron coupled to one nuclear degree of freedom in one dimension, the equation for the electronic wavefunction is (in 1D we replace by z and R respectively):

 (3) 
and that for the nuclear conditional wavefunction is:

[Ĥ_{n}(z, R, t) − ε_{e}(z, t)]χ_{z}(R, t) = iħ∂_{t}χ_{z}(R, t),  (4) 
with the nuclear Hamiltonian

 (5) 
where
_{n} is the nuclear kinetic energy operator,
Ŵ_{en}(
z,
R) (
Ŵ_{nn}(
R)) is the electron–nuclear (nuclear–nuclear) interaction, and
^{laser}_{n}(
R,
t) is the timedependent external potential acting on the nuclei. Here,
is the exact timedependent vector potential

 (6) 
and
ε_{e}(
z,
t) = 〈
χ_{z}(
t)
Ĥ_{n} −
iħ∂
_{t}
χ_{z}(
t)〉
_{R} is the exact timedependent electronic potential for an electron (see the next section). In one dimensional models,
can be set to zero as a choice of gauge and we adopt this gauge for the rest of this paper. In this case,
ε_{e}(
z,
t) is the sole potential that drives the electronic motion, which can be compared with the other traditional potentials that are used to study electronic dynamics.
A. Decomposition of the exact timedependent potential energy surface
The exact timedependent potential energy surface for electrons (eTDPES) contains the effects of the coupling to moving quantum nuclei as well as the external laser field. It can be written as 
 (7) 
which consists of: the approximate potential 
ε^{app}_{e}(z, t) = ^{laser}_{e}(z, t) + 〈χ_{z}(R, t)Ŵ_{en}(z, R) + Ŵ_{nn}(R)χ_{z}(R, t)〉_{R},  (8) 
the nuclearkinetic term 
 (9) 
the gauge dependent part of the potential
ε^{gd}_{e}(z, t) = ħ〈χ_{z}(R, t) − i∂_{t}χ_{z}(R, t)〉_{R}, 
and finally the electronickineticlike contribution
In ref. 1, we had found that this exact potential ε_{e}(z, t) has significant features that are missing in the traditional potentials based upon the quasistatic picture described above. Neglecting these features led to qualitatively incorrect predictions of ionization dynamics in the H_{2}^{+} molecule in strong fields. In ref. 1 it was found that, in general, all the four terms above are needed to reasonably reproduce the ionization in CREI processes. i.e. that propagation on any combination other than the full sum of the four contributions in eqn (7) gave qualitatively poor results. This means that in eventually developing approximations, all of the four terms must be considered.
Alternatively, one may instead decompose the eTDPES by exactly inverting the TDSE for the electronic wavefunction. The electronic wavefunction may be written in polar representation,

 (10) 
with
n_{e}(
z,
t) = 
Φ(
z,
t)
^{2} and
where

 (11) 
is the electronic current–density (which is equivalent to the true electronic current–density calculated from the full electron–nuclear wavefunction, when the vector potential can be set to zero. This is the case for our onedimensional model system. Here,
∇_{e} =
∂/∂
_{z}). Inserting this form into
eqn (3) and setting the timedependent vector potential to zero give

 (12) 
The first term is what survives in the absence of any dynamics, and it depends only on the instantaneous electron density; we denote it “adiabatic” in the spirit of timedependent density functional theory. The second term, denoted the “velocity term”, depends only on the electron velocity, namely on j_{e}(z, t)/n_{e}(z, t), while the last term, denoted the “acceleration term” depends on the spatial integral of the acceleration. In Section III B and C, we consider propagation on different contributions of the exact potential defined by the decomposition of eqn (12), again with a view of eventually developing approximations, possibly density functionals (see ref. 20), for the exact eTDPES.
B. Approximate electronic potentials based on conventional approximations
The standard approximation for the electronic potential assumes the socalled quasistatic approximation (qs) that treats the nuclei as classical point particles with positions that are either considered fixed, _{0}, as in the majority of studies of CREI,^{3,5–7,9,21} or move classically with classical trajectories (t) that are often described by mixed quantumclassical algorithms such as Ehrenfest or surfacehopping algorithms.^{22} In these methods, electrons on the other hand, regardless of whether the nuclei are frozen or move classically, follow the combined potential from the laser field and the electrostatic attraction of the nuclei, i.e., 
ε^{qs}(z, t(t)) = Ŵ_{en}(z, (t)) + ^{laser}_{e}(z, t).  (13) 
Considering the exact electronic potential eqn (7), we see that such approaches completely miss the dynamic electron–nuclear correlation effects contained in , , and ε^{gd}, and can be viewed as an approximation to ε^{app} alone: ε^{app} reduces to the qs approximation when the conditional nuclear wavefunction is approximated classically as a zindependent deltafunction at (t), i.e. n_{z}(R, t) ≈ δ(R(t), (t)) (in this limit Ŵ_{nn}((t)) becomes purely a timedependent constant and hence is dropped hereafter).
A step beyond the qs approximation for the electronic potential that can, in principle, account for the width and splitting of the nuclear wavefunction is the Hartree approximation^{23}

 (14) 
where
n(
R,
t) is the nuclear density obtained from nuclear wavepacket dynamics.
^{24} It can be easily seen that if the
zdependence in the conditional nuclear wavefunction is neglected,
i.e.,
n_{z}(
R,
t) ≈
n(
R,
t), the approximate potential simplifies to the Hartree approximation.
To provide a detailed analysis of the exact electronic potential, we compare the electronic dynamics resulting from different approximations. We will consider combinations of the terms of the exact potential decomposed according to eqn (12) to complement the analysis in ref. 1 given in terms of the decomposition in eqn (7), as well as the electronic dynamics resulting from the following approximations: (i) The qs approximation for which (t) = 〈R〉(t) is the average timedependent internuclear separation obtained from the exact calculations; (ii) The Hartree approximation for which the nuclear density in (14) is replaced by the exact timedependent nuclear density, indicated by ε^{exH}; (iii) The normalized Hartree, defined as

 (15) 
indicated as
ε^{nexH} when the exact nuclear density, obtained from the complete electron–nuclear wavefunction
Ψ(
z,
R),
i.e. 
 (16) 
is inserted in
eqn (15). This is to account for the loss of norm, as we use a mask function, where
Ψ(
z,
R,
t) is absorbed by a complex potential that is added to the full Hamiltonian;
^{25–27} (iv) The “selfconsistent” Hartree approximation (SCH) in which the full wavefunction is approximated as an uncorrelated product of electronic wavefunction and nuclear wavefunction,
, that treats the electrons and nuclei on the same footing: the electronic dynamics is described with the potential in
eqn (14) that is coupled to the nuclei that move under the influence of the analogous potential,
i.e.,

 (17) 
where
n_{e}(
z,
t) is the electron density obtained from the electronic dynamics. This approach does not involve the complications of dealing with the conditional nuclear wavefunction as in
ε^{app} but it fails to capture major effects of the correlated electron–nuclear dynamics due to its meanfield nature. Finally, we point out that
ε^{app} as well as all the approximations (i)–(iii) represent the electron–nuclear interaction in a Coulombic way only.
III. CREI in H_{2}^{+} isotopologues
We utilize a popular onedimensional model of the symmetric as well as asymmetric isotopologues of H_{2}^{+} subject to a linearly polarized laser field. As the motion of the nuclei and the electron in the true molecule is assumed to be restricted to the direction of the polarization axis of the laser field, the essential physics can be captured by a 1D Hamiltonian featuring “softCoulomb” interactions:^{28,29} 
 (18) 
where R and z are the internuclear distance and the electronic coordinate as measured from the nuclear centerofmass, respectively. The nuclear effective mass is denoted as while is the electronic reduced mass with M_{n} = M_{1} + M_{2}. The laser field, within the dipole approximation, is represented by 
^{laser}(R, z, t) = eE(t)(q_{e}z − ζR),  (19) 
where E(t) denotes the electric field amplitude is the reduced charge, and ζ = (M_{2} − M_{1})/M_{n} is the massasymmetry parameter. Such reduceddimensional models have been shown to qualitatively reproduce the experimental results (see ref. 30 for example). The parameters employed for the numerical simulations that are the same for all of the cases studied here are given in Table 1.
Table 1 Parameters employed for the numerical simulation that is performed on the realspace grid. The parameters are in atomic units. The timestep used for the integration of TDSE is given by Δt
z
_{min}

z
_{maz}

R
_{min}

R
_{max}

Δz 
ΔR 
Δt 
−204.6 
204.6 
0 
51.1 
0.4 
0.1 
0.0115026 
We study the symmetric isotopologues i.e., H_{2}^{+}, A_{2}^{+}, D_{2}^{+}, x_{2}^{+} as well as the asymmetric isotopologues HD^{+}, HT^{+}, Hx^{+}, HX^{+}. Here x(X) stands for the fictitious isotope of hydrogen that is 10(100) times heavier than that of H, while A_{2}^{+} is another fictitious isotopologue with the same effective nuclear mass as HT^{+} (see Table 2 in which the nuclear effective mass values of H_{2}^{+} and isotopologues are given).
Table 2 Nuclear effective mass corresponding to H_{2}^{+} and its symmetric and antisymmetric isotopologues in atomic units. Here, x(X) refers to the 10(100) times heavier fictitious isotope of hydrogen and A_{2}^{+} is another fictitious isotopologue with the same effective nuclear mass as HT^{+}
H_{2}^{+} 
HD^{+} 
HT^{+}(A_{2}^{+}) 
Hx^{+} 
HX^{+} 
D_{2}^{+} 
x_{2}^{+} 
918.076 
1223.742 
1376.228 
1669.229 
1817.973 
1834.533 
9180.76 
A. Comparison of ionization of isotopologues of H_{2}^{+} with different effective nuclear mass for a fixed field
We apply a field of duration 50cycles, wavelength λ = 800 nm (ω = 0.0569 a.u.) and intensity I = 2.02 × 10^{14} W cm^{−2}, with a sinesquared pulse envelope (Fig. 1), to each of the H_{2}^{+} isotopologues. First, we solve the electron–nuclear TDSE for Hamiltonian of eqn (18) numerically exactly, beginning in the initial groundstate of the molecule.

 Fig. 1 Laser field as a function of number of optical cycles t/T. The electric field amplitude is divided by the peak amplitude, .  
As the Hamiltonian (18) has been obtained after separating off the center of mass motion and the origin is set to be the nuclear center of mass the field couples directly to the nuclear motion only in the asymmetric cases; in the symmetric case, nuclear motion is driven purely by the electronic dynamics. This is also clear from eqn (19) where in the symmetric case, ζ = 0. In Fig. 2 we plot the ionization probability and the average internuclear distance, 〈R〉, as a function of number of cycles t/T where T denotes the duration of one cycle (T = 2.67 fs) for H_{2}^{+}, HD^{+}, HT^{+}, Hx^{+}, HX^{+}, A_{2}^{+}, D_{2}^{+} and x_{2}^{+}. Note that, the results are given in atomic units, e = m_{e} = h = 1, throughout the article, unless otherwise noted. The ionization probability is defined as,

 (20) 
with
z_{I} = 15 a.u. As can be seen in
Fig. 2, the ionization and average internuclear separation appear to be practically independent of the nuclear masssymmetry properties: they are identical for the asymmetric case of HT
^{+} and the symmetric case of A
_{2}^{+} systems with the same effective nuclear mass
μ_{n} and are very similar in the case of HX
^{+} and D
_{2}^{+} where the effective masses are only a little different. One can understand this from the fact that ionization probes the electron density in regions where
z ≫
R, while the nuclear massdependence in the Hamiltonian
eqn (18) to the lowest order in
R/
z is only via
μ_{n}, and the integrated nature of the observable reduces its sensitivity to the details of the distribution which is affected by the symmetry of the system.

 Fig. 2 Ionization probability and the average internuclear distance,〈R〉, as a function of number of cycles t/T (where T denotes the duration of one cycle (T = 2.67 fs)) of the H_{2}^{+}, HD^{+}, HT^{+}, Hx^{+} and HX^{+} molecules on the left panels and H_{2}^{+}, A_{2}^{+}, D_{2}^{+} and x_{2}^{+} on the right panels.  
It is also observed that as the nuclear effectivemass increases, the ionization decreases, with a dependence that appears to tend to 1/μ_{n} for intermediate masses. Furthermore, the growth over time of the average internuclear distance is less as the mass increases, suggesting that for larger masses the system reaches the critical internuclear distance for CREI at later times when the field has decreased from its peak value significantly, consequently resulting in lower ionization. To investigate whether the CREI process is still relevant or not and shed more light on the dependence of the ionization on the effective nuclear mass we utilize the concept of the timeresolved, Rresolved ionization probability defined as , with and z_{I} = 15 a.u.^{1,30} This quantity can be rewritten in terms of the concepts of exact nuclear wavefunction χ(R, t) and exact conditional electronic wavefunction Φ_{R}(z, t) introduced within the exact direct factorization framework, i.e.,

I(R, t) = χ(R, t)^{2}I_{cp}(R, t),  (21) 
where,
I_{cp}(
R,
t) follows the usual expression of the ionization probability but using the exact conditional electronic wavefunction
Φ_{R}(
z,
t),
i.e.,
that is coupled to the exact nuclear dynamics. Therefore,
I(
R,
t), which is the nucleardensityweighted conditional ionization probability, is analogous to the ionization probability calculated for a given nuclear configuration R in a quasistatic picture.
To analyze the dynamics of different cases, in Fig. 3 we have plotted I(R, t) along with the timedependent nuclear density for different isotopologues. As seen in the upper three rows of Fig. 3 with increasing effective nuclear mass the peak of I(R, t) moves slightly to larger times and its overall value decreases while remaining close to the CREI region as predicted by the internuclear separation of R_{c}. The lower set of panels of Fig. 3 shows the timedependent nuclear density in each case. For the case of x_{2}^{+} where only an exponentially small fraction of the nuclear density dissociates, the system does not reach the CREI regime and therefore the ionization is negligible as seen from the Rresolved, tresolved ionization probability. Notice the different scale of I(R, t) for the case of x_{2}^{+}.

 Fig. 3 The timeresolved, Rresolved ionization probability, I(R, t), (upper panels) and the nuclear density (lower panels) as a function of number of cycle for the asymmetric case in the left panels HD^{+}, HT^{+}, Hx^{+} and the symmetric case in the right panels H_{2}^{+}, D_{2}^{+} and x_{2}^{+}.  
Note that the ionization rate computed at any fixed R would be identical for all the isotopologues. Their different masses, however, lead to very different ionization rates when the full electron–nuclear dynamics is considered with the molecule beginning at equilibrium as we have shown. Still, there is some validation to the original quasistatic CREI prediction that ionization is enhanced for nuclear separations around R_{c}, as indicated by I(R, t), however, a modification to the statement is needed due to the spreading and splitting of the nuclear wavepacket: the enhancement occurs from electrons associated with the part of the nuclear density that is in the R_{c} region. Clearly, treating the nuclei as point particles will not work even for significant nuclear mass (except in the largemass limit) because of this. The question then arises whether accounting for the nuclear distribution is enough to capture CREI accurately, for example using one of the Coulombbased approximations of Section II, and, whether and how the errors decrease with the nuclear mass. To this end, we next consider adjusting the field strength so that the ionization is similar for the different isotopologues.
B. Comparison of potentials with similar ionization rates
As discussed in the previous section, the different isotopologues of H_{2}^{+} subject to the same field show significantly different degrees of ionization, since the differentmass systems reach the CREI region at different times with different probabilities. In particular, systems with larger nuclear masses hardly reach the internuclear distances for which CREI is expected to occur. Hence, to study the effect of the nuclear mass in the CREI regime for differentmass systems, we adjust the field intensity such that the ionization probability(rate) remains close to that of H_{2}^{+}. As the asymmetric and symmetric isotopologues with the same effective nuclear mass subject to the same field give almost the same ionization probability, average internuclear distance and I(R, t) (see Section III A), from here on we focus only on the symmetric isotopologues of H_{2}^{+}. The (crudely)optimized field intensity for differentmass systems is given in Table 3 while the other field parameters are kept unchanged.
Table 3 The (crudely)optimized field intensity for differentmass systems in the unit of W cm^{−2}
H_{2}^{+} 
A_{2}^{+} 
D_{2}^{+} 
x_{2}^{+} 
2.02 × 10^{14} 
2.14 × 10^{14} 
2.26 × 10^{14} 
3.1 × 10^{14} 
The ionization probability of the symmetric isotopologues subject to the (crudely)optimized field intensities is depicted in Fig. 4. For isotopologues heavier than H_{2}^{+} the ionization starts to set in about one optical cycle T later than the H_{2}^{+} case. In order to be able to compare the ionization yield (rate) of these systems with H_{2}^{+} visually better, in Fig. 4 we have shifted the timedependent ionization probability by one optical cycle, i.e. for isotopologues heavier than H_{2}^{+}, I((t + T)/T) has been plotted.

 Fig. 4 The field intensity is varied such that the ionization yield (rate) of different systems remains similar to that of H_{2}^{+} subject to a 50cycle pulse of wavelength λ = 800 nm (ω = 0.0569 a.u.) and intensity I = 2.02 × 10^{14} W cm^{−2}, with a sinesquared pulse envelope. For the isotopologues of H_{2}^{+}, to view better how the ionization yield (rate) compares to that of H_{2}^{+} we have shifted the ionization by one cycle (i.e. I((t + T)/T) is plotted).  
To analyze the mechanism of ionization in more detail, in Fig. 5 we plot the timedependent nuclear density for the various isotopologues. The nuclear density behaves similarly in all cases, except the heaviest case, namely x_{2}^{+}. That is before the field reaches its maximum the system becomes slightly ionized hence the nuclear density slightly spreads. As the field reaches its maximum a small fragment of the nuclear density starts to split off and dissociate from Coulomb explosion due to an increase in the ionization while the rest of the nuclear density remains bounded and oscillates around the equilibrium position. As an appreciable amount of dissociating fragment reaches the internuclear separation for CREI, the ionization gets enhanced. In the case of x_{2}^{+}, however, the nuclear dynamics exhibits a more classical behavior, i.e. during the first half of the pulse, the heavy nuclei only spread slightly around the equilibrium. The stronger field compared to the other cases enables the system to ionize initially (before the nuclear density reaches the critical R). This initial ionization will be followed then by a Coulomb explosion toward the middle of the field, after which the nuclear density spreads further but it hardly splits. The remaining part of the electronic density undergoes enhanced ionization as the nuclear density lies in the range of the internuclear separation associated with CREI. In this case, the average internuclear distance almost coincides with the peak of the nuclear wave packet and therefore for systems with very large effective nuclear mass the standard approximation to describe CREI, namely the qs approximation, is expected to perform better compared to the lighter isotopes. The different nature of ionization in large nuclearmass systems can also be seen from I(R, t) depicted in Fig. 6. While for not too heavy isotopologues namely A_{2}^{+} and D_{2}^{+}, I(R, t) has a similar structure to H_{2}^{+}, I(R, t) corresponding to x_{2}^{+} manifests a substantially distinctive structure compared to the lighter isotopologues: the internuclear separation for which ionization gets enhanced shifts to smaller values. In general, I(R, t) shifts to smaller R for heavier isotopologues which could be related to the stronger (crudely)optimized field intensities used.^{5,6}

 Fig. 5 The contour plot of the timedependent nuclear density for H_{2}^{+} (lower left panel), A_{2}^{+} (upper left panel), D_{2}^{+} (lower right panel), x_{2}^{+} (upper right panel) subject to the (crudely)optimized field intensities given in Table 3.  

 Fig. 6 The contour plot of I(R, t) for H_{2}^{+} (lower left panel), A_{2}^{+} (upper left panel), D_{2}^{+} (lower right panel), x_{2}^{+} (upper right panel) subject to the (crudely)optimized field intensities given in Table 3.  
Now, we investigate the performance of the various approximations introduced in Section II B for differentmass systems. In Fig. 7 the ionization probabilities calculated from propagating the electron on the exact, adiabatic(adiab), approximate(app), exactHartree(exH), normalized Hartree(nexH), quasistatic(qs) and selfconsistent Hartree (SCH) for various symmetric isotopologues of H_{2}^{+} are plotted.

 Fig. 7 Ionization probabilities calculated from propagating the electron on the exact, adiabatic, approximate, Hartree like, normalized Hartree and quasistatic PES for different isomers of ionic hydrogen.  
For all cases presented in Fig. 7, propagation on the qs potential gives rise to an underestimation of the ionization probability. The average 〈R〉 entering into the qs potential is always considerably less than the internuclear separation of the dissociating fragment, and so does not access the CREI region that long during the duration of the pulse. The exactHartree (which compared to the qs approximation accounts for the spreading and splitting of the nuclear wave packet) follows the qs results until the middle of the propagation time then overtakes the qs ionization and finally gives a rather large overestimation of the final ionization for all cases. This overestimation can be improved significantly by using the normalized Hartree approximation as defined in eqn (15). In fact, the nexH approximation performs clearly better than other conventional approximations (qs, SCH, exH) for H_{2}^{+}, A_{2}^{+} and x_{2}^{+}. Even the approximate potential that depends on the more rigorous and complicated concept of conditional nuclear density rather than the nuclear density that appears in the Hartreelike approximation is considerably worse than nexH. The SCH performs very poorly with a negligible ionization probability for not too heavy isomers, and only in the large mass limit shows some improvement in predicting the ionization probability, yielding a similar performance to the other conventional approximations. This very poor performance of the SCH stems from the fact that the SCH is an uncorrelated approach and the splitting of the nuclear wavepacket in the CREI regime cannot be accounted for. Hence, for the cases where the CREI occurs due to the splitting of a dissociating fragment of the nuclear wave packet the SCH cannot capture the right physics. In the large mass limit, however, the CREI mechanism relies on the spreading of the nuclear wave packet rather than splitting (see Fig. 5), therefore the SCH approximation is not as poor. In general, the approximations do not perform better with increasing nuclear mass: dynamic correlation effects that are missing in all the potentials shown (apart from adiab and adiab + acc) depend more critically on the nuclear velocities relative to the electronic velocities, and adjusting the field to get similar ionization results in these being similar (See also Section IV A).
The two curves in Fig. 7 remaining to be discussed are adiab and adiab + acc, which are components arising from the marginal decomposition (eqn (12)). The propagation of the electron on the adiabatic potential alone (adiab) gives rise to the complete ionization of the system rather early for all isotopologues, while the propagation on the potential composed of the adiabatic plus acceleration (adiab + acc) term yields an ionization probability in a very good agreement with the exact results with the velocity term adding only a small correction, except in the case of x_{2}^{+}, i.e. the large mass limit. In the following section by studying the structures of the terms of the marginal decomposition we try to shed some light on the reason behind the performance of these potentials.
C. The structure of the dynamic electron–nuclear terms in the eTDPES: marginal decomposition
We concluded the previous section by briefly discussing the electron dynamics on different terms/combinationsofterms of the marginal decomposition of the eTDPES (eqn (12)). In order to better understand the outcomes, in this section we study structures of different terms of marginal decompositions of the eTDPES for the two radically different cases of H_{2}^{+} and x_{2}^{+}. We refer the readers to ref. 1 for a discussion on the components of the conditional decomposition of the eTDPES. For the sake of simplifying the discussion, here, we divide the electronic coordinate into two regions: the “innerregion” that refers to the region with z < 5 a.u. and the “outerregion” that describes the rest of the axis for which z > 5 a.u.
1. H_{2}^{+} case.
In Fig. 8, we present the terms/combinationsofterms of marginal decomposition of the eTDPES for the case of H_{2}^{+} at four different snapshots in time.

 Fig. 8 The exact eTDPES and its component from the marginal decompositions at various snapshots in time of H_{2}^{+}.  
Adiabatic term.
The adiabatic term (the first term in eqn (12)) is the main constituent of the eTDPES when it is decomposed according to eqn (12). In this case, as it is seen in Fig. 8 (the upper left panel), it initially follows the exact eTDPES in the innerregion and to some extent in the outerregion while it shows a different behavior in the asymptotic region (deep in the outerregion). However, as initially there is a negligibly small amount of electronic density far from the innerregion, this deviation does not influence the dynamics significantly. As the field intensity increases, the adiabatic potential starts to deviate from the exact potential, both in the innerregion and in the outerregion. It can be seen in Fig. 8 (the lower left panel) that the shape of the adiabatic potential in the innerregion (especially the upfield part) and its (average) slope in the outerregion differ significantly from the exact eTDPES. The (average) slope of the exact potential follows the slope of the field in the outerregion as expected. This feature together with the upfield part of the exact potential are mainly responsible for controlling the ionization from the upfield direction. Indeed the overestimation of ionization probability corresponding to the propagation on the adiabatic potential discussed in the previous section (see the lower left panel of Fig. 7), which starts around the 20th optical cycle, is associated with the lack of the asymptotic slope and the error in the shape of the upfield well. In particular, the wrong asymptotic behavior of the adiabatic potential in the outer region allows for ionization from both upfield and downfield sides in each half cycle, resulting in a huge overestimation of ionization probability due to the addition of the upfield contribution to the ionization. However, an important feature of the exact eTDPES is captured in the adiabatic potential: the development of the four wells representing the branching of the nuclear wavefunction in the innerregion of the potential which is associated with the correlation between the electronic and nuclear motions. Towards the end of the dynamics, where the field intensity is small again the adiabatic potential follows the exact eTDPES closely as can be seen in Fig. 8 (lowerpanel, right).
Velocity term.
The velocity term is the second term in eqn (12) and its overall contribution to the eTDPES, in this case, is small especially in the inner region. As can be seen in Fig. 8, it slightly corrects the adiabatic potential in the upfield (inner)region as well as the outerregion but not enough to capture the essential features appearing in the exact eTDPES.
Acceleration term.
The acceleration term is the last term in eqn (12) that is initially very small in the innerregion but as the ionization sets in, the addition of this term to the adiabatic potential significantly improves the shape of the potential, particularly when the field approaches the peakintensity as evident in Fig. 8 (upperright and lowerleft panels) in the innerregion as well as the outerregion (see the “adiab + acc”). The latter is due to the much better (average) asymptotic behavior of this term in the outerregion compared to the adiabatic term. On the other hand, it also improves the upfield/downfield well in the innerregion when added to the adiabatic potential. As a result propagating on the combination of adiabatic and acceleration terms leads to an ionization probability in a good agreement with the exact results as it is shown in Fig. 7 (the lower left panel).
As the structures of the adiabatic, velocity, and acceleration terms in the case of A_{2}^{+} and D_{2}^{+} are very similar to H_{2}^{+}, we do not address them here.
2. x_{2}^{+} case.
In Fig. 9, the terms/combinationsofterms of marginal decomposition of the eTDPES for the case of x_{2}^{+} are plotted at four different times.

 Fig. 9 The exact eTDPES and its component from the marginal decompositions at various snapshots in time of x_{2}^{+}.  
Adiabatic term.
The adiabatic term behaves similar to the case of H_{2}^{+}, i.e. initially and finally it agrees well with the exact eTDPES while when the field intensity is large it fails to follow the shape of the exact potential. Again, this failure, in particular the wrong average slope of the adiabatic potential in the asymptotic region, results in a huge overestimation of the ionization probability.
Velocity term.
The velocity term plays a crucial role in this case as is seen in Fig. 9 (lowerpanel, left). In particular, as the field approaches its maximum intensity, it exhibits more structure in the inner region and contributes more significantly to the overall shape of the eTDPES. Specifically, it exhibits a valley in the center that is completely absent in the case of H_{2}^{+} and lowers the interatomic barrier when added to the “adiab + acc” potential for most of the times between t ≈ 28 T and t ≈ 40 T in which most of the ionization occurs. The more dominant role of the velocity term in this case could be attributed to the increase of ponderomotive (wiggling) motion of the electron driven by a stronger external laser field.
Acceleration term.
This term has a correct asymptotic behavior, similar to the case of H_{2}^{+} but leads to a wrong estimate for the innerbarrier (see Fig. 9, the lowerleft panel) when added to the adiabatic term. The wrong estimation of the innerbarrier is exaggerated after the field reaches its maximum intensity (especially due to the higher interatomic barrier around zero between t ≈ 28 T and t ≈ 40 T), leading to an inaccurate ionization probability after t ≈ 30 T, as seen in Fig. 7 (the upper right panel). In ref. 15 the importance of the innerbarrier (peaks) and upfield well (steps) to achieve the correct electronlocalization has been shown.
IV. Quasiclassical analysis
The equations for the electronic wavefunction and conditional nuclear wavefunction cannot be solved exactly for systems of more than a few degrees of freedom, just as solving the full molecular Schrödinger equation exactly for those systems is not possible. In many cases a quasiclassical treatment of the nuclear dynamics should be sensible; by quasiclassical, we mean an ensemble of classical trajectories, weighted according to the initial distribution, is evolved for the nuclei, rather than a single trajectory. This would allow the possibility to capture spreading and branching of the nuclear distribution. Such a procedure within the exact factorization in its reverse flavor involves taking the classical limit of the conditional nuclear wavefunction which does not satisfy an equation of the Schrödinger form. Therefore, it differs from quasiclassical treatments of usual Schrödinger equations that have also been discussed in various forms for the marginal nuclear wavefunction within the “direct factorization” framework.^{31–36} Here we begin taking the classical limit of the conditional nuclear wavefunction by representing it in polar form: 
χ_{z}(R, t) = A_{z}(R, t)e^{iSz(R, t)/ħ}  (22) 
where A_{z}(R, t) and S_{z}(R, t) are both real functions. Inserting this into the equation of motion for χ_{z}(R, t), eqn (4), and sorting the terms in orders of ħ, we find O(ħ^{0}): 
 (23) 
where we define V(z, R, t) = W_{nn}(R) + W_{en}(z, R) + V^{laser}(z, R, t), and _{e}(z, t) = −iħ∂Φ(z, t)/∂z (see shortly for more on this term). Similarly, keeping only O(ħ^{0}) terms in eqn (7) for the eTDPES, we find 
 (24) 
The terms on the righthandside correspond to classically evaluating and ε^{gd}_{e}(z, t) respectively. Notice that the electron–nuclear coupling operator U^{coup}_{en} has a classical counterpart, as it contributes already at the zerothorder to ħ with the term in both eqn (23) and (24).
Eqn (23) would be a standard Hamilton–Jacobi equation of the form (where H(q,p,t) is the Hamiltonian function), for the action S_{z}(R, t) of two particles, one of mass μ_{n} and the other of mass μ_{e} in a potential V(z, R, t) − ε_{e}(z, t), if the secondlast term on the lefthandside was not present. It is perhaps not surprising that we do not retrieve a standard Hamilton–Jacobi equation, given that the equation for χ_{z} is not a TDSE. Although an ħ multiplies this term, it does in fact contribute in the classical limit, as will be discussed shortly. Still, classical Newtonlike equations can be derived from eqn (23) by defining the velocity fields,

 (25) 
Then, taking ∂/∂
R of
eqn (23) yields

 (26) 
and

 (27) 
where d/d
t = ∂/∂
t +
u^{e}_{z}∂/∂
z +
u^{n}_{z}∂/∂
R is the timederivative in the Lagrangian frame defined by the velocities.
A. Quasiclassical analysis of the terms in the eTDPES
Whether the equations above, together with the solution of the TDSE eqn (3), could form the basis of a mixed quantumclassical method remains for future work. Here, instead we consider how a quasiclassical analysis of the entire coupled electron–nuclear system can shed light on the structure and nuclearmassdependence of the terms in the eTDPES. Many aspects of electron dynamics in strong fields can be treated classically, especially when tunneling and quantization are not driving the primary physics. Indeed, ref. 37 showed that classical trajectory calculations reproduce the essential features of CREI for both cases of fixed and moving nuclei.
To this end, we first evaluate _{e}(z, t) in eqn (23) to its lowestorder in ħ. Inserting Φ(z, t) = a(z, t)exp(is(z, t))/ħ, analogous to eqn (22), into the TDSE eqn (3) for the marginal wavefunction Φ(z, t), and keeping only the term that is lowestorder in ħ, gives

 (28) 
This is the Hamilton–Jacobi equation for the action
s(
z,
t) of a classical electron evolving in Hamiltonian
_{e} +
ε_{e}(
z,
t). Here
is the electronic kinetic energy associated with the Hamiltonian in the TDSE for the electronic wavefunction,
eqn (3). It is important to note that this is not the same as the true electronic kinetic energy (see eqn (69) in
ref. 17). The term in
eqn (23), −
iħ∂
Φ(
z,
t)/∂
z, thus becomes
in the classical limit.
Next, we note that eqn (23) for S_{z} is consistent with the corresponding classical limits taken for the electronic wavefunction and the full molecular wavefunction: writing χ_{z}(R, t) = Ψ(z, R, t)/Φ(z, t), with Ψ(z, R, t) = A(z, R, t)exp(iS(z, R, t)) and Φ(z, t) = a(z, t)exp(is(z, t)), then

S_{z}(R, t) = S(z, R, t) − s(z, t)  (29) 
where, in the classical limit,
s(
z,
t) satisfies
eqn (28) and
S(
z,
R,
t), satisfies

 (30) 
It is straightforward to check by substitution that
eqn (23), (28)–(30) are consistent with each other.
Eqn (28) and (30) are both standard Hamilton–Jacobi equations, so, for the classical trajectories that they describe, we readily identify the following:
the nuclear momentum;
the electronic momentum;
the classical energy; and, as in the earlier discussion,
∂s/∂z = _{e}, 
and
∂s/∂t = _{e} + ε_{e}. 
Now we imagine an ensemble of classical particles with the initial position and momenta distributed according to the initial wavefunction. We consider evaluating the various contributions to the classical eTDPES of
eqn (24) from these classical trajectories. The integral over
R in that equation then becomes a sum over classical trajectories which have arrived at
z at time
t, but with differing
R,
i.e. one sums over all trajectories that have reached
z at time
t.
For the first term in eqn (24), ,

 (31) 
where we noted that ∂
S_{z}(
R,
t)/∂
R = ∂
S(
z,
R,
t)/∂
R (from
eqn (29)), and replaced the integral with a sum over all
N^{traj}_{z} trajectories
I_{z} that reach
z(
t) =
z at time
t;
Ṙ_{Iz} is the nuclear velocity along the
I_{z}th trajectory. First, let us see what this says about the overall structure of this term. During the dynamics, we observed that part of the nuclear density oscillates around equilibrium, moving slower than the dissociating part of the density. This means that the trajectories for small
z are mostly associated with more slowlymoving nuclear dynamics near equilibrium, while those with larger
z are in the process of ionizing, and so are associated with faster nuclear speeds due to the net ionic charge they have left behind. Hence, for small
z, there are more trajectories with smaller
Ṙ's than for large
z. This quasiclassical analysis suggests that the kinetic term in
eqn (31) rises as
z gets larger which is indeed what we see in
Fig. 10. Consider first the black curve, H
_{2}^{+}. At early times, the middle region of
, where the bulk of the electron density is, is rather flat and takes its smallest value: initially, the conditional nuclear probability has very little
zdependence in the region where the electron density is appreciable,
i.e., at small times the classical positions of nuclear trajectories for
z in the region of appreciable electronic density are identical, reflecting the broadness of the electron distribution relative to the nuclei. The small nuclear kinetic energy at the early times near equilibrium is essentially from zeropoint motion. As the electrons begin to gain energy from the field and move out from the tails of the electronic distribution,
grows accordingly since these tails are associated with trajectories where the nuclei are beginning to move apart under Coulomb repulsion. The small middle flat region of
gets narrower as ionization takes away more of the electron density. As the final distribution sets in,
, one can identify three regions, reflecting the electronic distribution and the associated nuclear kinetic energies: an inner region associated with the density remaining near equilibrium, rising to the intermediate region where the shoulders of the electron density lie associated with the dissociated nuclear wavepacket, and then rising further out to the tails of the electron density which continue to oscillate in the field.

 Fig. 10 The kinetic term, , at various snapshots in time for H_{2}^{+} and its isotopologues.  
Regarding the μ_{n}dependence: it might appear from eqn (31) that grows as the mass μ_{n} increases, in the cases where the field is adjusted so that the average internuclear distance and speed remain the same. This is in fact not the case; in fact, the term remains about the same in size, as the nuclear mass increases as is evident in Fig. 10. This is because, for larger μ_{n}, the nuclear density tends to split less (see Fig. 6), i.e. for larger μ_{n}, a larger proportion of the nuclear density moves out to larger separations rather than remaining near the origin with small Ṙ. Yet about the same 〈Ṙ〉 is maintained by design by the different field strengths, which means that the fastest trajectories in the distribution are slower for larger masses μ_{n} than for smaller μ_{n}. This would contribute to a smaller rise of as one moves to larger z, as mass increases, compensating the larger μ_{n}. In Fig. 10 we plot for the different isotopologues at various snapshots in time. We see that the overall size of this term does not vary much with mass except for x_{2}^{+}. The exception is for the large mass x_{2}^{+} where there is comparatively weak zdependence of . This is because the nuclear distribution moves largely as a whole, showing less difference between nuclear speeds at different parts of the distribution.
The corresponding analysis for the other three terms of eqn (24) is unfortunately not as straightforward, because the evolution of the relevant classical trajectories themselves depends on the potential ε_{e}(z, t), which is the very object we wish to analyse! For example, the second term, K^{cond}_{e}(z, t), requires but _{e} is the kinetic energy of an electron in the potential ε_{e}. Future work along these lines will likely require further, possibly iterative, approximations to the potential, with the dual goal of analysing the structure of the terms and developing mixed quantumclassical schemes.
V. Conclusions
The exact potential driving the electron dynamics within the exact factorization of the full electron–nuclear wavefunction formally exactly accounts for the coupling to the nuclear subsystem as well as coupling to external fields. In order to develop adequate approximations to treat electronic nonadiabatic processes, it is important to understand the structure of the exact potential and pinpoint and analyse its important features in various situations. In a recent work^{1} we had presented the exact electronic potential driving the electron dynamics of H_{2}^{+} undergoing CREI. In this work, with the aim of developing a practical scheme to study the CREI for larger systems, we have presented a detailed investigation of the exact potential driving the electron dynamics in isotopologues of H_{2}^{+} undergoing CREI and studied the dependence of the correlated electron–nuclear dynamics on the nuclear mass.
We have shown that the elaborate concept of timeresolved, Rresolved ionization probability, I(R, t), provides an extremely useful tool, indicating the internuclear separations associated with ionization. The concept is analogous to the ionization probability calculated within the quasistatic picture (with clamped nuclei) but can be used unambiguously for the fully dynamic case. For the laser parameters and different isotopologues of H_{2}^{+} discussed in this work, I(R, t) exhibits only one peak along the Raxis. This is in agreement with most of the experimental findings^{38,39} and is in contrast to the predictions of CREI based on the standard clampednuclei quasistatic picture. We suggest that I(R, t) would be a very useful tool in future studies of CREI, for example to resolve the laser parameters for which the doublepeak structure could in fact be observed as in ref. 40.
We found that for fixed laser parameters, the ionization yields rapidly decrease as a function of the mass of the isotopologue because less of the nuclear density makes it to the CREI region during the time the laser is on. For all the isotopologues, I(R, t) indicated that the ionization is nevertheless dominated by the fraction of electrons associated with the nuclear density in the CREI region defined qualitatively by the original quasistatic argument. This implies that treating the nuclei as classical point particles will not work; one needs to account for their distribution, which, away from the largemass limit, displays a branched structure with part of the distribution oscillating near equilibrium separation while part of it, associated with the CREI electrons, dissociates.
Still, going beyond the quasistatic point nuclei picture and accounting for the nuclear distribution in a Coulombic way (as in Hartreetype approximations) is not adequate in capturing the dynamics accurately. The importance of going beyond a Coulombic description of electron–nuclear correlation is evident from our studies of how different approximate electronic potentials perform in describing the CREI for isotopologues of H_{2}^{+}. We have shown that, for the laser parameters used in this work, one must go beyond any purely Coulombic treatment of electron–nuclear correlation, and include truly dynamic aspects of the nuclear distribution and its coupling to the electronic system to get a good prediction of the ionization. In determining errors from conventional approximations and deviations of approximate potentials from the exact potential, nuclear velocities in the wavepacket are more important than the nucleartoelectronic mass ratios.
There are many different approaches to developing approximate methods based on the exact factorization. For example, one may consider the relative importance of different components of the exact potential when decomposed in terms of the conditional wavefunction as in ref. 1 (there, it was found that generally all terms were important). One may consider also approximations based on the marginal decomposition presented here, where we found that the “adiab + acc” component of the marginal decomposition describes the electronic dynamics accurately in all cases we studied here with the exception of a fictitious isotopologue with an effective nuclear mass 10 times larger than H_{2}^{+}. Further investigation of this potential may lead to development of an adequate approximation for practical purposes capable of describing the ionization dynamics accurately. Another approach is to develop quasiclassical approximations, and here we have sketched out a semiclassical derivation of the electronic and conditional nuclear equations of the exact factorization in its reverse form and analysed the structure of the nuclear kinetic term of the exact potential semiclassically.
This work has highlighted the effect of the complex interplay between electronic and nuclear dynamics in strong field enhanced ionization processes by demonstrating the large differences in the conventional potentials and the dynamics they cause with the exact potential driving the electron for systems of varying nuclear mass. The explorations of the details of the potential and approximation methods lay the groundwork for future development of accurate methods for coupled electron–ion dynamics in nonperturbative fields. Whether the dynamic electron–nuclear effects play a crucial a role in polyatomic molecules, and the scaling of these effects with respect to the number of electrons and the number of nuclei, is also an important avenue for future research.
Acknowledgements
We acknowledge support from the European Research Council (ERC2015AdG694097), Grupos Consolidados (IT57813), and the European Union's Horizon 2020 Research and Innovation programme under grant agreement no. 676580. A. K. and A. A. acknowledge funding from the European Union's Horizon 2020 research and innovation programme under the Marie SklodowskaCurie grant agreement no. 704218 and 702406, respectively. N. T. M. thanks the National Science Foundation, grant CHE1566197, for support. Open Access funding provided by the Max Planck Society.
References
 E. Khosravi, A. Abedi and N. T. Maitra, Phys. Rev. Lett., 2015, 115, 263002 CrossRef PubMed.
 T. Zuo, S. Chelkowski and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1993, 48, 3837 CrossRef CAS.
 T. Zuo and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1995, 52, R2511 CrossRef CAS.
 T. Seideman, M. Y. Ivanov and P. B. Corkum, Phys. Rev. Lett., 1995, 75, 2819 CrossRef CAS PubMed.
 S. Chelkowski, C. Foisy and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1998, 57, 1176 CrossRef CAS.
 S. Chelkowski and A. Bandrauk, J. Phys. B: At., Mol. Opt. Phys., 1995, 28, L723 CrossRef CAS.
 S. Chelkowski, A. Conjusteau, T. Zuo and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 54, 3235 CrossRef CAS.
 H. Yu, T. Zuo and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 54, 3290 CrossRef CAS.

A. D. Bandrauk and F. Légaré, Progress in Ultrafast Intense Laser Science VIII, Springer, 2012, pp. 29–46 Search PubMed.
 N. Takemoto and A. Becker, Phys. Rev. Lett., 2010, 105, 203004 CrossRef PubMed.
 N. Takemoto and A. Becker, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 84, 023401 CrossRef.
 C. Beylerian, S. Saugout and C. Cornaggia, J. Phys. B: At., Mol. Opt. Phys., 2006, 39, L105 CrossRef CAS.
 I. Bocharova, R. Karimi, E. F. Penka, J.P. Brichta, P. Lassonde, X. Fu, J.C. Kieffer, A. D. Bandrauk, I. Litvinyuk and J. Sanderson,
et al.
, Phys. Rev. Lett., 2011, 107, 063201 CrossRef PubMed.
 F. Légaré, I. V. Litvinyuk, P. W. Dooley, F. Quéré, A. D. Bandrauk, D. M. Villeneuve and P. B. Corkum, Phys. Rev. Lett., 2003, 91, 093002 CrossRef PubMed.
 Y. Suzuki, A. Abedi, N. T. Maitra, K. Yamashita and E. K. U. Gross, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 89, 040501 CrossRef.
 A. Abedi, N. T. Maitra and E. K. U. Gross, Phys. Rev. Lett., 2010, 105, 123002 CrossRef PubMed.
 A. Abedi, N. T. Maitra and E. K. U. Gross, J. Chem. Phys., 2012, 137, 22A530 CrossRef PubMed.
 A. Abedi, N. T. Maitra and E. K. U. Gross, J. Chem. Phys., 2013, 139, 087102 CrossRef PubMed.
 A. Abedi, F. Agostini, Y. Suzuki and E. K. U. Gross, Phys. Rev. Lett., 2013, 110, 263001 CrossRef PubMed.
 R. Requist and E. K. U. Gross, Phys. Rev. Lett., 2016, 117, 193001 CrossRef PubMed.
 H. Yu, T. Zuo and A. D. Bandrauk, J. Phys. B: At., Mol. Opt. Phys., 1998, 31, 1533 CrossRef CAS.
 M. Uhlmann, F. Großmann, T. Kunert and R. Schmidt, Phys. Lett. A, 2007, 364, 417 CrossRef CAS.
 J. C. Tully, Faraday Discuss., 1998, 110, 407 RSC.
 The Hartree potential (14) reduces to the qs expression (13) in the limit of very localized wavepackets centered around R.
 J. L. Krause, K. J. Schafer and K. C. Kulander, Phys. Rev. A: At., Mol., Opt. Phys., 1992, 45, 4998 CrossRef CAS.
 J. Muga, J. Palao, B. Navarro and I. Egusquiza, Phys. Rep., 2004, 395, 357 CrossRef CAS.
 R. Santra and L. S. Cederbaum, Phys. Rep., 2002, 368, 1 CrossRef CAS.
 J. Javanainen, J. H. Eberly and Q. Su, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3430 CrossRef CAS.
 T. Kreibich, M. Lein, V. Engel and E. K. U. Gross, Phys. Rev. Lett., 2001, 87, 103901 CrossRef CAS PubMed.
 K. C. Kulander, F. H. Mies and K. J. Schafer, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 53, 2562 CrossRef CAS.
 A. Abedi, F. Agostini and E. K. U. Gross, Europhys. Lett., 2014, 106, 33001 CrossRef.
 S. K. Min, F. Agostini and E. K. U. Gross, Phys. Rev. Lett., 2015, 115, 073001 CrossRef PubMed.
 F. Agostini, S. K. Min, A. Abedi and E. K. Gross, J. Chem. Theory Comput., 2016, 12, 2127 CrossRef CAS PubMed.
 Y. Suzuki, A. Abedi, N. Maitra and E. K. U. Gross, Phys. Chem. Chem. Phys., 2015, 17, 29271 RSC.
 F. Agostini, A. Abedi, Y. Suzuki, S. K. Min, N. T. Maitra and E. K. U. Gross, J. Chem. Phys., 2015, 142, 084303 CrossRef PubMed.
 Y. Suzuki and K. Watanabe, Phys. Rev. A, 2016, 94, 032517 CrossRef.
 D. M. Villeneuve, M. Y. Ivanov and P. B. Corkum, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 54, 736 CrossRef CAS.
 T. Ergler, A. Rudenko, B. Feuerstein, K. Zrost, C. D. Schröter, R. Moshammer and J. Ullrich, Phys. Rev. Lett., 2005, 95, 093001 CrossRef PubMed.
 I. BenItzhak, P. Q. Wang, A. M. Sayler, K. D. Carnes, M. Leonard, B. D. Esry, A. S. Alnaser, B. Ulrich, X. M. Tong, I. V. Litvinyuk, C. M. Maharjan, P. Ranitovic, T. Osipov, S. Ghimire, Z. Chang and C. L. Cocke, Phys. Rev. A: At., Mol., Opt. Phys., 2008, 78, 063419 CrossRef.
 H. Xu, F. He, D. Kielpinski, R. Sang and I. Litvinyuk, Sci. Rep., 2015, 5, 13527 CrossRef CAS PubMed.

This journal is © the Owner Societies 2017 