Electronic Non-adiabatic Dynamics in Enhanced Ionization of Isotopologues of H$_2^+$ from the Exact Factorization Perspective

It was recently shown that the exact potential driving the electron's dynamics in enhanced ionization of H$_2^+$ can have large contributions arising from dynamical electron-nuclear correlation, going beyond what any electrostatics-based model can provide[1]. This potential is defined via the exact factorization of the molecular wavefunction that allows the construction of a Schr\"odinger equation for the electronic system, in which the potential contains exactly the effect of coupling to the nuclear system and any external fields. Here we study enhanced ionization in isotopologues of H$_2^+$ in order to investigate 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-electrostatic components to the potential driving the ionizing electron.

It was recently shown that the exact potential driving the electron's dynamics in enhanced ionization of H + 2 can have large contributions arising from dynamical electron-nuclear correlation, going beyond what any electrostatics-based model can provide [1]. This potential is defined via the exact factorization of the molecular wavefunction that allows the construction of a Schrödinger equation for the electronic system, in which the potential contains exactly the effect of coupling to the nuclear system and any external fields. Here we study enhanced ionization in isotopologues of H + 2 in order to investigate 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-electrostatic components to the potential driving the ionizing electron.

I. INTRODUCTION
The phenomenon of charge-resonance enhanced ionization (CREI), predicted more than twenty years ago [2][3][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 Refs. [3][4][5][6][7][8][9] for which the time-dependent Schrödinger equation (TDSE) is solved for various clamped nuclear (cn) configurations R 0 , i.e., wherê H el (r, R 0 ) =T e +Ŵ ee (r)+Ŵ en (r, R 0 )+V laser e (r, t). (2) Here r and R are used to collectively denote the electronic and nuclear coordinates,T e is the electronic kinetic energy operator, andŴ ee is the electronelectron repulsion. Furthermore,Ŵ en (r, R 0 ) contains the electron-nuclear interaction, which, for a diatomic molecular ion with one electron, as will be considered here, has the form −Z 1 /|r − R 1 | − Z 2 /|r − R 2 |, i.e. a * elham.etn@gmail.com; Corresponding author † aliabedik@gmail.com; Corresponding author ‡ angel.rubio@mpsd.mpg.de; Corresponding author § nmaitra@hunter.cuny.edu; Corresponding author 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 Starkshift upwards while the inner barrier from the internuclear Coulombic potential also grows ( Fig. 2 in Ref. [3]). Provided 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 happen, the premise of the quasi-static picture can become a little shaky: in particular, the molecule must dissociate to the CREI region, which occurs predominantly by Coulomb explosion following ionization, but to actually observe CREI rapidly enough 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][13][14]. Typically a large fraction of the nuclear density remains arXiv:1612.00507v1 [physics.atm-clus] 1 Dec 2016 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 electrostatic nuclear-electron interactions. In Ref. [1], it was shown that the nuclear dynamics can contribute significant components to the exact potential that arise from dynamical electron-nuclear correlated motion. Neglect of 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 Sec.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 Sec. III we study a one-dimensional model for the (a)symmetric isotopologues of H + 2 subject to a linearly polarized laser field for two different situations, comparing with dynamics on different approximate potentials. In Sec. 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 Sec. V.

