Lucien Dupuy

Francesco Talotta

Federica Agostini
email: federica.agostini@universite-paris-saclay.fr

David Lauvergnat

Bill Poirier

Yohann Scribano
email: yohann.scribano@umontpellier.fr

Adiabatic and Nonadiabatic Dynamics with Interacting Quantum Trajectories

published or not. The documents may come

Introduction

The predictive power of molecular dynamics simulations crucially relies on the accuracy of trajectory-based representations of molecular systems to determine their structural and dynamical properties at a reasonable computational cost. A most common use of trajectorybased methods is to mimic the evolution of the nuclei by driving them with classical, quasiclassical, semiclassical, or quantum forces that represent the effect of the electrons, either in the ground state (i.e., for adiabatic dynamics), or including the effect of the excited states (i.e., for nonadiabatic dynamics). On the least computationally demanding end, purely classical trajectory simulations offer access to complex systems consisting of hundreds of atoms.

The primary limitation of the classical approximation is, however, also clear: inherently quantum effects such as delocalization, zero-point energy, tunnelling, and/or interferences, are completely missed.

Even for large systems, there are plenty of situations where quantum dynamical effects are known to be important-and therein lies the challenge. For example, the study of hydrogen transfer in proteins is important for understanding the multitude of biological functions supported by this fundamental reaction-which in some cases even requires consideration of deep proton tunneling. [START_REF] Benabbas | Deep proton tunneling in the electronically adiabatic and non-adiabatic limits: Comparison of the quantum and classical treatment of donor-acceptor motion in a protein environment[END_REF] Similarly, hydrogen tunneling dynamics was experimentally revealed in phenol-ammonia clusters activated by UV-photon absorption. [START_REF] Rodríguez | Direct evidence of hydrogen-atom tunneling dynamics in the excited state hydrogen transfer (ESHT) reaction of phenol-ammonia clusters[END_REF] Other aromatic biomolecules (e.g., indole and pyrole 3 in the presence of water and ammonia solvents) manifest a behavior that attests to quantum nuclear effects governing the dynamics in excited electronic states. These examples serve to highlight the importance of considering nuclear quantum effects, combined with multiple electronic states when necessary, that cannot be taken into account using only classical trajectories.

There is, however, no longer a need to restrict oneself to purely classical simulations. In the last decades, progress has been made towards pushing the limits imposed by the computational cost of exact or nearly exact quantum mechanical approaches (like time-dependent wavepacket methods, e.g., the multiconfiguration time-dependent Hartree, [START_REF] Meyer | The multi-configurational time-dependent Hartree approach[END_REF][START_REF] Meyer | Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF][START_REF] Gatti | Applications of Quantum Dynamics in Chemistry[END_REF] and even exact spectroscopic methods [START_REF] Sarka | Hitting the Trifecta: How to Simultaneously Push the Limits of Schroedinger Solution with Respect to System Size, Convergence Accuracy, and Number of Computed States[END_REF] ). Less expensive but more approximate computational methods are also available, such as semiclassical initial value representation, 8 approximated path-integral methods, [START_REF] Althorpe | Path-integral approximations to quantum dynamics[END_REF] nonadiabatic ring-polymer molecular dynamics [START_REF] Ananth | Mapping variable ring polymer molecular dynamics: A path-integral based method for nonadiabatic processes[END_REF][START_REF] Richardson | Communication: Nonadiabatic ring-polymer molecular dynamics[END_REF][START_REF] Bossion | Non-adiabatic ring polymer molecular dynamics with spin mapping variables[END_REF][START_REF] Weight | Ab initio symmetric quasi-classical approach to investigate molecular Tully models[END_REF] and nonadiabatic quantum instanton theory [START_REF] Lawrence | A general non-adiabatic quantum instanton approximation[END_REF][START_REF] Richardson | Ring-polymer instanton theory of electron transfer in the nonadiabatic limit[END_REF] -all of which have been applied to nonadiabatic processes.

An alternative route towards completely recovering nuclear quantum effects, while still maintaining a trajectory-based representation of nuclear dynamics, is provided by the socalled quantum trajectory methods (QTMs). In past years, such QTMs have been used to describe the quantum nuclear behaviour of different model [START_REF] Wyatt | Quantum dynamics with trajectories: introduction to quantum hydrodynamics[END_REF][START_REF] Babyuk | Multidimensional reactive scattering with quantum trajectories: Dynamics with 50-200 vibrational modes[END_REF][START_REF] Poirier | Reconciling Semiclassical and Bohmian Mechanics: IV. Multisurface Dynamics[END_REF][START_REF] Lopreore | Quantum wave packet dynamics with trajectories[END_REF][START_REF] Kendrick | Time-dependent wave packet propagation using quantum hydrodynamics[END_REF][START_REF] Garashchuk | Energy conserving approximations to the quantum potential: Dynamics with linearized quantum force[END_REF][START_REF] Parlant | Classical-like trajectory simulations for accurate computation of quantum reactive scattering probabilities[END_REF][23][START_REF] Talotta | Quantum trajectories for the dynamics in the exact factorization framework: a proof-of-principle test[END_REF] and molecular systems. [START_REF] Poirier | Reconciling Semiclassical and Bohmian Mechanics: VI. Multidimensional dynamics[END_REF][START_REF] Mazzuca | Description of proton transfer in soybean lipoxygenase-1 employing approximate quantum trajectory dynamics[END_REF][START_REF] Scribano | Communication: Adiabatic quantum trajectory capture for cold and ultra-cold chemical reactions[END_REF][START_REF] Talotta | Describing the photo-isomerization of a retinal chromophore model with coupled and quantum trajectories[END_REF] In this work, our aim is to focus on the theoretical and numerical aspects of one particular quantum trajectory formalism, and its application to adiabatic and nonadiabatic problems.

Therefore, we introduce here the theories that will be illustrated in the remainder of the paper as well as their position with respect to the state of the art.

Numerical QTMs that are referred to as synthetic are those that propagate the wavefunction's hydrodynamic fields (density and phase) along the quantum trajectories, and deduce trajectory dynamics at each time step from those field values. We refer the reader to the book of Wyatt [START_REF] Wyatt | Quantum dynamics with trajectories: introduction to quantum hydrodynamics[END_REF] for an excellent review. A major numerical benefit of such a trajectorybased representation in a Lagrangian frame is the fact that one essentially ends up with a moving grid (i.e., the trajectories themselves) that follows the probability flow. Using the moving grid to "sample" the nuclear configuration space in this manner thus allows one to substantially reduce the number of grid points needed to accurately describe, for instance, scattering processes. In many practical instances, use of QTMs rather than wavefunctionbased methods [START_REF] Babyuk | Multidimensional reactive scattering with quantum trajectories: Dynamics with 50-200 vibrational modes[END_REF] may be reasonably expected to provide better-than-exponential scaling with system dimensionality-whereas fixed-grid methods, in contrast, are always characterized by exponential scaling.

However, there are also severe numerical challenges introduced by the QTM approach.

The main challenge resides in the accurate and numerically stable computation of the quan-tum force that guides the quantum trajectory evolution-which is typically defined according to the Bohmian prescription. [START_REF] Bohm | A Suggested Interpretation of the Quantum Theory in Terms of "Hidden" Variables[END_REF][START_REF] Bohm | A Suggested Interpretation of the Quantum Theory in Terms of "Hidden" Variables[END_REF] This is a non-linear quantity involving third-order derivatives, which in a QTM context, must be computed on an unstructured moving grid.

