Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Electronic non-adiabatic dynamics in enhanced ionization of isotopologues of hydrogen molecular ions from the exact factorization perspective

Elham Khosravi *a, Ali Abedi *a, Angel Rubio *ba and Neepa T. Maitra *c
aNano-Bio Spectroscopy Group and ETSF, Universidad del País Vasco, CFM CSIC-UPV/EHU, 20018 San Sebastián, Spain. E-mail:;
bMax Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany. E-mail:
cDepartment of Physics and Astronomy, Hunter College and the Graduate Center of the City University of New York, 695 Park Avenue, New York, New York 10065, USA. E-mail:

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 H2+ can have large contributions arising from dynamic electron–nuclear correlation, going beyond what any Coulombic-based 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 H2+ in order to investigate the nuclear-mass-dependence 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 quasi-static or Hartree approximation with the exact propagation. A quasiclassical analysis is presented to help analyse the structure of different non-Coulombic components of the potential driving the ionizing electron.

I. Introduction

The phenomenon of charge-resonance 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 quasi-static argument involving instantaneously frozen nuclei in the pioneering works of ref. 3–9 for which the time-dependent Schrödinger equation (TDSE) is solved for various clamped nuclear (cn) configurations image file: c6cp08539c-t1.tif, i.e.,
image file: c6cp08539c-t2.tif(1)
image file: c6cp08539c-t3.tif(2)
Here image file: c6cp08539c-t4.tif and image file: c6cp08539c-t5.tif are used to collectively denote the electronic and nuclear coordinates, [T with combining circumflex]e is the electronic kinetic energy operator, and Ŵee is the electron–electron repulsion. Furthermore, image file: c6cp08539c-t6.tif contains the electron–nuclear interaction, which, for a diatomic molecular ion with one electron, as will be considered here, has the form −Z1/|rR1| − Z2/|rR2|, i.e. a Coulombic double-well. 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 up-field atom Stark-shift 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 up-field level (LUMO) does not significantly tunnel back to the down-field 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 Stark-shifted LUMO level exceeds the top of both the inner and outer field-modified Coulombic barriers, one finds Rc = 4.07/Ip 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 sub-cycle 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 quasi-static 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 un-ionized. 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 double-well 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 electron-nuclear 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 mass-dependence, 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 time-dependent 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 one-dimensional model for the (a)symmetric isotopologues of H2+ 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 time-dependent 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 image file: c6cp08539c-t7.tif, or image file: c6cp08539c-t8.tif, where partial normalization conditions on the conditional wavefunctions image file: c6cp08539c-t9.tif and image file: c6cp08539c-t10.tif, respectively, render each factorization unique up to a gauge transformation. In the first product the equation for the nuclear wavefunction, image file: c6cp08539c-t11.tif, is a TDSE, while in the second product the equation for the electronic factor image file: c6cp08539c-t12.tif 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 many-body nuclear density and current density of the complete electron–nuclear system while in the latter case it leads to the exact many-body 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 image file: c6cp08539c-t13.tif.

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 image file: c6cp08539c-t14.tif by z and R respectively):

image file: c6cp08539c-t15.tif(3)
and that for the nuclear conditional wavefunction is:
[Ĥn(z, R, t) − εe(z, t)]χz(R, t) = [thin space (1/6-em)]tχz(R, t),(4)
with the nuclear Hamiltonian
image file: c6cp08539c-t16.tif(5)
where [T with combining circumflex]n is the nuclear kinetic energy operator, Ŵen(z, R) (Ŵnn(R)) is the electron–nuclear (nuclear–nuclear) interaction, and [V with combining circumflex]lasern(R, t) is the time-dependent external potential acting on the nuclei. Here, image file: c6cp08539c-t17.tif is the exact time-dependent vector potential
image file: c6cp08539c-t18.tif(6)
and εe(z, t) = 〈χz(t)|Ĥn[thin space (1/6-em)]t|χz(t)〉R is the exact time-dependent electronic potential for an electron (see the next section). In one dimensional models, image file: c6cp08539c-t19.tif 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 time-dependent potential energy surface