II. EXACT FACTORIZATION APPROACH
The exact factorization of the time-dependent electron-nuclear wavefunction introduced in Refs. [15][16][17][18][19] enables the rigorous definitions of exact potentials acting on the nuclear subsystem and electronic subsystem in coupled electron-ion dynamics. These potentials follow from writing the full molecular wavefunction as a single product Ψ(r, R, t) = Φ R (r, t)χ(R, t), or Ψ(r, R, t) = χ r (R, t)Φ(r, t), where partial normalization conditions on the conditional wavefunctions Φ R (r, t) and χ r (R, t) respectively, render each factorization unique up to a gauge transformation. In the first product the equation for the nuclear wavefunction, χ(R, t), follows a TDSE, while in the second product the equation for the electronic factor Φ(r, t) is a TDSE. When the electronic dynamics is particularly of interest, as in our present study of ionization, we focus on the the second factorization, and investigate the potential that appears in the TDSE for Φ(r, t).
If we consider the case of one electron coupled to one nuclear degree of freedom in one dimension, the equa-tion for the electronic wavefunction is (in 1D we replace r and R by z and R respectively.) : (3) and that for the nuclear conditional wavefunction is: with the nuclear Hamiltonian whereT n is the nuclear kinetic energy operator, W en (z, R) (Ŵ nn (R)) is the electron-nuclear (nuclearnuclear) interaction, andv n ext (R, t) is time-dependent external potentials acting on the nuclei. Here, S(z, t) is the exact time-dependent vector potential and e (z, t) = χ z (t)|Ĥ n − i ∂ t |χ z (t) R is the exact timedependent electronic potential for electron (see next section). In one dimensional models, S(z, t) 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 electron (e-TDPES) contains the effects of the coupling to moving quantum nuclei as well as the external laser field. It can be written as the nuclear-kinetic term the 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 + 2 molecule in strong fields. Ref. [1] 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 Eq. 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)| 2 and α(z, t) = me z dz je(z ,t) is the electronic current-density. Inserting this form into Eq. 3 and setting the time-dependent vector potential to zero, gives 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 "velocity term", depends only on the electron velocity, namely on j e (z, t)/n e (z, t), while the last term, denoted "acceleration term" depends on the spatial integral of the acceleration. In Section III B and III C , we consider propagation on different contributions of the exact potential defined by the decomposition of Eq. 12, again with a view to 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 0 , as in the majority of studies of CREI [3, 5-7, 9, 21], or move classically with classical trajectoriesR(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(t)) =Ŵ en (z,R(t)) +V l e (z, t).
Considering the exact electronic potential Eq. 7, we see that such approaches completely miss the dynamical electron-nuclear correlation effects contained in T n (z, t), K cond e (z, t), 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 deltafunction atR(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 electrostatic or Hartree approximation [23] Hartree e (z, t) =V l e (z, t) + dRŴ en (z, R)n(R, t), (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., 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 Eq. 12 to complement the analysis in Ref. [1] given in terms of the decomposition in Eq. 7, as well as the electronic dynamics resulting from the following approximations: i) The qs approximation for whichR(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 timedependent nuclear density. Hence we write the corresponding electronic potential as ex−H to indicate that the exact nuclear density is substituted into the Hartree expression. iii) Due to a considerable loss of norm in the CREI regime, we normalize the electrostatic part of the potential in Eq. 14 to obtain normalized Hartree as indicated as n−ex−H when the exact nuclear density is inserted in Eq. 15. iv) the "self-consistent" Hartree approximation (SCH) in which the full wavefunction is approximated as an uncorrelated product of electronic wavefunction and nuclear wavefunction, Ψ H (r, R, t) = φ(r, t)χ(R, t) that treats the electrons and nuclei on the same footing: the electronic dynamics is described with the potential in Eq. 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 an electrostatic way only.

III. CREI IN H + 2 ISOTOPOLOGUES
We utilize a popular one-dimensional model of the symmetric as well as asymmetric isotopologues of H + 2 subject to a linearly polarized laser field. As the motion of the nuclei and the electron in the true molecule is assumed to be restricted to the direction of the polarization axis of the laser field, the essential physics can be captured by a 1D Hamiltonian featuring "soft-Coulomb" interactions [25,26]: 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 Here, x(X) refers to the 10(100) times heavier fictitious isotope of hydrogen and A + 2 is another fictitious isotopologues with the same effective nuclear mass as HT + .
is denoted as µ n = M1M2 Mn while µ e = Mn Mn+me is the electronic reduced mass with M n = M 1 + M 2 . The laser field, within the dipole approximation, is represented bŷ where E(t) denotes the electric field amplitude and q e = 1 + me Mn+me is the reduced charge and ζ = (M 2 − M 1 )/M n is the mass-asymmetry parameter. Such reduced-dimensional models have been shown to qualitatively reproduce experimental results (see Ref. [27] for example).
We study the symmetric isotopologues, i.e., H + 2 , A + 2 , D + 2 , x + 2 as well as the asymmetric isotopologues HD + , HT + , Hx + , HX + . Here x(X) stands for the fictitious isotope of hydrogen that is 10(100) times heavier than that of H, while A + 2 is another fictitious isotopologue with the same effective nuclear mass as HT + (See Table. I in which the nuclear effective mass of H + 2 and isotopologues are given.)