Moreover, since the "quantum potential" (whose gradient is the negative quantum force) has an inverse dependence on density, the numerical derivatives become unstable in the vicinity of wavefunction nodes, leading to a breakdown of the simulation-the so-called "node problem." [START_REF] Wyatt | Quantum dynamics with trajectories: introduction to quantum hydrodynamics[END_REF] Various treatments for these numerical issues have been proposed in the literature-e.g., the moving least squares (MLS) approach, [START_REF] Lopreore | Quantum wave packet dynamics with trajectories[END_REF] finite differences for constrained trajectories in an arbitrary Lagrangian-Eulerian (ALE) moving frame, [START_REF] Kendrick | Time-dependent wave packet propagation using quantum hydrodynamics[END_REF] artificial viscosity forces, [START_REF] Kendrick | Time-dependent wave packet propagation using quantum hydrodynamics[END_REF] "bipolar" decompositions of the wavefunction, [START_REF] Poirier | Reconciling Semiclassical and Bohmian Mechanics: IV. Multisurface Dynamics[END_REF][START_REF] Poirier | Reconciling Semiclassical and Bohmian Mechanics: VI. Multidimensional dynamics[END_REF][START_REF] Poirier | Reconciling Semiclassical and Bohmian Mechanics: I. Stationary states[END_REF][START_REF] Trahan | Reconciling Semiclassical and Bohmian Mechanics: II. Scattering States for Continuous Potentials[END_REF][START_REF] Poirier | Reconciling Semiclassical and Bohmian Mechanics: V. Wavepacket dynamics[END_REF] and semiclassical treatments based on linearized quantum forces. [START_REF] Garashchuk | Energy conserving approximations to the quantum potential: Dynamics with linearized quantum force[END_REF][START_REF] Mazzuca | Description of proton transfer in soybean lipoxygenase-1 employing approximate quantum trajectory dynamics[END_REF][START_REF] Garashchuk | Quantum dynamics with Bohmian trajectories: energy conserving approximation to the quantum potential[END_REF][START_REF] Rassolov | Bohmian dynamics on subspaces using linearized quantum force[END_REF] More recently, a fundamentally different QTM-type approach has been developed, based on a complete reformulation of quantum theory in terms of trajectories. [START_REF] Poirier | Bohmian mechanics without pilot waves[END_REF][START_REF] Poirier | CCP6: Daresbury Laboratory[END_REF][START_REF] Schiff | Communication: Quantum mechanics without wavefunctions[END_REF] In this approach one still works with the same ensemble of quantum trajectories as in the original Bohmian QTM formulation described above. However, instead of using those trajectories to propagate the wavefunction-based hydrodynamic fields, the quantum force is computed directly from the trajectories themselves-without making any reference to the wavefunction itself. That this is even possible is no trivial development. In any event, from a numerical standpoint, the spatial variable x is replaced with a "trajectory-labelling coordinate" (usually denoted C), in terms of which the numerical grid becomes structured and stationary. This, in turn, leads to enormous numerical advantages in terms of evaluating the derivatives needed to compute the quantum force. Indeed, the trajectory-based form of the quantum force is generally found to overcome the node problem, although in the time-dependent case under-sampling of trajectories in areas of low density impacts the accuracy of finite differences schemes. This is discussed in more details below.

The aforementioned trajectory-based reformulation can be derived either from the stationary or from the time-dependent Schrödinger equation, resulting in differing properties for the trajectories and numerical considerations. In a previous work, some of the authors applied the stationary formulation to the study of quantum scattering processes on a chemical abstraction reaction model [START_REF] Parlant | Classical-like trajectory simulations for accurate computation of quantum reactive scattering probabilities[END_REF] as well as in an adiabatic quantum capture model of cold and ultra-cold chemistry, [START_REF] Scribano | Communication: Adiabatic quantum trajectory capture for cold and ultra-cold chemical reactions[END_REF] in both cases within the Born-Oppenheimer (adiabatic) approximation. In such contexts, the numerical propagation is extremely stable, essentially because the trajectories are solutions of an ordinary differential equation (ODE).

In contrast, time-dependent applications necessarily involve a partial differential equation (PDE). Since, as mentioned, x(C, t) (e.g., for the one-dimensional case) is the sought-for PDE solution, it is in principle necessary to impose both initial conditions and boundary conditions [i.e., x(t) at both C-grid edges] in order to solve the requisite PDE. The problem is that the boundary conditions are not known a priori. Here, we find a noteworthy disadvantage, in comparison with wave-based PDE solutions for which zero Dirichlet boundary conditions may almost always be presumed. In practice, we find that trajectory dynamics in the interior are not sensitively dependent on the choice of boundary conditions, provided the grid interval is sufficiently large, as might be expected. On the other hand, the quantum trajectory dynamics can in many cases become numerically unstable, if there is a scattering potential present. This tends to manifest as errors propagating in from the grid edges, however, rather than originating from nodes-i.e., this is not the node problem.

Various techniques have been developed that provide some improvement to the numerical stability in the time-dependent case, [START_REF] Schiff | Communication: Quantum mechanics without wavefunctions[END_REF][START_REF] Hall | Quantum Phenomena Modeled by Interactions between Many Classical Worlds[END_REF][START_REF] Cruz-Rodríguez | Quantum dynamics modeled by interacting trajectories[END_REF] without the need to invoke additional, sometimes costly, computational smoothing procedures such as those mentioned above, nor supplementary approximations for the quantum potential and for the quantum force. Rather, the "tricks" used here are in the vein of a judicious choice of boundary conditions and finite difference discretizations, and are therefore numerically exact. Curiously, the particular choices made here can be interpreted as determining the precise form of the "interworld potential"-according to the "discrete" many-interacting-worlds (MIW) interpretation of quantum mechanics that has sprung up from the trajectory-based reformulation. [START_REF] Hall | Quantum Phenomena Modeled by Interactions between Many Classical Worlds[END_REF] In the original, "continuous" MIW interpretation, however, [START_REF] Poirier | Bohmian mechanics without pilot waves[END_REF][START_REF] Schiff | Communication: Quantum mechanics without wavefunctions[END_REF] these are merely choices for the numerical discretization.

In any event, in the present work, we reexamine the time-dependent QTM framework for the trajectory-based reformulation, and build upon previous successes in two important ways.

First, in addition to exploiting the numerical techniques described above, we also introduce new methods that appear to greatly improve the numerical stability of the time-dependent QTM calculations. This is an extremely important development from the perspective of practical, widespread adoption of the approach as a generic robust computational tool. Second, we aim to extend the time-dependent formalism to nonadiabatic processes in the same manner as for adiabatic processes. However, in a "standard" formulation of nonadiabatic dynamics, the effect of several electronic states is accounted for via the inclusion of multiple potential energy surfaces (PESs), and couplings among them; 41 this picture is fundamentally different from adiabatic dynamics, where a single electronic state, and thus a single PES, contributes to the "classical force" guiding nuclear dynamics (in addition to the quantum force already discussed). Note that wavefunction-based nonadiabatic quantum trajectory dynamics have certainly been considered in the past. [START_REF] Poirier | Reconciling Semiclassical and Bohmian Mechanics: IV. Multisurface Dynamics[END_REF][42][START_REF] Hussein | Quantum Trajectory Description of the Time-Independent (Inverse) Fermi Accelerators[END_REF][START_REF] Curchod | Trajectory-based solution of the nonadiabatic quantum dynamics equations: an on-the-fly approach for molecular dynamics simulations[END_REF][START_REF] Curchod | On trajectory-based nonadiabatic dynamics: Bohmian dynamics versus trajectory surface hopping[END_REF][START_REF] Curchod | Excited State Dynamics with Quantum Trajectories[END_REF] However, these were done in a standard Born-Huang (or related) representation of nonadiabatic dynamics involving multiple PES components-which, for technical reasons, poses a source of fundamental difficulties in the interacting trajectory-based context. Accordingly, in this work, we instead invoke the exact factorization formalism in combination with the interacting quantum trajectories approach presented above [ 36,38]. The exact factorization yields nuclear dynamics under the effect of a single time-dependent classical force accounting for the electronic excited states, thereby avoiding the aforementioned technical difficulty.

Note that some of the present authors have already attempted to combine the standard time-dependent Bohmian QTM formulation with exact factorization to include quantum nuclear effects within a trajectory description of coupled electron-nuclear dynamics. 23,[START_REF] Talotta | Quantum trajectories for the dynamics in the exact factorization framework: a proof-of-principle test[END_REF] Nonetheless, numerical issues associated with that study prevented the authors from en-visioning extensions to actual practical applications. Accordingly, we now present in this paper-for the first time-a proof-of-principle illustrative study that demonstrates the potential of the new scheme, and in particular, the value of combining the time-dependent interacting-trajectory-based QTM approach with exact factorization. Although only onedimensional examples are considered here, we note further that the "continuous" MIW theory adopted here (unlike the "discrete" version) readily generalizes for many-dimensional applications as well. [START_REF] Schiff | Communication: Quantum mechanics without wavefunctions[END_REF] The paper is organized as follows. Section 2 provides a brief description of the employed methodology for our quantum dynamics calculations, starting with the exact factorization formalism in Section 2.1 and continuing with the interacting quantum trajectories formalism in Section 2.2. In Section 3 we present our numerical results for various types of scattering potentials, focusing on the adiabatic case in Section 3.1 and on the nonadiabatic case in Section 3.2. We present our conclusions and perspectives in Section 4.

