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

Abstract

It was recently shown that the exact potential driving the electron's dynamics in enhanced ionization of H₂⁺ 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 H₂⁺ 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 , i.e.,whereHere and are used to collectively denote the electronic and nuclear coordinates, T̂_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₁/|r − R₁| − Z₂/|r − R₂|, 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 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 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 H₂⁺ 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 , 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 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 .

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):

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 Hamiltonianwhere T̂_n is the nuclear kinetic energy operator, Ŵ_en(z, R) (Ŵ_nn(R)) is the electron–nuclear (nuclear–nuclear) interaction, and V̂^laser_n(R, t) is the time-dependent external potential acting on the nuclei. Here, is the exact time-dependent vector potentialand ε_e(z, t) = 〈χ_z(t)|Ĥ_n − iħ ∂_t|χ_z(t)〉_R is the exact time-dependent 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 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 aswhich consists of: the approximate potential

ε^app_e(z, t) = V̂^laser_e(z, t) + 〈χ_z(R, t)|Ŵ_en(z, R) + Ŵ_nn(R)|χ_z(R, t)〉_R, (8)

the nuclear-kinetic termthe gauge dependent part of the potential

ε^gd_e(z, t) = ħ〈χ_z(R, t)| − i∂_tχ_z(R, t)〉_R,

and finally the electronic-kinetic-like 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₂⁺ 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,

