Exact Factorization Adventures: A Promising Approach for Non-bound States

Modeling the dynamics of non-bound states in molecules requires an accurate description of how electronic motion affects nuclear motion and vice-versa. The exact factorization (XF) approach offers a unique perspective, in that it provides potentials that act on the nuclear subsystem or electronic subsystem, which contain the effects of the coupling to the other subsystem in an exact way. We briefly review the various applications of the XF idea in different realms, and how features of these potentials aid in the interpretation of two different laser-driven dissociation mechanisms. We present a detailed study of the different ways the coupling terms in recently-developed XF-based mixed quantum-classical approximations are evaluated, where either truly coupled trajectories, or auxiliary trajectories that mimic the coupling are used, and discuss their effect in both a surface-hopping framework as well as the rigorously-derived coupled-trajectory mixed quantum-classical approach.


I. INTRODUCTION
A challenge in the theoretical simulation of molecules driven out of their ground-state by a strong field is how to account for, and interpret, the coupling between the electrons and ions.
If a molecule was only gently perturbed at frequencies below electronic excitations, the Born-Oppenheimer (BO) picture can still be useful: the disparity in the time-scales of nuclear and electronic motion allows for the nuclei to evolve on a single BO surface for much of the time, and regions where this breaks down tend to be somewhat localized in space and involving the interaction of only a small number of BO states. In those regions, a fully quantum treatment of the molecule based on the BO picture would involve propagating nuclear wavefunctions projected on the different surfaces, including coupling terms between them. If the nuclear degrees of freedom were to be treated classically, the effective force driving their motion would deviate from that of a single surface, or even a mean surface: different parts of the underlying nuclear wavepacket would be associated with different electronic surfaces, even within the same spatial region. The problems with capturing this electron-nuclear correlation in such mixed quantum-classical methods are well-known [1,2]. But with strong fields, or with high-frequency fields, the problem is even more challenging: the nuclear wavepacket straddles many BO states, a continuum beyond the ionization threshold. Fielddressed surfaces such as quasi-static [3][4][5] or Floquet [6][7][8] can be instructive in this regard; although both serve equally well as a basis for the nuclear dynamics, as does the BO basis, approximations inherent in mixed quantum-classical approaches may perform better in one or the other, for physical reasons or computational reasons, but for general field parameters, the jury is out. Further, the observables of interest depend on the application: for example, we may want to correlate electron ionization and nuclear dissociation channels to get a full understanding of the dynamics [9], or, we may be primarily interested in the behavior of the ionizing electron itself, with a time-dependent potential driving the electron helping to guide our understanding. As intuitive and instrumental the BO picture has been, it does not provide such a potential, nor is it convenient in strong fields when many surfaces get involved.
An alternative, but equally fundamental, picture is provided by the exact factorization (XF) approach [10][11][12][13][14][15][16][17][18]. XF re-casts the dynamics of interacting quantum subsystems in a way that provides unique potentials that drive the motion of each subsystem and account for complete coupling to others. Applied to the electron-nuclear case, the exact molecular wavefunction Ψ(r, R, t) is factored into a single correlated product of a marginal factor and a conditional factor. Choosing the marginal to be a function of R, one obtains a time-dependent Schrödinger equation (TDSE) for the nuclear wavefunction, in which the potentials contain the complete effect of coupling to the electrons [16]. Instead, choosing the marginal to be a function of r only, a TDSE is obtained for the electronic wavefunction, in which the potentials contain the complete effect of coupling to the nuclei [19]. These potentials drive the nuclear and electronic subsystem respectively, and have been instructive in interpreting correlated electron-ion dynamics, including for dissociation [16] and ionization processes [20] in strong fields. To be a practical approach, approximations must be made, and recently two mixed quantum-classical (MQC) approaches have been derived, in which a key ingredient in the electron-nuclear coupling terms is the so-called nuclear quantum momentum. The quantum momentum measures the spatial variation of the nuclear density during non-adiabatic dynamics, and couples strongly to the electronic evolution giving large corrections to mean-field approximations. Different schemes have proposed different ways to calculate this term, using either the support of coupled trajectories or of auxiliary trajectories to track the delocalization of the nuclear density, but until now a detailed comparison has not been made. Here we compare the use of auxiliary trajectories within an XF-based surface-hopping framework (SHXF) [21,22] to using coupled trajectories within the surfacehopping procedure (CTSH) [23] as well as within the original coupled-trajectory mixed quantum-classical approach (CTMQC) [24][25][26] that was derived directly from the XF. We further study the impact of a modified definition of the quantum momentum in the coupledtrajectories calculations CTSH and CTMQC, that ensures the physical condition of zero population transfer in regions of no non-adiabatic coupling.
The paper is organized as follows. In Sec. II, we review the formalism of the XF approach, and briefly meander down various avenues in which the XF idea has been extended and explored. Sec. III demonstrates its usefulness for the dynamics of non-bound systems, briefly reviewing the results of Refs. [16,17] for dissociation and Refs. [20,27] for ionization, in a one-dimensional model of H + 2 . The most practical impact of the XF approach is the development of rigorous MQC methods, and in Sec. IV we outline the CTMQC, CTSH, and SHXF schemes, before delving into the different ways the quantum momentum is computed in different implementations of these schemes in Sec. V. We offer some conclusions in Sec. VI.