Theoretical methods

Exact factorization of the electron-nuclear wavefunction

The non-relativistic molecular Hamiltonian

Ĥ(q, x) = 3Nn ν=1 -h 2 ∂ 2 ν 2M ν + Ĥel (q, x) (1)
describes a system of interacting electrons and nuclei, whose positions are collectively indicated as q and x, respectively. The nuclear kinetic energy operator is expressed in Cartesian coordinates and contains a sum over the 3N n nuclear degrees of freedom, each labeled with the index ν, with spatial derivatives ∂ ν with respect to nuclear positions; M ν are the nuclear masses. The electronic Hamiltonian Ĥel (q, x) is the sum of the electronic kinetic energy and of all interactions. The time evolution of the electron-nuclear system is dictated by the molecular time-dependent Schrödinger equation (TDSE)

ih∂ t Ψ(q, x, t) = Ĥ(q, x)Ψ(q, x, t) (2)
whose solution yields to the time-dependent molecular wavefunction Ψ(q, x, t).

In the exact factorization theory, [START_REF] Abedi | Exact factorization of the time-dependent electron-nuclear wave function[END_REF][START_REF] Agostini | Ultrafast dynamics with the exact factorization[END_REF] the molecular wavefunction is written as a product of two time-dependent functions:

Ψ(q, x, t) = ψ(x, t)ϕ(q, t; x) (3)
where the nuclear wavefunction ψ(x, t) evolves according to the nuclear TDSE

ih∂ t ψ(x, t) = 3Nn ν=1 [-ih∂ ν + A ν (x, t)] 2 2M ν + ε(x, t) ψ(x, t) (4)
and the electronic conditional factor ϕ(q, t; x), that depends parametrically on x, evolves according to the electronic equation ih∂ t ϕ(q, t; x) = Ĥel (q, x) + Ûen [ϕ, ψ]ε(x, t) ϕ(q, t; x)

Both evolution equations ( 4) and ( 5) contain the time-dependent vector potential (TDVP)

A ν (x, t) and the time-dependent potential energy surface (TDPES) ε(x, t) defined as

A ν (x, t) = ⟨ϕ(t; x)| -ih∂ ν ϕ(t; x)⟩ q (6) ε(x, t) = ⟨ϕ(t; x)| Ĥel (x) + Ûen [ϕ, ψ] -ih∂ t |ϕ(t; x)⟩ q (7)
The symbol ⟨•⟩ q stands for an integration over the electronic degrees of freedom, and we removed all dependencies on q within ⟨•⟩ q to imply that this variable is integrated out.

The TDVP and the TDPES are electronic quantities that evolve in time as effect of electronic dynamics, which is coupled to the nuclear dynamics in nonadiabatic conditions, and completely determine the nuclear evolution via the nuclear TDSE (4).

The electronic equation ( 5) contains the electron-nuclear coupling operator

Ûen [ϕ, ψ] = 3Nn ν=1 1 M ν [-ih∂ ν -A ν (x, t)] 2 2 + -ih ∂ ν ψ(x, t) ψ(x, t) + A ν (x, t) (-ih∂ ν -A ν (x, t)) (8)
which depends explicitly on the nuclear wavefunction, implicitly on the electronic wavefunction via the presence of the TDVP, and acts on the parametric dependence of ϕ(q, t; x) as a spatial derivative.

The ambiguity of the product form of the molecular wavefunction in Eq. ( 3) is partially removed by imposing the normalization condition |ϕ(q, t; x)| 2 dq = 1 ∀ x, t, which allows us to identify |ϕ(q, t; x)| 2 as a conditional probability density and |ψ(x, t)| 2 = |Ψ(q, x, t)| 2 dq as a marginal probability density that yields the nuclear density from the full wavefunction.

Equation ( 3) is thus unique up to a gauge encoded in a phase factor e (i/h)θ(x,t) (with θ(x, t) a real function): multiplying the nuclear wavefunction by this factor and the electronic term by its complex conjugate, Eq. (3) remains unaffected and the time-dependent potentials transform as standard gauge potentials. Therefore, this ambiguity has to be eliminated by imposing a choice of gauge.

In Section 2.2 we will work in a gauge where the TDVP is identically zero, which is a possible choice of gauge only for one-dimensional problems (in nuclear space), as those presented in Section 3. Another possible choice of gauge in high-dimensional situations, is to put to zero the TDVP along a "selected" direction in nuclear configuration space. Then, one could combine a classical and a quantum trajectory-based representation of nuclear dynamics, adopting the refined quantum description only along the "selected" direction. An alternative strategy is to reformulate Section 2.2 including the effect of the TDVP, along with the TDPES. Note that quantum trajectories' velocity field is a gauge-invariant quantity, thus they are the same for any choice of gauge. We are aware that one-dimensional models may not fully unravel the variety of quantum effects associated with nonadiabatic dynamics processes, for example, related to conical intersections in the Born-Huang picture. However, we do not expect fundamental or numerical complications in moving towards higher-dimensional problems, in particular from the perspective of the exact factorization.

Specifically, as shown in previous work, [START_REF] Ibele | A photochemical reaction in different theoretical representations[END_REF][START_REF] Agostini | When the Exact Factorization Meets Conical Intersections[END_REF][START_REF] Curchod | On the Dynamics through a Conical Intersection[END_REF] the time-dependent potentials arising in the exact factorization are capable of capturing effects related to conical intersections, with the advantage of not manifesting singularities, as is instead the case, for instance, of the nonadiabatic coupling vectors arising in the Born-Huang picture. In any case, it is clear that many interesting routes, currently under investigation, can be undertaken to combine the interacting quantum trajectories formalism with exact factorization.

One-dimensional quantum trajectories method

In this preliminary study, we apply the interacting quantum trajectories approach [START_REF] Poirier | Bohmian mechanics without pilot waves[END_REF][START_REF] Schiff | Communication: Quantum mechanics without wavefunctions[END_REF] for a single nuclear degree of freedom in the exact factorization framework; thus the vectorial notation x reduces to the single nuclear coordinate x. As stated in Section 2.1, we choose to work in the gauge where only the TDPES affects nuclear dynamics. With this gauge choice, the nuclear current density reduces to

J(x, t) = h M Im ψ * (x, t) ∂ψ(x, t) ∂x (9)
which becomes

J(x, t) = ρ(x, t) 1 M ∂S(x, t) ∂x (10)
when considering the polar form ψ(x, t) = ρ(x, t)e iS(x,t)/h , with ρ(x, t) the nuclear probability density and S(x, t) the phase of the nuclear wavefunction.

From the above, we identify M ẋ = ∂S(x, t)/∂x as the velocity field for the quantum trajectories (since we deal with one nuclear degree of freedom we remove all dependencies on the index ν). Inserting the polar form of the nuclear wavefunction in the nuclear TDSE (4), and separating real and imaginary parts, one obtains quantum hydrodynamic equations for the fields ρ(x, t) and S(x, t) in an Eulerian frame of reference,

∂ρ(x, t) ∂t = - ∂J(x, t) ∂x (11)
- ∂S(x, t) ∂t = 1 2M ∂S(x, t) ∂x 2 + ε(x, t) + Q ψ (x, t) (12)
where the quantum potential Q ψ (x, t) is defined as

Q ψ (x, t) = - h2 4M 1 ρ(x, t) ∂ 2 ρ(x, t) ∂x 2 - 1 2 1 ρ(x, t) ∂ρ(x, t) ∂x 2 (13)
Switching to a Lagrangian frame by using the relation

d dt = ∂ ∂t + ẋ ∂ ∂x (14)
and the definition of the velocity field, the coupled evolution equations along the quantum trajectories are then found to be as follows:

ρ(x, t) = -ρ(x, t) ∂ ∂x ẋ (15) Ṡ(x, t) = 1 2M ∂S(x, t) ∂x 2 -ε(x, t) -Q ψ (x, t) (16)
Note that in Eq. ( 15) the spatial derivative with respect to x acts on the nuclear velocity field in square brackets. Finally, taking the x-derivative of the phase equation above, one finds that the local shape of the quantum hydrodynamic fields drives the quantum trajectories via the sum of both classical and quantum potentials as follows:

M ẍ = - ∂ε(x, t) ∂x - ∂Q Ψ (x, t) ∂x (17)
The same basic equations apply in the trajectory-based formulation, except that as stated, the quantum potential is not computed from Eq. ( 13). To understand how Q is defined, it is first necessary to discuss the trajectory labelling coordinate C, in terms of which the PDE solution for the quantum trajectory ensemble is expressed as x(C, t). In principle, there is complete freedom in terms of how C is defined. One of the simplest choices is to take C to be the initial value of a given trajectory at time t = 0-i.e., C = x 0 = x(t = 0). [START_REF] Poirier | Bohmian mechanics without pilot waves[END_REF] Through probability conservation [i.e., Eq. ( 15)], one then obtains the following relation for the density at any time t:

ρ(x, t) = ρ 0 (x 0 ) x ′ (x 0 , t) , (18)
where ρ 0 (x 0 ) = ρ(x 0 , t = 0), x = (x 0 , t), and x ′ (x 0 , t) = ∂x(x 0 , t) ∂x 0 | t . Note that the "spatial" derivative of x is taken with respect to the labelling coordinate C, which in this case is just the initial value x 0 . The dimensionless quantity x ′ thus becomes a measure of the relative spacing of nearby trajectories over time, as compared to the initial spacing at t = 0. 16

Inserting the expressions for ρ(x, t) [from Eq. ( 18)] and for x ′ (x 0 , t) into the above wavefunction-based quantum hydrodynamic equations, all reference to the time-evolved wavefunction are now entirely removed. The resultant trajectory-based dynamical PDE

[used to solve for x(x 0 , t)] then becomes

M ẍ = - ∂ε(x, t) ∂x - 1 x ′ ∂Q(x ′ , x ′′ , x ′′′ , x 0 ) ∂x 0 , (19)
where the trajectory-based form of the quantum potential Q is now given as follows:

Q(x ′ , x ′′ , x ′′′ , x 0 ) = - h2 2M x ′2 ρ 1/2 0 (x 0 ) ∂ 2 ∂x 2 0 ρ 0 (x 0 ) 1/2 + h2 M ρ 0 (x 0 ) 1/2 x ′′ x ′3 ∂ ∂x 0 ρ 0 (x 0 ) 1/2 + h2 4M x ′′′ x ′3 - 5 2 x ′′2 x ′4 (20)
The trajectory ensemble dynamical PDE of Eq. ( 19) is fourth-order in "space" (i.e., x 0 ) and second order in time. Note that all trajectory interactions (which arise from the x 0 derivatives) are due to quantum forces; otherwise, the trajectories would not interact, and Newton's classical ODE would result.

We note that Eq. ( 19) is written out explicitly in a form that depends only on x(x 0 , t), constants, and the initial density ρ 0 (x 0 ). The presence of the latter quantity is simply an artifact associated with the coordinate choice C = x 0 . Indeed, the ρ 0 dependence can be easily removed by transforming to a "uniformizing" choice for C, in terms of which [START_REF] Schiff | Communication: Quantum mechanics without wavefunctions[END_REF] The explicit transformation from x 0 to the uniformizing C can be defined as follows:

ρ 0 (C, t = 0) = ρ(C, t) = 1.
C(x 0 ) = x 0 -∞ ρ 0 (x ′ )dx ′ (21)
Note that any two suitable labelling coordinates must be related to via a bijective function (i.e., monotonic, in the one-dimensional case). In addition, we stress that the labelling coordinate (which could also be termed the "Lagrangian coordinate") is necessarily timeindependent along a given trajectory: C(x(x 0 , t)) = C(x(x 0 , 0)) ∀t. Finally, for each trajectory in the ensemble, we have

x = x(C, t) (22)
As discussed, the uniformizing C of Eq. ( 21) implies trajectories x(C, t) that all bear the same probability density, in contrast to the x(x 0 , t) ensemble. From Eq. ( 18), however, the shape of the density profile ρ(x, t) can be retrieved at any time t as follows:

x ′ = ∂x ∂C = 1 ρ(x, t) (23)
This leads to a simplified trajectory-based dynamical PDE that depends only on x(C, t):

M ẍ = F ε (x, t) + F Q (x ′ , x ′′ , x ′′′ , x ′′′′ ), (24)
with classical force

F ε (x, t) = - ∂ε(x, t) ∂x , (25)
and quantum force

F Q (x ′ , x ′′ , x ′′′ , x ′′′′ ) = - h2 4M x ′′′′ x ′4 -8 x ′′′ x ′′ x ′5 + 10 x ′′3 x ′6 . ( 26
)
The classical force is indicated as the spatial (x) derivative of the TDPES of Eq. ( 7), which reduces to the adiabatic potential if the Born-Oppenheimer approximation is valid, as was shown in Ref. [ 52]. An explicit uniformizing-C expression for Q is straightforwardly derived by inserting Eq. ( 23) into Eq. ( 13), and making use of the chain rule:

Q(x ′ , x ′′ , x ′′′ ) = h2 4M x ′′′ x ′3 - 5 2 x ′′2 x ′4 (27)
Finally, we note that a more rigorous derivation of the above quantum trajectory expressions is also possible, based on a Lagrangian/action extremization procedure. [START_REF] Poirier | Bohmian mechanics without pilot waves[END_REF][START_REF] Poirier | CCP6: Daresbury Laboratory[END_REF][START_REF] Schiff | Communication: Quantum mechanics without wavefunctions[END_REF] In this approach, the Lagrangian is equal to the usual classical one, with the addition of a quantum contribution-i.e., L = T -V -L Q , where T is the kinetic energy. It is important to note that a gauge freedom exists in the definition of L Q . Thus in general, L Q need not be Qalthough this particular choice is allowed, and has the advantage that the resultant action, obtained by integrating the Lagrangian over time, or equivalently

Ṡ(x, t) = L(x, ẋ, x ′ , x ′′ , x ′′′ , t) = 1 2 M ẋ2 -ε(x, t) -Q(x ′ , x ′′ , x ′′′ ) (28)
yields the phase of the wavefunction. In this manner [i.e., together also with Eq. ( 23)], the wavefunction can be recovered at any time from the positions and actions of the quantum trajectories-although the ensemble of trajectories alone is sufficient to completely predict the dynamical evolution of any quantum system.

Numerical evaluation of the quantum force

In any numerical implementation, the continuous trajectory labelling coordinate C must be discretized-with specific, discrete values C i corresponding to the individual discrete quantum trajectories used in the calculation. Note that these values do not change over time-the C i grid is thus fixed, rather than moving. It can therefore also be structuredusually such that the grid-point spacing, (C i+1 -C i ) = ∆C, is uniform across the grid,

1 ≤ i ≤ n.
If, in addition, the uniformizing choice of C is presumed (as will be the case throughout the rest of this article), then each of the n discrete quantum trajectories carries the same probability, 1/n.

Having defined the discrete quantum trajectories, x i = x(C i , t), some kind of finitedifference scheme is needed to evaluate C derivatives numerically, in order to compute the quantum forces F Q i , acting on each x i trajectory. In previous numerical applications of this method, [START_REF] Schiff | Communication: Quantum mechanics without wavefunctions[END_REF][START_REF] Hall | Quantum Phenomena Modeled by Interactions between Many Classical Worlds[END_REF][START_REF] Cruz-Rodríguez | Quantum dynamics modeled by interacting trajectories[END_REF] the following expression for the quantum force was used, which can yield a reasonably stable numerical propagation (especially for free particles), and also accounts correctly for quantum effects:

F Q i ≃ h2 4M 1 (x i+1 -x i ) 2 1 x i+2 -x i+1 - 2 x i+1 -x i + 1 x i -x i-1 - 1 (x i -x i-1 ) 2 1 x i+1 -x i - 2 x i -x i-1 + 1 x i-1 -x i-2 (29)
Interestingly, Eq. ( 29) does not depend at all on ∆C-a consequence of the scaling laws that apply in the trajectory-based formulation. [START_REF] Tsai | Exploring the propagation of relativistic quantum wavepackets in the trajectory-based formulation[END_REF] Because the dynamical PDE is fourth-order in C, evaluation of the quantum force requires interaction between a given trajectory and its four nearest symmetrical neighbors.