The exact time-dependent potential energy surface for electrons (e-TDPES) contains the effects of the coupling to moving quantum nuclei as well as the external laser field. It can be written as
image file: c6cp08539c-t20.tif(7)
which consists of: the approximate potential
εappe(z, t) = [V with combining circumflex]lasere(z, t) + 〈χz(R, t)|Ŵen(z, R) + Ŵnn(R)|χz(R, t)〉R,(8)
the nuclear-kinetic term
image file: c6cp08539c-t21.tif(9)
the gauge dependent part of the potential
εgde(z, t) = ħχz(R, t)| − itχz(R, t)〉R,
and finally the electronic-kinetic-like contribution
image file: c6cp08539c-t22.tif
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 H2+ 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 e-TDPES by exactly inverting the TDSE for the electronic wavefunction. The electronic wavefunction may be written in polar representation,

image file: c6cp08539c-t23.tif(10)
with ne(z, t) = |Φ(z, t)|2 and image file: c6cp08539c-t24.tifwhere
image file: c6cp08539c-t25.tif(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 one-dimensional model system. Here, e = [k with combining circumflex]∂/∂z). Inserting this form into eqn (3) and setting the time-dependent vector potential to zero give
image file: c6cp08539c-t26.tif(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 time-dependent density functional theory. The second term, denoted the “velocity term”, depends only on the electron velocity, namely on je(z, t)/ne(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 e-TDPES.

B. Approximate electronic potentials based on conventional approximations

The standard approximation for the electronic potential assumes the so-called quasistatic approximation (qs) that treats the nuclei as classical point particles with positions that are either considered fixed, [R with combining macron]0, as in the majority of studies of CREI,3,5–7,9,21 or move classically with classical trajectories [R with combining macron](t) that are often described by mixed quantum-classical algorithms such as Ehrenfest or surface-hopping 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|[R with combining macron](t)) = Ŵen(z, [R with combining macron](t)) + [V with combining circumflex]lasere(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 image file: c6cp08539c-t27.tif, image file: c6cp08539c-t28.tif, 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 z-independent delta-function at [R with combining macron](t), i.e. nz(R, t) ≈ δ(R(t), [R with combining macron](t)) (in this limit Ŵnn([R with combining macron](t)) becomes purely a time-dependent 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 approximation23

image file: c6cp08539c-t29.tif(14)
where n(R, t) is the nuclear density obtained from nuclear wavepacket dynamics.24 It can be easily seen that if the z-dependence in the conditional nuclear wavefunction is neglected, i.e., nz(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 [R with combining macron](t) = 〈R〉(t) is the average time-dependent internuclear separation obtained from the exact calculations; (ii) The Hartree approximation for which the nuclear density in (14) is replaced by the exact time-dependent nuclear density, indicated by εex-H; (iii) The normalized Hartree, defined as

image file: c6cp08539c-t30.tif(15)
indicated as εn-ex-H when the exact nuclear density, obtained from the complete electron–nuclear wavefunction Ψ(z, R), i.e.
image file: c6cp08539c-t31.tif(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 “self-consistent” Hartree approximation (SCH) in which the full wavefunction is approximated as an uncorrelated product of electronic wavefunction and nuclear wavefunction, image file: c6cp08539c-t32.tif, 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.,
image file: c6cp08539c-t33.tif(17)
where ne(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 mean-field 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 H2+ isotopologues

We utilize a popular one-dimensional model of the symmetric as well as asymmetric isotopologues of H2+ 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 “soft-Coulomb” interactions:28,29
image file: c6cp08539c-t34.tif(18)
where R and z are the internuclear distance and the electronic coordinate as measured from the nuclear center-of-mass, respectively. The nuclear effective mass is denoted as image file: c6cp08539c-t35.tif while image file: c6cp08539c-t36.tif is the electronic reduced mass with Mn = M1 + M2. The laser field, within the dipole approximation, is represented by
[V with combining circumflex]laser(R, z, t) = eE(t)(qezζR),(19)
where E(t) denotes the electric field amplitude image file: c6cp08539c-t37.tif is the reduced charge, and ζ = (M2M1)/Mn is the mass-asymmetry parameter. Such reduced-dimensional 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 real-space grid. The parameters are in atomic units. The time-step 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., H2+, A2+, D2+, x2+ 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 A2+ is another fictitious isotopologue with the same effective nuclear mass as HT+ (see Table 2 in which the nuclear effective mass values of H2+ and isotopologues are given).

Table 2 Nuclear effective mass corresponding to H2+ and its symmetric and antisymmetric isotopologues in atomic units. Here, x(X) refers to the 10(100) times heavier fictitious isotope of hydrogen and A2+ is another fictitious isotopologue with the same effective nuclear mass as HT+
H2+ HD+ HT+(A2+) Hx+ HX+ D2+ x2+
918.076 1223.742 1376.228 1669.229 1817.973 1834.533 9180.76

A. Comparison of ionization of isotopologues of H2+ with different effective nuclear mass for a fixed field

We apply a field of duration 50-cycles, wavelength λ = 800 nm (ω = 0.0569 a.u.) and intensity I = 2.02 × 1014 W cm−2, with a sine-squared pulse envelope (Fig. 1), to each of the H2+ isotopologues. First, we solve the electron–nuclear TDSE for Hamiltonian of eqn (18) numerically exactly, beginning in the initial ground-state of the molecule.
image file: c6cp08539c-f1.tif
Fig. 1 Laser field as a function of number of optical cycles t/T. The electric field amplitude is divided by the peak amplitude, image file: c6cp08539c-t65.tif.

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 H2+, HD+, HT+, Hx+, HX+, A2+, D2+ and x2+. Note that, the results are given in atomic units, e = me = h = 1, throughout the article, unless otherwise noted. The ionization probability is defined as,

image file: c6cp08539c-t38.tif(20)
with zI = 15 a.u. As can be seen in Fig. 2, the ionization and average internuclear separation appear to be practically independent of the nuclear mass-symmetry properties: they are identical for the asymmetric case of HT+ and the symmetric case of A2+ systems with the same effective nuclear mass μn and are very similar in the case of HX+ and D2+ 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 zR, while the nuclear mass-dependence 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.

image file: c6cp08539c-f2.tif
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 H2+, HD+, HT+, Hx+ and HX+ molecules on the left panels and H2+, A2+, D2+ and x2+ on the right panels.

It is also observed that as the nuclear effective-mass 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 time-resolved, R-resolved ionization probability defined as image file: c6cp08539c-t39.tif, with image file: c6cp08539c-t40.tif and zI = 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)|2Icp(R, t),(21)
where, Icp(R, t) follows the usual expression of the ionization probability but using the exact conditional electronic wavefunction ΦR(z, t), i.e., image file: c6cp08539c-t41.tif that is coupled to the exact nuclear dynamics. Therefore, I(R, t), which is the nuclear-density-weighted conditional ionization probability, is analogous to the ionization probability calculated for a given nuclear configuration R in a quasi-static picture.

To analyze the dynamics of different cases, in Fig. 3 we have plotted I(R, t) along with the time-dependent 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 Rc. The lower set of panels of Fig. 3 shows the time-dependent nuclear density in each case. For the case of x2+ 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 R-resolved, t-resolved ionization probability. Notice the different scale of I(R, t) for the case of x2+.

image file: c6cp08539c-f3.tif
Fig. 3 The time-resolved, R-resolved 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 H2+, D2+ and x2+.

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 Rc, 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 Rc region. Clearly, treating the nuclei as point particles will not work even for significant nuclear mass (except in the large-mass 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 Coulomb-based 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 H2+ subject to the same field show significantly different degrees of ionization, since the different-mass 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 different-mass systems, we adjust the field intensity such that the ionization probability(rate) remains close to that of H2+. 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 H2+. The (crudely-)optimized field intensity for different-mass systems is given in Table 3 while the other field parameters are kept unchanged.
Table 3 The (crudely-)optimized field intensity for different-mass systems in the unit of W cm−2
H2+ A2+ D2+ x2+
2.02 × 1014 2.14 × 1014 2.26 × 1014 3.1 × 1014

The ionization probability of the symmetric isotopologues subject to the (crudely-)optimized field intensities is depicted in Fig. 4. For isotopologues heavier than H2+ the ionization starts to set in about one optical cycle T later than the H2+ case. In order to be able to compare the ionization yield (rate) of these systems with H2+ visually better, in Fig. 4 we have shifted the time-dependent ionization probability by one optical cycle, i.e. for isotopologues heavier than H2+, I((t + T)/T) has been plotted.

image file: c6cp08539c-f4.tif
Fig. 4 The field intensity is varied such that the ionization yield (rate) of different systems remains similar to that of H2+ subject to a 50-cycle pulse of wavelength λ = 800 nm (ω = 0.0569 a.u.) and intensity I = 2.02 × 1014 W cm−2, with a sine-squared pulse envelope. For the isotopologues of H2+, to view better how the ionization yield (rate) compares to that of H2+ 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 time-dependent nuclear density for the various isotopologues. The nuclear density behaves similarly in all cases, except the heaviest case, namely x2+. 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 x2+, 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 nuclear-mass systems can also be seen from I(R, t) depicted in Fig. 6. While for not too heavy isotopologues namely A2+ and D2+, I(R, t) has a similar structure to H2+, I(R, t) corresponding to x2+ 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

image file: c6cp08539c-f5.tif
Fig. 5 The contour plot of the time-dependent nuclear density for H2+ (lower left panel), A2+ (upper left panel), D2+ (lower right panel), x2+ (upper right panel) subject to the (crudely-)optimized field intensities given in Table 3.

image file: c6cp08539c-f6.tif
Fig. 6 The contour plot of I(R, t) for H2+ (lower left panel), A2+ (upper left panel), D2+ (lower right panel), x2+ (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 different-mass systems. In Fig. 7 the ionization probabilities calculated from propagating the electron on the exact, adiabatic(adiab), approximate(app), exact-Hartree(ex-H), normalized Hartree(n-ex-H), quasi-static(qs) and self-consistent Hartree (SCH) for various symmetric isotopologues of H2+ are plotted.

image file: c6cp08539c-f7.tif
Fig. 7 Ionization probabilities calculated from propagating the electron on the exact, adiabatic, approximate, Hartree like, normalized Hartree and quasi-static 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 exact-Hartree (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 n-ex-H approximation performs clearly better than other conventional approximations (qs, SCH, ex-H) for H2+, A2+ and x2+. 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 Hartree-like approximation is considerably worse than n-ex-H. 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 wave-packet 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 x2+, 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 e-TDPES: marginal decomposition

We concluded the previous section by briefly discussing the electron dynamics on different terms/combinations-of-terms of the marginal decomposition of the e-TDPES (eqn (12)). In order to better understand the outcomes, in this section we study structures of different terms of marginal decompositions of the e-TDPES for the two radically different cases of H2+ and x2+. We refer the readers to ref. 1 for a discussion on the components of the conditional decomposition of the e-TDPES. For the sake of simplifying the discussion, here, we divide the electronic coordinate into two regions: the “inner-region” that refers to the region with |z| < 5 a.u. and the “outer-region” that describes the rest of the axis for which |z| > 5 a.u.
1. H2+ case. In Fig. 8, we present the terms/combinations-of-terms of marginal decomposition of the e-TDPES for the case of H2+ at four different snapshots in time.
image file: c6cp08539c-f8.tif
Fig. 8 The exact e-TDPES and its component from the marginal decompositions at various snapshots in time of H2+.

Adiabatic term. The adiabatic term (the first term in eqn (12)) is the main constituent of the e-TDPES 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 e-TDPES in the inner-region and to some extent in the outer-region while it shows a different behavior in the asymptotic region (deep in the outer-region). However, as initially there is a negligibly small amount of electronic density far from the inner-region, 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 inner-region and in the outer-region. It can be seen in Fig. 8 (the lower left panel) that the shape of the adiabatic potential in the inner-region (especially the up-field part) and its (average) slope in the outer-region differ significantly from the exact e-TDPES. The (average) slope of the exact potential follows the slope of the field in the outer-region as expected. This feature together with the up-field part of the exact potential are mainly responsible for controlling the ionization from the up-field 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 up-field well. In particular, the wrong asymptotic behavior of the adiabatic potential in the outer region allows for ionization from both up-field and down-field sides in each half cycle, resulting in a huge overestimation of ionization probability due to the addition of the up-field contribution to the ionization. However, an important feature of the exact e-TDPES is captured in the adiabatic potential: the development of the four wells representing the branching of the nuclear wavefunction in the inner-region 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 e-TDPES closely as can be seen in Fig. 8 (lower-panel, right).
Velocity term. The velocity term is the second term in eqn (12) and its overall contribution to the e-TDPES, 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 up-field (inner-)region as well as the outer-region but not enough to capture the essential features appearing in the exact e-TDPES.
Acceleration term. The acceleration term is the last term in eqn (12) that is initially very small in the inner-region 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 peak-intensity as evident in Fig. 8 (upper-right and lower-left panels) in the inner-region as well as the outer-region (see the “adiab + acc”). The latter is due to the much better (average) asymptotic behavior of this term in the outer-region compared to the adiabatic term. On the other hand, it also improves the up-field/down-field well in the inner-region 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 A2+ and D2+ are very similar to H2+, we do not address them here.

2. x2+ case. In Fig. 9, the terms/combinations-of-terms of marginal decomposition of the e-TDPES for the case of x2+ are plotted at four different times.
image file: c6cp08539c-f9.tif
Fig. 9 The exact e-TDPES and its component from the marginal decompositions at various snapshots in time of x2+.

Adiabatic term. The adiabatic term behaves similar to the case of H2+, i.e. initially and finally it agrees well with the exact e-TDPES 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 (lower-panel, 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 e-TDPES. Specifically, it exhibits a valley in the center that is completely absent in the case of H2+ 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 H2+ but leads to a wrong estimate for the inner-barrier (see Fig. 9, the lower-left panel) when added to the adiabatic term. The wrong estimation of the inner-barrier 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 inner-barrier (peaks) and up-field well (steps) to achieve the correct electron-localization has been shown.

IV. Quasi-classical 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) = Az(R, t)eiSz(R, t)/ħ(22)
where Az(R, t) and Sz(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):
image file: c6cp08539c-t42.tif(23)
where we define V(z, R, t) = Wnn(R) + Wen(z, R) + Vlaser(z, R, t), and [p with combining tilde]e(z, t) = −Φ(z, t)/∂z (see shortly for more on this term). Similarly, keeping only O(ħ0) terms in eqn (7) for the e-TDPES, we find
image file: c6cp08539c-t43.tif(24)
The terms on the right-hand-side correspond to classically evaluating image file: c6cp08539c-t44.tif and εgde(z, t) respectively. Notice that the electron–nuclear coupling operator Ucoupen has a classical counterpart, as it contributes already at the zeroth-order to ħ with the term image file: c6cp08539c-t45.tif in both eqn (23) and (24).

Eqn (23) would be a standard Hamilton–Jacobi equation of the form image file: c6cp08539c-t46.tif (where H(q,p,t) is the Hamiltonian function), for the action Sz(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 second-last term on the left-hand-side 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 Newton-like equations can be derived from eqn (23) by defining the velocity fields,

image file: c6cp08539c-t47.tif(25)
Then, taking ∂/∂R of eqn (23) yields
image file: c6cp08539c-t48.tif(26)
image file: c6cp08539c-t49.tif(27)
where d/dt = ∂/∂t + uez∂/∂z + unz∂/∂R is the time-derivative in the Lagrangian frame defined by the velocities.

A. Quasiclassical analysis of the terms in the e-TDPES

Whether the equations above, together with the solution of the TDSE eqn (3), could form the basis of a mixed quantum-classical 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 nuclear-mass-dependence of the terms in the e-TDPES. 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 [p with combining tilde]e(z, t) in eqn (23) to its lowest-order 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 lowest-order in ħ, gives

image file: c6cp08539c-t50.tif(28)
This is the Hamilton–Jacobi equation for the action s(z, t) of a classical electron evolving in Hamiltonian [T with combining tilde]e + εe(z, t). Here image file: c6cp08539c-t51.tif 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), −Φ(z, t)/∂z, thus becomes image file: c6cp08539c-t52.tif in the classical limit.

Next, we note that eqn (23) for Sz 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)[thin space (1/6-em)]exp(is(z, t)), then

Sz(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
image file: c6cp08539c-t53.tif(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:
S/∂R = Pn = μn,
the nuclear momentum;
S/∂z = pe,
the electronic momentum;
S/∂t = E(t),
the classical energy; and, as in the earlier discussion,
s/∂z = [p with combining tilde]e,
s/∂t = [T with combining tilde]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 e-TDPES 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), image file: c6cp08539c-t54.tif,

image file: c6cp08539c-t55.tif(31)
where we noted that ∂Sz(R, t)/∂R = ∂S(z, R, t)/∂R (from eqn (29)), and replaced the integral with a sum over all Ntrajz trajectories Iz that reach z(t) = z at time t; Iz is the nuclear velocity along the Iz-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 slowly-moving 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 quasi-classical 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, H2+. At early times, the middle region of image file: c6cp08539c-t56.tif, where the bulk of the electron density is, is rather flat and takes its smallest value: initially, the conditional nuclear probability has very little z-dependence 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 zero-point motion. As the electrons begin to gain energy from the field and move out from the tails of the electronic distribution, image file: c6cp08539c-t57.tif 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 image file: c6cp08539c-t58.tif gets narrower as ionization takes away more of the electron density. As the final distribution sets in, image file: c6cp08539c-t59.tif, 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.

image file: c6cp08539c-f10.tif
Fig. 10 The kinetic term, image file: c6cp08539c-t66.tif, at various snapshots in time for H2+ and its isotopologues.

Regarding the μn-dependence: it might appear from eqn (31) that image file: c6cp08539c-t60.tif 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 image file: c6cp08539c-t61.tif as one moves to larger z, as mass increases, compensating the larger μn. In Fig. 10 we plot image file: c6cp08539c-t62.tif 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 x2+. The exception is for the large mass x2+ where there is comparatively weak z-dependence of image file: c6cp08539c-t63.tif. 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, Kconde(z, t), requires image file: c6cp08539c-t64.tif but [T with combining tilde]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 quantum-classical 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 non-adiabatic 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 work1 we had presented the exact electronic potential driving the electron dynamics of H2+ 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 H2+ undergoing CREI and studied the dependence of the correlated electron–nuclear dynamics on the nuclear mass.

We have shown that the elaborate concept of time-resolved, R-resolved 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 H2+ discussed in this work, I(R, t) exhibits only one peak along the R-axis. This is in agreement with most of the experimental findings38,39 and is in contrast to the predictions of CREI based on the standard clamped-nuclei 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 double-peak 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 large-mass 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 Hartree-type 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 H2+. 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 nuclear-to-electronic 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 H2+. 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 ground-work for future development of accurate methods for coupled electron–ion dynamics in non-perturbative 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.


We acknowledge support from the European Research Council (ERC-2015-AdG-694097), Grupos Consolidados (IT578-13), 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 Sklodowska-Curie grant agreement no. 704218 and 702406, respectively. N. T. M. thanks the National Science Foundation, grant CHE-1566197, for support. Open Access funding provided by the Max Planck Society.


  1. E. Khosravi, A. Abedi and N. T. Maitra, Phys. Rev. Lett., 2015, 115, 263002 CrossRef PubMed.
  2. T. Zuo, S. Chelkowski and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1993, 48, 3837 CrossRef CAS.
  3. T. Zuo and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1995, 52, R2511 CrossRef CAS.
  4. T. Seideman, M. Y. Ivanov and P. B. Corkum, Phys. Rev. Lett., 1995, 75, 2819 CrossRef CAS PubMed.
  5. S. Chelkowski, C. Foisy and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1998, 57, 1176 CrossRef CAS.
  6. S. Chelkowski and A. Bandrauk, J. Phys. B: At., Mol. Opt. Phys., 1995, 28, L723 CrossRef CAS.
  7. S. Chelkowski, A. Conjusteau, T. Zuo and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 54, 3235 CrossRef CAS.
  8. H. Yu, T. Zuo and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 54, 3290 CrossRef CAS.
  9. A. D. Bandrauk and F. Légaré, Progress in Ultrafast Intense Laser Science VIII, Springer, 2012, pp. 29–46 Search PubMed.
  10. N. Takemoto and A. Becker, Phys. Rev. Lett., 2010, 105, 203004 CrossRef PubMed.
  11. N. Takemoto and A. Becker, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 84, 023401 CrossRef.
  12. C. Beylerian, S. Saugout and C. Cornaggia, J. Phys. B: At., Mol. Opt. Phys., 2006, 39, L105 CrossRef CAS.
  13. 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.
  14. 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.
  15. 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.
  16. A. Abedi, N. T. Maitra and E. K. U. Gross, Phys. Rev. Lett., 2010, 105, 123002 CrossRef PubMed.
  17. A. Abedi, N. T. Maitra and E. K. U. Gross, J. Chem. Phys., 2012, 137, 22A530 CrossRef PubMed.
  18. A. Abedi, N. T. Maitra and E. K. U. Gross, J. Chem. Phys., 2013, 139, 087102 CrossRef PubMed.
  19. A. Abedi, F. Agostini, Y. Suzuki and E. K. U. Gross, Phys. Rev. Lett., 2013, 110, 263001 CrossRef PubMed.
  20. R. Requist and E. K. U. Gross, Phys. Rev. Lett., 2016, 117, 193001 CrossRef PubMed.
  21. H. Yu, T. Zuo and A. D. Bandrauk, J. Phys. B: At., Mol. Opt. Phys., 1998, 31, 1533 CrossRef CAS.
  22. M. Uhlmann, F. Großmann, T. Kunert and R. Schmidt, Phys. Lett. A, 2007, 364, 417 CrossRef CAS.
  23. J. C. Tully, Faraday Discuss., 1998, 110, 407 RSC.
  24. The Hartree potential (14) reduces to the qs expression (13) in the limit of very localized wave-packets centered around R.
  25. J. L. Krause, K. J. Schafer and K. C. Kulander, Phys. Rev. A: At., Mol., Opt. Phys., 1992, 45, 4998 CrossRef CAS.
  26. J. Muga, J. Palao, B. Navarro and I. Egusquiza, Phys. Rep., 2004, 395, 357 CrossRef CAS.
  27. R. Santra and L. S. Cederbaum, Phys. Rep., 2002, 368, 1 CrossRef CAS.
  28. J. Javanainen, J. H. Eberly and Q. Su, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3430 CrossRef CAS.
  29. T. Kreibich, M. Lein, V. Engel and E. K. U. Gross, Phys. Rev. Lett., 2001, 87, 103901 CrossRef CAS PubMed.
  30. K. C. Kulander, F. H. Mies and K. J. Schafer, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 53, 2562 CrossRef CAS.
  31. A. Abedi, F. Agostini and E. K. U. Gross, Europhys. Lett., 2014, 106, 33001 CrossRef.
  32. S. K. Min, F. Agostini and E. K. U. Gross, Phys. Rev. Lett., 2015, 115, 073001 CrossRef PubMed.
  33. F. Agostini, S. K. Min, A. Abedi and E. K. Gross, J. Chem. Theory Comput., 2016, 12, 2127 CrossRef CAS PubMed.
  34. Y. Suzuki, A. Abedi, N. Maitra and E. K. U. Gross, Phys. Chem. Chem. Phys., 2015, 17, 29271 RSC.
  35. 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.
  36. Y. Suzuki and K. Watanabe, Phys. Rev. A, 2016, 94, 032517 CrossRef.
  37. D. M. Villeneuve, M. Y. Ivanov and P. B. Corkum, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 54, 736 CrossRef CAS.
  38. 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.
  39. I. Ben-Itzhak, 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.
  40. 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