Nonadiabatic Molecular Dynamics Based on Trajectories

Performing molecular dynamics in electronically excited states requires the inclusion of nonadiabatic effects to properly describe phenomena beyond the Born-Oppenheimer approximation. This article provides a survey of selected nonadiabatic methods based on quantum or classical trajectories. Among these techniques, trajectory surface hopping constitutes an interesting compromise between accuracy and efficiency for the simulation of mediumto large-scale molecular systems. This approach is, however, based on non-rigorous approximations that could compromise, in some cases, the correct description of the nonadiabatic effects under consideration and hamper a systematic improvement of the theory. With the help of an in principle exact description of nonadiabatic dynamics based on Bohmian quantum trajectories, we will investigate the origin of the main approximations in trajectory surface hopping and illustrate some of the limits of this approach by means of a few simple examples.


Introduction
Traditionally, ab-initio molecular dynamics (AIMD) is described within the so-called Born-Oppenheimer approximation, which assumes that the electronic and nuclear dynamics can be adiabatically separated [1], due to a large difference in mass between nuclei and electrons.Within this approximation, one usually solves the time-independent electronic Schrödinger equation for a given nuclear configuration [2] and then computes the quantum mechanical forces acting on the nuclei from the gradient of the corresponding eigenvalues, which depend parametrically on the nuclear coordinates and form the so-called potential energy surfaces (PES).However, in the description of most photophysical and photochemical processes, the electronic and nuclear dynamics become entangled, and therefore, more accurate nonadiabatic molecular dynamics schemes that go beyond the Born-Oppenheimer (BO) approximation are required.The most commonly used ab initio nonadiabatic molecular dynamics schemes are those based on the mixed quantum/classical propagation of an ensemble of (quasi-) classical trajectories [3][4][5][6], which, to some extent, reproduce the quantum dynamics of the nuclei.These mixed quantum/classical methods are especially popular, because they only require that the necessary electronic structure properties be computed on-the-fly, i.e., only at the points in the configuration space visited during the dynamics, therefore making the calculation of the full potential energy surfaces unnecessary.These approaches can be implemented numerically using electronic structure methods, such as Kohn-Sham Density Functional Theory (DFT) [7,8] and its time-dependent version (TDDFT) [9][10][11][12] or wavefunction-based approaches, such as Complete Active Space Self-Consistent Field (CASSCF), Multireference Configuration Interaction (MRCI) and Second-Order Approximate Coupled-Cluster (CC2) [13].
Among all nonadiabatic AIMD schemes, Tully's fewest switches trajectory surface hopping [14,15] (TSH) and its extensions to mixed quantum/classical dynamics [16] are probably the most widely used.
In the framework of TSH, the nuclear wave packets on the different PESs are described as a swarm of independent classical trajectories, while the nonadiabatic couplings induce hops of the trajectories from one electronic state to another; the occurrence of a trajectory hop is governed by the evaluation of a hopping probability, which depends on the temporal evolution of state amplitudes (Tully's coefficients) and on the value of the nonadiabatic couplings.
Despite the success of the nonadiabatic trajectory-based approaches, there are many quantum mechanical phenomena that cannot be entirely captured within this framework, namely nuclear quantum effects, like wave packets interference [22] and decoherence [26][27][28] and tunneling.Quantum dynamics methods based on a quantum mechanical representation of both electronic and nuclear degrees of freedom have also become available (see, for example, [29]).However, their high computational cost and the need for a numerical fit of the relevant PESs prior to propagation have limited their application to just a few nuclear degrees of freedom, and therefore, they are not yet suited for the simulation of complex molecular systems.
One possible way to account for quantum nuclear effects within a trajectory-based framework consists in the use of quantum (or Bohmian) trajectories [30][31][32].This approach emerges from a transformation of the time-dependent Schrödinger equation using a polar representation of the complex nuclear wavefunction (see [33] and Section 2 below).Robert.E. Wyatt and coworkers have recently introduced a numerical formulation of Bohmian dynamics using a trajectory-based solution of the so-called quantum hydrodynamics equations [34], named the quantum trajectory method (QTM).In their approach, the spatial support of the nuclear wave packet is split into fluid elements (FEs) that represent volume elements in configuration space carrying quantum information (amplitude and phase).Each of them is propagated according to a Newton-like equation of motion augmented by a nonlocal quantum potential.The latter supplies correlation between the FEs and is, therefore, responsible for most quantum nuclear effects.The QTM approach has been employed to address challenging quantum dynamics problems in low dimensional model systems (see [35][36][37] for an extended presentation of quantum trajectory methods).Generalizations of QTM for multi-electronic states have also been proposed [38][39][40][41][42].These are, however, based on a diabatic representation of the PESs.In an attempt to extend this type of dynamics to the investigation of molecular systems, we have recently developed an in principle exact QTM approach, named NABDY (nonadiabatic Bohmian dynamics), which solves the non-relativistic quantum dynamics of nuclei and electrons within the framework of quantum hydrodynamics, using the adiabatic representation of the electronic states [43,44].
In this article, we review a number of trajectory-based nonadiabatic molecular dynamics schemes together with our recent work on nonadiabatic Bohmian dynamics.Our aim is to provide a unified picture of the field by trying to "derive" the different approaches starting from a common framework, namely the quantum hydrodynamics reformulation of the molecular time-dependent Schrödinger equation.
In particular, we propose a classification of the different trajectory-based approaches based on the choice of the initial expansion of the molecular wavefunction (that depends on both the nuclear and the electron degrees of freedom) into a sum or a single product of electronic and nuclear wavefunctions.Finally, we propose a rationalization of the TSH equation of motion based on our exact nonadiabatic Bohmian dynamics scheme, showing by means of tests on two simple model systems the origin of some typical failures of TSH.

Nonadiabatic Dynamics with Classical and Quantum Trajectories
In this Section, we briefly review the theoretical background of the different nonadiabatic molecular dynamics schemes that we have selected for this study.The selection is based on the fact that all these trajectory-based approaches can be classified according to the way the molecular wavefunction is represented in terms of the electronic and nuclear components.
Starting from the Born-Huang representation of the total molecular wavefunction, we first introduce the Born-Oppenheimer molecular dynamics (BO-MD), which is based on the adiabatic separation of the electronic and nuclear dynamics, the latter being described by a single classical trajectory.Nonadiabaticity is then reintroduced following different strategies.In trajectory surface hopping (TSH), when the classical trajectories enter a region of strong coupling between different PESs, they are allowed to hop from one surface to another according to a hopping algorithm designed by Tully [15].An interesting improvement of this scheme consists in adding Gaussian-expanded nuclear wavefunctions to the propagating trajectories; this approach is named Full Multiple Spawning [45][46][47][48] and is characterized by an interesting balance between accuracy and numerical efficiency.Finally, we will describe a trajectory-based approach in which classical trajectories are replaced by Bohmian quantum trajectories that evolve under the influence of quantum adiabatic and nonadiabatic potentials.All these methods make use of the computationally advantageous adiabatic representation of all involved electronic states.
In the second part of this review, we discuss nonadiabatic AIMD approaches that can be derived from a single product ansatz for the total molecular wavefunction.Two of these methods will be investigated, namely the approximated Ehrenfest dynamics and the exact solution, named "Exact Factorization", which has recently been proposed by Gross and coworkers.
We begin by introducing the time-dependent Schrödinger equation (TDSE) for a molecular system, which, neglecting the nuclear and electronic spins, is given by: where Ψ(r, R, t) is the total wavefunction of the system, r = (r 1 , . . ., r k , . . ., r N el ) is the collective position vector for the N el electrons and R = (R 1 , . . ., R γ . . ., R Nn ) the corresponding one for the N n nuclei of mass M γ .The molecular Hamiltonian can be expressed in the following form: where Ĥel (r; R) is the electronic Hamiltonian, which is parametrically dependent on the nuclear coordinates.In Equation ( 2) and in the ones that follow, atomic units are used, except for the reduced Planck's constant, , that will be kept for clarity.

Methods Based on the Born-Huang Expansion
The Born-Huang expansion gives an exact expression for the total wavefunction [49,50]: The total wavefunction, Ψ(r, R, t), is expanded in the complete basis set of of electronic eigenfunctions of Ĥel (r; R), which depend parametrically on the nuclear positions, R. The expansion "coefficients", Ω i (R, t), are functions of the nuclear coordinates, R, and are explicitly dependent on time.Inserting Equation (3) into the TDSE, multiplying by Φ * j (r; R) and then integrating over r gives the equation of motion for the amplitudes, Ω j (R, t): where the F ji (R) are the elements of the nonadiabatic coupling matrix: These terms induce nonadiabatic coupling between different electronic states (for i = j) due to nuclear motion.Equation ( 4) can be interpreted as a Schrödinger-like equation for "nuclear" wavefunctions Ω j (R, t) augmented by nonadiabatic coupling terms.In fact, the amplitudes Ω j (R, t) can be interpreted as nuclear wavefunctions in state j only when the coupling terms vanish.

The Born-Oppenheimer Approximation and Adiabatic Dynamics
The BO approximation consists in neglecting all off-diagonal terms, F ji (R), in Equation (4) (i.e., neglecting inter-state couplings, but keeping intra-state electronic-nuclear couplings).The molecular wavefunction on each PES is therefore represented by the simple product Ψ(r, R, t) = Ω j (R, t)Φ j (r; R).If the diagonal terms, F jj (R), are also neglected, then we obtain what is usually called the adiabatic BO approximation [51].Introducing the polar representation of Ω j (R, t), we obtain: where both the amplitude, A j (R, t), and the phase, S j (R, t), are real.By inserting Equation ( 6) into Equation ( 4) and separating the real and imaginary parts, we obtain, within the adiabatic BO approximation, two separate, but coupled, equations for the amplitude and the phases: Taking the so-called classical limit → 0 in Equation ( 7) leads to something akin to the classical Hamilton-Jacobi equation of motion: where S j (R, t) can now be interpreted as the classical Hamilton's principal function.In this case, ∇ γ S j (R, t) corresponds to the nuclear momentum, p γ j (t), and by rearranging Equation ( 9), we finally obtain the Newtonian equation of motion for the nuclei: The nuclei, therefore, evolve on a given potential energy surface, E el j (R(t)) (selected by the initial conditions), while the electrons adiabatically follow the nuclei along their classical trajectories, R(t).Equation (8) represents a continuity equation for the nuclear amplitudes, A j (R, t), in an arbitrary state, j, which, in the classical limit, is trivially fulfilled, because of the conservation of the number of trajectories.The BO-MD method therefore consists in solving the time-independent electronic Schrödinger equation to get the potential and the forces acting on the nuclei; these are then used to propagate the nuclei for time step dt using Equation (10), and the process is iterated until the desired propagation time is reached.

Tully's Trajectory Surface Hopping
One of the most successful methods for nonadiabatic dynamics is Tully's trajectory surface hopping [14,15].In this method, the nuclei are treated classically, and the only nuclear quantum effect that is accounted for is the nonadiabatic transfer of "amplitude" between electronic states.This is achieved classically through hops of trajectories from one electronic state to another according to a hopping probability determined by the strength of the nonadiabatic couplings and the values of the state amplitudes, C [α] i (t), defined here below.A swarm of trajectories needs to be propagated in order to reproduce the probability distribution associated with corresponding nuclear quantum wave packet.
In this section, we only give a brief introduction to TSH, while a more detailed description of the method is given in Section 3, where we attempt a "derivation" of TSH, starting from the nonadiabatic Bohmian dynamics equations of motion.
The main ansatz in TSH is given by the following description of the molecular wavefunction [5,15]: which, in a way, constitutes a simplified version of the original Born-Huang expansion.When we introduce Equation ( 11) into the electronic time-dependent Schrödinger equation, we get a set of coupled equations of motion for the complex nuclear state amplitudes, C [α] where d γ ji (R) = dr Φ * j (r; R)∇ γ Φ i (r; R) are the first-order nonadiabatic couplings (see Equation ( 5)).These coupled equations will be solved along a classical trajectory, α, evolving adiabatically in a given electronic state, j.The probability, g ji (t, t + dt), for the trajectory α to jump from state j to state i during the time interval [t, t + dt] is given by: where A surface hop between two PESs, j and i (j → i), occurs "stochastically" when, for a randomly generated number, ζ ∈ [0, 1], we get: This algorithm guarantees that a minimum number of hops is performed along each trajectory; for this reason, the method is also referred to as the "fewest switches algorithm".

Full Multiple Spawning
Full Multiple Spawning (FMS) [45][46][47][48] proposes an interesting compromise between accuracy and efficiency by representing nuclear wavefunctions as sums of time-dependent Gaussian basis functions, whose width is frozen and whose center evolves adiabatically according to classical mechanics.This ansatz on the classical evolution of the Gaussian centers is consistently applied throughout the full derivation of the FMS equations of motion.
In the FMS method, the nuclear wavefunction in electronic state i, Ω i (R, t) is represented by a linear combination of multidimensional Gaussian wave packets , products of one-dimensional Gaussian functions [52]: In Equation ( 15), the multidimensional Gaussian basis functions are labeled with index J, their time-independent width by α i J and their time-dependent position, momentum and nuclear phase by R i J (t), p i J (t) and γ i J (t), respectively.N i (t) gives the number of Gaussian basis functions used in order to describe the nuclear wavefunction in electronic state i, and its time dependence comes from the possible "spawning" of new basis functions (as further discussed here below).The nuclear phases are propagated semi-classically, whereas the positions and momenta of the center of the Gaussians obey classical equations of motion in a given electronic state [52].
The time-evolution of the expansion coefficients, C i (t), is obtained through the solution of the following differential equation: This equation is derived by plugging the Born-Huang expansion (Equation ( 3)) and the ansatz for the nuclear wavefunctions (Equation ( 15)) into the time-dependent Schrödinger equation.In Equation ( 16), the bold symbols emphasize that, for each electronic state, i, there is a time-dependent coefficient per each Gaussian basis function, and S ii and Ṡii represent different overlap matrices of the Gaussian functions (see [52] for the more details).The matrix elements of H ij are given by: where Ĥel is the electronic Hamiltonian and TR the kinetic energy operator for the nuclei.
In Equation (17), The spawning procedure takes place when a region of nonadiabaticity is detected along a trajectory (by monitoring the strength of nonadiabatic couplings in the adiabatic representation) and allows for the generation of new Gaussian basis functions (children), placed in the newly populated electronic state, according to physical rules (like position or momentum conservation [52]) maximizing the coupling between the parent and children Gaussian basis functions [52], until the system leaves the region of strong nonadiabatic coupling [53] The spawning procedure, therefore, limits the number of Gaussian basis function used in the calculation by defining precisely where and when they are needed.Moreover, the FMS method offers a numerically exact [48] solution when all matrix elements are computed exactly, and a complete Gaussian basis set is used.
While keeping a trajectory-based formalism, FMS fully incorporates nuclear quantum effects that are missing in methods like TSH.Furthermore, the nuclear propagation can be performed on-the-fly, by computing any electronic structure property needed, like electronic energies (H ii el in Equation ( 17)) or nonadiabatic couplings ( Φ i |∇ R |Φ j r and Φ i |∇ 2 R |Φ j r ; note that the G iKjK terms are normally small and usually neglected) with either an ab initio electronic structure or semiempirical methods (Ab Initio Multiple Spawning, AIMS [54]).AIMS, therefore, overcomes the limitations in accuracy of TSH, preserving efficiency all the while.For additional information about the derivation and numerical procedure of this method, the interested reader is referred to [52].

Nonadiabatic Bohmian Dynamics
Just as for the previous three methods, nonadiabatic Bohmian dynamics (NABDY) is also based on the propagation of trajectories.However, this time, the trajectories evolve under the action of additional quantum potentials (adiabatic and nonadiabatic), which make the dynamics exact in principle.In other words, this approach is able to capture all adiabatic and nonadiabatic nuclear quantum effects through the propagation of a sufficiently large (i.e., converged) number of trajectories.
The derivation of the NABDY equations of motion starts from the insertion of the polar representation of the nuclear wavefunction in Equation ( 6) into Equation (4).After separation of the real and imaginary parts, we obtain: and: where φ ij (R, t) = 1 (S i (R, t) − S j (R, t)) and D γ ji (R) are the second-order nonadiabatic couplings.Equation ( 18) is equivalent to the classical Hamilton-Jacobi equation augmented by terms that are O( ) and O( 2 ).The third term, Q(R, t), on the right-hand side of Equation ( 18) is called the quantum potential, and it includes all adiabatic quantum effects (adiabatic in the sense that the potential Q(R, t) acts on a single PES and does not include contributions from other surfaces).Unlike "classical" potentials, it is non-local in space, in the sense that it depends on the gradient in configuration space [55].The last three terms on the right-hand side of Equation ( 18) describe inter-state nonadiabatic quantum effects and, like the quantum potential, Q(R, t), do not have a classical equivalent.
Equation ( 19) represents a continuity equation the for probability density, |A j (R, t)| 2 , with corresponding probability density flux J (R, t) [33,35,44].The first two terms on the right-hand side describe the "adiabatic" probability density flow within a given state, j, while the remaining terms that depend on the first-order and second-order nonadiabatic couplings induce probability density exchanges across different electronic states.Of course, the overall nuclear amplitude (summed up over all states) is conserved.
The two equations for the phases and the amplitudes are coupled, and they therefore need to be solved simultaneously.Instead of solving complex differential equations for the two fields, (A j (R, t) and S j (R, t)), we reintroduce trajectories in configuration space that drive the dynamics of "infinitesimal" volume elements called "fluid elements" [35].The derivation of the equations of motion is similar to that described in Section 2.1.1.for the BO approach, with the important difference that in NABDY, new fluid elements can be created at any time on any other PES according to the size of the nonadiabatic terms in Equation (19).The details of the numerical implementation of NABDY are given in [44], while a possible extension of NABDY to large dimensions (in the adiabatic case) is proposed in [56].

Ehrenfest Dynamics
The equation of motion that drives Ehrenfest dynamics (EHD) is derived from a simpler ansatz for the total wavefunction than the Born-Huang expansion (Equation (3)) used for the methods presented in Section 2.1.
In EHD, the molecular wavefunction is described by the simple product: where Φ(r, t) is the electronic wavefunction and Ω(R, t) is the nuclear wavefunction.Note that in this case, both amplitudes, Φ(r, t) and Ω(R, t), depend explicitly on time.In addition, they also have a parametric dependence on the other set of coordinates (Φ(r, t) on R and Ω(R, t) on r), which is not explicitly shown, so as to simplify the notation.
The exponential in Equation ( 20) is named the phase term and is defined as: It guarantees that the product wavefunction, Ψ(r, R, t), fulfills the corresponding time-dependent Schrödinger equation.
Following the derivation proposed by Tully [57], we can substitute Equation ( 20) into the TDSE, multiply by Ω * (R, t) and integrate over R, and assuming that Ω(R, t) is normalized, we finally obtain: dR where V (r, R) is the sum of all the potential energy terms in the molecular Hamiltonian.
Applying an analogous procedure, we can also derive the equation of motion for Ω(R, t): Using the relations [57]: and: where E is the expectation value of the molecular Hamiltonian for the wavefunction appearing in Equation (20), and the fact that E = E el (t) + T R ( T R is the expectation value of the nuclear kinetic energy), we can further simplify Equations ( 22) and ( 23) and obtain the following differential equations for the two amplitudes: and: These are mean field coupled equations, in which the electrons move in a field generated by the nuclei (second term on the right hand side (r.h.s.) of Equation ( 26)) and the nuclei move in a field generated by the electrons (second term on the r.h.s. of Equation ( 27)).Strictly speaking, these are not yet the EHD equations of motion, but, rather, a version of the time-dependent self-consistent field equations.EHD implies the passage to the classical limit for the nuclear amplitudes, which is again accomplished through the use of the polar representation of the nuclear wavefunction (Equation ( 6)) in Equation (27).
Once again, we obtain a classical Hamilton-Jacobi equation, which can be transformed into a Newton equation of motion given by: Notice that the potential acting on the nuclei is now given by the expectation value of the electronic Hamiltonian computed using the time-dependent electronic "wavefunction" Φ(r, t), which is not necessarily an eigenstate of Ĥel (r, R(t)), but which can be expressed as a linear combination of the static solutions of the corresponding time-independent electronic Schrödinger equation for the same nuclear position, R, at time t .For this reason, EHD is called a mean-field solution of the time-dependent molecular Schrödinger equation.The equation of motion for the electronic amplitudes, Equation ( 26), also depends on the nuclear amplitudes, Ω(R, t).However, since the nuclei are treated as classical particles, we can set: that is to say, we induce localization of the nuclear densities at a fixed position, R γ (t).By plugging Equation ( 29) into Equation ( 26) we obtain a TDSE for the electronic amplitude: where the Hamiltonian and the wavefunction both depend parametrically on the nuclear positions, which induces the coupling with the nuclear equation of motion (Equation ( 28)).As we mentioned before, in EHD, the nuclei will evolve on a single time-dependent PES, which can be expressed at any instant of time as a linear combination of all adiabatic PESs.This implies that in EHD nonadiabatic effects are taken into account through the propagation of the electronic wavefunction [58]; a perspective that is indeed very different from what is observed in the approaches derived from the Born-Huang expansion (Section 2.1).

The Exact Factorization-Based Dynamics
Recently, Gross et al. [59,60] have shown that the (exact) solution of the molecular TDSE can be factorized into the product of an electronic and a nuclear wavefunction [61] (even when the Hamiltonian includes coupling to external electromagnetic fields): Equation ( 31) might seem counter-intuitive at first, because the molecular Hamiltonian is not separable.In fact, while the factorization in Equation ( 31) can be made at any time, t, and at any position, r or R, the persistence of this kind of solution along the time propagation of the two wavefunctions is less obvious (as can be seen from the resulting evolution equations [59]).As discussed in [60], the factorization in Equation ( 31) can be justified using multivariate statistics, according to which any probability distribution (here, the square of the molecular wavefunction) can be factored into a marginal probability and a conditional probability.In this respect, it is also important to notice that Φ(r; R, t) depends (parametrically) on the nuclear coordinates, R, and there is, therefore, no loss of generality in applying Equation (31).
The factorization of Ψ(r, R, t) does not simplify, per se, the task of solving the molecular TDSE.Nonetheless, this approach has a great interpretive power, since Φ(r; R, t) and Ω(R, t) have both a clear physical meaning: they are the exact electronic and nuclear time-dependent wavefunctions, respectively.One crucial requirement for this to be true is the so-called partial normalization condition: This condition allows for the interpretation of dr|Φ(r; R, t)| 2 as the conditional probability of finding an electron in volume element, dr, at position r given a nuclear configuration, R; that is to say, |Φ(r; R, t)| 2 is an electronic probability density function.According to the standard interpretation of quantum mechanics, Φ(r; R, t) is then the corresponding electronic wavefunction.Similarly, Ω(R, t) is the marginal probability density for the nuclear position, R (marginal, and not conditional, because r is unknown), and Ω(R, t) is the corresponding nuclear wavefunction.Interestingly, just as in EHD, this factorization leads to the definition of a single time-dependent potential energy surface (because of the time-dependence of the electronic wavefunction), which, this time, is, however, exact and unique.What is lost is the picture of a time-dependent nuclear wave packet (or corresponding trajectories) evolving on an ensemble of static PESs; a picture that has provided important insights for the understanding of many photophysical and photochemical processes.The time evolution of the wavefunctions, Φ(r; R, t) and Ω(R, t), is described by two connected differential equations, which contain, besides the usual interaction terms, additional scalar and vector potential terms [59,60,62,63].

Trajectory Surface Hopping from the Nonadiabatic Bohmian Dynamics Equations
When it comes to nonadiabatic molecular dynamics, TSH is probably the most popular simulation scheme.As stated in Section 2, it relies on the description of nuclear wave packets by means of a swarm of classical trajectories.A complex coefficient, C [α] j , for each electronic state, j, is propagated along a given classical trajectory, α, according to Equation ( 12).The classical trajectory may "hop" from its current electronic state, i, to another at any point in time, and the probability that a hop to state j occurs is given by Equation ( 13) [15,[64][65][66].
In this Section, starting from Equations ( 18) and ( 19), we will present a "rationalization" of the TSH equations of motion based on the nonadiabatic Bohmian dynamics equations.
The following steps were reported in [44] and can be summarized as follows: (a) The nuclear wave packet dynamics is discretized into a swarm of classical trajectories.Within the independent trajectory approximation, the quantum potential is set to zero, Q j (R, t) = 0, ∀j and so is the divergence of the current, ∇ R • J [α] j (R, t) = 0, ∀j, ∀α.The independent trajectory approximation arises from the assumption that all adiabatic non-local terms (involving a single electronic state, j) are set to zero.This will ensure that there is no adiabatic quantum transfer of amplitude among trajectories.In the adiabatic regime (d ji (R) and D ji (R) → 0), the trajectories evolve independently from each other, and their equation of motion is obtained from the solution by characteristics of Equation (18) in the classical limit, which corresponds to the classical Newton equation of motion with classical forces −∇ R E el j (R).
(b) For the description of the nonadiabatic components of the dynamics (the three last terms in Equations ( 18) and ( 19)), we first move to a reference frame evolving along the classical trajectories for which: and: Note that due to the independent trajectory approximation, we assume that there is no amplitude exchange among the FEs propagated along the different trajectories.
Neglecting the second-order nonadiabatic couplings, D ji (R), due to their usually small size [67], we are left with an equation of motion for the phases and the amplitudes, which is equivalent to the following nuclear wavefunction time-evolution equation: where we have used the definition of the momentum operator pγ = −i ∇ γ .
(c) In the derivation of the equation of motion for the nuclear amplitude coefficients, we start by assigning delta-like wave packets (denoted as the TSH wave packet in the following) to each trajectory, α, defined as: where Ã[α] j (t) and S[α] j (t)/ are real functions representing the amplitude and the phase of the TSH nuclear wave packet at R [α] (t) in electronic state j.The function ) , localized at the position of the classical trajectory, α, is normalized to dR g λ (R − R [α] (t)) = 1 and becomes a δ-function in the limit lim λ→∞ g The total probability density of the nuclear wave packet in state , where N traj is the total number of trajectories.The independent trajectory approximation invoked in point (a) also has an important impact on the nonadiabatic component of the nuclear dynamics (due to their nonlocality; see Equations ( 18) and ( 19)).Indeed, it has the consequence that, for a given trajectory, α, the complex amplitudes, Ωλ,[α] j (R, t), for each and every electronic state, j, share the same support (localized around the instantaneous nuclear position, R, in configuration space).Said otherwise, the TSH nuclear wave packet component, g λ (R − R [α] (t)), will be the same for all electronic states of a trajectory, α, at any time, t (this is why g λ does not have an electronic state index).This is indubitably the strongest approximation made in the "derivation", since it induces "overcoherence" in the dynamics of the amplitudes, C [α] j (t), and suppresses all (nonadiabatic) decoherence effects that could occur, for example, at and after the branching of nuclear wave packets.
(d) Since we are working in the Lagrangian frame, we need only consider the explicit time-dependence of the amplitudes and phases.As a consequence, the TSH nuclear wave packet evolving in electronic state j follows the classical trajectory, α, on the support of the function, g λ (R − R [α] ) (where R [α] is the position vector in the Lagrangian frame), and is described by If we substitute Ω j (R, t) in Equation ( 35) by the form given in Equation ( 36) and then apply points (a), (b) and (d), we obtain [44]: Here, are the nuclear velocities at time t along trajectory α.Notice that Equations (37) and (38) are equivalent to the TSH equations for the complex coefficients, C [α] j (t) (Equation ( 12)), which is obtained using a polar representation of the complex coefficients, C [α] t) .We have described until now the dynamics of TSH nuclear wave packets following a single classical trajectory, α.At this point, we have to account for the fact that the nuclear dynamics in TSH is described by a "swarm" of classical trajectories that evolve according to the adiabatic and nonadiabatic components of the equation of motion (points (a) to (d)).In order to to this, we have to require that the following be maintained: at all times, for a sufficiently large number of trajectories.In Equation ( 39), (A T SH j (R, t)) 2 is computed as the density (histogram) of configuration space points in the volume element, dR, at R(t) in state j that are sampled by the ensemble of N traj trajectories: while the right-hand side is the corresponding nuclear density obtained from the corresponding quantum mechanically propagated nuclear wave packets.Note that Equation ( 39) is only valid when correlated quantum (Bohmian) trajectories are used [36].
In TSH, the balance described in Equation ( 39) is maintained (in an approximative way) through the use of the switching algorithm given in Equations ( 13) and ( 14), which can be motivated by the following considerations: (e) In the independent trajectory approximation, the nonadiabatic terms in the time-evolution of the TSH nuclear wave packets (Equations ( 37) and ( 38)) induce trajectory surface hop transitions between different states according to the probability, g ji (t, t + dt), which is a function of local variables computed along the propagation of a single trajectory.
(f ) The switching probability is obtained from quantum mechanical arguments [15] under the assumption that: where the |C j (t)| av are the norms of the coefficients defined in Equation ( 12) averaged over the ensemble of trajectories.Equation ( 41) is the internal consistency criterion described in [68].However, in practice, one replaces the |C j (t)| 2 av with the corresponding amplitudes computed along a single trajectory, |C [α] j (t)| 2 , which are the coefficients that appear in Equation (13).The reason for this modification is that in the independent trajectory approximation, one computes single trajectories, and therefore, the average over the ensemble is not available during propagation.This replacement of |C j (t)| 2 av by |C [α] j (t)| 2 is, in our opinion, more of an assumption than an approximation and remains without formal justification.
In summary, starting from the exact formulation of the nonadiabatic dynamics within the nonadiabatic Bohmian dynamics framework, we proposed a series of approximations/assumptions (points (a) to (f )) that help rationalize Tully's TSH equations of motion for the nuclear trajectories and amplitudes.In particular, the independent trajectory approximation (point (a)) implies that the amplitudes and phases associated with the classical trajectories are uncorrelated (which is also evident from the fact that the trajectories are propagated separately) and that quantum nonlocality is, therefore, lost.The assumption made in point (f ) is particularly strong, as it states that the averaged TSH population amplitude (on a given electronic state, j) taken over the ensemble of trajectories can be replaced by the corresponding amplitude, C [α] j , computed along a single trajectory, α.Furthermore, the nuclear amplitudes associated to each electronic state are evaluated strictly at the same position in space, at any time t, even though the different curvature of the PESs involved in the dynamics may drive the nuclear wave packets towards different regions in configuration space.This implies that TSH is strictly local in space and time or, equivalently, that equal-time corresponds to equal-space events, which leads to the loss of quantum mechanical nonlocality [37].This is the case, even if we allow for retardation (causality), since the TSH equations have no memory.In other words, Equation ( 12) is obtained from: with the kernel, F (t − t ), replaced by a delta function, δ(t − t ).Some implication of these approximations will be described in Section 4 for simple one-dimensional model systems.