Equation ( 29) is easily derived by expressing the quantum force as the C-derivative of a field quantity P , rather than as the x-derivative of Q, i.e.

F Q (x ′ , x ′′ , x ′′′ , x ′′′′ ) = ∂ ∂C P (x ′ , x ′′ , x ′′′ ) ( 30
)
with

P (x ′ , x ′′ , x ′′′ ) = 1 x ′2 ∂ 2 ∂C 2 1 x ′ (31)
Using central finite differences to re-express C-derivatives, along with the approximation

x ′ i ≃ (x i+1/2 -x i-1/2
)/∆C, Eq. ( 29) is recovered.

A numerical issue arises with the (two) first and (two) last trajectories at the edges of the discrete C i -grid ensemble, since these lack a sufficient number of neighbors to evaluate spatial derivatives using the finite-difference scheme. In reality, this is nothing but the aforementioned boundary condition difficulty. This problem is mitigated by introducing additional fixed virtual trajectories x 0 and x -1 set at -∞, and x n+1 and x n+2 at +∞. Note that (x 0x -1 ) = (x n+2x n+1 ) = +∞ is also presumed. With these choices, reasonable (i.e., non-singular) values for F Q are obtained using Eq. ( 29)-although as discussed, errors can still propagate in from the edges, particularly if there is a scattering PES present.

To improve numerical stability still further, some additional measures are also implemented here. Consider, as an example, a Gaussian wavepacket, for which the quantum force is linear in x. Clearly, Eq. ( 29) yields more accurate F Q values in the interior of the wavepacket than in the periphery, since the trajectories are less densely distributed in regions of low density. Numerical errors near the edges do not necessarily lead to a significant deterioration of the overall propagation, however-again, because probabilities are small in the periphery. If the edge errors are stable at least (even if somewhat large) then as discussed, comparatively large errors at the periphery have little effect in the interior. In any event, errors tend to manifest as trajectories that oscillate around their true course. Numerical instability is signalled by oscillations that grow in magnitude over time, and/or propagate into the interior region.

One very simple strategy that we now use to mitigate the onset of unstable numerical oscillations is to perform the calculation using an extremely small integration time step (in the present study, such small integration step is of the order 10 -4 fs for an ensemble of 400 coupled trajectories). Though surprisingly effective, this approach is not a panacea, and in any event rather expensive to implement. We therefore also introduce an initial grid-point relaxation procedure, prior to the numerical propagation, which operates as follows. To begin with, for every calculation performed here, the initial wavepacket is Gaussian, and therefore the initial quantum force is known analytically everywhere. As discussed, grid points are distributed uniformly in C, and remain uniform over time. However, even for a Gaussian wavepacket, such a distribution does not lead to numerically exact F Q i values via Eq. ( 29). The purpose of the relaxation procedure, therefore, is to shift the trajectories starting positions slightly (especially in the periphery) such that the numerically computed

F Q i errors are minimized.
More specifically, the initial relaxation step is performed by using Newton's algorithm [START_REF] Press | Numerical Recipes in FORTRAN[END_REF] to displace the initial positions x i (t = 0) in order to find the root of the n functions,

f i (x i-2 , x i-1 , x i , x i+1 , x i+2 ) = F Q i (x i-2 , x i-1 , x i , x i+1 , x i+2 ) -F exact Q (x i ), (32)
representing the differences between exact and numerical initial quantum forces. Since the extremal virtual trajectories are fixed, the n functions f 1≤i≤n depend only on the n positions

x 1≤j≤n , and so the inverse of the Jacobian matrix J ij = ∂f i ∂x j used in the Newton iterative root search is well defined. Note that the exact quantum force values, F exact Q (x i ) can be obtained analytically from the known Gaussian initial density through Eq. ( 13). In particular, the Gaussian form

ρ 0 (x) = e -(x-xc) 2 /(2γ 2 0 ) √ 2πγ 0 , (33)
where γ 0 and x c are the Gaussian width and center, respectively, leads to the following linear quantum force expression:

F exact Q (x) = h2 4mγ 4 0 x -x c (34)
Following the relaxation procedure described above, the trajectories' displacements are small enough that their representation of the initial density is not significantly deteriorated.

Yet, initial quantum force errors are substantially diminished, to the extent that subsequent numerical propagation in general becomes much more stable. The initial grid-point relaxation procedure thus strikes a nice balance between dynamical fidelity with (nearly) uniform wavefunction representation. In any event, how we recover the density from the trajectory distribution is detailed below.

Density synthesis

In the trajectory-based approach, the relation between x and C is known at a given time only through the discrete values, x i (C i , t). Hence, interpolation is needed to estimate the density at arbitrary x. Our interpolation procedure is as follows: we first generate a monotonic interpolation of C(x, t) by the means of a monotone quintic splines interpolation algorithm. [START_REF] Lux | An Algorithm for Constructing Monotone Quintic Interpolating Splines[END_REF] Then, taking the derivative of the resulting interpolating function yields the density at arbitrary points, thanks to the definition of C as ρ(x, t) = ∂C/∂x. The resulting density synthesis strictly obeys norm conservation, which is implicitly guaranteed by the definition of quantum trajectories and use of the trajectory labelling coordinate C. This method is favored over using finite differences to estimate the density from the spacing of trajectories [i.e., taking ρ(x i , t) ≈ 2∆C/(x i+1 -x i-1 ), then interpolating ln [ρ(x, t)] using cubic splines] as that procedure might violate the conservation law because of spline oscillations or "ringing".

3 Numerical results

Adiabatic case: Dynamics on time-independent PESs

The method presented in Section 2 is illustrated here by applying it to the scattering of an initially Gaussian wavepacket through one-dimensional time-independent PES profiles.

Assuming an initial Gaussian density of the form of Eq. ( 33), the initial phase is given by S(x, t = 0) = p 0 (xx c ), where x c < 0, and p 0 > 0 is the initial momentum, so that the wavepacket is incident from the left side of the PES barrier. Unless explicitly stated, we use atomic units henceforth.

The numerical procedure based on the propagation of quantum trajectories will be henceforth referred to as the "time dependent quantum trajectories" (TDQT) approach. Given that the trajectories are discretized and coupled through the quantum force as evaluated in Eq. ( 29), the dynamical PDE ( 24) can be thought of as having been replaced by the following set of coupled ODEs in terms of positions {x i } and momenta {p i = M ẋi } for the discrete ensemble of trajectories:

ẋi = p i /M ṗi = F Q i (x i-2 , x i-1 , x i , x i+1 , x i+2 ) + F ε (x, t) (35)
This set of coupled equations may thus be propagated using an ODE integrator. The time integration is performed using the well-known adaptive time-step Bulirsch-Stoer scheme [START_REF] Press | Numerical Recipes in FORTRAN[END_REF] with error tolerance set to 10 -9 . Moreover, an upper bound is set on the integration step so that no neighboring trajectories spacing should decrease by a factor above 40% from one time step to the next. For comparison, we also performed a wavefunction-based calculation of the time evolution of the corresponding TDSE, using the Crank-Nicholson algorithm.

Table 1 lists the parameters used for this calculation, which were sufficient to converge the results to an accuracy substantially beyond that obtained using TDQT. The TDSE results

will therefore serve as benchmark data.

We are interested in determining the transmission probability over time, which is obtained by integrating the density in the product region defined as [x P ; +∞], for some sufficiently large x P . For the TDQT calculation, this is straightforward, as every trajectory carries the same probability of 1/n. Accordingly, at any given time t, one need only count the number of trajectories, n trans (t), for which x(C i , t) > x P , and compare with n. Note that for this reason, the dynamical simulation cannot be expected to provide a transmission probability resolution finer than 1/n. Indeed, the exact result is only bound to be between n trans (t)/n and (n trans (t) + 1)/n. Hence, a sensible estimate of the transmission probability is obtained as

P T (t) = n trans (t) + 1/2 n (36)
with an uncertainty of 1 2n .

In this section, three PES profiles are considered:

1. The asymmetric Eckart potential:

V (x) = V 0 1 -α 1 + e -2ax + 1 + √ α 2 cosh(ax) 2 (37)
with parameters chosen to mimic the H + PH 3 → H 2 + PH 2 energy profile along the minimum energy path (see Table 1, second column).

2. The potential ramp:

V (x) = λ 1 + e -b (x-d) (38)
with parameters given in Table 1 (first column).

3. The symmetric Eckart potential:

V (x) = V 0 cosh(αx) (39)
parameterized to reproduce the one-dimensional hydrogen exchange reaction H+H 2 → H 2 + H (see Table 1, third column).

For the symmetric and asymmetric Eckart potentials, several calculations were performed by increasing the initial kinetic energy as indicated in Fig. 1; for the potential ramp, four calculations were performed by increasing the value of λ as indicated in Fig. 1 with an initial kinetic energy of E 0 = 3.42 × 10 -2 E H . In all cases, the value x P = 0 was used to compute the transmission probability.

E = V 0 /10 E = V 0 × 2/5 E = V 0 × 7/5 E = V 0 × 13/5
Transmission probability Table 1: List of parameters defining: the initial nuclear wavepacket, via γ 0 and x c ; the system and the potentials, via the mass m and a, α, V 0 for the Eckart barriers or b, d for the ramp; the numerical procedure for the TDQT propagation, via the ODE error tolerance and the number of trajectories; the numerical procedure for the Crank-Nicholson integration via the grid boundaries, the spacing dx and the time step dt. All quantities are given in atomic units. In Fig. 1 the panels on the left indicate the PES profiles, whereas the panels on the right present a comparison between the TDQT and benchmark Crank-Nicholson TDSE results.

Time [fs] TDQT C-N (b) 0 0.2 0.4 0.6 0.8 1 -8 -7 -6 -5 -4 -3 -2 -1 0 V (x) [λ] x [a 0 ] (c)

Parameters

Near perfect agreement is achieved in every case. In particular, absolute differences between the two calculations-which yield an estimate of the TDQT error-is always found to be less than the maximum expected value of ∆P max T = 1/n. For the Eckart PES problems, this value is ∆P max T = 2.5 × 10 -3 ; for the potential ramp problem, it is 3.3 × 10 -3 .

In Fig. 2 we study the robustness of TDQT results with respect to the number of trajectories. This is done by calculating the transmission probabilities as functions of time using

E = V 0 /10 E = V 0 × 2/5 E = V 0 × 7/5 E = V 0 × 13/5
Transmission probability Time [fs] n=400 n=300 n=200 n=100 mission probabilities become less smooth, since this quantity is estimated by counting the number of transmitted trajectories [see Eq. ( 36)], owing to their 1/n sampling capability. We conclude, then, that small ensembles of trajectories remain very accurate, under-sampling being the only (arguably mild) limiting factor to their accuracy.

Table 2: CPU execution time for propagating a wavepacket through the asymetrical Eckart barrier with initial kinetic energy E = V 0 /10. The Crank-Nicolson (C-N) method (results shown in Fig. 1) is compared to TDQT propagations varying number of trajectories (results in Fig. 2). Concerning the computational cost of the method, its execution time is similar to the Crank-Nicolson scheme for the numerical parameters of Table 1. In terms of memory use, storing information on a few hundreds of trajectories is preferable to storing the values of the complex wavefunction values on ten thousand grid points together with the evolution matrix.

In addition, the numerical results of Fig. 2 show that 400 trajectories is a reasonable number of trajectories to obtain accurate results. For a quantitative comparison between the Crank-Nicolson algorithm and TDQT, Table 2 compares their execution times in the case of the scattering through the asymmetric Eckart barrier with E = V 0 /10. The numerical cost of the TDQT approach does not evolve linearly with the number of trajectories, as the typical time step decreases with the minimal nearest neighbor spacing. The execution time goes from the order of two hours to two minutes when the number of trajectories decreases from 400 to 100, but we have to emphasize that the current implementation of the TDQT scheme is not completely optimised, and our Crank-Nicolson scheme is in its simplest form. Hence, there is more value in their relative efficiency than in absolute CPU time. These results clearly illustrate the advantage of a trajectory-based quantum dynamical propagation scheme over one using a fixed grid, and in particular the combined efficiency and accuracy of the TDQT method.

Nonadiabatic case: Dynamics on time-dependent PESs