A. Comparison of ionization of isotopologues of H + 2 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 14 W/cm 2 , with a sine-squared pulse envelope ( Fig. 1), to each of the H + 2 isotopologues. First, we solve the electron-nuclear TDSE for Hamiltonian of Eq. 17 numerically exactly, beginning in the initial ground-state of the molecule.
As the Hamiltonian (17) 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 Eq. 18 where in the symmetric case, ζ = 0. In Figs. 2 we plot the ionization probability and average internuclear distance , R , as a function of number of cycles t/T where T denotes duration of one cycle (T = 2.67 fs) for H + 2 , HD + , HT + , Hx + , HX + , A + 2 , D + 2 and x + 2 . Note that, the results are given in atomic units, e = m e = h = 1, throughout the article, unless otherwise noted. The ionization probability is defined as a.u. As it can be seen in Figs. 2, the ionization and aver- age internuclear separation appear to be practically independent of nuclear mass-symmetry properties: they are identical for asymmetric case of HT + and symmetric case of A + 2 systems with the same effective nuclear mass µ n and very similar in the case of HX + and D + 2 where the effective masses are only a little different. One can understand this from that fact that the ionization probes the electron density in regions where z R, where the nuclear mass-dependence in the Hamiltonian Eq. 17 to 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.
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 I( [1,27]. 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 factorization framework, i. e., where, I cp follows the usual expression of the ionization probability but using the exact conditional electronic wavefunction Φ R (z, t) , i.e., I cp (R, t) = 1 − z I −z I dz|Φ R (z, t)| 2 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 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 the 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 + 2 where only an exponentially small fraction of the nuclear density dissociates, the system does not reach the CREI regime and therefore the ionization is negligible as seen from the R−resolved, t−resolved ionization probability. Notice the different scale of I(R, t) for the case of x + 2 .  Note that the ionization rate computed at any fixed R would be identical for all the isotopologues. Their different masses, however, lead to very different ionization rates when the full electron-nuclear dynamics is considered with the molecule beginning at equilibrium as we have shown. Still, there is some validation to the original quasistatic CREI prediction that ionization is enhanced for nuclear separations around R c , as indicated by I (R, t), however, a modification to the statement is needed due to the spreading and splitting of the nuclear wavepacket: the enhancement occurs from electrons associated with the part of the nuclear density that is in the R c region. Clearly, treating the nuclei as point particles will not work even for significant nuclear mass (except in the 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 electrostatically-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 + 2 subject to the same field show significantly different degrees of ionization, since the differentmass systems reach the CREI region at different times with different probabilities. In particular, systems with larger nuclear masses hardly reach the internuclear distances for which CREI is expected to occur. Hence, to study the effect of the nuclear mass in the CREI regime for different-mass systems, we adjust the field intensity such that the ionization probability(rate) remains close to that of H + 2 . As the asymmetric and symmetric isotopologues with the same effective nuclear mass subject to the same field give almost the same ionization probability, average internuclear distance and I(R, t) (see section III A), from here on we focus only on the symmetric isotopologues of H + 2 . The optimized field intensity for different-mass systems is given in Table. II while the other field parameters are kept unchanged.
The ionization probability of the symmetric isotopologues subject to the optimized field is depicted in Fig. 4. For isotopologues heavier than H + 2 the ionization starts to set in about one optical cycle T later than the H + 2 case. In order to be able to compare the ionization yield (rate) of these systems with H + 2 visually better, in Fig. 4 we have shifted the time-dependent ionization probability  table II. by one optical cycle, i.e. for isotopologues heavier than H + 2 , I(t − T ) has been 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 + 2 . That is before the field reaches its maximum the system becomes slightly ionized hence the nuclear density slightly spreads. As the field reaches its maximum a small fragment of the nuclear density starts to split off and dissociate from Coulomb explosion due to an increase in the ionization while the rest of the nuclear density remains bounded and oscillates around the equilibrium position. As an appreciable amount of dissociating fragment reaches the internuclear separation for CREI, the ionization gets enhanced. In the case of x + 2 , however, the nuclear dynamics exhibits a more classical behavior, i. e. during the first half of the pulse, the heavy nuclei only spreads 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 to 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 + 2 and D + 2 , the I(R, t) has a similar structure to H + 2 , the I(R, t) corresponding to x + 2 manifests a substantially distinctive structure compared to the lighter isotopologues: the internuclear separation for which ionization gets enhanced, shifts to smaller values. In general, the I(R, t) shifts to smaller R for heavier isotopologues which could be related to the stronger optimized field used [5,6]. Now, we investigate the performance of the various approximations introduced in Sec. II B for differentmass 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 selfconsistent Hartree (SCH) for various symmetric isotopologues of H + 2 are plotted. 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 doesn't 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 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 Eq. 15. In fact, the n-ex-H approximation performs clearly better than other conventional approximations (qs, SCH, ex-H) for H + 2 , A + 2 and x + 2 . Even the approximate potential that depends on the more rigorous and complicated concept of conditional nuclear density rather than the nuclear density that appears in the 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 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 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: dynamical 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, and adjusting the field to get similar ionization results in these being similar (See also Sec 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 (Eq. 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, expect in 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 dynamical electron-nuclear terms in 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 e-TDPES (Eq. 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 + 2 and x + 2 . We refer the readers to Ref. [1] for a discussion on the components of the conditional decomposition of e-TDPES. For the sake of simplifying the discussion, here, we divide the electronic coordinate into two regions: "inner-region" that refers to the region with |z| < 5 a.u. and "outer-region" that describes the rest of the axis for which |z| > 5 a.u.

H + 2 case
In Fig. 8, we present the terms/combinations-of-terms of marginal decomposition of the e-TDPES for the case of H + 2 at four different snapshots in time. Adiabatic term: the adiabatic term (first term in Eq. 12) is the main constituent of the e-TDPES when it is decomposed according to Eq. 12. In this case, as it is seen in Fig. 8 (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 innerregion and outer-region. It can be seen in Fig. 8 (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 the exact potential are mainly responsible for controlling the ionization from the up-field direction. Indeed the over-ionization 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 to 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 sides (up-field and down-field) in each half cycle, resulting in a huge overestimation of ionization probability. 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 to 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: is the second term in Eq. 12 and its overall contribution to the e-TDPES, in this case, is small especially in the inner region. As it 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: is the last term in Eq. 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 peakintensity 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 (lower left panel).
As the structures of the adiabatic, velocity, and acceleration terms in the case of A + 2 and D + 2 are very similar to H + 2 , we do not address them here.

x + 2 case
In Fig. 9, the terms/combinations-of-terms of marginal decomposition of the e-TDPES for the case of x + 2 are plotted at four different times. Adiabatic term: behaves similarly to the case of H + 2 , 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 adia- batic 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 + 2 and lowers the interatomic barrier when added to the adiab+acc potential for most of the times between t ≈ 28 T and t ≈ 40 T in which most of the ionization happens. The more dominant role of the velocity term in this case, could be attributed to the increase of ponderomotive (wiggling) motion of the electron driven by a stronger external laser field.
Acceleration term: this term has a correct asymptotic behavior, similar to the case of H + 2 but leads to a wrong estimate for the inner-barrier (see Fig. 9, lowerleft 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  (uper 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 dynam- ics 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 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 [28][29][30][31][32][33]. Here we begin taking the classical limit of the conditional nuclear wavefunction by representing it in polar form: 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), Eq. 4, and sorting the terms in orders of , we find to O( 0 ): where we define V (z, R, t) = W nn (R) + W en (z, R) + V l (z, R, t), andp e (z, t) = −i ∂ Φ(z, t)/∂z (see shortly for more on this term). Similarly, keeping only O( 0 ) terms in Eq. 7 for the e-TDPES, we find cl e (z, t) = dR|χ z (R, t)| 2 1 2µ n The terms on the right-hand-side correspond to classically evaluating T n (z, t), K cond e (z, t), app e (z, t) 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 zeroth-order in with the term 1 2µe ∂Sz ∂z 2 in both Eq (21) and (22).
Eq. 21 would be a standard Hamilton-Jacobi equation of the form H(q, ∂Sz ∂q , t) + ∂S z /∂t = 0 (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 Eq. 21 by defining the velocity fields, and where d/dt = ∂/∂t + u e z ∂/∂z + u n z ∂/∂R is the timederivative in the Lagrangian frame defined by the velocities.

A. Quasiclassical analysis of the terms in the e-TDPES
Whether the equations above, together with the solution of the TDSE Eq. 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-massdependence 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. [34] showed that classical trajectory calculations reproduce the essential features of CREI for both cases of fixed and moving nuclei.
To this end, we first evaluatep e (z, t) in Eq. 21 to its lowest-order in .
Inserting Φ(z, t) = a(z, t) exp(is(z, t))/ , analogous to Eq. 20, into the TDSE Eq. 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 HamiltonianT e + e (z, t). HereT e = 1 2µe ∂s ∂z 2 is the electronic kinetic energy associated with the Hamiltonian in the TDSE for the electronic wavefunction, Eq. 3. It is important to note that this is not the same as the true electronic kinetic energy (see Eq. 69 in Ref. [17]). The term in Eq. 21, −i ∂Φ(z, t)/∂z thus becomesp e ≡ 2µ eTe in the classical limit. Next, we note that Eq. 21 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 where, in the classical limit, s(z, t) satisfies Eq. 26 and S(z, R, t) satisfies For the first term in Eq. 22, T n (z, t), where we noted that ∂S z (R, t)/∂R = ∂S(z, R, t)/∂R (from Eq. 27), 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, and so we expect a potential that rises as z gets larger. This is indeed what we see in Fig. 10. Consider first the black curve, H + 2 . At early times, the middle region of T n (z, t), 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 zeropoint motion. As the electrons begin to gain energy from the field and move out from the tails of the electronic distribution, the T n (z, t) 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 T n (z, t) gets narrower as ionization take away more of the electron density. As the final distribution sets in, T n (z, t), 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.
Regarding now the µ n -dependence: it might appear from Eq. 29 that T n (z, t) grows as the mass µ n increases, in the cases where the field is adjusted so that the average internuclear distance and speed remains the same. This is in fact not the case; in fact, the term remains about the same 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˙ R 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 T n (z, t) as one moves to larger z, as mass increases, compensating the larger µ n . In Fig. 10 we plot T n (z, t) for the different isotopologues at various snapshots in time. We see that the overall size of this term does not vary much with mass except for x + 2 . The exception is for the large mass x + 2 where there is comparatively weak z-dependence of T n (z, t). 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 Eq. 22 is unfortunately not as straightforward, because the evolution of the relevant classical trajectories themselves depend on the potential e (z, t), which is the very object we wish to analyse! For example, the second term, K cond 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 exter-nal 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 this work, we have presented a detailed investigation of the exact potential driving the electron dynamics in isotopologues of H + 2 undergoing CREI and studied the dependence of the correlated electron-nuclear dynamics on the nuclear mass.
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 dynamical case. For the laser parameters and different isotopologues of H + 2 discussed in this work, I(R, t), exhibits only one peak. This is in agreement with most of the experimental findings [35,36] and 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 double-peak structure could in fact be observed as in Ref. [37].
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 an electrostatic way (as in Hartree-type approximations) is not adequate in capturing the dynamics accurately. The importance of going beyond an electrostatic description of electron-nuclear correlation is evident from our studies of how different approximate electronic potentials perform in describing the CREI for isotopologues of H + 2 . We have shown that, for the laser parameters used in this work, one must go beyond any purely electrostatic treatment of electron-nuclear correlation, and include truly dynamical 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 approx-imate potentials from the exact potential, what are more important than the nuclear-to-electronic mass ratios, are the nuclear velocities in the wavepacket.
There are many different approaches to developing approximate methods based on the exact factorization. For example, one may consider the relative importance of different components of the exact potential when decomposed in terms of the conditional wavefunction as in Ref. [1] (there, it was found that generally all terms were important). One may consider also approximations based on the marginal decomposition presented here, where we found that the "adiab+acc" component of the marginal decomposition describes the electronic dynamics accurately in all cases we studied here with the exception of a fictitious isotopologue with an effective nuclear mass 10 times larger than H + 2 . Further investigation of this potential may lead to development of an adequate approximation for practical purposes capable of describing the ionization dynamics accurately. Another approach is to develop quasiclassical approximations, and here we have sketched out a semiclassical derivation of the electronic and conditional nuclear equations of the exact factorization in its reverse form and analysed the structure of the nuclear kinetic term of the exact potential semiclassically. This work has highlighted the effect of the complex interplay of 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 dynamical electron-nuclear effects play as crucial a role for polyatomic molecules, and the scaling of these effects with respect to the number of electrons and number of nuclei, is also an important avenue for future research