with n_e(z, t) = |Φ(z, t)|² and whereis 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̂∂/∂_z). Inserting this form into eqn (3) and setting the time-dependent vector potential to zero give

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 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 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̄₀, as in the majority of studies of CREI,^3,5–7,9,21 or move classically with classical trajectories R̄(t) that are often described by mixed quantum-classical algorithms such as Ehrenfest or surface-hopping algorithms.²² 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̄(t)) = Ŵ_en(z, R̄(t)) + V̂^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 z-independent delta-function at R̄(t), i.e. n_z(R, t) ≈ δ(R(t), R̄(t)) (in this limit Ŵ_nn(R̄(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 approximation²³

where n(R, t) is the nuclear density obtained from nuclear wavepacket dynamics.²⁴ It can be easily seen that if the z-dependence 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 R̄(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

indicated as ε^n-ex-H when the exact nuclear density, obtained from the complete electron–nuclear wavefunction Ψ(z, R), i.e.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, , 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.,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 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 H₂⁺ isotopologues

We utilize a popular one-dimensional model of the symmetric as well as asymmetric isotopologues of H₂⁺ 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,29where 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 while is the electronic reduced mass with M_n = M₁ + M₂. The laser field, within the dipole approximation, is represented by

V̂^laser(R, z, t) = eE(t)(q_ez − ζR), (19)

where E(t) denotes the electric field amplitude is the reduced charge, and ζ = (M₂ − M₁)/M_n 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., H₂⁺, A₂⁺, D₂⁺, x₂⁺ 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₂⁺ is another fictitious isotopologue with the same effective nuclear mass as HT⁺ (see Table 2 in which the nuclear effective mass values of H₂⁺ and isotopologues are given).

Table 2 Nuclear effective mass corresponding to H₂⁺ 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₂⁺ 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 H₂⁺ 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 × 10¹⁴ W cm⁻², with a sine-squared pulse envelope (Fig. 1), to each of the H₂⁺ isotopologues. First, we solve the electron–nuclear TDSE for Hamiltonian of eqn (18) numerically exactly, beginning in the initial ground-state 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₂⁺, HD⁺, HT⁺, Hx⁺, HX⁺, A₂⁺, D₂⁺ and x₂⁺. 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,

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 mass-symmetry properties: they are identical for the asymmetric case of HT⁺ and the symmetric case of A₂⁺ systems with the same effective nuclear mass μ_n and are very similar in the case of HX⁺ and D₂⁺ 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 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.

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₂⁺, HD⁺, HT⁺, Hx⁺ and HX⁺ molecules on the left panels and H₂⁺, A₂⁺, D₂⁺ and x₂⁺ 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 , 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)|²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 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 R_c. The lower set of panels of Fig. 3 shows the time-dependent nuclear density in each case. For the case of x₂⁺ 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 x₂⁺.

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 H₂⁺, D₂⁺ and x₂⁺.

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 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 H₂⁺ 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 H₂⁺. 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₂⁺. 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⁻²

H2^+	A2^+	D2^+	x2^+
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₂⁺ the ionization starts to set in about one optical cycle T later than the H₂⁺ case. In order to be able to compare the ionization yield (rate) of these systems with H₂⁺ visually better, in Fig. 4 we have shifted the time-dependent ionization probability by one optical cycle, i.e. for isotopologues heavier than H₂⁺, 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₂⁺ subject to a 50-cycle pulse of wavelength λ = 800 nm (ω = 0.0569 a.u.) and intensity I = 2.02 × 10¹⁴ W cm⁻², with a sine-squared pulse envelope. For the isotopologues of H₂⁺, to view better how the ionization yield (rate) compares to that of H₂⁺ 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 x₂⁺. 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₂⁺, 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 A₂⁺ and D₂⁺, I(R, t) has a similar structure to H₂⁺, I(R, t) corresponding to x₂⁺ 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 time-dependent nuclear density for H₂⁺ (lower left panel), A₂⁺ (upper left panel), D₂⁺ (lower right panel), x₂⁺ (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₂⁺ (lower left panel), A₂⁺ (upper left panel), D₂⁺ (lower right panel), x₂⁺ (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 H₂⁺ are plotted.

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 H₂⁺, A₂⁺ and x₂⁺. 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 x₂⁺, 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 H₂⁺ and x₂⁺. 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. H₂⁺ case. In Fig. 8, we present the terms/combinations-of-terms of marginal decomposition of the e-TDPES for the case of H₂⁺ at four different snapshots in time.

Fig. 8 The exact e-TDPES and its component from the marginal decompositions at various snapshots in time of H₂⁺.

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 A₂⁺ and D₂⁺ are very similar to H₂⁺, we do not address them here.

2. x₂⁺ case. In Fig. 9, the terms/combinations-of-terms of marginal decomposition of the e-TDPES for the case of x₂⁺ are plotted at four different times.

Fig. 9 The exact e-TDPES and its component from the marginal decompositions at various snapshots in time of x₂⁺.

Adiabatic term. The adiabatic term behaves similar to the case of H₂⁺, 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 H₂⁺ 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₂⁺ 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) = A_z(R, t)e^{iS_z(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(ħ⁰):where we define V(z, R, t) = W_nn(R) + W_en(z, R) + V^laser(z, R, t), and p̃_e(z, t) = −iħ∂Φ(z, t)/∂z (see shortly for more on this term). Similarly, keeping only O(ħ⁰) terms in eqn (7) for the e-TDPES, we findThe terms on the right-hand-side 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 zeroth-order 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 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,

Then, taking ∂/∂R of eqn (23) yieldsandwhere d/dt = ∂/∂t + u^e_z∂/∂z + uⁿ_z∂/∂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̃_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

This is the Hamilton–Jacobi equation for the action s(z, t) of a classical electron evolving in Hamiltonian T̃_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), satisfiesIt 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 = P_n = μ_nṘ,

the nuclear momentum;

∂S/∂z = p_e,

the electronic momentum;

∂S/∂t = E(t),

the classical energy; and, as in the earlier discussion,

∂s/∂z = p̃_e,

and

∂s/∂t = 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 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), ,

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 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, H₂⁺. 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 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, 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₂⁺ 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₂⁺. The exception is for the large mass x₂⁺ where there is comparatively weak z-dependence 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 T̃_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 work¹ we had presented the exact electronic potential driving the electron dynamics of H₂⁺ 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₂⁺ 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 H₂⁺ discussed in this work, I(R, t) exhibits only one peak along the R-axis. 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 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 H₂⁺. 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 H₂⁺. 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.

Acknowledgements

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.

References

1. E. Khosravi, A. Abedi and N. T. Maitra, Phys. Rev. Lett., 2015, 115, 263002.
2. T. Zuo, S. Chelkowski and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1993, 48, 3837.
3. T. Zuo and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1995, 52, R2511.
4. T. Seideman, M. Y. Ivanov and P. B. Corkum, Phys. Rev. Lett., 1995, 75, 2819.
5. S. Chelkowski, C. Foisy and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1998, 57, 1176.
6. S. Chelkowski and A. Bandrauk, J. Phys. B: At., Mol. Opt. Phys., 1995, 28, L723.
7. S. Chelkowski, A. Conjusteau, T. Zuo and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 54, 3235.
8. H. Yu, T. Zuo and A. D. Bandrauk, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 54, 3290.
9. A. D. Bandrauk and F. Légaré, Progress in Ultrafast Intense Laser Science VIII, Springer, 2012, pp. 29–46.
10. N. Takemoto and A. Becker, Phys. Rev. Lett., 2010, 105, 203004.
11. N. Takemoto and A. Becker, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 84, 023401.
12. C. Beylerian, S. Saugout and C. Cornaggia, J. Phys. B: At., Mol. Opt. Phys., 2006, 39, L105.
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.
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.
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.
16. A. Abedi, N. T. Maitra and E. K. U. Gross, Phys. Rev. Lett., 2010, 105, 123002.
17. A. Abedi, N. T. Maitra and E. K. U. Gross, J. Chem. Phys., 2012, 137, 22A530.
18. A. Abedi, N. T. Maitra and E. K. U. Gross, J. Chem. Phys., 2013, 139, 087102.
19. A. Abedi, F. Agostini, Y. Suzuki and E. K. U. Gross, Phys. Rev. Lett., 2013, 110, 263001.
20. R. Requist and E. K. U. Gross, Phys. Rev. Lett., 2016, 117, 193001.
21. H. Yu, T. Zuo and A. D. Bandrauk, J. Phys. B: At., Mol. Opt. Phys., 1998, 31, 1533.
22. M. Uhlmann, F. Großmann, T. Kunert and R. Schmidt, Phys. Lett. A, 2007, 364, 417.
23. J. C. Tully, Faraday Discuss., 1998, 110, 407.
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.
26. J. Muga, J. Palao, B. Navarro and I. Egusquiza, Phys. Rep., 2004, 395, 357.
27. R. Santra and L. S. Cederbaum, Phys. Rep., 2002, 368, 1.
28. J. Javanainen, J. H. Eberly and Q. Su, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3430.
29. T. Kreibich, M. Lein, V. Engel and E. K. U. Gross, Phys. Rev. Lett., 2001, 87, 103901.
30. K. C. Kulander, F. H. Mies and K. J. Schafer, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 53, 2562.
31. A. Abedi, F. Agostini and E. K. U. Gross, Europhys. Lett., 2014, 106, 33001.
32. S. K. Min, F. Agostini and E. K. U. Gross, Phys. Rev. Lett., 2015, 115, 073001.
33. F. Agostini, S. K. Min, A. Abedi and E. K. Gross, J. Chem. Theory Comput., 2016, 12, 2127.
34. Y. Suzuki, A. Abedi, N. Maitra and E. K. U. Gross, Phys. Chem. Chem. Phys., 2015, 17, 29271.
35. F. Agostini, A. Abedi, Y. Suzuki, S. K. Min, N. T. Maitra and E. K. U. Gross, J. Chem. Phys., 2015, 142, 084303.
36. Y. Suzuki and K. Watanabe, Phys. Rev. A, 2016, 94, 032517.
37. D. M. Villeneuve, M. Y. Ivanov and P. B. Corkum, Phys. Rev. A: At., Mol., Opt. Phys., 1996, 54, 736.
38. T. Ergler, A. Rudenko, B. Feuerstein, K. Zrost, C. D. Schröter, R. Moshammer and J. Ullrich, Phys. Rev. Lett., 2005, 95, 093001.
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.
40. H. Xu, F. He, D. Kielpinski, R. Sang and I. Litvinyuk, Sci. Rep., 2015, 5, 13527.

---