As an illustrative study of the combination of the TDQT approach with exact factorization for nonadiabatic dynamics using the TDPES (7), we simulate the well-known Tully models. [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] These are one-dimensional models in nuclear space (as are those discussed in Section 3. The numerical calculation of the TDPES has been carried out as follows: (1) exact vibronic wavepacket dynamics calculations have been performed with the ElVibRot program 57 on the Tully models; (2) the output of ElVibRot has been analyzed by the Exact Factorization Analysis Code (EFAC) 58 to produce the TDPES. The TDPES can be easily expressed in terms of "standard" (a)diabatic nuclear amplitudes (and their spatial and time derivatives) arising from the Born-Huang representation of the molecular wavefunction. We refer the interested reader to Refs. [ 49-51] for a detailed discussion on those expressions. In any case, ElVibRot outputs the necessary information in the diabatic representation at various time steps throughout the propagation, and EFAC reads this information as input in order to reconstruct the TDPES. Additionally, EFAC imposes the gauge condition, which in the present case is simply taken to be A(x, t) = 0 ∀ t. Following from its definition given in Eq. ( 6), and expressing the electronic wavefunction as the ratio of the molecular and nuclear wavefunctions, from Eq. ( 3), the TDVP reads

A(x, t) = h Im[⟨Ψ(x, t)|∂ x Ψ(x, t)⟩ q ] |ψ(x, t)| 2 -∂ x S(x, t) (40)
with S(x, t) the phase of the nuclear wavefunction. It is easy to see that imposing the gauge condition A(x, t) = 0 yields an integral equation that defines S(x, t) in terms of the nuclear momentum field 59 -i.e., the first term on the right-hand side of Eq. ( 40).

ElVibRot calculations have been performed in the diabatic basis by initializing a Gaussian wavepacket in the lowest-energy electronic state in the negative x region. The Gaussian is centered at x c = -8.0 a 0 , and we considered three different values of the initial momentum p 0 = hk 0 , defined by the values k 0 = 10, 15, 20 a -1 0 , with an initial width γ 0 = 10/k 0 . Exact propagation has been carried out with the Chebychev scheme for which the evolution operator is expanded on Chebychev polynomials 60 with a time step 0.1 fs. The Hamiltonian has been normalized in order to ensure its spectral range lies in the interval [-1, 1]. The grid size and the number of grid points have been carefully optimized for each Tully model depending on the value of the initial momentum k 0 , in order to avoid reflections at the boundaries and to accurately capture the fine features of the nuclear wavepacket everywhere in space, even in very delocalized situations. The computational parameters for each calculation are listed in Table 3. ElVibRot provides as output the nuclear wavepackets in the diabatic basis, as well as their spatial and time derivatives. This information is required by EFAC to reconstruct the TDPES (7) in the gauge where the TDVP is zero.

Table 3: Grid parameters used in the exact quantum vibronic wavepacket propagations. In the remainder of this section, we show results obtained using the TDQT method, where the classical force is obtained as (minus) the gradient of the numerical TDPES. In addition to benchmarking the performance of TDQT against quantum wavepacket propagations, we also compare TDQT against a time-dependent classical trajectory (TDCT) approach, using only the classical force in the integration procedure. Note that our TDQT formalism does not involve any theoretical approximation, thus our "benchmarking" against exact wavepacket simulations only concerns the accuracy of the proposed numerical scheme, e.g., the numerical evaluation of the quantum force. On the other hand, TDCT are inherently approximate as they completely neglect the nuclear quantum potential and force. The purpose of the "benchmark" based on the comparison between TDQT and TDCT is, thus, to quantify the increase in accuracy that a quantum treatment of nuclear motion by the means of TDQT can provide. For each calculation, 5000 TDCT classical trajectories have been propagated using the same initial condition parameters used for the exact vibronic calculations, sampling the positions and momenta according to the Wigner distribution. Concerning the TDQT simulations, 400 trajectories were used in all cases presented in this section, with an ODE error tolerance of 10 -9 . Supplementary calculations were realized for models 1 and 3 to assess the robustness of TDQT results by varying the number of trajectories used. As quantum trajectory propagation is performed using an adaptive time step method, it was necessary to access the TDPES and associated forces at arbitrary times. A time interpolation procedure was therefore introduced, assuming linear evolution of the TDPES between two adjacent Chebychev time steps. Similarly, a linear spatial interpolation was performed to obtain the TDPES-derived force at arbitrary positions between grid points.

As observed in the Introduction, the results presented in this section are intended as a proof-of-principle illustration of the possibility to combine the exact factorization with a TDQT method. We find this illustration a crucial step towards the development of a trajectory-based nonadiabatic simulation scheme accounting for nuclear quantum effects.

Previous attempts have been made in this direction by some of the authors; 23,24 however, the outcome of those studies were not extremely promising for general applications, as already discussed in the Introduction. In contrast-and despite some small residual errors related to numerics, that will be documented below-the nonadiabatic dynamics produced by TDQT, and presented here as illustration, appears to be a promising route for further investigation.

We analyze in detail below the three models using k 0 = 10 a -1 0 . An overall assessment of the performance of TDQT for k 0 = 10, 15, 20 a -1 0 is provided at the end of the section.

In addition to that, we present as Supporting Information the movies for the three models and for the three values of k 0 representing the full dynamics and qualitatively comparing the vibronic wavepacket evolution with the TDQT. In all movies, we report: in the upper panel, the adiabatic, or Born-Oppenheimer (BO), PESs along with the TDPES at the positions of the trajectories; in the middle panel, the exact nuclear density and its estimate using TDQT;

in the lower panel, the phase of the nuclear wavefunction given by TDQT. Note that in some cases, the TDPES appears rigidly shifted along the y-axis and is not superimposed to any of the BOPESs. This feature has been already observed in previous work [START_REF] Abedi | Dynamical steps that bridge piecewise adiabatic shapes in the exact time-dependent potential energy surface[END_REF][START_REF] Agostini | Mixed quantum-classical dynamics on the exact time-dependent potential energy surfaces: A novel perspective on nonadiabatic processes[END_REF][START_REF] Agostini | The exact forces on classical nuclei in non-adiabatic charge transfer[END_REF][START_REF] Curchod | An exact factorization perspective on quantum interferences in nonadiabatic dynamics[END_REF] and does not affect in anyway the nuclear dynamics with TDQT, which depends only on the gradient of the TDPES. Upper panels: BOPES (black) and TDPES (gold). Lower panels: Nuclear densities from exact calculations (black), TDCT (dashed blue) and TDQT (red).

In Fig. 3 we present some results for the single avoided crossing model where the initial momentum of the incoming wavepacket is k 0 = 10 a -1 0 . Two snapshots are shown in the figure, at t = 40 fs (left panels) and at t = 120 fs (right panels): in the upper panels the (static) adiabatic, BOPESs are reported in black, together with the TDPES in gold at the two snapshots; in the lower panels the exact nuclear density is shown in black, and is compared to the density reconstructed from TDQT (red) and from TDCT (blue). The evolution of the TDPES manifests the nonadiabatic event by forming a pronounced peak between 0 and 2 a 0 that deforms the Gaussian shape of the nuclear density at t = 40 fs, and also by developing a small bump between 10 a 0 and 15 a 0 that splits the nuclear density into two portions.

In particular, the right most portion of the nuclear density is associated with the lower electronic states, whereas the left portion is what has been "transferred" to the upper state at the avoided crossing. For this simple model, TDQT is in extremely good agreement with the benchmark calculations, even though the TDPES is known only numerically on a spatial grid. For instance, evaluating the integrated density beyond x = 12 a 0 at t = 120 fs, one gets a transmitted probability value of 0.84625 using TDQT, which is very close to the exact value of 0.84518 (TDCT yields 0.88584).

It is worth mentioning that the TDCT density is produced based on the classical distribu- The dual avoided crossing model, with initial wavepacket launched towards the avoided crossings with initial momentum of k 0 = 10 a -1 0 , presents very quantum-mechanical behavior.

Specifically, after the first crossing, where the nuclear density branches into a ground-state and an excited-state portion, the nuclear wavepackets meet again at the second avoided crossing and interfere. As shown in Fig. 5, at short time, TDQT is capable of reproducing the recombination of the two portions of the nuclear wavepacket (t = 40 fs). Furthermore, TDQT captures quite well the portion of density that remains localized in the lower-state well at around 0, even though the fine details manifesting interferences are missed. On the other hand, TDCT nearly completely misses this part of the density, and already at t = 40 fs we observe some deviations from the reference. In general, we observed that for this model the numerical TDPES is quite noisy, which is the most likely reason for deviations of TDQT results from the quantum density at t = 120 fs. Note that in all nonadiabatic calculations, the TDPES is used directly as provided by EFAC and no smoothing procedure is applied to cope with the numerical inaccuracies which are due to the lack of resolution in time and space. In addition, it is important to underline, that errors arising from grid periphery are essentially made negligible by suitable choices of the grid boundaries.

In this case, we do not compare calculations performed using different number of TDQT trajectories because, as the density trapping phenomenon at around 0 involves only a tiny fraction of the probability/trajectories, we expect that decreasing the number of trajectories does not yield accurate resolution of the density in this region.

Figure 6 shows numerical results for the Tully model with extended coupling region and reflection, using k 0 = 10 a -1 0 as the initial momentum. When the nuclear wavepacket travelling in the lower state passes through the coupling region, it transfers population to the upper state, before reaching the branching portions of the BOPESs. Afterwards, for such a low-energy, the lower-state wavepacket is transmitted towards positive values of x, and decoheres from the upper-state wavepacket. Furthermore, the latter is reflected and crosses the coupling region again: the oscillations appearing in the TDPES and in the nuclear density at t = 35 fs attest to the recoherence of this portion of the wavepacket. These subtle quantum effects are all captured well by TDQT, which is in extremely good agreement with the exact benchmark results, whereas TDCT misses completely the quantum oscillations. Later in time, at t = 60 fs, we observe that the TDQT results deviate from the benchmark (as does TDCT), which is probably due to the fact that the TDPES shows numerical instabilities before developing smooth behavior, as evident at t = 60 fs. Pseudo-nodes become severe and cannot be resolved with the number of trajectories employed. The distribution becomes noisy, but still maintains qualitative agreement with exact results. It is worth noting that when using exact factorization and the TDPES, it is not possible to recover information related to individual electronic states-neither in the adiabatic nor in the diabatic basis representations. A single PES produces the evolution of the nuclear wavepacket, which is not resolved according to its adiabatic components. Therefore, information about the occupation of the electronic states is not accessible. However, in order to present a more in-depth analysis of the models and to circumvent this feature of exact factorization, we show in Fig. 7 the transmission probability through a dividing surface placed at x d = 2 a 0 , for the third model discussed in this section. This observable provides indirect information about the occupation of the electronic states, because the transmitted density, at the low initial momenta used here, is only found in the lower electronic state. The results shown in Fig. 7 attest to the importance of nuclear quantum effects to recover the correct dynamics, as TDCT (short-dashed lines) does not reproduce exact results (continuous lines).

On the other hand, TDQT (long-dashed lines) reproduces very well the benchmark calculations, despite the numerical instabilities described above. The step observed in TDCT results on the transmission probability for the case k 0 = 10 a Finally, we evaluate the errors in reproducing the whole quantum dynamics using the TDQT and TDCT methods, by defining and calculating a distance between the exact nuclear density and the density reconstructed using the trajectory-based methods, namely

ϵ(t) = ∞ -∞ |ρ traj (x, t) -ρ exact (x, t)| dx (41)
ǫ(t) k 0 = 10 a -1 0 k 0 = 15 a -1 0 k 0 = 20 a -1 0 ǫ(t) ǫ(t) t/t f

Conclusion

The TDQT approach presented in this work for propagating quantum trajectories in both adiabatic and nonadiabatic conditions shines by its very simplicity, at least in the onedimensional test studies proposed in this work. In particular, no fitting procedure is needed to recover the quantum force (e.g., via moving weighted least squares), nor is any modification of the trajectory ensemble itself required during the ensemble propagation (as would be necessary, e.g., in the ALE approach). These and other complicated refinements of conventional QTM (as described in the Introduction) can provide modest benefit, but often only at the expense of substantial and highly problem-dependent "parameter-tweaking". At heart, these are all means of addressing the inherent limitations associated with the node problem-and especially, with numerical errors caused by the use of moving, unstructured grids in x space.

In contrast, the simplicity and accuracy of the TDQT approach stems from its use of fixed, structured grids in C space, that never change over time. This fortuitous state of affairs is ultimately due to the trajectory-based reformulation, and its replacement of x with C as the requisite "spatial" coordinate-a change that also allows for much more natural comparison and integration with classical theories. As a result, the node problem does not lead to a breakdown of the simulation, quasi-nodes inducing a much milder under-sampling problem in our case, rather than a fatal numerical blow up of errors as in standard QTM approaches.

Differentiation errors are greatly reduced and can be evaluated without much fanfare. While this promise was recognized early on in the development of the TDQT theory, [START_REF] Poirier | Bohmian mechanics without pilot waves[END_REF][START_REF] Schiff | Communication: Quantum mechanics without wavefunctions[END_REF] other numerical issues have prevented the approach from reaching its fullest potential. Until now, that is-at least according to what the present results seem to suggest.

In keeping with our theme of simplicity, integration of the TDQT approach with exact factorization, as a means of addressing the nonadiabatic regime, also appears to provide the "perfect marriage" of methodologies. In particular, the collection of multiple BOPES surfaces (and their couplings) that characterize the standard approach, is replaced with just a single TDPES (ignoring the TDVP for now)-that is treated in TDQT in exactly the same manner as a single-PES adiabatic calculation. This is key, because it appears to be extremely difficult to extend TDQT theory to multiple components. Conversely, within the exact factorization framework, TDQT appears to offer a much cleaner and more effective quantum trajectory methodology than do conventional wavefunction-based QTMs. In particular, the behavior of the TDQT trajectory ensemble evolution-even when computed from a numerically-determined, and thus noisy, TDPES-demonstrates its stability without introducing ad hoc smoothing procedures (as for example, in the form of viscosity forces also used by more traditional QTMs).

The interacting trajectory-based reformulation, and ensuing TDQT methodology, hence allow one to recover key nuclear quantum effects necessary to describe the correct quantum dynamics of both adiabatic and nonadiabatic processes. In future, we envisage several clear developments of TDQT, to extend calculations to higher dimensions (for which the theoretical equations have already been derived), and to combine it with TDCT (which is straightforward using trajectory-Lagrangian-based action extremization). Treatment of nonadiabatic processes with several nuclear degrees of freedom will only allow for setting the TDVP to zero by gauge choice along a single dimension, but this does not undermine in any way the relevance of what is presented here as quantum trajectories are gauge-invariant. In addition, the possibility of combining the TDQT approach with CT-MQC, the coupled-trajectory mixed quantum-classical algorithm derived from the exact factorization, is currently being explored.

E = 5 ×E = 3 × 10 - 2

 53102 10 -5 E H E = 1.562 × 10 -3 E H E = 6.25 × 10 -3 E H E = 1.25 × 10 -2 E H E = 2 × 10 -2 E H

λFigure 1 :

 1 Figure 1: Left panels: PES profiles for the symmetric (a) and asymmetric (b) Eckart barriers, and the uphill ramp (c). Right panels: Transmission probability as function of time (in fs) through the symmetric (a) and asymmetric (b) Eckart barriers, and the uphill ramp (c), indicating TDQT results in red and Crank-Nicholson results in blue.

n

  = 100, 200 and 300 trajectories (they are shown as violet, green and blue curves, respectively, in the figure) and by comparing them to previous results (red curves in the figure) using an ensemble of size n = 400. Only results for the asymmetric Eckart barrier are shown,as the behavior in all systems is similar. Figure2shows an excellent agreement among all calculations. As the number of trajectories decreases, the curves representing the trans-

Figure 2 :

 2 Figure 2: Transmission probability through the asymmetric Eckart barrier obtained by propagating n = 100, 200 and 300 trajectories (represented as violet, green and blue curves, respectively) and compared to the previous results using 400 trajectories (red curves).

  1), and include two coupled electronic states: Tully model 1 presents a singled avoided crossing between the potential energy curves; Tully model 2 is a dual avoided crossing model; Tully model 3 shows an extended coupling region and manifests possible reflections of the incoming wavepacket. In this work we use the original models, 56 whose electronic Hamiltonian Ĥel is given in the diabatic basis representation.

