Electronic Non-adiabatic Dynamics in Enhanced Ionization of Isotopologues of H + 2 from the Exact Factorization Perspective
EElectronic Non-adiabatic Dynamics in Enhanced Ionization of Isotopologues of H +2 from theExact Factorization Perspective Elham Khosravi, ∗ Ali Abedi, † Angel Rubio,
2, 1, ‡ and Neepa T. Maitra § Nano-Bio Spectroscopy Group and ETSF, Universidad del Pa´ıs Vasco, CFM CSIC-UPV/EHU, 20018 San Sebasti´an, Spain Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany Department of Physics and Astronomy, Hunter College and the Graduate Center ofthe City University of New York, 695 Park Avenue, New York, New York 10065, USA (Dated: November 11, 2018)It was recently shown that the exact potential driving the electron’s dynamics in enhanced ioniza-tion of H +2 can have large contributions arising from dynamical electron-nuclear correlation, goingbeyond what any electrostatics-based model can provide [1]. This potential is defined via the exactfactorization of the molecular wavefunction that allows the construction of a Schr¨odinger equationfor the electronic system, in which the potential contains exactly the effect of coupling to the nuclearsystem and any external fields. Here we study enhanced ionization in isotopologues of H +2 in order toinvestigate nuclear-mass-dependence of these terms for this process. We decompose the exact poten-tial into components that naturally arise from the conditional wavefunction, and also into componentsarising from the marginal electronic wavefunction, and compare the performance of propagation onthese different components as well as approximate potentials based on the quasi-static or Hartreeapproximation with the exact propagation. A quasiclassical analysis is presented to help analyse thestructure of different non-electrostatic components to the potential driving the ionizing electron. I. INTRODUCTION
The phenomenon of charge-resonance enhanced ion-ization (CREI), predicted more than twenty yearsago [2–4], is a prominent example of the complex cou-pling of electronic motion, ionic motion, and stronglaser fields. At a critical range of internuclear sep-arations, the ionization rate of a molecule in a laserfield can be orders of magnitude greater than the ratefrom the constituent atoms. The ionization rate hasbeen explained by a quasi-static argument involving in-stantaneously frozen nuclei in the pioneering works ofRefs. [3–9] for which the time-dependent Schr ¨odingerequation (TDSE) is solved for various clamped nuclear(cn) configurations R , i.e., ˆ H el ( r , R )Φ cn R ( r , t ) = i (cid:126) ∂ t Φ cn R ( r , t ) , (1)where ˆ H el ( r , R ) = ˆ T e + ˆ W ee ( r )+ ˆ W en ( r , R )+ ˆ V lasere ( r , t ) . (2)Here r and R are used to collectively denote theelectronic and nuclear coordinates, ˆ T e is the elec-tronic kinetic energy operator, and ˆ W ee is the electron-electron repulsion. Furthermore, ˆ W en ( r , R ) containsthe electron-nuclear interaction, which, for a diatomicmolecular ion with one electron, as will be consideredhere, has the form − Z / | r − R | − Z / | r − R | , i.e. a ∗ [email protected]; Corresponding author † [email protected]; Corresponding author ‡ [email protected]; Corresponding author § [email protected]; Corresponding author Coulombic double-well. The critical internuclear sep-aration for CREI can be qualitatively explained by thefollowing argument: at stretched geometries in a givenstatic field, the energy levels in the up-field atom Stark-shift upwards while the inner barrier from the internu-clear Coulombic potential also grows (Fig. 2 in Ref. [3]).Provided the field is turned on fast enough such thatany population in the up-field level (LUMO) does notsignificantly tunnel back to the down-field atomic level(HOMO), the molecule can rapidly ionize over both theinner and outer barriers, which gives rise to the ioniza-tion rate observed to be enhanced by orders of magni-tude compared to the atomic rate. By requiring that theStark-shifted LUMO level exceeds the top of both theinner and outer field-modified Coulombic barriers, onefinds R c = 4 . /I p for the critical internuclear separa-tion for CREI. The analysis can be generalized to thecase of a laser field, where the field’s period is shorterthan the tunneling time [2, 3, 10, 11].It has been pointed out that, in actuality, the under-lying assumption in this picture of the electron adiabat-ically following the field is not quite adequate, as theionization tends to occur in multiple sub-cycle ioniza-tion bursts, not at the peak of the field cycles [10, 11].Moreover, when applied to an experiment where themolecule’s initial geometry is far from where the CREIis expected to happen, the premise of the quasi-staticpicture can become a little shaky: in particular, themolecule must dissociate to the CREI region, which oc-curs predominantly by Coulomb explosion followingionization, but to actually observe CREI rapidly enoughthat appreciable electron density still remains largelyun-ionized. In fact, in many experiments, CREI is notobserved because too little electron density reaches theCREI region over the course of the applied field [12–14].Typically a large fraction of the nuclear density remains a r X i v : . [ phy s i c s . a t m - c l u s ] D ec ear equilibrium, while a part of it dissociates. So, (i)representing its potential on the electron as a Coulombicdouble-well is not appropriate, and (ii) the fragment ve-locities can be comparable to the electronic velocities insome dissociating channels questioning the very notionof electrostatic nuclear-electron interactions. In Ref. [1],it was shown that the nuclear dynamics can contributesignificant components to the exact potential that arisefrom dynamical electron-nuclear correlated motion. Ne-glect of these contributions leads to severe errors in theprediction of the electronic motion.In this paper, we analyse the exact potential in moredetail by providing a complementary decomposition tothe one that was studied in Ref. [1], comparing with dif-ferent approximations, studying the mass-dependence,and investigating a quasiclassical treatment. The rest ofthe paper is organized as follows: First, a brief reviewof the exact factorization approach is presented in Sec.II,with a focus on the reverse factorization introduced inRef. [15]. This provides us with a TDSE for electronsevolving on a single time-dependent potential that ac-counts for the coupling to the nuclear system and anyexternal field in an exact way. We present two differentways of decomposing this potential and some approxi-mate 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 linearlypolarized laser field for two different situations, com-paring with dynamics on different approximate poten-tials. In Sec. IV we take a first step to analyse quasiclas-sically the structure of various components of the poten-tial driving the ionizing electron. Finally, some conclud-ing remarks are presented in Sec. V. II. EXACT FACTORIZATION APPROACH
The exact factorization of the time-dependentelectron-nuclear wavefunction introduced in Refs. [15–19] enables the rigorous definitions of exact potentialsacting on the nuclear subsystem and electronic subsys-tem in coupled electron-ion dynamics. These potentialsfollow from writing the full molecular wavefunctionas a single product Ψ( r , R , t ) = Φ R ( r , t ) χ ( R , t ) , or Ψ( r , R , t ) = χ r ( R , t )Φ( r , t ) , where partial normalizationconditions on the conditional wavefunctions Φ R ( r , t ) and χ r ( R , t ) respectively, render each factorizationunique up to a gauge transformation. In the firstproduct the equation for the nuclear wavefunction, χ ( R , t ) , follows a TDSE, while in the second productthe 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 thethe second factorization, and investigate the potentialthat appears in the TDSE for Φ( r , t ) .If we consider the case of one electron coupled to onenuclear 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.) : (cid:34) ( − i (cid:126) ∂/∂z + S ( z, t )) m e + (cid:15) e ( z, t ) (cid:35) Φ( z, t ) = i (cid:126) ∂ t Φ( z, t ) . (3)and that for the nuclear conditional wavefunction is: (cid:104) ˆ H n ( z, R, t ) − (cid:15) e ( z, t ) (cid:105) χ z ( R, t ) = i (cid:126) ∂ t χ z ( R, t ) , (4)with the nuclear Hamiltonian ˆ H n ( z, R, t ) = ˆ T n ( R ) + ˆ W en ( z, R ) + ˆ W nn ( R ) + ˆ v n ext ( R, t )+ 1 m e (cid:20) [ − i (cid:126) ∂/∂ z − S ( z, t )]
2+ ( − i (cid:126) ∂ Φ /∂ z Φ + S ( z, t ))( − i∂/∂ z − S ( z, t )) (cid:21) (5)where ˆ T n is the nuclear kinetic energy operator, ˆ W en ( z, R ) ( ˆ W nn ( R ) ) is the electron-nuclear (nuclear-nuclear) interaction, and ˆ v n ext ( R, t ) is time-dependent ex-ternal potentials acting on the nuclei. Here, S ( z, t ) is theexact time-dependent vector potential S ( z, t ) = (cid:104) χ z ( R, t ) | − i (cid:126) ∂/∂ z χ z ( R, t ) (cid:105) R (6)and (cid:15) e ( z, t ) = (cid:104) χ z ( t ) | ˆ H n − i (cid:126) ∂ t | χ z ( t ) (cid:105) R is the exact time-dependent electronic potential for electron (see next sec-tion). In one dimensional models, S ( z, t ) can be set tozero as a choice of gauge and we adopt this gauge forthe rest of this paper. In this case, (cid:15) e ( z, t ) is the sole po-tential that drives the electronic motion, which can becompared with the other traditional potentials that areused to study electronic dynamics. A. Decomposition of the exact time-dependent potentialenergy surface
The exact time-dependent potential energy surface forelectron ( e -TDPES) contains the effects of the couplingto moving quantum nuclei as well as the external laserfield. It can be written as (cid:15) e ( z, t ) = (cid:15) app e ( z, t ) + T n ( z, t ) + K cond e ( z, t ) + (cid:15) gd e ( z, t ) , (7)which consists of: the approximate potential (cid:15) app e ( z, t ) = (cid:104) χ z ( R, t ) | ˆ v e ext ( z, t ) | χ z ( R, t ) (cid:105) R + (cid:104) χ z ( R, t ) | ˆ W en ( z, R ) + ˆ W nn ( R ) | χ z ( R, t ) (cid:105) R , (8)the nuclear-kinetic term T n ( z, t ) = − (cid:126) m n (cid:104) χ z ( R, t ) | ∂ ∂R | χ z ( R, t ) (cid:105) R , (9)the gauge dependent part of the potential (cid:15) gd e ( z, t ) = (cid:126) (cid:104) χ z ( R, t ) | − i∂ t χ z ( R, t ) (cid:105) R , K cond e ( z, t ) = (cid:126) m e (cid:28) ∂∂z χ z ( R, t ) | ∂∂z χ z ( R, t ) (cid:29) R . In Ref. [1], we had found that this exact potential (cid:15) e ( z, t ) ,has significant features that are missing in the tradi-tional potentials based upon the quasistatic picture de-scribed above. Neglecting these features led to quali-tatively incorrect predictions of ionization dynamics inthe H +2 molecule in strong fields. Ref. [1] found that,in general, all the four terms above are needed to rea-sonably reproduce the ionization in CREI processes. i.e.that propagation on any combination other than the fullsum of the four contributions in Eq. 7 gave qualitativelypoor results. This means that in eventually developingapproximations, all of the four terms must be consid-ered.Alternatively, one may instead decompose the e -TDPES by exactly inverting the TDSE for the electronicwavefunction. The electronic wavefunction may bewritten in polar representation, Φ( z, t ) = (cid:112) n e ( z, t ) exp ( iα ( z, t )) (10)with n e ( z, t ) = | Φ( z, t ) | and α ( z, t ) = m e (cid:126) (cid:82) z dz (cid:48) j e ( z (cid:48) ,t ) n e ( z (cid:48) ,t ) where j e ( z, t ) = (cid:126) m e (cid:61) (Φ ∗ ( z, t ) ∇ e Φ( z, t )) , (11)is the electronic current-density. Inserting this form intoEq. 3 and setting the time-dependent vector potential tozero, gives (cid:15) e ( z, t ) = (cid:126) m e (cid:34) ∇ e (cid:112) n e ( z, t ) (cid:112) n e ( z, t ) (cid:35) − m e (cid:18) j e ( z, t ) n e ( z, t ) (cid:19) − ∂ t α ( z, t ) (12)The first term is what survives in the absence of anydynamics, and it depends only on the instantaneouselectron density; we denote it “adiabatic” in the spiritof time-dependent density functional theory. The sec-ond term, denoted “velocity term”, depends only onthe electron velocity, namely on j e ( z, t ) /n e ( z, t ) , whilethe last term, denoted “acceleration term” depends onthe spatial integral of the acceleration. In Section III Band III C , we consider propagation on different contri-butions of the exact potential defined by the decompo-sition of Eq. 12, again with a view to eventually devel-oping approximations, possibly density functionals (seeRef. [20]), for the exact e -TDPES. B. Approximate electronic potentials based onconventional approximations
The standard approximation for the electronic poten-tial assumes the so-called quasistatic approximation (qs) that treats the nuclei as classical point particles with po-sitions that are either considered fixed , ¯ R , as in themajority of studies of CREI [3, 5–7, 9, 21], or move clas-sically with classical trajectories ¯ R ( t ) that are often de-scribed by mixed quantum-classical algorithms such asEhrenfest or surface-hopping algorithms [22]. In thesemethods, electrons on the other hand, regardless ofwhether the nuclei are frozen or move classically, fol-low the combined potential from the laser field and theelectrostatic attraction of the nuclei, i. e., (cid:15) qs ( z, t | ¯ R ( t )) = ˆ W en ( z, ¯ R ( t )) + ˆ V le ( z, t ) . (13)Considering the exact electronic potential Eq. 7, we seethat such approaches completely miss the dynamicalelectron-nuclear correlation effects contained in T n ( z, t ) , K cond e ( z, t ) , and (cid:15) gd , and can be viewed as an approx-imation to (cid:15) app alone: (cid:15) app reduces to the qs approx-imation when the conditional nuclear wavefunctionis approximated classically as a z -independent delta-function at ¯ R ( t ) , i.e. n z ( R, t ) ≈ δ ( R ( t ) , ¯ R ( t )) (in thislimit ˆ W nn ( ¯ R ( t )) becomes purely a time-dependent con-stant and hence is dropped hereafter).A step beyond the qs approximation for the electronicpotential that can, in principle, account for the widthand splitting of the nuclear wavefunction is the electro-static or Hartree approximation [23] (cid:15) Hartree e ( z, t ) = ˆ V le ( z, t ) + (cid:90) dR ˆ W en ( z, R ) n ( R, t ) , (14)where n ( R, t ) is the nuclear density obtained from nu-clear wavepacket dynamics [24]. It can be easily seenthat if the z -dependence in the conditional nuclearwavefunction is neglected, i. e., n z ( R, t ) ≈ n ( R, t ) , theapproximate potential simplifies to the Hartree approx-imation.To provide a detailed analysis of the exact electronicpotential, we compare the electronic dynamics result-ing from different approximations. We will considercombinations of the terms of the exact potential decom-posed according to Eq. 12 to complement the analy-sis in Ref. [1] given in terms of the decomposition inEq. 7, as well as the electronic dynamics resulting fromthe following approximations: i) The qs approximationfor which ¯ R ( t ) = (cid:104) R (cid:105) ( t ) is the average time-dependentinternuclear separation obtained from the exact calcu-lations. ii) The Hartree approximation for which thenuclear density in (14) is replaced by the exact time-dependent nuclear density. Hence we write the corre-sponding electronic potential as (cid:15) ex − H to indicate thatthe exact nuclear density is substituted into the Hartreeexpression. iii) Due to a considerable loss of norm in theCREI regime, we normalize the electrostatic part of thepotential in Eq. 14 to obtain normalized Hartree as (cid:15) n − Hartree ( z, t ) = ˆ V le ( z, t )+ (cid:82) dR ˆ W en ( z, R ) n ( R, t ) (cid:82) dR n ( R ) , (15)3
10 20 30 40 50t / T-1-0.500.51 E ( t ) / E FIG. 1. Laser field as a function of number of optical cycles t/T .The electric field amplitude is divided by the peak amplitude, E = √ I . indicated as (cid:15) n − ex − H when the exact nuclear density isinserted in Eq. 15. iv) the “self-consistent” Hartree ap-proximation (SCH) in which the full wavefunction isapproximated as an uncorrelated product of electronicwavefunction and nuclear wavefunction, Ψ H ( r , R , t ) = φ ( r , t ) χ ( R , t ) that treats the electrons and nuclei on thesame footing: the electronic dynamics is described withthe potential in Eq. 14 that is coupled to the nuclei thatmove under the influence of the analogous potential,i.e., (cid:15) Hartree n ( R, t ) = ˆ V ln ( R, t ) + (cid:90) dz ˆ W en ( z, R ) n e ( z, t ) , (16)where n e ( z, t ) is the electron density obtained from theelectronic dynamics. This approach does not involvethe complications of dealing with the conditional nu-clear wavefunction as in (cid:15) app but it fails to capture majoreffects of the correlated electron-nuclear dynamics dueto its mean-field nature. Finally, we point out that (cid:15) app as well as all the approximations (i)–(iii) represent theelectron-nuclear interaction in an electrostatic way only. III. CREI IN H +2 ISOTOPOLOGUES
We utilize a popular one-dimensional model of thesymmetric as well as asymmetric isotopologues of H +2 subject to a linearly polarized laser field. As the mo-tion of the nuclei and the electron in the true moleculeis assumed to be restricted to the direction of the po-larization axis of the laser field, the essential physicscan be captured by a 1D Hamiltonian featuring “soft-Coulomb” interactions [25, 26]: ˆ H ( t ) = − (cid:126) µ e ∂ ∂z − (cid:126) µ n ∂ ∂R − (cid:113) z − M M n R ) − (cid:113) z + M M n R ) + 1 √ .
03 + R + ˆ V l ( R, z, t ) (17)where R and z are the internuclear distance and theelectronic coordinate as measured from the nuclearcenter-of-mass, respectively. The nuclear effective mass H +2 HD + HT + (A +2 ) Hx + HX + D +2 x +2 +2 and itssymmetric and antisymmetric isotopologues in atomic units.Here, x(X) refers to the 10(100) times heavier fictitious isotopeof hydrogen and A +2 is another fictitious isotopologues withthe same effective nuclear mass as HT + . is denoted as µ n = M M M n while µ e = M n M n + m e is the elec-tronic reduced mass with M n = M + M . The laserfield, within the dipole approximation, is represented by ˆ V l ( R, z, t ) = eE ( t )( q e z − ζR ) , (18)where E ( t ) denotes the electric field amplitude and q e = 1 + m e M n + m e is the reduced charge and ζ =( M − M ) /M n is the mass-asymmetry parameter. Suchreduced-dimensional models have been shown to quali-tatively reproduce experimental results (see Ref. [27] forexample).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 iso-tope of hydrogen that is 10(100) times heavier than thatof H, while A +2 is another fictitious isotopologue withthe same effective nuclear mass as HT + (See Table. Iin which the nuclear effective mass of H +2 and isotopo-logues are given.) A. Comparison of ionization of isotopologues of H +2 withdifferent effective nuclear mass for a fixed field We apply a field of duration 50-cycles, wavelength λ = 800 nm ( ω = 0 . a.u.) and intensity I =2 . × W/cm , with a sine-squared pulse envelope(Fig. 1), to each of the H +2 isotopologues. First, we solvethe electron-nuclear TDSE for Hamiltonian of Eq. 17 nu-merically exactly, beginning in the initial ground-stateof the molecule.As the Hamiltonian (17) has been obtained after sep-arating off the center of mass motion and the originis set to be the nuclear center of mass the field cou-ples directly to the nuclear motion only in the asym-metric cases; in the symmetric case, nuclear motion isdriven purely by the electronic dynamics. This is alsoclear from Eq. 18 where in the symmetric case, ζ = 0 .In Figs. 2 we plot the ionization probability and aver-age internuclear distance , (cid:104) R (cid:105) , as a function of num-ber of cycles t/T where T denotes duration of one cy-cle ( T = 2 . fs) for H +2 , HD + , HT + , Hx + , HX + , A +2 ,D +2 and x +2 . Note that, the results are given in atomicunits, e = m e = h = 1 , throughout the article, unlessotherwise noted. The ionization probability is definedas I ( t ) = 1 − (cid:82) ∞−∞ dR (cid:82) z I − z I dz | Ψ( z, R, t ) | , with z I = 15 a.u. As it can be seen in Figs. 2, the ionization and aver-4
10 20 30 40 5000.10.20.30.40.5 I on i z a t i on H HD + HT + Hx + HX + < R > ( a . u . ) H D x FIG. 2. Ionization probability and average internuclear dis-tance , (cid:104) R (cid:105) , as a function of number of cycles t/T where T de-notes duration of one cycle ( T = 2 . fs) of the H +2 , HD + ,HT + , Hx + and HX + molecule on the left panels and A +2 , D +2 and x +2 on the right panels.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 asymmetric case in the leftpanels HD + , HT + , Hx + and symmetric case in the right panelsH +2 , D +2 and x +2 . age internuclear separation appear to be practically in-dependent of nuclear mass-symmetry properties: theyare identical for asymmetric case of HT + and symmetriccase of A +2 systems with the same effective nuclear mass µ n and very similar in the case of HX + and D +2 wherethe effective masses are only a little different. One canunderstand this from that fact that the ionization probesthe electron density in regions where z (cid:29) R , where thenuclear mass-dependence in the Hamiltonian Eq. 17 tolowest order in R/z is only via µ n , and the integratednature of the observable reduces its sensitivity to the de-tails of the distribution which is affected by the symme-try of the system.It is also observed that as the nuclear effective-massincreases, the ionization decreases, with a dependencethat appears to tend to /µ n for intermediate masses.Furthermore, the growth over time of the average inter-nuclear distance is less as the mass increases, suggest-ing that for larger masses the system reaches the criti-cal internuclear distance for CREI at later times whenthe field has decreased from its peak value significantly,consequently resulting in lower ionization. To investi-gate whether the CREI process is still relevant or notand shed more light on the dependence of the ioniza-tion on the effective nuclear mass we utilize the conceptof the time-resolved, R − resolved ionization probabil-ity defined as I ( R, t ) = (cid:82) z (cid:48) I dz | Ψ( z, R, t ) | , with (cid:82) z (cid:48) I = (cid:82) − z I −∞ + (cid:82) ∞ z I and z I = 15 a.u. [1, 27]. This quantity can berewritten in terms of the concepts of exact nuclear wave-function χ ( R, t ) and exact conditional electronic wave-function Φ R ( z, t ) introduced within the exact factoriza-tion framework, i. e., I ( R, t ) = | χ ( R, t ) | I cp , (19)where, I cp follows the usual expression of the ioniza-tion probability but using the exact conditional elec-tronic wavefunction Φ R ( z, t ) , i.e., I cp ( R, t ) = 1 − (cid:82) z I − z I dz | Φ R ( z, t ) | that is coupled to the exact nucleardynamics. Therefore, I ( R, t ) , which is the nuclear densityweighted conditional ionization probability is analogous to theionization probability calculated for a given nuclear configu-ration R in quasi-static picture. To analyze the dynamics of different cases, in Fig. 3 wehave plotted I ( R, t ) along with the time-dependent nu-clear density for different isotopologues. As seen in theupper three rows of Fig. 3 with increasing the effectivenuclear mass the peak of I ( R, t ) moves slightly to largertimes and its overall value decreases while remainingclose to the CREI region as predicted by the internuclearseparation of R c . The lower set of panels of Fig. 3 showsthe time-dependent nuclear density in each case. Forthe case of x +2 where only an exponentially small frac-tion of the nuclear density dissociates, the system doesnot reach the CREI regime and therefore the ionization isnegligible as seen from the R − resolved, t − resolved ion-ization probability. Notice the different scale of I ( R, t ) for the case of x +2 .5 +2 A +2 D +2 x +2 . × . × . × . × TABLE II. The optimized field intensity for different-mass sys-tems in the unit of W/cm . Note that the ionization rate computed at any fixed R would be identical for all the isotopologues. Their dif-ferent masses, however, lead to very different ionizationrates when the full electron-nuclear dynamics is consid-ered with the molecule beginning at equilibrium as wehave shown. Still, there is some validation to the orig-inal quasistatic CREI prediction that ionization is en-hanced for nuclear separations around R c , as indicatedby I ( R, t ) , however, a modification to the statement isneeded due to the spreading and splitting of the nuclearwavepacket: the enhancement occurs from electrons as-sociated with the part of the nuclear density that is in the R c region. Clearly, treating the nuclei as point particleswill not work even for significant nuclear mass (exceptin the large-mass limit) because of this. The questionthen arises whether accounting for the nuclear distribu-tion is enough to capture CREI accurately, for exampleusing one of the electrostatically-based approximationsof Section II , and, whether and how the errors decreasewith the nuclear mass. To this end, we next consider ad-justing the field strength so that the ionization is similarfor the different isotopologues. B. Comparison of potentials with similar ionization rates
As discussed in the previous section, the different iso-topologues of H +2 subject to the same field show signifi-cantly different degrees of ionization, since the different-mass systems reach the CREI region at different timeswith different probabilities. In particular, systems withlarger nuclear masses hardly reach the internuclear dis-tances for which CREI is expected to occur. Hence, tostudy the effect of the nuclear mass in the CREI regimefor different-mass systems, we adjust the field intensitysuch that the ionization probability(rate) remains closeto that of H +2 . As the asymmetric and symmetric iso-topologues with the same effective nuclear mass subjectto the same field give almost the same ionization proba-bility, average internuclear distance and I ( R, t ) (see sec-tion III A), from here on we focus only on the symmet-ric isotopologues of H +2 . The optimized field intensityfor different-mass systems is given in Table. II while theother field parameters are kept unchanged.The ionization probability of the symmetric isotopo-logues subject to the optimized field is depicted in Fig. 4.For isotopologues heavier than H +2 the ionization startsto 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 wehave shifted the time-dependent ionization probability
20 25 30 35 40t/T00.10.20.30.40.5 I on i za ti on H A D x FIG. 4. The field intensity is varied such that the ionizationyield(rate) of different systems remains similar to that of H +2 subject to a -cycle pulse of wavelength λ = 800 nm ( ω =0 . a.u.) and intensity I = 2 . × W/cm , with a sine-squared pulse envelope. For the isotopologues of H +2 , to viewbetter how the ionization yield(rate) compares to that of H +2 we have shifted the ionization by one cycle ( i.e. I((t-T)/T) isplotted).FIG. 5. The contour plot of the time-dependent nuclear densityfor H +2 (lower left panel), A +2 (upper left panel) , D +2 (lowerright panel), x +2 (upper right panel) subject to the optimizedfield given in table II. by one optical cycle, i.e. for isotopologues heavier thanH +2 , I ( t − T ) has been plotted.To analyze the mechanism of ionization in more de-tail, in Fig. 5 we plot the time-dependent nuclear den-sity for the various isotopologues. The nuclear densitybehaves similarly in all cases, except the heaviest case,namely x +2 . That is before the field reaches its maximumthe system becomes slightly ionized hence the nucleardensity slightly spreads. As the field reaches its maxi-mum a small fragment of the nuclear density starts tosplit off and dissociate from Coulomb explosion due toan increase in the ionization while the rest of the nu-clear density remains bounded and oscillates around theequilibrium position. As an appreciable amount of dis-sociating fragment reaches the internuclear separationfor CREI, the ionization gets enhanced. In the case ofx +2 , however, the nuclear dynamics exhibits a more clas-sical behavior, i. e. during the first half of the pulse,the heavy nuclei only spreads slightly around the equi-librium. The stronger field compared to the other casesenables the system to ionize initially (before the nucleardensity reaches the critical R ). This initial ionization6 IG. 6. The contour plot of I ( R, t ) for H +2 (lower left panel),A +2 (upper left panel) , D +2 (lower right panel), x +2 (upper rightpanel) subject to the optimized field given in table II. will be followed then by a Coulomb explosion towardthe middle of the field, after which the nuclear den-sity spreads further but it hardly splits. The remain-ing part of the electronic density undergoes enhancedionization as the nuclear density lies in the range of theinternuclear separation associated to CREI. In this case,the average internuclear distance almost coincides withthe peak of the nuclear wave packet and therefore forsystems with very large effective nuclear mass the stan-dard approximation to describe CREI, namely the qs ap-proximation is expected to perform better compared tothe lighter isotopes. The different nature of ionizationin large nuclear-mass systems, can also be seen from I ( R, t ) depicted in Fig. 6. While for not too heavy iso-topologues namely A +2 and D +2 , the I ( R, t ) has a similarstructure to H +2 , the I ( R, t ) corresponding to x +2 man-ifests a substantially distinctive structure compared tothe lighter isotopologues: the internuclear separationfor which ionization gets enhanced, shifts to smaller val-ues. In general, the I ( R, t ) shifts to smaller R for heavierisotopologues which could be related to the stronger op-timized field used [5, 6].Now, we investigate the performance of the vari-ous approximations introduced in Sec. II B for different-mass systems. In Fig. 7 the ionization probabilities cal-culated from propagating the electron on the exact, adi-abatic(adiab), approximate(app), exact-Hartree(ex-H),normalized Hartree(n-ex-H), quasi-static(qs) and self-consistent Hartree (SCH) for various symmetric isotopo-logues of H +2 are plotted.For all cases presented in Fig. 7, propagation on the qspotential gives rise to an underestimation of the ioniza-tion probability. The average (cid:104) R (cid:105) entering into the qs po-tential is always considerably less than the internuclearseparation of the dissociating fragment, and so doesn’taccess the CREI region that long during the duration ofthe pulse. The exact-Hartree (which compared to the qsapproximation accounts for the spreading and splittingof the nuclear wave packet) follows the qs results un-til the middle of propagation time then overtakes the qs
20 30 40t/T 00.20.40.6 I on i z a t i on I on i z a t i on
20 30 40t/Texactadiabappex-Hn-ex-Hqsadiab+accSCH D x H A FIG. 7. Ionization probabilities calculated from propagatingthe electron on the exact, adiabatic, approximate, Hartree like,normalized Hartree and quasi-static PES for different isomerof ionic Hydrogen. ionization and finally gives a rather large overestima-tion of the final ionization for all cases. This overestima-tion can be improved significantly by using the normal-ized Hartree approximation as defined in Eq. 15. In fact,the n-ex-H approximation performs clearly better thanother conventional approximations (qs, SCH, ex-H) forH +2 , A +2 and x +2 . Even the approximate potential thatdepends on the more rigorous and complicated conceptof conditional nuclear density rather than the nucleardensity that appears in the Hartree-like approximation,is considerably worse than n-ex-H. The SCH performsvery poorly with a negligible ionization probability fornot too heavy isomers, and only in the large mass limitshows some improvement in predicting the ionizationprobability, yielding a similar performance to the otherconventional approximations. This very poor perfor-mance of the SCH stems from the fact that the SCH isan uncorrelated approach and the splitting of nuclearwave-packet in the CREI regime cannot be accountedfor. Hence, for the cases where the CREI occurs due tothe splitting of a dissociating fragment of nuclear wavepacket the SCH cannot capture the right physics. In thelarge mass limit, however, the CREI mechanism relieson the spreading of the nuclear wave packet rather thansplitting (see Fig. 5) , therefore the SCH approximationis not as poor. In general, the approximations do not perform better with increasing nuclear mass: dynami-cal correlation effects that are missing in all the poten-tials shown (apart from adiab and adiab+acc) dependmore critically on the nuclear velocities relative to theelectronic, and adjusting the field to get similar ioniza-tion results in these being similar (See also Sec IV A).The two curves in Fig. 7 remaining to be discussed7re adiab and adiab+acc, which are components arisingfrom the marginal decomposition (Eq. 12). The prop-agation of the electron on the adiabatic potential alone(adiab), gives rise to the complete ionization of the sys-tem rather early for all isotopologues, while the prop-agation on the potential composed of the adiabatic plusacceleration (adiab + acc) term yields an ionization prob-ability in a very good agreement with the exact resultswith the velocity term adding only a small correction,expect in the large mass limit. In the following sectionby studying the structures of the terms of the marginaldecomposition we try to shed some light on the reasonbehind the performance of these potentials. C. The structure of the dynamical electron-nuclear termsin e -TDPES: marginal decomposition We concluded the previous section by brieflydiscussing the electron dynamics on differentterms/combinations-of-terms of the marginal de-composition of e -TDPES (Eq. 12). In order to betterunderstand the outcomes, in this section we studystructures of different terms of marginal decomposi-tions of the e -TDPES for the two radically differentcases of H +2 and x +2 . We refer the readers to Ref. [1]for a discussion on the components of the conditionaldecomposition of e -TDPES. For the sake of simplifyingthe discussion, here, we divide the electronic coordinateinto two regions: ”inner-region” that refers to the regionwith | z | < a.u. and ”outer-region” that describes therest of the axis for which | z | > a.u.
1. H +2 case In Fig. 8, we present the terms/combinations-of-termsof marginal decomposition of the e -TDPES for the caseof H +2 at four different snapshots in time. Adiabatic term : the adiabatic term (first term inEq. 12) is the main constituent of the e -TDPES when itis decomposed according to Eq. 12. In this case, as it isseen in Fig. 8 (upper left panel), it initially follows theexact e -TDPES in the inner-region and to some extent inthe outer-region while it shows a different behavior inthe asymptotic region (deep in the outer-region). How-ever, as initially there is a negligibly small amount ofelectronic density far from the inner-region, this devi-ation does not influence the dynamics significantly. Asthe field intensity increases, the adiabatic potential startsto deviate from the exact potential, both in the inner-region and outer-region. It can be seen in Fig. 8 (lowerleft panel) that the shape of the adiabatic potential in theinner-region (especially the up-field part) and its (aver-age) slope in the outer-region differ significantly fromthe exact e -TDPES. The (average) slope of the exact po-tential follows the slope of the field in the outer-regionas expected. This feature together with the up-field part the exact potential are mainly responsible for controllingthe ionization from the up-field direction. Indeed theover-ionization corresponding to the propagation on theadiabatic potential discussed in the previous section (seethe lower left panel of Fig. 7), which starts around the th optical cycle is associated to the lack of the asymp-totic slope, and the error in the shape of the up-fieldwell. In particular, the wrong asymptotic behavior of theadiabatic potential in the outer region, allows for ioniza-tion from both sides (up-field and down-field) in eachhalf cycle, resulting in a huge overestimation of ioniza-tion probability. However, an important feature of theexact e -TDPES is captured in the adiabatic potential: thedevelopment of the four wells representing the branch-ing of the nuclear wavefunction in the inner-region ofthe potential which is associated to the correlation be-tween the electronic and nuclear motions. Towards theend of the dynamics, where the field intensity is smallagain the adiabatic potential follows the exact e -TDPESclosely as can be seen in Fig. 8 (lower-panel, right). Velocity term : is the second term in Eq. 12 and itsoverall contribution to the e -TDPES, in this case, is smallespecially 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 enoughto capture the essential features appearing in the exact e -TDPES. Acceleration term : is the last term in Eq. 12 that isinitially very small in the inner-region but as the ion-ization sets in, the addition of this term to the adiabaticpotential significantly improves the shape of the poten-tial, particularly when the field approaches the peak-intensity as evident in Fig. 8 (upper-right and lower-leftpanels) in the inner-region as well as the outer-region(see the ”adiab+acc”). The latter is due to the muchbetter (average) asymptotic behavior of this term in theouter-region compared to the adiabatic term. On theother hand, it also improves the up-field/down-fieldwell in the inner-region when added to the adiabatic po-tential. As a result propagating on the combination ofadiabatic and acceleration terms leads to an ionizationprobability in a good agreement with the exact resultsas it is shown in Fig. 7 (lower left panel).As the structures of the adiabatic, velocity, and accel-eration terms in the case of A +2 and D +2 are very similarto H +2 , we do not address them here.
2. x +2 case In Fig. 9, the terms/combinations-of-terms ofmarginal decomposition of the e -TDPES for the case ofx +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 tofollow the shape of the exact potential. Again, this fail-ure, in particular the wrong average slope of the adia-8 IG. 8. The exact e -TDPES and its component from themarginal decompositions at various snapshots in time of H +2 . batic potential in the asymptotic region, results in a hugeoverestimation of the ionization probability. Velocity term : the velocity term plays a crucial role inthis case as is seen in Fig. 9 (lower-panel, left). In par-ticular, as the field approaches its maximum intensity,it exhibits more structure in the inner region and con-tributes more significantly to the overall shape of the e -TDPES. Specifically, it exhibits a valley in the center thatis completely absent in the case of H +2 and lowers theinteratomic barrier when added to the adiab+acc poten-tial for most of the times between t ≈ T and t ≈ T in which most of the ionization happens. The moredominant role of the velocity term in this case, could beattributed to the increase of ponderomotive (wiggling)motion of the electron driven by a stronger external laserfield. Acceleration term : this term has a correct asymp-totic behavior, similar to the case of H +2 but leads to awrong estimate for the inner-barrier (see Fig. 9, lower-left panel,) when added to the adiabatic term. Thewrong estimation of the inner-barrier is exaggerated af-ter the field reaches its maximum intensity (especiallydue to the higher interatomic barrier around zero be-tween t ≈ T and t ≈ T ) , leading to an inaccurateionization probability after t ≈ T , as seen in Fig. 7(uper right panel). In ref. [15] the importance of theinner-barrier (peaks) and up-field well (steps) to achievethe correct electron-localization has been shown. IV. QUASI-CLASSICAL ANALYSIS
The equations for the electronic wavefunction andconditional nuclear wavefunction cannot be solved ex-actly for systems of more than a few degrees of freedom,just as solving the full molecular Schr ¨odinger equa-tion exactly for those systems is not possible. In manycases a quasiclassical treatment of the nuclear dynam-
FIG. 9. The exact e -TDPES and its component from themarginal decompositions at various snapshots in time of x +2 . ics should be sensible; by quasiclassical, we mean an ensemble of classical trajectories, weighted according tothe initial distribution, is evolved for the nuclei, ratherthan a single trajectory. This would allow the possibil-ity to capture spreading and branching of the nucleardistribution. Such a procedure within the exact factor-ization in its reverse flavor involves taking the classi-cal limit of the conditional nuclear wavefunction whichdoes not satisfy an equation of Schr ¨odinger form. There-fore, it differs from quasiclassical treatments of usualSchr ¨odinger equations that have also been discussedin various forms for the marginal nuclear wavefunc-tion within the ”direct factorization” framework [28–33]. Here we begin taking the classical limit of the condi-tional nuclear wavefunction by representing it in polarform: χ z ( R, t ) = A z ( R, t ) e iS z ( R,t ) / (cid:126) (20)where A z ( R, t ) and S z ( R, t ) are both real functions. In-serting this into the equation of motion for χ z ( R, t ) ,Eq. 4, and sorting the terms in orders of (cid:126) , we find to O ( (cid:126) ) : µ n (cid:18) ∂S z ∂R (cid:19) + 12 µ e (cid:18) ∂S z ∂z (cid:19) + V ( z, R, t ) − (cid:15) e ( z, t )+ ˜ p e ( z, t ) µ e ∂S z ∂z + (cid:18) ∂S z ∂t (cid:19) = 0 (21)where we define V ( z, R, t ) = W nn ( R ) + W en ( z, R ) + V l ( z, R, t ) , and ˜ p e ( z, t ) = − i (cid:126) ∂ Φ( z, t ) /∂z (see shortly formore on this term). Similarly, keeping only O ( (cid:126) ) termsin Eq. 7 for the e -TDPES, we find (cid:15) cle ( z, t ) = (cid:90) dR | χ z ( R, t ) | (cid:16) µ n (cid:18) ∂S z ∂R (cid:19) + 12 µ e (cid:18) ∂S z ∂z (cid:19) + V ( z, R, t ) + ∂S z ∂t (cid:17) . (22)9he terms on the right-hand-side correspond toclassically evaluating T n ( z, t ) , K cond e ( z, t ) , (cid:15) app e ( z, t ) and (cid:15) gd e ( z, t ) respectively. Notice that the electron-nuclearcoupling operator U coupen has a classical counterpart, asit contributes already at zeroth-order in (cid:126) with the term µ e (cid:0) ∂S z ∂z (cid:1) in both Eq (21) and (22).Eq. 21 would be a standard Hamilton-Jacobi equationof the form H ( q , ∂S z ∂ q , t ) + ∂S z /∂t = 0 (where H ( q , p , t ) is the Hamiltonian function), for the action S z ( R, t ) oftwo particles, one of mass µ n and the other of mass µ e in a potential V ( z, R, t ) − (cid:15) e ( z, t ) , if the second-last termon the left-hand-side was not present. It is perhaps notsurprising that we do not retrieve a standard Hamilton-Jacobi equation, given that the equation for χ z is nota TDSE. Although an (cid:126) multiplies this term, it does infact contribute in the classical limit, as will be discussedshortly. Still, classical Newton-like equations can be de-rived from Eq. 21 by defining the velocity fields, u nz ( R, t ) = 1 µ n ∂S z ( R, t ) ∂R , and u ez ( R, t ) = 1 µ e ∂S z ( R, t ) ∂z (23)Then, taking ∂/∂R of Eq. 21 yields µ n du nz ( R, t ) dt = − ∂∂R ( V ( z, R, t ) + ˜ p e ( z, t ) u ez ( R, t )) (24)and µ e du ez ( R, t ) dt = − ∂∂z ( V ( z, R, t ) − (cid:15) e ( z, t ) + ˜ p e ( z, t ) u ez ( R, t )) (25)where d/dt = ∂/∂t + u ez ∂/∂z + u nz ∂/∂R is the time-derivative in the Lagrangian frame defined by the ve-locities. A. Quasiclassical analysis of the terms in the e -TDPES Whether the equations above, together with the so-lution of the TDSE Eq. 3 could form the basis of amixed quantum-classical method remains for futurework. Here, instead we consider how a quasiclassi-cal analysis of the entire coupled electron-nuclear sys-tem can shed light on the structure and nuclear-mass-dependence of the terms in the e -TDPES. Many aspectsof electron dynamics in strong fields can be treatedclassically, especially when tunneling and quantizationare not driving the primary physics. Indeed, Ref. [34]showed that classical trajectory calculations reproducethe essential features of CREI for both cases of fixed andmoving nuclei.To this end, we first evaluate ˜ p e ( z, t ) in Eq. 21to its lowest-order in (cid:126) . Inserting Φ( z, t ) = a ( z, t ) exp( is ( z, t )) / (cid:126) , analogous to Eq. 20, into the TDSEEq. 3 for the marginal wavefunction Φ( z, t ) , and keepingonly the term that is lowest-order in (cid:126) , gives µ e (cid:18) ∂s ( z, t ) ∂z (cid:19) + (cid:15) e ( z, t ) + ∂s∂t = 0 . (26) This is the Hamilton-Jacobi equation for the action s ( z, t ) of a classical electron evolving in Hamiltonian ˜ T e + (cid:15) e ( z, t ) . Here ˜ T e = µ e (cid:0) ∂s∂z (cid:1) is the electronic kinetic en-ergy associated with the Hamiltonian in the TDSE forthe electronic wavefunction, Eq. 3. It is important tonote that this is not the same as the true electronic ki-netic energy (see Eq. 69 in Ref. [17]). The term in Eq. 21, − i (cid:126) ∂ Φ( z, t ) /∂z thus becomes ˜ p e ≡ (cid:113) µ e ˜ T e in the classi-cal limit.Next, we note that Eq. 21 for S z is consistentwith the corresponding classical limits taken for theelectronic wavefunction and the full molecular wave-function: 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 ) (27)where, in the classical limit, s ( z, t ) satisfies Eq. 26 and S ( z, R, t ) satisfies µ n (cid:18) ∂S∂R (cid:19) + 12 µ e (cid:18) ∂S∂z (cid:19) + V ( z, R, t ) + (cid:18) ∂S∂t (cid:19) = 0 (28)It is straightforward to check by substitution thatEqs. 26,28,21 and 27 are consistent with each other.Eqs. (26) and (28) are both standard Hamilton-Jacobiequations, so, for the classical trajectories that they de-scribe, we readily identify the following: ∂S/∂R = P n = µ n ˙ R, 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 + (cid:15) e . Now we imagine an ensemble of classical particles withinitial position and momenta distributed according tothe initial wavefunction. We consider evaluating thevarious contributions to the classical e -TDPES of Eq. 22from these classical trajectories. The integral over R inthat equation then becomes a sum over classical trajecto-ries 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 .10 IG. 10. The kinetic term, T n ( z, t ) , at various snapshots intime for H +2 and its isotopologues. For the first term in Eq. 22, T n ( z, t ) , T cln ( z, t ) = (cid:90) dR | χ z ( R, t ) | µ n (cid:18) ∂S z ∂R (cid:19) = (cid:90) dR | χ z ( R, t ) | µ n R = µ n N trajz N trajz (cid:88) I z ˙ R I z (29)where we noted that ∂S z ( R, t ) /∂R = ∂S ( z, R, t ) /∂R (from Eq. 27), and replaced the integral with a sum overall N trajz trajectories I z that reach z ( t ) = z at time t ; ˙ R I z is the nuclear velocity along the I z -th trajectory. First, letus see what this says about the overall structure of thisterm. During the dynamics, we observed that part ofthe nuclear density oscillates around equilibrium, mov-ing slower than the dissociating part of the density. Thismeans that the trajectories for small z are mostly associ-ated with more slowly-moving nuclear dynamics nearequilibrium, while those with larger z are in the pro-cess of ionizing, and so are associated with faster nu-clear speeds due to the net ionic charge they have leftbehind. Hence, for small z , there are more trajectorieswith smaller ˙ R ’s than for large z , and so we expect apotential that rises as z gets larger. This is indeed whatwe see in Fig. 10. Consider first the black curve, H +2 .At early times, the middle region of T n ( z, t ) , where thebulk of the electron density is, is rather flat and takes itssmallest value: Initially, the conditional nuclear proba-bility has very little z -dependence in the region wherethe electron density is appreciable, i. e., at small timesthe classical positions of nuclear trajectories for z in theregion of appreciable electronic density are identical, re-flecting the broadness of the electron distribution rela-tive to the nuclei. The small nuclear kinetic energy atthe early times near equilibrium is essentially from zero-point motion. As the electrons begin to gain energy fromthe field and move out from the tails of the electronic distribution, the T n ( z, t ) grows accordingly since thesetails are associated with trajectories where the nuclei arebeginning to move apart under Coulomb repulsion. Thesmall middle flat region of T n ( z, t ) gets narrower as ion-ization take away more of the electron density. As thefinal distribution sets in, T n ( z, t ) , one can identify threeregions, reflecting the electronic distribution and the as-sociated nuclear kinetic energies: an inner region asso-ciated with the density remaining near equilibrium, ris-ing to the intermediate region where the shoulders ofthe electron density lie associated with the dissociatednuclear wavepacket, and then rising further out to thetails of the electron density which continue to oscillatein the field.Regarding now the µ n -dependence: it might appearfrom Eq. 29 that T n ( z, t ) grows as the mass µ n increases,in the cases where the field is adjusted so that the aver-age internuclear distance and speed remains the same.This is in fact not the case; in fact, the term remainsabout the same size, as the nuclear mass increases as isevident in Fig. 10. This is because, for larger µ n , the nu-clear density tends to split less (see Fig. 6), i.e. for larger µ n , a larger proportion of the nuclear density moves outto larger separations rather than remaining near the ori-gin with small ˙ R . Yet about the same ˙ (cid:104) R (cid:105) is maintainedby design by the different field strengths, which meansthat the fastest trajectories in the distribution are slowerfor larger masses µ n than for smaller µ n . This wouldcontribute to a smaller rise of T n ( z, t ) as one moves tolarger z , as mass increases, compensating the larger µ n .In Fig. 10 we plot T n ( z, t ) for the different isotopologuesat various snapshots in time. We see that the overall sizeof this term does not vary much with mass except forx +2 . The exception is for the large mass x +2 where thereis comparatively weak z -dependence of T n ( z, t ) . Thisis because the nuclear distribution moves largely as awhole, showing less difference between nuclear speedsat different parts of the distribution.The corresponding analysis for the other three termsof Eq. 22 is unfortunately not as straightforward, be-cause the evolution of the relevant classical trajectoriesthemselves depend on the potential (cid:15) e ( z, t ) , which is thevery object we wish to analyse! For example, the secondterm, K cond e ( z, t ) , requires ∂S z /∂z = ( √ T e − (cid:112) ˜ T e ) but ˜ T e is the kinetic energy of an electron in the potential (cid:15) e .Future work along these lines will likely require further,possibly iterative, approximations to the potential, withthe dual goal of analysing the structure of the terms anddeveloping mixed quantum-classical schemes. V. CONCLUSIONS
The exact potential driving the electron dynamicswithin the exact factorization of the full electron-nuclearwavefunction formally exactly accounts for the couplingto the nuclear subsystem as well as coupling to exter-11al fields. In order to develop adequate approxima-tions to treat electronic non-adiabatic processes, it is im-portant to understand the structure of the exact poten-tial and pinpoint and analyse its important features invarious situations. In this work, we have presenteda detailed investigation of the exact potential drivingthe electron dynamics in isotopologues of H +2 undergo-ing CREI and studied the dependence of the correlatedelectron-nuclear dynamics on the nuclear mass.The elaborate concept of time-resolved, R -resolvedionization probability, I ( R, t ) , provides an extremelyuseful tool, indicating the internuclear separations as-sociated with ionization. The concept is analogous tothe ionization probability calculated within the qua-sistatic picture (with clamped nuclei) but can be usedunambiguously for the fully dynamical case. For thelaser parameters and different isotopologues of H +2 dis-cussed in this work, I ( R, t ) , exhibits only one peak.This is in agreement with most of the experimental find-ings [35, 36] and in contrast to the predictions of CREIbased on the standard clamped-nuclei quasistatic pic-ture. We suggest that I ( R, t ) would be a very usefultool in future studies of CREI, for example to resolve thelaser parameters for which double-peak structure couldin fact be observed as in Ref. [37].We found that for fixed laser parameters, the ioniza-tion yields rapidly decrease as a function of the massof the isotopologue because less of the nuclear densitymakes it to the CREI region during the time the laser ison. For all the isotopologues, I ( R, t ) indicated that theionization is nevertheless dominated by the fraction ofelectrons associated with the nuclear density in the CREIregion defined qualitatively by the original quasistaticargument. This implies that treating the nuclei as clas-sical point particles will not work; one needs to accountfor their distribution, which, away from the large-masslimit, displays a branched structure with part of the dis-tribution oscillating near equilibrium separation whilepart of it, associated with the CREI electrons, dissoci-ates.Still, going beyond the quasistatic point nuclei pictureand accounting for the nuclear distribution in an elec-trostatic way (as in Hartree-type approximations) is notadequate in capturing the dynamics accurately. The im-portance of going beyond an electrostatic description ofelectron-nuclear correlation is evident from our studiesof how different approximate electronic potentials per-form in describing the CREI for isotopologues of H +2 .We have shown that, for the laser parameters used inthis work, one must go beyond any purely electrostatictreatment of electron-nuclear correlation, and includetruly dynamical aspects of the nuclear distribution andits coupling to the electronic system to get a good pre-diction of the ionization. In determining errors fromconventional approximations and deviations of approx- imate potentials from the exact potential, what are moreimportant than the nuclear-to-electronic mass ratios, arethe nuclear velocities in the wavepacket.There are many different approaches to developingapproximate methods based on the exact factorization.For example, one may consider the relative importanceof different components of the exact potential when de-composed in terms of the conditional wavefunction asin Ref. [1] (there, it was found that generally all termswere important). One may consider also approxima-tions based on the marginal decomposition presentedhere, where we found that the ”adiab+acc” componentof the marginal decomposition describes the electronicdynamics accurately in all cases we studied here withthe exception of a fictitious isotopologue with an effec-tive nuclear mass 10 times larger than H +2 . Further in-vestigation of this potential may lead to developmentof an adequate approximation for practical purposes ca-pable of describing the ionization dynamics accurately.Another approach is to develop quasiclassical approx-imations, and here we have sketched out a semiclassi-cal derivation of the electronic and conditional nuclearequations of the exact factorization in its reverse formand analysed the structure of the nuclear kinetic term ofthe exact potential semiclassically.This work has highlighted the effect of the complexinterplay of electronic and nuclear dynamics in strongfield enhanced ionization processes by demonstratingthe large differences in the conventional potentials andthe dynamics they cause with the exact potential drivingthe electron for systems of varying nuclear mass. Theexplorations of the details of the potential and approxi-mation methods lay the ground-work for future devel-opment of accurate methods for coupled electron-iondynamics in non-perturbative fields. Whether the dy-namical electron-nuclear effects play as crucial a role forpolyatomic molecules, and the scaling of these effectswith respect to the number of electrons and number ofnuclei, is also an important avenue for future research ACKNOWLEDGMENTS
We acknowledge support from the European Re-search Council(ERC-2015-AdG-694097), Grupos Con-solidados (IT578-13), and the European Union’s Hori-zon 2020 Research and Innovation programme undergrant agreement no. 676580. A.K. and A.A acknowl-edge funding from the European Uninos Horizon 2020research and innovation programme under the MarieSklodowska-Curie grant agreement no. 704218 and702406, respectively. N.T.M thanks the National ScienceFoundation, grant CHE-1566197 for support.12
1] E. Khosravi, A. Abedi, and N. T. Maitra, Phys. Rev. Lett. , 263002 (2015).[2] T. Zuo, S. Chelkowski, and A. D. Bandrauk, Phys. Rev. A , 3837 (1993).[3] T. Zuo and A. D. Bandrauk, Phys. Rev. A , R2511 (1995).[4] T. Seideman, M. Y. Ivanov, and P. B. Corkum, Phys. Rev.Lett. , 2819 (1995).[5] S. Chelkowski, C. Foisy, and A. D. Bandrauk, Phys. Rev.A , 1176 (1998).[6] S. Chelkowski and A. Bandrauk, J. Phys. B: Atomic,Molecular and Optical Physics , L723 (1995).[7] S. Chelkowski, A. Conjusteau, T. Zuo, and A. D. Ban-drauk, Phys. Rev. A , 3235 (1996).[8] H. Yu, T. Zuo, and A. D. Bandrauk, Phys. Rev. A , 3290(1996).[9] A. D. Bandrauk and F. L´egar´e, inProgress in Ultrafast Intense Laser Science VIII (Springer,2012) pp. 29–46.[10] N. Takemoto and A. Becker, Phys. Rev. Lett. , 203004(2010).[11] N. Takemoto and A. Becker, Phys. Rev. A , 023401(2011).[12] C. Beylerian, S. Saugout, and C. Cornaggia, Journal ofPhysics B: Atomic, Molecular and Optical Physics ,L105 (2006).[13] I. Bocharova, R. Karimi, E. F. Penka, J.-P. Brichta, P. Las-sonde, X. Fu, J.-C. Kieffer, A. D. Bandrauk, I. Litvinyuk,J. Sanderson, et al., Phys. Rev. Lett. , 063201 (2011).[14] F. L´egar´e, I. V. Litvinyuk, P. W. Dooley, F. Qu´er´e, A. D.Bandrauk, D. M. Villeneuve, and P. B. Corkum, Phys. Rev.Lett. , 093002 (2003).[15] Y. Suzuki, A. Abedi, N. T. Maitra, K. Yamashita, andE. K. U. Gross, Phys. Rev. A , 040501 (2014).[16] A. Abedi, N. T. Maitra, and E. K. U. Gross, Phys. Rev.Lett. , 123002 (2010).[17] A. Abedi, N. T. Maitra, and E. K. U. Gross, J. Chem. Phys. , 22A530 (2012).[18] A. Abedi, N. T. Maitra, and E. K. U. Gross, J. Chem. Phys , 087102 (2013).[19] A. Abedi, F. Agostini, Y. Suzuki, and E. K. U. Gross, Phys.Rev. Lett. , 263001 (2013). [20] R. Requist and E. K. U. Gross, Phys. Rev. Lett. , 193001(2016).[21] H. Yu, T. Zuo, and A. D. Bandrauk, J. Phys. B: Atomic,Molecular and Optical Physics , 1533 (1998).[22] M. Uhlmann, F. Großmann, T. Kunert, and R. Schmidt,Phys. Lett. A , 417 (2007).[23] J. C. Tully, Faraday Discuss. , 407 (1998).[24] The Hartree potential (14) reduces to the qs expres-sion (13) in the limit of very localized wave-packets cen-tered around ¯ R .[25] J. Javanainen, J. H. Eberly, and Q. Su, Phys. Rev. A ,3430 (1988).[26] T. Kreibich, M. Lein, V. Engel, and E. K. U. Gross, Phys.Rev. Lett. , 103901 (2001).[27] K. C. Kulander, F. H. Mies, and K. J. Schafer, Phys. Rev. A , 2562 (1996).[28] A. Abedi, F. Agostini, and E. K. U. Gross, Europhys. Lett. , 33001 (2014).[29] S. K. Min, F. Agostini, and E. K. U. Gross, Phys. Rev. Lett. , 073001 (2015).[30] F. Agostini, S. K. Min, A. Abedi, and E. K. Gross, Journalof chemical theory and computation (2016).[31] Y. Suzuki, A. Abedi, N. Maitra, and E. K. U. Gross, Phys.Chem. Chem. Phys. , (2015).[32] F. Agostini, A. Abedi, Y. Suzuki, S. K. Min, N. T. Maitra,and E. K. U. Gross, J. Chem. Phys. , 084303 (2015).[33] Y. Suzuki and K. Watanabe, Phys. Rev. A , 032517(2016).[34] D. M. Villeneuve, M. Y. Ivanov, and P. B. Corkum, Phys.Rev. A , 736 (1996).[35] T. Ergler, A. Rudenko, B. Feuerstein, K. Zrost, C. D.Schr¨oter, R. Moshammer, and J. Ullrich, Phys. Rev. Lett. , 093001 (2005).[36] 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. Os-ipov, S. Ghimire, Z. Chang, and C. L. Cocke, Phys. Rev. A , 063419 (2008).[37] H. Xu, F. He, D. Kielpinski, R. Sang, and I. Litvinyuk,Scientific reports (2015).(2015).