II. THE XF APPROACH IN A NUTSHELL
Consider the time evolution of a system of coupled quantum subsystems, possibly subject to some external fields:Ĥ Ψ(q, Q, . . . , t) = i∂ t Ψ(q, Q, . . . , t) (Atomic units, e 2 = = m e = 1, are used throughout this article unless otherwise stated).
These could be different types of particles, for example with q representing electronic coordinates and Q nuclear coordinates, and perhaps another set of coordinates for photons.
Or, for example, we could have a system of the same type of particle, say electrons, but with q being occupation numbers of some set of orbitals, e.g. weakly-correlated and Q being occupation numbers of strongly correlated orbitals. The problem is quite general, and in the XF the full wavefunction is factored into a marginal term and a conditional term [10,[15][16][17]: Here ...|... q represents an inner product over all q variables only. The factorization Eq. (2) is exact and unique up to a Q-and t-dependent phase, where e iθ(Q,t) multiplies χ while e −iθ(Q,t) multiplies Φ Q , yielding the same product. For Hamiltonians that have the form where V (q, Q, t) is a purely multiplicative operator in q, Q, applying the Dirac-Frenkel variational principle shows that the marginal factor χ(Q, t) satisfies a TDSE with a scalar and vector potential that depend on Φ Q (q, t) while the conditional factor satisfies a more complicated equation, with coupling terms dependent on χ(Q, t) [16,17,28]. In Eq. (3) we used the symbols M I and m i to indicate the masses -in atomic units -of the two sets of particles with positions Q and q, respectively. Specifically, for the electron-nuclear problem, with Q = R the nuclear coordinates, and q = r the electronic coordinates, and the partial normalization condition dr|Φ R (r, t)| 2 = 1, the wavefunctions Φ R (r, t) and χ(R, t) satisfy: with the electronic HamiltonianĤ el defined aŝ and with electron-nuclear coupling term Nuclear masses in atomic units are indicated with the symbol M ν , with ν labelling the nuclei.
HereĤ BO is the usual BO Hamiltonian (the sum of the electron kinetic energy, electronelectron, nuclear-nuclear, and electron-nuclear interaction operators). TheV ext n (R, t) or V ext e (r, t) potential is any externally applied potential, such as a laser field, acting on the nuclei or the electrons, respectively. We denote the time-dependent scalar potential appearing in Eqs. (4) and (5) as the time-dependent potential energy surface (TDPES). Together with the time-dependent vector potential appearing in Eqs. (5) and (7) and the electron-nuclear coupling termÛ en , these three terms embody the exact electronnuclear coupling. Refs. [16,17] showed Eqs. (4-9) are form-invariant under a gauge-like transformation that multiplies Φ R by an R, t-dependent phase e iθ(R,t) , and χ by its inverse; the potentials transform correspondingly as A → A + ∇θ(R, t) and (R, t) → (R, t) + ∂ t θ(R, t). Further, it was shown that we can identify χ(R, t) as the nuclear wavefunction, in the sense that its modulus-square gives the N n -body nuclear density, and the exact N n -body nuclear current-density can be extracted in the usual way from its phase together with the vector potential [16].
These equations emphasize the distinction from the far simpler BO approximation, that also involves a single correlated product, but with a time-independent electronic wavefunction and much simpler potential terms. The TDSE-form of the nuclear equation means that the potentials appearing in it can be directly interpreted as driving the nuclear motion, fully incorporating the effect of non-adiabatic and time-dependent coupling to the electronic system.
Notice that the potentials appearing in the nuclear equation, (R, t) and A ν (R, t), involve the electronic wavefunction Φ R (r, t), while the electron-nuclear coupling term appearing in the electronic equation,Û en , involves the nuclear wavefunction χ: thus the evolution of each factor is intimately correlated with the other.

A. XF Adventures in Different Realms
Although originally presented for coupled electron-nuclear dynamics, the XF approach applies quite generally to interacting quantum subsystems. It has been explored and extended along myriad different paths over the past decade. We summarize some of these here.
Exact potentials driving the nuclear motion The first applications [16,17] showed the usefulness of the XF approach in providing potentials that directly correlate with nuclear motion in laser-driven dynamics, aiding in the interpretation of the physical processes. We defer further discussion of this to Sec. III. The exact surfaces have been instructive for analysis of basis-dependence of MQC surface-hopping schemes for laser-driven bond-softening and bond-hardening processes [29], and control of laser-induced electron localization [30].
For field-free dynamics following a photo-excitation, the exact potential was found to track the BO surface initially occupied, develop a diabatic-like character as the wavepacket passed through a non-adiabatic coupling region and then form steps that bridged the different BO surfaces associated to the local nuclear wavepacket [31][32][33][34]. The steps appeared in a gaugedependent term which, when gauge-transformed to a vector potential, could be interpreted as a momentum-rescaling. Interestingly, in situations where nuclear wavepackets associated with different electronic states meet, the exact TDPES develops features such that independent classical nuclear trajectories evolving on it mimic the true nuclear distribution of interfering nuclear wavepackets [35]. In all these studies, the exact potentials were obtained from an inversion of the numerical solution of the full molecular TDSE.
Mixed quantum-classical and quantum trajectory methods The most practical impact that the XF has had so far has been in the development of rigorous MQC methods [24-26, 36, 37].
Similarly, using a hydrodynamic description of the XF equations, a trajectory-based representation for the electrons had been introduced in Ref. [60] coupled to quantum nuclei.
Inclusion of spin-orbit coupling Non-adiabatic processes mediated by spin-orbit coupling (SOC), i.e. intersystem crossings, play a key role in unraveling many excited-state processes such as photo-induced dynamics in organometallic complexes [61][62][63][64][65][66] and even collision reactions involving systems of light elements [67][68][69]. Refs. [70,71] extended the XF approach to derive a MQC approach that addresses both spin-allowed electronic transitions mediated by non-adiabatic coupling (NAC) as well as spin-forbidden transitions mediated by the SOC in the same theoretical construction.
Geometric phase The concept of a geometric phase arises naturally in XF and can be defined as the circulation integral of the vector potential along a closed path, both in the stationary formulation of XF [72][73][74] and in its time-dependent version [75][76][77][78]. It is worth underlining, though, that a non-vanishing geometric phase (not topological, as it is not quantized) in the framework of XF does not require invoking the adiabatic separation of electronic and nuclear motion in a molecule, thus it is a more general concept that does not rely on the particular representation (adiabatic vs. diabatic, for instance) used to describe the electronic subsystem.
Density functionalization The XF has been used to extend density functional theory (DFT) to non-adiabatic situations, where the conditional electron density n R (r) = N e |Φ R (r, r 2 ...r N )| 2 d 3 r 2 ...d 3 r N replaces the N e -body wavefunction as basic variable [79] of the theory. The full nuclear wavefunction χ(R) is retained, and the exchange-correlation functional of the usual adiabatic DFT is generalized to include nonadiabatic contributions arising from theÛ en coupling term; it has functional dependence on the conditional electronic density and paramagnetic current-density, nuclear wavefunction, vector potential, and the quantum geometric tensor. Ref. [80] proposed a local conditional density approximation, and applied the XF-based DFT to study dissociation in a Hubbard-model of the LiF molecule, showing that non-adiabatic electron-nuclear coupling results in a 0.5 Bohr elongation of the bond-length at which electron-transfer is onset.
Exact electron factorization The XF formalism can be applied within a purely electronic system, where the marginal factor is a function of one electronic coordinate r 1 , with the conditional factor a function of the rest, conditionally dependent on r 1 [81]. From the partial normalization condition and the antisymmetry of the full electronic wavefunction, it follows that the marginal density gives the exact one-body electronic density and current-density.
The exact XF potentials for an electron ionized by a field have been studied [81][82][83], and approximated, and they provide a new perspective for one-electron approaches commonly used in strong-field physics, such as the single active electron approximation.
Electronic embedding: strong correlation The XF idea can also be applied to an electronic wavefunction in Fock space [84][85][86], where some orbitals (e.g. the more strongly correlated ones) are selected to be those of the marginal (the "fragment") while all others are conditionally dependent on the occupation numbers of the marginal. XF then yields an embedding Hamiltonian on the fragment, which is approximated by solving the full system at a low-level; the embedded fragment is then solved with a high-level method. This has been demonstrated to be accurate for ground-state energies over a full range of weak to strong correlation in Hubbard model systems [85]. A generalized Kohn-Sham approach has also been used to find the conditional factor [86], using an orbital-dependent functional approximation, which accurately captured the topological phase diagram of a multiband Hubbard model.
Reverse factorization There is nothing in the factorization presented above that requires us to choose the marginal as the nuclear coordinate, and the conditional as the electronic.
It would be formally equally possible to choose the marginal factor to be the wavefunction for the electronic coordinate, and then the nuclear wavefunction is conditionally dependent on the electronic one. Although less intuitive, this is particularly useful for problems where we are primarily interested in electronic dynamics, since it yields a TDSE for the electrons, in which the potentials contain the effect of coupling to the nuclei. This was first studied in Ref. [19] with the example of laser-control of electron localization: the time-delay between the excitation of a nuclear wavepacket from the ground to the first excited state of H + 2 and an applied infrared pulse controls the localization asymmetry, i.e. the probability the electron ends up on the "left" or "right" nucleus. We will return to this briefly in Sec. III [19,20,27].
Mathematical analysis Formally, the XF expression of the time-dependent molecular wavefunction closely resembles the BO approximation. However, in the latter case, the elec-tronic conditional term is constrained to be an eigenstate of the BO HamiltonianĤ BO (r, R) all along the dynamics and, thus, it does not depend on time. Ref. [87] showed that the BO approximation can be recovered from XF as a limiting case for the electron-nuclear mass ratio µ tending to zero. Not only is this procedure an elegant scheme to connect XF to the well-known BO approximation, but it has been employed as well to develop algorithms that solve the coupled electronic (4) and nuclear (5) equations of XF. More specifically, the expansion of those equations in powers of µ helped in identifying leading non-adiabatic contributions in the electronic equation (4) to be retained in the approximate CTMQC algorithm (see Sec. IV), as well as in suggesting strategies for a perturbative treatment of non-adiabatic effects [88][89][90][91] (discussed below).
It is of fundamental interest whether the set of equations, Eqs. (4)(5)(6)(7)(8)(9), are numerically stable to propagate, but it is also of practical interest, since in making approximations, we would like to know which terms might need care in developing methods that converge stably.
In fact the mathematical form of the equations is unprecedented, and the usual numerical methods for TDSE's fail when applied in a straightforward way to the exact equations [92]; a deeper mathematical analysis can be found in Ref. [93].
Polaritonic chemistry The application of XF to systems of coupled electrons, nuclei, and photons [94], has provided new perspectives in the burgeoning field of polaritonic chemistry [95]. Confining a molecule to an optical or plasmonic cavity enhances the light-matter coupling strength that scales as the square-root of the volume, and distorts potential energy surfaces according to how close the energy gaps at a given nuclear configuration are to the resonant frequencies of the cavity, thus altering reaction pathways. Choosing the marginal factor for the nuclear coordinate provides the exact TDPES that drives the nuclear motion, directly correlating with the nuclear dynamics in contrast to the polaritonic surfaces [96] which are generalizations of the BO surfaces; Ref. [97] showed how its features result in the suppression of proton-coupled electron-transfer in a model system, unlike the polaritonic surfaces. To properly model the complex systems used in the experiments, approximations such as MQC would be required, and running classical dynamics on the exact TDPES [98] gives a big improvement than when using simply the weighted polaritonic surface. Work is in progress to extend the self-consistent XF-based MQC methods to the polaritonic case for practical applications. To account for many cavity modes, a classical Ehrenfest treatment for the dynamics of the photonic system has been explored [99], and using the XF with the marginal as the photonic displacement coordinate could explain the general underestimation of photon numbers in this approach [100] through comparing the exact potential driving the photons with that underlying the Ehrenfest approach. A further useful factorization could be Ψ(r, R, q, t) = χ(R, q, t)Φ R,q (r, t) since this would yield a TDSE for both the nuclear and photonic systems in which the potential "exactifies" the concept of cavity-BO surfaces introduced in Ref. [101]. Also, in general, when several identical, or non-identical particles compose the subsystems, nested commutators along the lines of the formalism of Ref. [102] could be considered.
Perturbative limit The XF has been employed to address diverse classes of static and dynamic problems in molecules [88][89][90] and solids [91] manifesting essential, though not strong, non-adiabatic effects. For instance, it has been observed that, in determining vibrational circular dichroism of chiral molecules [90,103,104] or in computing electronic current densities [88,105,106], the BO approximation fails to fully account for electronic effects during the dynamics. Even though non-adiabaticity is not strong in these examples because the system does not evolve close to avoided crossings or conical intersections, including excited-state effects as a perturbation to the BO (unperturbed) state is essentialand sufficient -to recover the correct behavior. In the context of XF, the electron-nuclear coupling term of Eq. (7) has been treated as a perturbation to the BO Hamiltonian [90], where the nuclear velocity is the small parameter. This results in a vector potential that can be seen as a position-dependent mass correction to the bare nuclear masses and that can be used to correct vibrational harmonic frequencies of hydrogen-based molecules [89] or phonon frequencies in solids [91].
Quantum-mechanical clock The separation of degrees of freedom at the basis of the stationary version of the XF was used in Ref. [107] to introduce the concept of time in quantum mechanics via the classical limit of the quantum-mechanical clock's degrees of freedom R represented by the marginal wavefunction. A clock-dependent continuity equation for the conditional wavefunction was then derived.

III. EXACT POTENTIALS DRIVING DISSOCIATION AND IONIZATION
Returning to the theme of dynamics of non-bound states, the exact TDPES of the XF has proved to be a useful interpretative tool. Consider a laser field causing dissociation and ionization in a diatomic molecule. If we are primarily interested in the dissociation aspect, the original factorization of Eqs. (4) -(9) might be most helpful for analysis since it provides a rigorous definition of the potentials that drive the nuclear motion. If we are interested in the ionization aspect, we may instead want to utilize the reverse factorization, to study the exact potential driving the electronic system. In this section, we discuss examples of these two potentials, previously studied in Refs. [16,17] and Refs. [20,27] respectively, for the case of a one-dimensional (1D) model of H + 2 , subject to a laser field. The model used in Ref. [16,17] is from Refs. [108][109][110][111], and involves one electronic coordinate and one internuclear (relative) coordinate with soft-Coulomb-interactions. We note that the laser does not couple directly to the relative nuclear coordinate: the dissociation is driven solely by the electronic motion. In the XF picture, from the nuclear viewpoint this is represented entirely by the TDPES and the time-dependent vector potential. For onedimensional problems, we can always pick a gauge in which the vector potential vanishes everywhere, and then the TDPES represents the complete potential driving the nuclear system. This is the gauge chosen in this study [16,17], and the TDPES is obtained by inversion of Eq. (5). The two-dimensional TDSE for the electron-nuclear wavefunction is first solved, starting in the ground state of the molecule, and then χ is constructed: its magnitude is the nuclear density obtained by integrating the full wavefunction over electronic coordinates and its phase is settled by the condition that A(R, t) = 0 [17].
The laser field has a wavelength of 228nm, which corresponds to a frequency of about twice the dissociation energy. A sketch is shown in the top middle panel of Fig. 1. We see from the black curves in the upper left panel (b), that for an intensity of I 1 = 10 14 W/cm 2 , the system dissociates, accompanied by significant ionization. This suggests a Coulomb explosion mechanism. On the other hand, for the weaker intensity of I 2 = 2.5 × 10 13 W/cm 2 , over the same time period there is dissociation with R reaching about half the distance as for the stronger field, but very little ionization (0.008 of an electron by the end of the time shown, as opposed to a little more than 0.5 for the stronger field). The other curves in these plots show the predictions from approximate methods: Ehrenfest, time-dependent Hartree, and exact-Ehrenfest. In the Ehrenfest approach [112,113], a single nuclear classical trajectory, initially located at the R of the initial nuclear probability density and with zero momentum, is propagated by the classical Hamilton's equations with a force given by the average of the soft-Coulomb nuclear-nuclear interaction and electron-density-weighted electron-nuclear interaction. In exact-Ehrenfest, on the other hand, a classical nuclear trajectory is propagated under the gradient of the exact TDPES [16,17]. These methods both treat the nuclear dynamics classically, while time-dependent Hartree (time-dependent self-consistent field) involves a quantum propagation of a nuclear wavefunction based on an uncorrelated factorization of the molecular wavefunction [113] resulting in the potential driving the nuclei determined self-consistently through the electronic density. We observe that while all three of the approximate methods predict dissociation in the stronger field I 1 , with exact-Ehrenfest being closest to the exact result, none of them is able to predict the dissociation in the weaker field I 2 case. The exact TDPES can help us understand why.
The TDPES, shown in the lower set of panels in Fig. 1b and 1c, shows a strong deviation from the BO surface, with large oscillations in time in the region from near the equilibrium distance on out. The tail of the TDPES alternately falls sharply and returns in correspondence with the field, which softens the bond, releasing the nuclear density. In this way the TDPES mediates the transfer of energy from the field-accelerated electronic motion to the nuclei. The features are similar for both the stronger and weaker field cases, and why a classical particle evolving on the TDPES does not therefore dissociate for the case of I 2 is perhaps not immediately obvious. The answer is revealed when we take a closer look at the structure of the surfaces close to the equilibrium distance and plot the classical energy of the particle evolving in the TDPES (the exact-Ehrenfest calculation) as a function of its position. These are black dots shown in the time snapshots. While for the stronger field I 1 , the particle quickly becomes loose of the binding potential, being kicked up by the TDPES oscillations, in the case of the weaker field it gets trapped in a narrow potential well not so far from the equilibrium bond length. Without tunneling, it cannot escape, and therefore dissociation does not occur in the classical calculation. This shows that the mechanism for dissociation of the quantum system in the weaker field is tunneling: the TDPES transfers field energy to the nuclei via tunneling, and although the classical exact-Ehrenfest calculation shows a larger amplitude of oscillation of R than the other approximation methods, it ultimately cannot tunnel through the barrier. Despite treating the nuclei quantum mechanically, so in principle being able to capture effects like tunneling, the time-dependent Hartree method still fails. This is because of its uncorrelated nature: the electronic density determining the potential for the nuclear motion is the same for every R, and does not deviate far from the exact electronic density near the equilibrium bond-length, due to the weak field strength. The reader is referred to Ref. [17] where the Hartree surfaces are explicitly shown. This example emphasizes that the mechanism of dissociation via tunneling requires not only an quantum mechanical description of the nuclei but also an adequate description of electron-nuclear correlation. We now briefly turn to the electron's viewpoint, and ask what is the potential that is driving its motion? In this case, we choose the reverse factorization where the marginal has the electronic coordinate: Ψ(r, R, t) = Φ r (R, t)χ(r, t), so that the equations are exactly the same as Eq. (4) -(9) but with r ↔ R. The electronic wavefunction satisfies a TDSE, and, again picking a gauge where A = 0, the TDPES appearing in the equation for Φ(r, t) (e-TDPES) together with the external laser field account completely for the electronic system's motion.
Although formally the same as the direct factorization equations Eq. (4) -(9), the structure of the TDPES is significantly different, with different terms in Eq. (8) contributing in different ways. In particular, the term in the e-TDPES that comes from the electronnuclear coupling operatorÛ en , 1 2me Φ r |∇ 2 j Φ r R (where ∇ j is the gradient with respect to the jth electronic coordinate) is significantly larger in the reverse factorization than in the direct factorization, due to the division by the much smaller electron mass. This term tends to give peak structures in the e-TDPES, which can have a significant impact on the electron dynamics. For example, in the case of laser-controlled electron localization [114][115][116], mentioned in Sec. II A, the peak led to an interatomic barrier that enhanced the one in the traditional pictures that comes from Coulomb interaction, resulting in an earlier localization time and smaller localization asymmetry than predicted by traditional methods [19]. The traditional methods use a "quasistatic" picture where the electron-nuclear electrostatic potential is added to the laser field to determine the electronic wavefunction's evolution. This work showed that the quasistatic approach misses important effects coming from dynamical electron-nuclear correlation, that can have a notable influence on the predicted observables.
Dynamical electron-nuclear correlation was also shown to be important in understanding the process of charge-resonance enhanced ionization (CREI). CREI refers to the phenomenon that the ionization rate from a molecule can be several orders of magnitude higher than the rate from the constituent atoms at a critical range of internuclear separations [109,[117][118][119][120][121].
This was originally explained by a quasistatic argument, where nuclei are instantaneously fixed point particles. The combined electrostatic electron-nuclear attraction, together with the field, leads to Stark shifted atomic levels, and the potential acting on the electrons involves both an outer barrier on the down-field side, and an internuclear barrier that grows as internuclear separation increases. Provided that the field is turned on fast enough such that population in the up-field level remains and does not tunnel back to the down-field atomic level, the molecule can rapidly ionize over both the inner and outer barriers, at an optimal internuclear separation. This gives rise to the enhanced ionization rate. By requiring that the Stark-shifted up-field atomic level exceeds the top of both the inner and outer field modified Coulombic barriers, one finds the critical internuclear separation for CREI as 4.07 divided by the ionization potential. The analysis can be generalized to the case of a laser field whose period is shorter than the tunneling time [117,[122][123][124].
This physical picture is very instructive, but in modeling the experiment, the premise of a quasi-static picture is questionable. A time-resolved study showed that even with static nuclei, the ionization tends to occur in multiple sub-cycle bursts [123,124], not at the peak of the field cycles, which was assumed in the prior analysis. Furthermore, to observe CREI, the molecule needs to be stretched to the CREI region rapidly enough that appreciable electron density remains un-ionized. In fact, in some experiments, CREI is not observed because there is too much ionization before the critical distance is reached [125][126][127]. Typically, a fraction of the nuclear density remains near equilibrium separations, while a fraction dissociates, and it is the electron density associated with the dissociating fragment that is responsible for CREI. Thus not only the point particle picture for the nuclei is not accurate, but also the nuclear dynamics means the quasistatic picture is not appropriate, especially as in some dissociating channels the fragment velocities can be comparable to the electronic velocities.
Ref. [20,27] showed that the dynamical electron-nuclear correlation terms in the exact e-TDPES are needed to accurately predict the ionization characteristics when modeling a real CREI experiment. Going beyond the quasistatic treatment by only accounting for the width and splitting of the nuclear wave packet is generally not enough to get the correct dynamics of CREI.
These examples show that the exact TDPES for nuclei or e-TDPES for electrons can provide an illuminating picture on the nature of dynamical electron-nuclear correlation in laser-driven dissociation and ionization, but such studies can only be done when the full molecular TDSE can be somehow solved exactly or highly accurately. In the next section, we discuss approximations based on the XF, which move the XF to be a useful approach for real systems.

IV. XF-BASED MIXED QUANTUM-CLASSICAL METHODS
Being an exact reformulation of the molecular TDSE, the XF equations are no easier to solve than the original molecular TDSE; in fact, due to numerical stability issues, they appear to be harder [92]! Although the form of the exact coupling terms can give us insight and help us interpret phenomena related to electron-nuclear correlation, to have a practical impact, the XF requires approximations. MQC approximations are a natural direction for this, justified by the large mass of the nuclei [87].
By taking a classical limit for the nuclear motion, and discarding terms in the equations justified by the exact studies to have a small effect [33], Refs. [24,25] derived the coupledtrajectory MQC (CTMQC) approach. In CTMQC the electrons are propagated quantummechanically whereas the nuclear dynamics consists of classical trajectories coupled through the nuclear quantum momentum When writing the nuclear wavefunction in terms of its modulus |χ(R, t)| and phase S(R, t), the term in parenthesis in Eq. (7) that depends on χ(R, t) can be rewritten as The first two terms on the right-hand side sum up to the nuclear momentum field, which has a clear classical interpretation. The last term, i.e. the quantum momentum, appears as an imaginary correction to the classical momentum, with no classical counterpart.
The conditional electronic wavefunction associated with each trajectory α is expanded Then the CTMQC equations for the electronic and nuclear evolution arė where we note the shorthand f (R (α) (t), t) = f (α) (t). The first terms in the electronic and nuclear equations are Ehrenfest-like termṡ with d (α) ν,lk being the non-adiabatic coupling vector (NAC) along the ν'th nuclear coordinate between BO states l and k evaluated at the coordinate R α (t), i.e. Φ l R (r)|∇ ν Φ k R (r) R α and E (α) k the k-th BO electronic surface. The second terms in Eqs. (12) and (13) are the corrections coming from XF : The term f (α) ν,l is a momentum that equals the time-integrated adiabatic force accumulated on the lth surface along the trajectory α: and Q In a trajectory-based scheme, at any given time, information of the position of all trajectories is required to reconstruct the nuclear spatial distribution and compute the quantum momentum. Hence, this is a coupled-trajectory scheme, not an independent-trajectory one.
CTMQC has been successfully employed to simulate several excited-state processes, such as the ring-opening process of oxirane [26,128] and the photo-isomerization of 2-cis-penta-2,4dienimiun cation, a retinal photoreceptor [129]. Key to its success in these problems is its ability to capture decoherence from first principles, unlike traditional MQC methods as we will shortly discuss.
Another approach is to use the electronic equation derived from XF, i.e. Eq. (12), in a surface-hopping scheme [21,22] where the nuclei evolve on a single BO surface at any time, making hops between them according to a stochastic algorithm. This XF-based surfacehopping (SHXF) has been applied to a range of complex systems, from organic molecules to molecular motors [130][131][132][133][134][135]. Part of what has made SHXF able to treat such large systems, is that the quantum momentum term is computed via auxiliary trajectories thus enabling an independent trajectory scheme; these auxiliary trajectories approximately mimic the local coupling of a trajectory to nearby ones. Instead, computing the quantum momentum using coupled trajectories in the same way as in CTMQC, but within a surface-hopping scheme, yields the coupled-trajectory surface-hopping (CTSH) method [23]. More details on the different ways to compute the quantum momentum are given in Sec. IV A, and an investigation of their effect on dynamics in Sec. V.
In traditional MQC schemes such as Ehrenfest and surface-hopping [1,2,113,136], the electronic system suffers from "overcoherence": after passing through a region of nonadiabatic interactions, the nuclear wavepacket splits with different parts of the wavepacket being correlated with different electronic surfaces. However, the Ehrenfest nuclear distribution follows a mean-field surface and does not split (even if it approximates averaged observables adequately), while, although the surface-hopping nuclear distribution is able to split with different branches evolving on different electronic surfaces, the electronic coefficients remain coherent. In a sense, surface-hopping is more unsettling than Ehrenfest, because of the disconnect between electronic and nuclear systems. Various ad hoc decoherence corrections have been imposed on the surface-hopping algorithm [1,137,138], most often involving an empirical parameter; they are often adequate to reproduce experimental results, but not always, and not reliably.
The quantum momentum term in Eq. (12) has been shown to give a first-principles description of decoherence in a number of complex systems as well as in detailed studies on model systems [21,25,26,129,134,139,140]. But its role goes beyond just decoherence. It was recently shown to be crucial to describe processes where multiple electronic states are simultaneously coupled in some regions of the nuclear configuration space via non-adiabatic coupling [135].

A. Computation of the Quantum Momentum
Currently, available codes to perform CTMQC simulations are G-CTMQC [23,141] and PyUNIxMD [22]. G-CTMQC includes the coupled-trajectory schemes CTMQC and CTSH, with hopping calculated either via fewest-switches hopping, or Landau-Zener, and it is interfaced with the potential library ModelLib [142] Furthermore, it has been extended to treat spin-forbidden electronic transitions mediated by spin-orbit coupling, both in the spin-diabatic and spin-adiabatic representations [70,71]. PyUNIxMD has CTMQC and SHXF (with fewest switches hopping probability) capabilities, and is interfaced with a number of electronic structure codes. Recently, the simulation of photo-isomerization dynamics of a protonated Schiff base was calculated using both SHXF and CTMQC, achieving results that are in good agreement with each other [143].
In the CTMQC algorithm there are two main ways to compute the quantum momentum.
The first is based on the idea of reconstructing the nuclear density using a sum of gaussians centered at the position of each classical trajectory Note that we have indicated with the symbol g σ (β) ν (t)) a normalized threedimensional gaussian (for each nucleus ν) centered at R  ν . Then, the quantum momentum is computed analytically applying Eq. (19). In the original CTMQC implementation [24][25][26] the variance of the gaussians was time-dependent and determined from the distribution of classical trajectories along the dynamics. However, to stabilize the dynamics and for practical reasons, a frozen gaussian approach with time-independent widths is currently used. This results in the following expression for the quantum momentum: where the slope of the quantum momentum is computed as Γ and the center of the quantum momentum is defined as Note that if, for each nucleus ν, all the gaussians carry the same time-independent width σ ν , given e.g. by the widths of the initial nuclear wavepacket, then the slope becomes independent of the trajectory index: where ρ v,lk = 0, the second term on the right-hand-side of Eq. (25) is required to be zero, and this gives an expression for the quantum momentum [26].
Note that this means the overall population transfer of the whole swarm of trajectories comes solely from the NAC term in Eq. (25) and the quantum momentum affects the electronic population dynamics only indirectly through the evolution of the coefficients appearing in the first term. The condition is imposed pairwise in the electronic states and for each degree of freedom, which results in a different quantum momentum center for each pair of electronic states which is required to be independent of the trajectory index, R R R Imposing the condition, the numerical expression satisfies where now the center of the quantum momentum is computed as A alternative approach is to choose the y-intercept of Eq.
The variance of the gaussian centered on the auxiliary trajectory σ ν in each degree of freedom ν is an input parameter, and, to eliminate empiricism is chosen to be that of the initial distribution. This has given good results for the real molecules studied. Recently, a new implementation of SHXF [144] was proposed in which a time-dependent width σ ν (t) is computed by imposing maximum overlap between nuclear wavefunctions on different BO states in phase space.
Both SHXF and CTSH follow the same equations, except for the computation of the quantum momentum, and, in particular, Eq. (29) where R The main advantages of SHXF are computational efficiency and numerical stability.
However, both SHXF and CTSH suffer from deficiencies inherited from the ad hoc nature of the surface-hopping procedure itself: the ambiguity of velocity-adjustments after a hop [134,[145][146][147], and the incorrect forces in the vicinity of a NAC region where the exact force tends to have a diabatic or averaged character rather than that from one surface. Further, although the disconnect between the single-surface evolution of the nuclei and the coherent electronic evolution is partially ameliorated with the XF term, reducing the difference between the electronic population at a given time and the fraction of trajectories evolving on a given surface, this so-called "internal consistency" is not guaranteed.
Finally, from the XF side, we note that the modified computation of the quantum momen- We choose our model system as the Tully model with an extended coupling region with reflection (ECR) [136], where an incoming nuclear wavepacket enters a region of extended non-adiabatic coupling between two BO surfaces that are asymptotically parallel. Along with other Tully models, the performance of CTMQC was studied in detail in Ref. [25].
The analytical form of the electronic Hamiltonian matrix elements, the BO surfaces, and the NAC are shown in Fig. 2. We will study two cases, a Gaussian nuclear wavepacket incoming on the lower surface from the left with a higher (k 0 = 30 a.u.) and a lower (k 0 = 10 a.u.) initial momentum. In the NAC region, some population is transferred to the upper surface and in both cases, the offdiagonal elements of the density matrix in the BO basis (coherences) rise and then eventually vanish after passing the non-adiabatic region: For the higher initial momentum, the nuclear wavepacket associated with the lower surface moves faster, separating in nuclear space from the part that has transferred to the upper surface, while for the low initial momentum case one branch of the nuclear wavepacket is reflected and the other is transmitted. The decay of the coherence cannot be accounted for by uncorrected surface hopping since the electronic equations are propagated fully coherently along each independent trajectory [136], and the method becomes "internally inconsistent" in that BO population determined from the modulus-square of the associated electronic coefficient differs from that determined by the fraction of nuclear trajectories evolving on that surface [1,2,113,[136][137][138]. the original definition. Lower panels: electronic coherences,    (Fig. 7). This is responsible for the rise in the CTSH populations between 3000 and 4000 a.u. Arguably the trend of the fraction of trajectories measure of populations in Fig. 5 is closer to the exact than that from the electronic populations, just as the running state population is closer to that of CTMQC in Fig. 7. The reason why the CTSH trajectories stay longer in the NAC region and reflect earlier than the CTMQC trajectories is related to the effective potential energies they experience: the ones which have hopped to the upper surface evolve on the more sharply rising upper BO surface, with a velocity penalty to conserve energy from making the hop, while the CTMQC ones evolve on an approximate TDPES which tends to be bounded by the upper BO surface but can be below it. The CTMQC trajectories (black dots moving on the thin red lines representing the BO energy curves in the left panels) leave the NAC region earlier than the CTSH ones (black dots in the right panels), and move further to the right, before the reflection occurs. While at the beginning of the non-adiabatic event (shown at time t = 2875 a.u. in Fig. 7) the distributions of CTMQC and CTSH trajectories are similar, the issues just described can be identified at later times, for instance at t = 3215 a.u. and t = 3395 a.u. in Fig. 7. The distribution generated by the CTSH evolution at t = 3795 a.u. and t = 4675 a.u. has a splitting that agrees better with CTMQC than at previous times.
Movies describing the full dynamics as documented in Fig. 7 are provided as Supplementary Material.
In Fig. 7, comparing the CTSH electronic populations (green dots) in the top right-hand plots at each time with the running state (red dots, or, compare with the bottom plots) attests to the breakdown of the internal consistency of the algorithm, discussed above.
Specifically, we note that in some regions of space and already at early times, the running state for each trajectory (the red dots are at zero for the ground electronic state and one for the excited state), differs from the value of the excited-state population for those same trajectories. When a trajectory is on the ground-state potential energy curve, the corresponding population associated to the excited state (ρ α 22 (t)) should be zero if internal consistency holds. However, this is not always the case, as evident from Fig. 7. For CTMQC, the issue of internal consistency does not exist. We show for completeness, in the same figure but in the middle plots of each panel, the contribution of the quantum-momentum XF term in the electronic evolution equation (in particular, we show as blue dots the second term in the right-hand side of Eq. (25) but without the trajectory-sum). We note that at times t = 3215, 3395, 3795 a.u. the effect of such term is non-zero, in contrast to what is predicted by CTMQC.
We turn now to the CTMQC and CTSH results using the original definition of the quantum momentum Eq. (23). With this definition we see from the lower left panel of  Fig. 8) can be larger in the transmitted part than it is in the NAC region, which is a big distinction with the modified quantum momentum result where it was zero in the transmitted region (compare Fig. 7 and 8). Driven by the quantum momentum (which is predominantly negative in this leading edge [139]) populations in this region start transferring back to the ground state after previously having transferred to the excited state, earlier than in when the modified quantum momentum definition is used. The effective TDPES shown by the electronic energies associated with each CTMQC trajectory shows a significant difference from that when using the modified definition of the quantum momentum in Fig. 7   The most practical impact of the XF approach so far has been in the development of rigorous MQC methods, and CTMQC and SHXF have both successfully predicted dynamics in a number of complex situations [26,[128][129][130][131][132][133][134], giving improved results from first-principles compared to traditional MQC methods [135]. These methods involve computing a term that depends on the nuclear quantum momentum, for which three different implementations exist.
In particular, with coupled trajectories, there is the original definition coming directly from The effects of these different methods were studied here on two cases of dynamics through Tully's model of an extended coupling region between two BO surfaces. In one case, the incoming momentum of the nuclear wavepacket is high enough that there is a single pass through the NAC region with significant transfer of population to the other BO surface, and we showed that while CTMQC with the modified definition yields accurate population transfer and coherence properties, CTMQC with the original definition of the quantum momentum suffers from the violation of the condition. The original definition however yields accurate population transfer and coherences when used within surface-hopping framework, provided the fraction of trajectories measure is used, but at the expense of violating internal consistency. This is because the hopping probability is only non-zero in regions of the NAC, so in this sense, takes care of the condition. This is true both for when the quantum momentum is computed using coupled trajectories, or from auxiliary trajectories.
In the second case, the incoming momentum is low enough that part of the wavepacket that transfers to the upper surface, reflects and recrosses the NAC region a second time. Here CTMQC with the modified definition does a good job, although underestimates the second rise in coherence, while the other methods are less accurate and also less in agreement with each other. Using the modified definition in a surface-hopping approach however yields poor electronic populations likely due to the incorrect forces on the nuclei from the BO surfaces in the region of interaction, and incorrect imposition of energy conservation on an individual trajectory level, that lead the trajectories to spend longer in these regions. Using the original definition of the quantum momentum leads to unphysical population transfer. Again the fraction of trajectories measure reduces the error and, temporarily, gives a reasonable trend.
The method that is most well-grounded from a theoretical point of view is CTMQC, and, with the modified definition of the quantum momentum, gives reliable results for both the populations and coherences. From a practical standpoint, it has the advantage of needing fewer trajectories for convergence, and no adjustable parameters. However, the method can be more numerically challenging to deal with than methods like surface hopping. In CTSH the forces that propagate the trajectories are purely adiabatic, and the NAC vectors, which are computational expensive for molecular systems, are not needed along the dynamics unless the velocity rescaling after a surface hop has occurred is performed along the direction of the corresponding NAC vector. But with coupled trajectories, patience is required to wait for all the trajectories at each time-step in order to compute the coupled-trajectory term.
The computational advantage of SHXF is evident due to its use of independent trajectories supplemented with auxiliary trajectories. So, in addition to the first-principles nature of the electronic equation, SHXF has an attraction due to its independent trajectories, but surface-hopping in itself is unsatisfactory due to the ad hoc aspects such as the very hopping procedure, choice of velocity adjustment, energy conservation on the independent trajectory level, and the choice of how to deal with forbidden hops. These problems contribute to the internal consistency, which is generally partly reduced by the XF-based quantum momentum term as seen in previous studies on decoherence, but we have seen here that it can also exacerbate it. Still, even then, we can sometimes temporarily take advantage of this internal inconsistency by using the fraction of trajectories measure, since it takes care of the condition that there should be zero net population transfer when there is zero NAC; this is only temporary since the evolving populations will continue to be wrong, and will give poor results when the system then meets another interaction region.
In conclusion, the XF offers a new playground for the development of non-adiabatic dynamics of electrons and ions, and we hope this work will be helpful to inform future adventures in this area.