Figure 3 :

 3 Figure 3: Snapshots at times t = 40 fs (left panels) and 120 fs (right panels) for the single avoided crossing model with k 0 = 10 a -10 . Upper panels: BOPES (black) and TDPES (gold). Lower panels: Nuclear densities from exact calculations (black), TDCT (dashed blue) and TDQT (red).

  tion of trajectories evolved on the TDPES. A histogram is generated from such distribution, and by centering Gaussians on each trajectory at each time step, the smooth function shown in the lower panels of Fig.3is reconstructed.

Figure 4 : 29 Fig. 4 Figure 5 :

 42945 Figure 4: Snapshots at times t = 40 fs (left panel) and 120 fs (right panel) for the single avoided crossing model with k 0 = 10 a -10 . The exact density recovered from the reference wavepacket propagation (black curve) is compared to TDQT results using 400 (red dots), 300 (blue dots), 200 (green dots) and 100 (violet dots) trajectories.

Figure 6 :

 6 Figure 6: Same as in Fig. 3, for the extended coupling region with reflection model at t = 35 fs (left) and 60 fs (right) with k 0 = 10 a -1 0 .

- 1 0Figure 7 : 1 0

 171 Figure 7: Transmission probability for the model of extended coupling region with reflection. Results of exact benchmark calculations are shown as continuous lines, of TDCT as shortdashed lines, and of TDQT as long-dashed lines. The different colors indicate different initial momenta of the incident nuclear wavepacket, namely red for k 0 = 10 a -1 0 , blue for k 0 = 15 a -1 0

Figure 8 :

 8 Figure 8: Transmission probability for the model of extended coupling region with reflection for the case k 0 = 10 a -1 0 using different numbers of trajectories in TDQT: n = 100 in violet, n = 200 in green, n = 300 in blue and n = 400 in red. The result of the exact calculation is shown in black.

Figure 9 :

 9 Figure 9: Distance between the exact density and the densities reconstructed using TDQT (continuous lines) and TDCT (dotted lines), as defined in Eq. (41). The colors indicate different initial momenta, namely red for k 0 = 10 a -1 0 , blue for k 0 = 15 a -1 0 and green for k 0 = 20 a -1 0 , used to initialize the dynamics in Tully model 1 (top), Tully model 2 (middle), and Tully model 3 (bottom). The time axis is normalized to the final time used in each simulation.

Figure 9

 9 Figure 9 presents ϵ(t) as function of time, where the time axis is scaled to the final times t f of the simulations, as given in Table4.

Table 4 :

 4 figure attests to the improvement in numerical accuracy that TDQT can provide over TDCT, even insofar as reproducing wavepacket dynamics is concerned. It should also be remembered that the TDCT calculations were performed using more than one order of magnitude more trajectories than the TDQT calculations. While decreasing the number of classical trajectories does not affect the qualitative behavior of the ensemble, we observed that n = 5000 is the minimum number of TDCT trajectories allowing us to reproduce a smooth nuclear density.

Acknowledgments

This work has been supported by the Agence Nationale de la Recherche (ANR-HYTRAJ), Contract No. ANR-19-CE30-0039-01. Y.S. and L.D. thank support from the High Performance Computing Platform MESO@LR at the University of Montpellier. BP acknowledges the Robert A. Welch Foundation for support through grant D-1523. Finally, all authors would like to dedicate this submission to the memory of Christoph Meier-a true devotee of QTM, particularly of the interacting-trajectory-based variety.

Supporting Information Available

The following files are provided as Supporting Information: tully1k10.mp4, tully1k15.mp4, tully1k20.mp4, tully2k10.mp4, tully2k15.mp4, tully3k20.mp4, tully3k10.mp4, tully3k15.mp4, tully3k20.mp4. The movies for the three Tully models and for the three values k 0 = 10, 15, 20 a -1 0 represent the full dynamics and qualitatively compare the vibronic wavepacket evolution with the TDQT.

TOC Graphic