Trajectory Surface Hopping at Work
While TSH is an elegant compromise between accuracy and efficiency for the simulation of nonadiabatic phenomena, its accuracy (either in its fewest-switches version or with additional corrections) has been challenged several times in the literature (see [22,[68][69][70][71][72] for an non-exhaustive list).Recently, a series of simple one-dimensional model systems were used to highlight potential failures of the standard TSH, even with high initial momenta [28,[73][74][75][76].The "double arch" model is composed of a couple of potential energy curves, whose shapes strongly differ in the region where they are not degenerate (−10 ≤ x ≤ 10 a.u. in Figure 1).In this model, a Gaussian wave packet launched from x = −20 a.u. with a positive initial momentum will first reach a region of strong nonadiabaticity (Figure 2, upper panel), leading to a population of both the ground state (GS) and the first excited state (S 1 ).Right after this nonadiabatic event occurs, the two potential energy curves will diverge, one exhibiting a strong positive slope (S 1 state), the other a negative one (GS).The wave packet contribution in each electronic state will therefore be spatially split and eventually recombined in a second nonadiabatic region at x = 10 a.u.(Figure 2, upper panel).The final population on S 1 after the second nonadiabatic region strongly depends on the spatial decoherence between the nuclear wave packets.However, such peculiar decoherence is hardly captured by TSH, due to the independent trajectory approximation (and other approximations discussed in Section 3).This is observed from its deviation with respect to an exact nuclear wave packet propagation in the lower panel of Figure 2. TSH in general fails qualitatively for all different initial momenta tested here, which correspond in all cases to a propagation with no back reflections.Changing the initial conditions of TSH strongly alters the final population of S 1 , but does not improve it substantially [75].On the other hand, the correlated quantum trajectories (NABDY) provide an accurate description of the nuclear wave packet propagation with only minor deviations from the exact propagation (full numerical details can be found in [44])., for different initial momenta ("TSH": initial conditions sampled from a Gaussian distribution for positions and momenta, 1,500 trajectories; "TSH * ": same initial conditions, momentum and position, for all 1,500 trajectories; "NABDY" is based on a maximum total number of 162 trajectories).The maximum total number of quantum trajectories used in NABDY is 162.We further investigate the effects of overcoherence on the TSH dynamics by means of a second model system consisting of two coupled harmonic potentials, as depicted in the upper inset of Figure 3.A swarm of trajectories (and a corresponding Gaussian wave packet for the exact propagation) is initialized in the excited state (S 1 ) at x = 0 a.u., with a positive initial momentum p 0 = 40 a.u..In this model system, a single nonadiabatic region is located at x = 10 a.u.; the initial conditions are chosen in such a way that the wave packets (and the classical trajectories) will reflect back shortly after the first transition through the nonadiabatic coupling region, recrossing, therefore, the same coupling region, a second time, with opposite velocity (for a total of two nonadiabatic events, see the lower inset of Figure 3).The GS (S 1 ) potential energy curve is represented with a continuous (dashed) black line and the nonadiabatic coupling with a blue dotted line.The initial nuclear wave packet is displayed in grey.(Lower inset) Time series of potential energies along a TSH trajectory.The trajectory is initially in S 1 , then jumps to the GS after the first coupling and, finally, hops back to S 1 after it reaches back to the coupling region.This representation highlights that the model describes two nonadiabatic events with a single nonadiabatic region.The S 1 wave packet enters the strong coupling region at x = 10 a.u.(for t < 1, 000 a.u.) and populates the GS (87%, Figure 3).This first nonadiabatic event is perfectly described by TSH (Figure 3, 1000 ≤ t ≤ 2000 a.u.).Due to the difference between the potential energy curves (slope of E el S 1 larger than the one of E el GS ), the wave packet component in the GS travels further towards positive x values, while the weak contribution in S 1 inverts the direction of its propagation and rapidly returns towards the nonadiabatic region at x = 10 a.u..In the exact propagation, there is no interference with the wave packet evolving in the GS, since the two wave packets (GS and S 1 ) are spatially separated.As for the first transition through the nonadiabatic region at x ∼ 10 a.u., the S 1 wave packet is transferred almost entirely to the other electronic state (now, the GS, Figure 3, t ≥ 3000 a.u.).On the other hand, in TSH, each independent trajectory carries a set of coherently coupled complex amplitudes (see point (c) of Section 3).When reaching the nonadiabatic coupling at x = 10 a.u.for the second time, the complex amplitudes, C S 1 (t), evolved along a given trajectory, α, in S 1 , couple coherently, because they share the same support (same position in space for any time t).This induces "overcoherence" in dynamics for the amplitudes, which leads to deviations from the exact propagation (Figure 3, 3000 ≤ t ≤ 5000 a.u.).The total population in S 1 increases back to ∼78% of the t = 0 value when the wave packet in the GS recrosses the nonadiabatic region, at t = 6000 a.u.. Some additional deviations of TSH with respect to the exact propagation are observed, and they are linked to subsequent recrossings of the nuclear wave packets.

Conclusions
The description of the nonadiabatic dynamics of molecular systems is a challenging task for theory, due to the difficulty of providing both the electronic structure of a system beyond its electronic ground state and the corresponding nuclear dynamics beyond the Born-Oppenheimer approximation.Moreover, nonadiabatic phenomena require a description of the nuclear degrees of freedom that goes beyond the classical approximation, and finally, a good compromise between accuracy and efficiency is also required when realistic molecular systems are investigated.In this article, we have summarized some of the main techniques for describing the nonadiabatic dynamics of molecular systems, namely Ehrenfest dynamics, nonadiabatic Bohmian dynamics, Multiple Spawning, the recently proposed Exact Factorization and trajectory surface hopping.We have also shown how the latter method can be rationalized starting from the "exact" nonadiabatic Bohmian dynamics equations.Trajectory surface hopping is indeed one of the most commonly applied on-the-fly trajectory-based methods to describe the dynamics of molecular systems beyond the Born-Oppenheimer approximation in the (unconstrained) configuration space.This is possible at the cost of describing the nuclear wave packet dynamics with a swarm of uncorrelated classical trajectories with the consequent banishing of all quantum (de)coherence effects.Understanding the underlying limitations of trajectory surface hopping is of foremost importance for the future improvement of the theory, and in our opinion, quantum Bohmian dynamics can give valuable contributions in this direction.

Figure 1 .
Figure1.The double arch model in the adiabatic representation.The ground state (GS) (S 1 ) potential energy curve is represented with a red (dashed) line and nonadiabatic coupling with a blue dotted line.The initial nuclear wave packet is displayed in grey.

Figure 2 .
Figure 2. Nonadiabatic dynamics for the double arch system.(Upper panel) Time series (gray scale) of the nuclear wave packet probability density, |A j (x, t)| 2 , and trajectory surface hopping (TSH) histograms for p 0 = 45 a.u.(lower panel = GS; upper panel = S 1 ).The adiabatic potential energy curves are given in red, while the nonadiabatic coupling vectors are shown in blue.(Lower panel) Deviation of the final population in S 1 from an exact nuclear wave packet propagation obtained with TSH and nonadiabatic Bohmian dynamics (NABDY), for different initial momenta ("TSH": initial conditions sampled from a Gaussian distribution for positions and momenta, 1,500 trajectories; "TSH * ": same initial conditions, momentum and position, for all 1,500 trajectories; "NABDY" is based on a maximum total number of 162 trajectories).The maximum total number of quantum trajectories used in NABDY is 162.

Figure 3 .
Figure 3. Nonadiabatic dynamics on two coupled harmonic potential energy curves.Population in the first excited state (S 1 ) along the dynamics for 3,444 TSH trajectories (green) and an exact propagation (red).(Upper inset) Schematic representation of the model.The GS (S 1 ) potential energy curve is represented with a continuous (dashed) black line and the nonadiabatic coupling with a blue dotted line.The initial nuclear wave packet is displayed in grey.(Lower inset) Time series of potential energies along a TSH trajectory.The trajectory is initially in S 1 , then jumps to the GS after the first coupling and, finally, hops back to S 1 after it reaches back to the coupling region.This representation highlights that the model describes two nonadiabatic events with a single nonadiabatic region.