Carlotta Pieroni

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

Nonadiabatic Dynamics with Coupled Trajectories

In this paper we discuss coupled-trajectory schemes for molecular-dynamics simulations of excited-state processes. New coupled-trajectory strategies to capture decoherence effects, revival of coherence and nonadiabatic interferences in long-time dynamics are proposed, and compared to independent-trajectory schemes. The working framework is provided by the exact factorization of the electron-nuclear wavefunction, and it exploites ideas emanating from various surface-hopping schemes. The new coupledtrajectory algorithms are tested on a one-dimensional two-state system using different model parameters allowing to induce different dynamics. Benchmark is provided by the numerically exact solution of the time-dependent Schrödinger equation.

Introduction

The development of excited-state molecular-dynamics schemes and codes [START_REF]Quantum chemistry and dynamics of excited states: Methods and applications[END_REF][START_REF] Agostini | Handbook of Materials Modeling[END_REF][START_REF] Curchod | Ab Initio Nonadiabatic Quantum Molecular Dynamics[END_REF][START_REF] Agostini | Different flavors of nonadiabatic molecular dynamics[END_REF][START_REF] Richter | SHARC: ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings[END_REF][START_REF] Westermayr | Machine Learning for Electronically Excited States of Molecules[END_REF][7][START_REF] Akimov | Advanced Capabilities of the PYXAID Program: Integration Schemes, Decoherence Effects, Multiexcitonic States, and Field-Matter Interaction[END_REF][START_REF] Barbatti | Newton-X: a surface-hopping program for nonadiabatic molecular dynamics[END_REF][START_REF] Freixas | Nonadiabatic Excited-State Molecular Dynamics Methodologies: Comparison and Convergence[END_REF] is a topic which has steadily evolved for few decades now, and raises continuous interest in the theoretical 1 chemistry and chemical physics communities. Novel challenges are often presented to theoreticians following, for instance, advances in spectroscopy techniques that propose new setups to study light-matter interactions. [START_REF] Nogly | Retinal isomerization in bacteriorhodopsin captured by a femtosecond x-ray laser[END_REF][START_REF] Johnson | Local vibrational coherences drive the primary photochemistry of vision[END_REF][START_REF] Schnedermann | Evidence for a vibrational phase-dependent isotope effect on the photochemistry of vision[END_REF][START_REF] Polli | Conical intersection dynamics of the primary photoisomerization event in vision[END_REF][START_REF] Lara-Astiaso | Attosecond Pump-Probe Spectroscopy of Charge Dynamics in Tryptophan[END_REF][START_REF] Yang | Diffractive imaging of coherent nuclear motion in isolated molecules[END_REF][START_REF] Hua | Monitoring conical intersections in the ring opening of furan by attosecond stimulated X-ray Raman spectroscopy[END_REF][START_REF] Shalabney | Enhanced Raman Scattering from Vibro-Polariton Hybrid States[END_REF][START_REF] Ebbesen | Hybrid light-matter states in a molecular and material science perspective[END_REF][START_REF] Flick | Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry[END_REF] Tackling those new challenges requires that algorithms/codes evolve and adapt, aiming to extend their range of applicability and to become more efficient or more accurate.

In this context the present paper proposes new advances focusing on coupled-trajectory methods for excited-state dynamics. The working framework is provided by the exact factorization of the time-dependent electron-nuclear wavefunction. [START_REF] Abedi | Exact factorization of the time-dependent electron-nuclear wave function[END_REF][START_REF] Agostini | Quantum chemistry and dynamics of excited states: Methods and applications[END_REF][START_REF] Agostini | Ultrafast dynamics with the exact factorization[END_REF] The exact-factorization Ansatz allows one to factor the molecular wavefunction as the product of a nuclear wavefunction and a conditional electronic factor, and is then proven to be an exact rewriting of the molecular time-dependent wavefunction. Evolution equations for the two terms of the factored wavefunction can be derived from the time-dependent Schrödinger equation, [START_REF] Abedi | Correlated electron-nuclear dynamics: Exact factorization of the molecular wave-function[END_REF] and have led, in previous work by Gross and collaborators, [START_REF] Min | Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes[END_REF][START_REF] Agostini | Quantum-classical non-adiabatic dynamics: Coupled-vs. independent-trajectory methods[END_REF][START_REF] Min | Ab Initio Nonadiabatic Dynamics with Coupled Trajectories: A Rigorous Approach to Quantum (De)Coherence[END_REF][START_REF] Curchod | CT-MQC -A Coupled-Trajectory Mixed Quantum/Classical method including nonadiabatic quantum coherence effects[END_REF] to derive a coupled-trajectory mixed quantum-classical (CT-MQC) scheme for excited-state nonadiabatic dynamics. As the name of the algorithm suggests, CT-MQC describes nuclear dynamics in terms of classicallike trajectories that retain some quantumness via their coupling, which enables to naturally account for decoherence effects. [START_REF] Agostini | An exact-factorization perspective on quantum-classical approaches to excited-state dynamics[END_REF] In CT-MQC, the electrons are, instead, propagated in time quantum mechanically. Since the original derivation, CT-MQC has evolved in two main directions: Min and coworkers exploited information provided by the exact factorization and by CT-MQC to propose a new trajectory surface hopping procedure [START_REF] Tully | Trajectory Surface Hopping Approach to Nonadiabatic Molecular Collisions: The Reaction of H + with D 2[END_REF][START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] that includes decoherence effects, i.e., the DISH-XF algorithm (decoherence-induced surface hopping based on the exact factorization); [START_REF] Ha | Surface hopping dynamics beyond nonadiabatic couplings for quantum coherence[END_REF] Agostini and coworkers extended the original CT-MQC aiming to include spin-orbit coupling [START_REF] Talotta | Spin-Orbit Interactions in Ultrafast Molecular Processes[END_REF] as well as the effect of an external time-dependent driving field employing the Floquet formalism. [START_REF] Schirò | Quantum-classical nonadiabatic dynamics of Floquet driven systems[END_REF] DISH-XF relies on the use of independent trajectories in order to boost the computational efficiency of the algorithm, thus providing access to molecular systems and to their photo-induced dynamics. [START_REF] Filatov | Direct Nonadiabatic Dynamics by Mixed Quantum-Classical Formalism Connected with Ensemble Density Functional Theory Method: Application to trans-Penta-2,4-dieniminium Cation[END_REF][START_REF] Vindel-Zandbergen | A study of the decoherence correction derived from the exact factorization approach for non-adiabatic dynamics[END_REF][START_REF] Filatov | Direct Nonadiabatic Dynamics by Mixed Quantum-Classical Formalism Connected with Ensemble Density Functional Theory Method: Application to trans-Penta-2,4-dieniminium Cation[END_REF][START_REF] Filatov | Non-adiabatic dynamics of ring opening in cyclohexa-1,3-diene described by an ensemble density-functional theory method[END_REF][START_REF] Filatov | Theoretical modelling of the dynamics of primary photoprocess of cyclopropanone[END_REF][START_REF] Filatov | Design and photoisomerization dynamics of a new family of synthetic 2-stroke light driven molecular rotary motors[END_REF] The generalized CT-MQC (G-CT-MQC) scheme, instead, has been designed to treat on the same footing [START_REF] Talotta | Spin-Orbit Interactions in Ultrafast Molecular Processes[END_REF][START_REF] Talotta | Internal conversion and intersystem crossing with the exact factorization[END_REF] spin-allowed nonadiabatic electronic transitions, i.e., internal conversion, and spin-forbidden electronic transitions mediated by spin-orbit coupling, i.e., intersystem crossing; the Floquet-based CT-MQC (F-CT-MQC) scheme has been derived by combining the exact factorization with the Floquet formalism [START_REF] Sambe | Steady States and Quasienergies of a Quantum-Mechanical System in an Oscillating Field[END_REF] to treat cases where systems are periodically driven by external classical fields. [START_REF] Schirò | Quantum-classical nonadiabatic dynamics of Floquet driven systems[END_REF] Motivated by the appeal of a surface-hopping-like scheme [START_REF] Tully | Trajectory Surface Hopping Approach to Nonadiabatic Molecular Collisions: The Reaction of H + with D 2[END_REF][START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] and by the success of DISH-XF, we propose and test in this work new strategies to improve the performance, in terms of computational efficiency and accuracy, of coupled-trajectory schemes derived from the exact factorization. In particular, we introduce:

(i) the independent bundle approximation (IBA) in CT-MQC, where bundles of coupled trajectories are evolved independently from each other;

(ii) the fewest-switches coupled-trajectory trajectory surface hopping (FS-CT-TSH) scheme, where coupled trajectories are evolved adiabatically and are subject to stochastic hops between electronic states, in the spirit of surface hopping, according to the fewestswitches procedure; [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] (iii) the Landau-Zener coupled-trajectory trajectory surface hopping (LZ-CT-TSH) scheme, where coupled trajectories are evolved adiabatically and are subject to stochastic hops between electronic states according to the Landau-Zener probability. [START_REF] Landau | Zur theorie der energieubertragung i[END_REF][START_REF] Landau | Zur theorie der energieubertragung ii[END_REF][START_REF] Zener | Non-adiabatic crossing of energy levels[END_REF][START_REF] Suchan | Pragmatic Approach to Photodynamics: Mixed Landau-Zener Surface Hopping with Intersystem Crossing[END_REF][START_REF] Belyaev | Landau-Zener type surface hopping algorithms[END_REF][START_REF] Jean | Vibrational coherence effects on electronic curve crossing[END_REF][START_REF] Belyaev | Nonadiabatic nuclear dynamics of atomic collisions based on branching classical trajectories[END_REF][START_REF]Approximate switching algorithms for trajectory surface hopping[END_REF] The key ingredient of all these strategies is the coupling among the trajectories, which is fundamental to capture decoherence effects in excited-state dynamics, beyond ad hoc correction schemes. [START_REF] Akimov | Advanced Capabilities of the PYXAID Program: Integration Schemes, Decoherence Effects, Multiexcitonic States, and Field-Matter Interaction[END_REF][START_REF] Subotnik | Understanding the surface hopping view of electronic transitions and decoherence[END_REF][START_REF] Bittner | Quantum decoherence in mixed quantum-classical systems: Nonadiabatic processes[END_REF][START_REF] Bittner | Decoherent histories and nonadiabatic quantum molecular dynamics simulations[END_REF][START_REF] Jaeger | Decoherence-induced surface hopping[END_REF][START_REF] Zhu | Non-Born-Oppenheimer trajectories with selfconsistent decay of mixing[END_REF][START_REF] Shenvi | Simultaneous-trajectory surface hopping: A parameter-free algorithm for implementing decoherence in nonadiabatic dynamics[END_REF][START_REF] Jain | An Efficient, Augmented Surface Hopping Algorithm That Includes Decoherence for Use in Large-Scale Simulations[END_REF] In addition, as the numerical case studies proposed here confirm, the coupling among the trajectories results to be crucial to capture additional quantum effects such as revival of coherence [START_REF] Jean | Vibrational coherence effects on electronic curve crossing[END_REF][START_REF] Ibele | A molecular perspective on Tully models for nonadiabatic dynamics[END_REF] and nonadiabatic interferences [START_REF] Granucci | Coherent excitation of wavepackets in two electronic states. Interference effects at an avoided crossing[END_REF][START_REF] Curchod | An exact factorization perspective on quantum interferences in nonadiabatic dynamics[END_REF][START_REF] Donoso | Simulation of nonadiabatic wave packet interferometry using classical trajectories[END_REF][START_REF] Granucci | Coherent excitation of wavepackets in two electronic states. Interference effects at an avoided crossing[END_REF][START_REF] Romstad | Nonadiabatic transitions and interference in photodissociation dynamics[END_REF][START_REF] Matsuzakia | Electronic and nuclear fluxes induced by quantum interference in the adiabatic and nonadiabatic dynamics in the Born-Huang representation[END_REF][START_REF] Gu | Generalized Theory for the Timescale of Molecular Electronic Decoherence in the Condensed Phase[END_REF] due, for instance, to multiple crossings of photo-excited species through nonadiabatic regions. [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] Those phenomena have been observed recently in the so-called molecular Tully models, namely in 4-N,N'-dimethylaminobenzonitrile (DMABN) and in fulvene. [START_REF] Ibele | A molecular perspective on Tully models for nonadiabatic dynamics[END_REF] The reason for testing differ-ent hopping schemes in CT-TSH is to purse flexibility of the algorithm in handling processes of diverse nature, as, for instance, fewest switches performs better when the nonadiabatic couplings are spatially localized, 66 e.g., at conical intersections, whereas Landau-Zener can be better exploited in the presence of delocalized couplings, 67 e.g., for spin-orbit coupling.

The paper presents in detail the procedures illustrated above and is organized as follows. In Section 2 we recall briefly the exact factorization and the main ingredients for the derivation of CT-MQC. Section 3 is devoted to the analysis of the coupling among the trajectories, first derived in the framework of CT-MQC, and then applied to FS/LZ-CT-TSH. In Section 4 numerical studies are presented on a one-dimensional two-states system using different model parameters in order to induce different dynamics. The performance of IBA in CT-MQC and of FS/LZ-CT-TSH are benchmarked against the exact solution of the time-dependent Schrödinger equation. In addition, in order to point out the importance of the coupling among the trajectories, results obtained with the standard (independenttrajectory) fewest-switches trajectory surface hopping algorithm [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] are presented, without (FS-TSH) and with (FS-TSH-EDC) energy-decoherence corrections. [START_REF] Granucci | Critical appraisal of the fewest switches algorithm for surface hopping[END_REF] All algorithms employed in this work are implemented in G-CTMQC 69 freely available on GitLab under GNU Lesser General Public License (LGPL). [START_REF]Quantum chemistry and dynamics of excited states: Methods and applications[END_REF]

Exact factorization and CT-MQC algorithm

A system of interacting electrons and nuclei evolving quantum mechanically is described by the time-dependent Schrödinger equation (TDSE)

ih∂ t Ψ(r, R, t) = Tn (R) + ĤBO (r, R) Ψ(r, R, t) (1)
where Ψ(r, R, t) is the time-dependent molecular wavefunction, Tn (R) is the nuclear kinetic energy, and ĤBO (r, R) is the electronic, Born-Oppenheimer (BO), Hamiltonian containing [START_REF]Quantum chemistry and dynamics of excited states: Methods and applications[END_REF] Input files and initial conditions for the calculations on Model I performed in this work are also available.

the electronic kinetic energy and the interaction potentials. The variables r and R indicate the positions of all electrons and all nuclei, respectively. According to the exact-factorization Ansatz, [START_REF] Agostini | Different flavors of nonadiabatic molecular dynamics[END_REF][START_REF] Abedi | Exact factorization of the time-dependent electron-nuclear wave function[END_REF] the molecular wavefunction can be factored as

Ψ(r, R, t) = χ(R, t)Φ R (r, t) (2)
with χ(R, t) being the nuclear wavefunction and Φ R (r, t) being the electronic conditional factor that parametrically depends on nuclear positions. Imposing that the electronic term is normalized to one ∀ R, t, the molecular TDSE yields 24,70,71

ih∂ t χ(R, t) = Nn ν=1 [-ih∇ ν + A ν (R, t)] 2 2M ν + (R, t) χ(R, t) (3)
ih∂ t Φ R (r, t) = ĤBO (r, R) + Ûen [Φ R , χ] -(R, t) Φ R (r, t) (4)
when Eq. ( 2) is inserted into Eq. ( 1). In Eq. ( 3) the sum over the index ν is performed over the N n nuclei, with masses M ν ; the spatial derivative with respect to the nuclear positions is indicated as ∇ ν . In the coupled nuclear and electronic equations, Eq. (3) and Eq. ( 4), the time-dependent vector potential A ν (R, t), [START_REF] Curchod | On the Dynamics through a Conical Intersection[END_REF] the time-dependent potential energy surface (TDPES) (R, t), [START_REF] Abedi | Dynamical steps that bridge piecewise adiabatic shapes in the exact time-dependent potential energy surface[END_REF][START_REF] Agostini | The exact forces on classical nuclei in non-adiabatic charge transfer[END_REF] and the electron-nuclear coupling operator Ûen [Φ R , χ] 75 mediate the dynamical coupling between electrons and nuclei beyond the adiabatic regime. These quantities are defined as

A ν (R, t) = Φ R (r, t)| -ih∇ ν Φ R (r, t) r (5) (R, t) = Φ R (r, t)| ĤBO (r, R) + Ûen [Φ R , χ] -ih∂ t |Φ R (r, t) r (6) Ûen [Φ R , χ] = Nn ν=1 [-ih∇ ν -A ν (R, t)] 2 2M ν + Nn ν=1 1 M ν -ih∇ ν χ(R, t) χ(R, t) + A ν (R, t) (-ih∇ ν -A ν (R, t)) (7)
and the symbol • r stands for an integration over electronic positions. The product form (2) of the molecular wavefunction is invariant under a gauge-like (R, t)-dependent phase transformation of the electronic and nuclear components. [START_REF] Agostini | The exact forces on classical nuclei in non-adiabatic charge transfer[END_REF] Fixing the gauge means to choose an expression for such a phase or imposing a condition on the gauge fields that indirectly defines the phase. The gauge fields are the vector potential and the TDPES, since they transform as standard gauge potentials if the electronic and the nuclear wavefunctions are modified by the gauge phase. Solving Eqs. ( 3) and ( 4) with any given choice of gauge leads to a unique solution of the TDSE, as expected. Note that the sum of the first two terms in Eq. ( 6)

is gauge-invariant 73 (this will be used in Section 4 to compare exact and trajectory-based calculations).

Adopting a classical-like, trajectory-based perspective for the nuclear evolution, the coupled-trajectory mixed quantum-classical (CT-MQC) algorithm has been derived in previous work. [START_REF] Agostini | Quantum-classical non-adiabatic dynamics: Coupled-vs. independent-trajectory methods[END_REF] The idea is to rewrite the nuclear equation (3) employing the polar representation of the nuclear wavefunction, i.e., expressing it in terms of its modulus/density and phase.

Two coupled equations are obtained by separating real and imaginary parts. The equation for the phase is a Hamilton-Jacobi equation [START_REF] Agostini | Nuclear Quantum Effects in Electronic (Non)Adiabatic Dynamics[END_REF] which can be solved with the method of characteristics by introducing a set Hamilton-like ordinary differential equations to evolve positions and momenta in time. The equation for the density is a continuity equation, which we do not solve explicitly, because we suppose that for "short enough" times the spatial distribution of classical-like trajectories remains localized in regions where the nuclear probability density is large [START_REF] Talotta | Quantum trajectories for the dynamics in the exact factorization framework: a proof-of-principle test[END_REF] (this hypothesis will be tested in Section 4).

In the following we recall the basic equations and main ingredients leading to the derivation of CT-MQC equations, while we refer to Refs. [ 26,27,41] for a detailed description of the underlying approximations.

Henceforth, the term "trajectory" refers to a collection of 3N n nuclear positions that evolve in time, and all quantities depending on R in the quantum formulation acquire a dependence on the trajectory R(t).

When the nuclear wavefunction is replaced by an ensemble of N traj trajectories, it is particularly interesting to analyze the electron-nuclear coupling operator. In fact, as shown in Eq. ( 7), the (second and) leading term [START_REF] Eich | The adiabatic limit of the exact factorization of the electronnuclear wave function[END_REF] in Eq. ( 7) depends explicitly on the nuclear wavefunction. Therefore, along a trajectory R(t), the electron-nuclear coupling operator becomes

Ûen Nn ν=1 Ṙν (t) + i P ν R(t), t M ν -ih∇ ν -A ν R(t), t (8)
where Ṙν (t) if the classical velocity and P ν R(t), t has been dubbed quantum momentum. [START_REF] Min | Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes[END_REF] This expression can be derived by expressing the nuclear wavefunction in terms of its modulus |χ(R(t), t)| and phase S(R(t), t), and by defining the velocity as 76

Ṙν (t) = ∇ ν S(R(t), t) + A ν (R(t), t) M ν (9)
and the quantum momentum as

P ν R(t), t = -h∇ ν |χ(R(t), t)| 2 2 |χ(R(t), t)| 2 (10)
The quantum momentum (divided by the nuclear mass) appears as a purely imaginary correction [START_REF] Garashchuk | Energy conserving approximations to the quantum potential: Dynamics with linearized quantum force[END_REF] to the classical (real) velocity in the electron-nuclear coupling operator in the electronic evolution equation (4). It has been argued that it is responsible for inducing quantum decoherence effects on electronic dynamics, 2,29 as will be shown in Section 4. Calculation of the quantum momentum requires to reconstruct the nuclear density at all times, which can be achieved from the knowledge of the spatial distribution of the trajectories. There-fore, CT-MQC trajectories cannot be propagated independently, thus they are "coupled".

Section 3 describes the quantum momentum is more details.

Since N traj trajectories are necessary to reconstruct the nuclear distribution, each trajectory will be labeled with the symbol α. At each nuclear position along the trajectory α, the electronic wavefunction is expressed as a linear combination of adiabatic states ϕ

(m)
R α (t) (r), i.e., the eigenstates of ĤBO , namely

Φ R α (t) (r, t) = m C m (R α (t), t) ϕ (m) R α (t) (r) (11)
The index m runs over the electronic states considered in the expansion. When Eq. ( 11)

is inserted into Eq. ( 4), the electronic evolution equation for the expansion coefficients is derived

Ċα m (t) = Ċα m,Eh (t) + Ċα m,qm (t) (12)
The implicit time dependence of the coefficients via the trajectory is indicated by the index α, and since the coefficients also depend explicitly on time, this time dependence is explicitly indicated. The first (Eh) "Ehrenfest-like" term in Eq. ( 12) is

Ċα m,Eh (t) = - i hE α m C α m (t) - l Nn ν=1 Ṙα ν (t) • d α ν,ml C α l (t) (13)
where E α m is the m-th eigenvalue of ĤBO along the trajectory α, Ṙα ν (t) is the velocity of the nucleus ν along trajectory α, and

d α ν,ml = ϕ (m) R α (t) |∇ ν ϕ (l)
R α (t) r is the nonadiabatic coupling vector between states m and l along the trajectory α (we suppose that d α ν,ml are real). Note that Eq. ( 13) is the electronic equation used in the trajectory surface hopping scheme as well. [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] The second term in Eq. (12), labeled "qm" as it depends on the quantum momentum

P α ν (t), is Ċα m,qm (t) = Nn ν=1 P α ν (t) hM ν • f α ν,m -A α ν (t) C α m (t) (14)
where the time-dependent vector potential along the trajectory α, i.e., A α ν (t), appears. The term f α ν,m is the m-th adiabatic force accumulated over time along the trajectory α, namely

f α ν,m = t 0 (-∇ ν E α m ) dτ (15)
The accumulated force appears in the approximated expression of the vector potential, namely

A α ν (t) = h m,l Cα m (t)C α l (t) d α ν,ml + m |C α m (t)| 2 f α ν,m m |C α m (t)| 2 f α ν,m (16)
with the symbol Cα m (t) indicating the complex conjugated of C α m (t). As indicated in the last equality, the contribution containing the nonadiabatic coupling is neglected, [START_REF] Min | Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes[END_REF][START_REF] Agostini | Quantum-classical non-adiabatic dynamics: Coupled-vs. independent-trajectory methods[END_REF] as it is localized in space (d α ν,ml are usually spatially localized) and, thus, expected to be negligible when compared to the second contribution, which is accumulated in time.

The set of electronic equations (13) for the expansion coefficients for each α has to be solved along the corresponding trajectory, and information about all other trajectories is encoded in the quantum momentum.

The nuclear trajectory R α ν (t) for the ν-th nucleus is generated by the CT-MQC force

F α ν (t) = F α ν,Eh (t) + F α ν,qm (t) (17)
The expression of the force contains two terms, similarly to Eq. ( 12): the first is a (standard)

Ehrenfest-like term

F α ν,Eh (t) = m |C α m (t)| 2 (-∇ ν E α m ) + m,l Cα m (t)C α l (t) (E α m -E α l ) d α ν,ml (18)
the second term in Eq. ( 17) depends on the quantum momentum and on the accumulated adiabatic force, namely

F α ν,qm (t) = 2 h m |C α m (t)| 2 Nn µ=1 P α µ (t) • f α µ,m f α ν,m -A α ν (t) (19)
From Eqs. ( 12) and ( 17) it is clear that CT-MQC equations are basically Ehrenfest-like with correction terms depending on the quantum momentum. The gauge has been chosen such

that α (t) + ν Ṙα ν (t) • A α ν (t) = 0.
As mentioned above, the quantum momentum couples the trajectories, thus allowing to include non-local effects into the electronic and nuclear evolution equations. This coupling appears to be essential to capture decoherence effects, therefore, it will be further analyzed in Section 3.

Coupling among the trajectories

The phenomenon of quantum decoherence is a well-studied critical issue in trajectory-based algorithms, 8,51-57,80-84 such as trajectory surface hopping [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] and Ehrenfest dynamics, [START_REF] Tully | Mixed quantum-classical dynamics[END_REF] which are based on the independent-trajectory approximation. Methods such as ab initio multiple spawning (AIMS) [START_REF] Curchod | Ab Initio Nonadiabatic Quantum Molecular Dynamics[END_REF] and the direct-dynamics variational multi-configurational Gaussian (DD-vMCG) approach [START_REF] Worth | A novel algorithm for non-adiabatic direct dynamics using variational Gaussian wavepackets[END_REF][START_REF] Lasorne | Direct quantum dynamics using variational multi-configuration Gaussian wavepackets[END_REF][START_REF] Lasorne | Direct quantum dynamics using variational multiconfiguration Gaussian wavepackets. Implementation details and test case[END_REF][START_REF] Richings | Quantum dynamics simulations using Gaussian wavepackets: the vMCG method[END_REF] are, instead, derived from the molecular TDSE, and, as such, are able to capture quantum decoherence. Similarly, CT-MQC is derived from the (exact) nuclear (3) and electronic (4) equations of the exact factorization (2), and via the quantum momentum, it is able to capture decoherence effects. However, it should be mentioned that CT-MQC equations are derived by performing the classical limit on Eqs. ( 3) and ( 4), which means that the effect of the quantum potential is neglected in the expression of the CT-MQC force and that the continuity equation for the nuclear density is not solved along the characteristics. [START_REF] Agostini | Nuclear Quantum Effects in Electronic (Non)Adiabatic Dynamics[END_REF][START_REF] Talotta | Quantum trajectories for the dynamics in the exact factorization framework: a proof-of-principle test[END_REF] Consequently, convergence to the exact quantum mechanical results should not be expected (not even in the limit of an infinite number of trajectories).

The procedure to determine the quantum momentum has been illustrated in previous work, [START_REF] Min | Ab Initio Nonadiabatic Dynamics with Coupled Trajectories: A Rigorous Approach to Quantum (De)Coherence[END_REF][START_REF] Agostini | Nuclear Quantum Effects in Electronic (Non)Adiabatic Dynamics[END_REF] but we find instructive to discuss it here as well. In addition, in Appendix A we prove that the expression used below can be derived based on a semiclassical analysis of the nuclear wavefunction. [START_REF] Heller | Frozen Gaussians: A very simple semiclassical approximation[END_REF][START_REF] Heller | Time-dependent approach to semiclassical dynamics[END_REF] As a first step, we notice that in the region where the nonadiabatic couplings between the state m and all other states l, i.e., d α ν,ml , are identically zero, the population of state m along a trajectory α changes in time as

d |C α m (t)| 2 dt = Nn ν=1 2P α ν (t) hM ν • f α ν,m - l |C α l (t)| 2 f α ν,l |C α m (t)| 2 (20)
The expression in parentheses is the time-dependent vector potential from Eq. ( 16). In the term in square brackets we multiply f α ν,m by the normalization condition of the electronic wavefunction, which is expressed in the adiabatic basis as l |C α l (t)| 2 = 1. In addition, we suppose that the quantum momentum can be decomposed into contributions involving only pairs of states (in analogy with the nonadiabatic couplings). We thus obtain

d |C α m (t)| 2 dt = Nn ν=1 2 hM ν l P α ν,ml (t) • f α ν,m -f α ν,l |C α l (t)| 2 |C α m (t)| 2 (21)
which is, in general, non-zero. However, this is not an issue along a single CT-MQC trajectory, since it is the average over all trajectories that is expected to yield zero population transfer from state m to l when the corresponding nonadiabatic couplings are zero. Therefore, we have to impose

N -1 traj N traj α=1 d |C α m (t)| 2 /dt = 0, or equivalently 1 N traj N traj α=1 Nn ν=1 2 hM ν l P α ν,ml (t) • f α ν,m -f α ν,l |C α l (t)| 2 |C α m (t)| 2 = 0 ∀ l (22)
Due to the approximate nature of CT-MQC equations and of the quantum momentum, this condition is not automatically satisfied, and has to be enforced by a suitable choice of the quantum momentum. In particular, we impose the following condition, for each ν and for each Cartesian component appearing in the scalar product above

N traj α=1 P α i,ml (t) f α i,m -f α i,l |C α l (t)| 2 |C α m (t)| 2 = 0 ∀ l, i (23)
with i = 1, . . . , 3N n .

In "perfect decoherence" conditions for the nonadiabatic electron-nuclear problem, the effect of decoherence obtained in CT-MQC would be to collapse the electronic wavefunction along the trajectory α onto an electronic state, for instance onto state m. Consequently, the selected state is associated to a population that is equal to one, i.e., |C α m (t)| 2 = 1. Norm conservation along that trajectory implies that |C α l (t)| 2 = 0 ∀ l = m. Therefore, Eq. ( 22), as well as Eq. ( 23), are naturally satisfied in perfect decoherence conditions. If, however, the condition |C α l (t)| 2 = 0 ∀ l = m is not exactly satisfied, Eq. ( 23) has to be enforced in order to avoid spurious population transfer between electronic states in the regions where the nonadiabatic couplings are zero.

Based on these observations, the quantum momentum can be computed such that it satisfies condition (23). This is done by reconstructing the nuclear density as a sum of frozen Gaussians centered at the trajectories positions, as shown in Appendix A. The linear-like expression of the quantum momentum between states m and l is

P α ν,ml (t) = Γ α ν (t)R α ν (t) -R α ν,ml (t) (24)
However, it may not naturally satisfy condition (23). Therefore, while we compute the slope Γ α ν (t) with the analytical expression Γ α ν (t) = h σν derived in Eq. ( 35), the expression of the

y-intercept R α ν,ml (t) is modified such that it is independent of the trajectory index α R i,ml (t) = N traj β=1 Γ β i (t)R β i (t) f β i,m -f β i,l C β l (t) 2 C β m (t) 2   N traj γ=1 f γ i,m -f γ i,l |C γ l (t)| 2 |C γ m (t)| 2   -1 (25)
and the quantum momentum between states m and l becomes

P α i,ml (t) = h σ i R α i (t) -R i,ml (t) (26)
In order to recover Eq. ( 21), the y-intercept in Eq. ( 25) is chosen independent of α so as to satisfy the condition given in Eq. ( 23). Note that when the denominator in Eq. ( 25) is zero because of decoherence2 , the quantum momentum itself is set to zero. If the denominator in Eq. ( 25) is small for reasons other than decoherence, then the analytical expression of the y-intercept is used (given in Appendix A in Eq. ( 36)). We monitor this divergence by introducing a cutoff-radius ε as input parameter, and if |R i,ml (t) -R α i (t)| > ε, then Eq. ( 36) is used.

It is clear from the above equations that computing the quantum momentum along CT-MQC trajectories requires knowledge of the positions of all trajectories at all times. The quantum momentum is a spatially non-local quantity, [START_REF] Garashchuk | Energy conserving approximations to the quantum potential: Dynamics with linearized quantum force[END_REF] thus encoding some quantumness in the, otherwise purely classical, dynamics. However, calculation of the quantum momentum, due to the coupling among the trajectories, has been identified as a major bottleneck of CT-MQC when aiming at applying it to the study of complex molecular systems. [START_REF] Min | Ab Initio Nonadiabatic Dynamics with Coupled Trajectories: A Rigorous Approach to Quantum (De)Coherence[END_REF] Furthermore, as it is clear from Eq. ( 17), calculation of all nonadiabatic coupling vectors along the trajectories are required to determine the force in CT-MQC. Such calculation as well is extremely time consuming when performed on molecular systems. To circumvent these issues, the DISH-XF algorithm, discussed in the Introduction, has been proposed. [START_REF] Ha | Surface hopping dynamics beyond nonadiabatic couplings for quantum coherence[END_REF] In DISH-XF, single trajectories are propagated independently from one another and, in the vicinity of avoided crossings and conical intersections, new trajectories are created in the coupled states, similarly to a spawning procedure, thus allowing to determine the quantum momentum via the knowledge of the dynamics of the coupled trajectory. Such a procedure has been proven successful for the study of the photo-dynamics of several molecular systems. [START_REF] Ha | Surface hopping dynamics beyond nonadiabatic couplings for quantum coherence[END_REF][START_REF] Filatov | Direct Nonadiabatic Dynamics by Mixed Quantum-Classical Formalism Connected with Ensemble Density Functional Theory Method: Application to trans-Penta-2,4-dieniminium Cation[END_REF][START_REF] Vindel-Zandbergen | A study of the decoherence correction derived from the exact factorization approach for non-adiabatic dynamics[END_REF][START_REF] Filatov | Direct Nonadiabatic Dynamics by Mixed Quantum-Classical Formalism Connected with Ensemble Density Functional Theory Method: Application to trans-Penta-2,4-dieniminium Cation[END_REF][START_REF] Filatov | Non-adiabatic dynamics of ring opening in cyclohexa-1,3-diene described by an ensemble density-functional theory method[END_REF][START_REF] Filatov | Theoretical modelling of the dynamics of primary photoprocess of cyclopropanone[END_REF][START_REF] Filatov | Design and photoisomerization dynamics of a new family of synthetic 2-stroke light driven molecular rotary motors[END_REF] In the present work, we follow different strategies to improve the performance of CT-MQC, with the goal to remain as close as possible to the original idea of coupled trajectories which is the strength of the method. Furthermore, we aim to avoid dealing with (i) analysis of trajectories having different duration due to the spawning procedure, (ii) growing number of trajectories when multiple passages through the nonadiabatic region(s) occur, (iii) introducing parameters to regulate the creation of new trajectories. [START_REF] Curchod | Ab Initio Nonadiabatic Quantum Molecular Dynamics[END_REF] Concerning point (ii), additional procedures would probably be necessary to limit/reduce the number of trajectories along a simulation, as it has been recently shown in the context of AIMS. [START_REF] Curchod | SSAIMS-Stochastic-Selection Ab Initio Multiple Spawning for Efficient Nonadiabatic Molecular Dynamics[END_REF][START_REF] Ibele | Comparing (stochasticselection) ab initio multiple spawning with trajectory surface hopping for the photodynamics of cyclopropanone, fulvene, and dithiane[END_REF] Keeping these points in mind, and aiming to propose an efficient and accurate simulation scheme employing coupled trajectories, we illustrate our ideas in Section 3.1.

New coupled-trajectory nonadiabatic schemes

In the framework of coupled-trajectory schemes to nonadiabatic dynamics based on the exact factorization and on CT-MQC, the new strategies briefly presented in the Introduction are described here, and tested in Section 4 aiming to improve on computational efficiency and/or accuracy.

The first idea is to adopt in CT-MQC an approximation that is reminiscent of the independent first generation approximation used in AIMS. [START_REF] Mignolet | A walk through the approximations of ab initio multiple spawning[END_REF] We introduce the independentbundle approximation (IBA): if we want to run a simulation with N traj trajectories, then we decompose the whole ensemble into n bundles with n traj = N traj /n trajectories each, and the final observables are computed by averaging over the bundles. If n = N traj we recover the independent-trajectory approximation underlying Ehrenfest dynamics, because the quantum momentum is identically zero for a run with a single trajectory. When n bundles of trajectories are created and evolved independently from each other, a natural question arises: Is there a best way to select their initial conditions? Different sampling schemes are proposed and tested in Section 4.

The second idea is to exploit a surface-hopping scheme to drive nuclear dynamics, while maintaining the full coupling among the trajectories. In this way, each trajectory is propagated using a purely adiabatic force, rather than the full CT-MQC force of Eq. ( 17), while electronic coefficients are evolved according to the CT-MQC equation (12). Therefore, only the scalar product between the nonadiabatic coupling vector and the nuclear velocity along a trajectory in Eq. ( 13) needs to be computed for the evolution. As in standard surface hopping, the active state, or the force state, is chosen stochastically at each time step. Note that the quantum-momentum term (14) in the electronic evolution equation (12) shall not be viewed as a decoherence correction to surface hopping, since it has been derived from the electronic equation ( 4) in the framework of the exact factorization. Therefore, the populations of the electronic states can be, and shall be, computed as the squared moduli of the electronic coefficients evolved according to Eq. (12). It is important to stress that with this procedure, which we refer to as CT-TSH, i.e., coupled-trajectory trajectory surface hopping, the trajectories are intended as a support to sample the nuclear configuration space and the electronic energy landscape (see discussion in Section 3.2). Thus, various strategies for the hopping scheme can be envisaged, and similar outcomes are, indeed, expected.

Using adiabatic forces to propagate the coupled trajectories in CT-TSH, rather than the full CT-MQC force (17), is expected to improve the stability of the numerical procedure at long times, especially when multiple nonadiabatic crossings are observed, as the situations tested in Section 4 or those observed in previous work. [START_REF] Pieroni | Relaxation dynamics through a conical intersection: Quantum and quantum-classical studies[END_REF][START_REF] Marsili | Quantum and quantum-classical studies of the photoisomerization of a retinal chromophore model[END_REF] Even though the trajectories do not carry phase information, the electronic coefficients of Eq. ( 12) do (approximately), and this affects the force (17). In long-time dynamics, and in the presence of interferences and revivals of coherence, however, phase information is difficult to capture accurately, and the accumulated errors in the phases produce noisy (vector and scalar) potentials, [START_REF] Curchod | An exact factorization perspective on quantum interferences in nonadiabatic dynamics[END_REF] and thus, a CT-MQC force which can be numerically unstable (as observed in Refs. [95,96]). Therefore, employing CT-TSH completely circumvents the issue.

In the following section we will test the performance of the (original) fewest-switches (FS) formula [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] and of the Landau-Zener (LZ) formula [START_REF] Suchan | Pragmatic Approach to Photodynamics: Mixed Landau-Zener Surface Hopping with Intersystem Crossing[END_REF][START_REF] Belyaev | Landau-Zener type surface hopping algorithms[END_REF] as hopping schemes, and we denote the two procedures as FS-CT-TSH and LZ-CT-TSH, respectively. Since information about electronic population is encoded in the electronic coefficients used to expand the (exact) electronic wavefunction in the adiabatic basis, Eq. ( 11), even in the Landau-Zener scheme, the electronic evolution equation (12) has to be solved along the trajectories. In the fewestswitches scheme, the hopping probability along the trajectory α from state m to state l at time t is given by 97

P FS α,m→l (t) = max   0, - 2 Cα m (t)C α l (t) ν d α ν,lm • Ṙα ν (t) |C α m (t)| 2 dt   ( 27
)
which depends on the time step dt, whereas in the (single-hopping) Landau-Zener scheme the hopping probability [START_REF] Suchan | Pragmatic Approach to Photodynamics: Mixed Landau-Zener Surface Hopping with Intersystem Crossing[END_REF] is

P LZ α,m→l = exp   - 2π h (H α ml ) 2 ν Ṙα ν (t) • ∇ ν (H α mm -H α ll )   (28)
where H α ml and H α mm (or H α ll ) are the off-diagonal and diagonal elements of the electronic Hamiltonian in diabatic representation at the position of the trajectory α. Equation (28) is given using diabatic quantities, which will be available for the model systems used in the numerical studies presented in Section 4. Landau-Zener schemes employing adiabatic quantities have been proposed as well, [START_REF] Suchan | Pragmatic Approach to Photodynamics: Mixed Landau-Zener Surface Hopping with Intersystem Crossing[END_REF][START_REF] Belyaev | Landau-Zener type surface hopping algorithms[END_REF] therefore, extending our approach to more general situations where the diabatic quantities are not available is indeed possible. It is worth mentioning that the Landau-Zener procedure can be extended to multiple dimensions. [START_REF] Belyaev | Nonadiabatic nuclear dynamics of the ammonia cation studied by surface hopping classical trajectory calculations[END_REF] In FS-CT-TSH and LZ-CT-CTSH, at each time step the hopping probability is calculated and compared to a random number ξ uniformly distributed between 0 and 1 to decide if a hop to a new state occurs, exactly as it is done in standard surface hopping. [START_REF] Barbatti | Nonadiabatic dynamics with trajectory surface hopping method[END_REF] Further computational details and information about the treatment of frustrated hops and velocity rescaling after a successful hop are given in Appendix B.

Rationalization of CT-TSH

In the new coupled-trajectory schemes proposed in Section 3.1 some elements are borrowed from a standard surface-hopping procedure while others emanate from the exact factorization. In this section, our aim is to clarify these points, justifying why we expect a better performance of CT-TSH over CT-MQC and over surface hopping, as will be documented in Section 4 based on numerical tests.

First, extensive analysis [START_REF] Agostini | An exact-factorization perspective on quantum-classical approaches to excited-state dynamics[END_REF][START_REF] Abedi | Dynamical steps that bridge piecewise adiabatic shapes in the exact time-dependent potential energy surface[END_REF][START_REF] Agostini | The exact forces on classical nuclei in non-adiabatic charge transfer[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] of the properties of the TDPES (6) of the exact factorization confirmed that nuclear dynamics is basically driven by adiabatic potentials, which is a common feature with surface hopping. The TDPES, in fact, (i) develops often adiabatic shapes, (ii) with steps appearing in regions of small nuclear density, and bridging between regions of different (adiabatic) slopes, and (iii) it becomes "more diabatic" in the regions where adiabatic surfaces are close or degenerate. Therefore, the dynamics produced by the TDPES can be replaced/approximated by purely adiabatic dynamics (see point (i)), as it is done in surface hopping, mimicking the effect of the diabatic shape of the TDPES (see point (iii)) via the stochastic hops in the regions of strong nonadiabaticity. The steps (see point (ii)) are missed within a surface-hopping-like procedure. However, being associated to regions of low nuclear probability density, the error does not strongly affect the final averages when a large number of trajectories are considered.

Second, the time-dependent vector potential (5) [START_REF] Curchod | On the Dynamics through a Conical Intersection[END_REF][START_REF] Agostini | The exact forces on classical nuclei in non-adiabatic charge transfer[END_REF]100 couples to the nuclear kinetic energy in Eq. (3). As it has been discussed in Ref. [ 74], its effect is to modify the nuclear momentum after the TDPES has developed the steps, namely when different portions of the nuclear wavepackets propagate "on" different adiabatic surfaces and are, thus, associated to different potential energies. This effect is reproduced as well in a surface-hopping procedure via the (discontinuous) velocity rescaling after a hop has occurred.

Third, the major weaknesses of surface hopping are: the electronic equation, which is not derived from the molecular TDSE (1), and the independent-trajectory approximation.

Both issues are overcome when employing the electronic evolution equation derived from (the exact) Eq. ( 4) rather than the usual surface hopping equation (13), which requires accounting for the coupling among the trajectories via the quantum momentum.

Fourth, since the hopping scheme in CT-TSH is merely a way to reproduce the effect of the diabatic shape of the TDPES at the nonadiabatic crossing, we dot not expect that different prescriptions for the hopping probability yield substantially different results. In particular because, as mentioned above, in CT-TSH the trajectoires become a support to sample the nuclear configuration space and the energy landscape. It follows that electronic populations shall not be estimated by counting the trajectories in each state, but rather from the average of the quantum electronic populations computed by solving the electronic equation (12). This last point is particularly important: in CT-TSH we witness a paradigm shift if compared to standard surface hopping. In fewest-switches surface hopping, it is the electronic equation that is merely a support for the propagation of the trajectories, because it is required to compute the hopping probability (and in Landau-Zener surface hopping, solving the electronic equation is not required at all); in CT-TSH, instead, the opposite is true. However, in some sense, according to the justifications presented above relating the TDPES and the vector potential to some of the ingredients of a surface-hopping algorithm, CT-TSH trajectories can be used for physical interpretation of the real nuclear dynamics.

Nonadiabatic dynamics with coupled trajectories

In this section we present test studies of nonadiabatic dynamics in one-dimesional, in nuclear space, two-level systems, [START_REF] Jean | Vibrational coherence effects on electronic curve crossing[END_REF]101,102 i.e., with two electronic states, employing the coupledtrajectory schemes presented in Section 2 and Section 3, namely CT-MQC (without and with IBA), FS-CT-TSH and LZ-CT-TSH. In addition, in order to evaluate the effect of the coupling among the trajectories, we include results based on surface hopping using the standard independent-trajectory approximation. For this, we use the fewest-switches algorithm as hopping scheme; results obtained without considering any decoherence correction are referred to as FS-TSH, whereas in FS-TSH-EDC we employ energy-decoherence corrections. [START_REF] Granucci | Critical appraisal of the fewest switches algorithm for surface hopping[END_REF] The BO Hamiltonian is given in the diabatic basis as In order to create different physical situations, we use different initial conditions in the three models. In Model I, the initial nuclear wavepacket starts in the electronic excited state centered at R 0 = R 2 , so as to reproduce the effect of a short laser pulse that photo-excites the system and, thus, initiates its (weak) nonadiabatic dynamics. The same initial condition is used for Model II, however, the resulting dynamics will be different since the model is strongly nonadiabatic. In both cases, the nuclear coordinate R can be interpreted as a reaction coordinate that guides the de-excitation of the system. In Model III, the nuclear wavepacket starts in the ground state and is centered in R 0 = -2.5 a 0 . Also in this case, 20 the resulting dynamics is nonadiabatic, and models, for instance, a proton-coupled electron transfer process 102,103 initiated from the non-equilibrium position of the proton. The initial nuclear wavefunction is a real normalized Gaussian centered in R 0 wit variance σ = 0.223 a 0

ĤBO (R) =    H 11 (R) H 12 (R) H 12 (R) H 22 (R)    =    1 2 k (R -R 1 ) 2 b exp [-a(R -R 3 ) 2 ] b exp [-a(R -R 3 ) 2 ] 1 2 k (R -R 2 ) 2 + ∆    (29
χ(R, t = 0) = 4 1 πσ 2 exp - (R -R 0 ) 2 2σ 2 (30)
The used nuclear mass is M = 20000 m e in all models. Reference calculations are obtained by solving the TDSE (1) in the diabatic basis. Quantum-classical results are obtained by using • CT-MQC: N traj coupled nuclear trajectories are evolved according to Eq. ( 17) and electronic coefficients are evolved according to Eq. ( 12);

• CT-MQC with IBA: n independent bundles of n traj coupled nuclear trajectories are evolved according to Eq. ( 17) and electronic coefficients are evolved according to Eq. ( 12);

• FS-TSH: N traj independent nuclear trajectories are evolved according to adiabatic forces and electronic coefficients are evolved according to Eq. ( 13); the hopping probability is given by Eq. ( 27);

• FS-TSH-EDC: N traj independent nuclear trajectories are evolved according to adiabatic forces and electronic coefficients are evolved according to Eq. ( 13); the hopping probability is given by Eq. ( 27); at each time step, the coefficients associated to the non-active states are damped exponentially over a characteristic decoherence time (see Appendix B) to induce decoherence;

• FS-CT-TSH: N traj coupled nuclear trajectories are evolved according to adiabatic forces and electronic coefficients are evolved according to Eq. ( 12); the hopping probability is given by Eq. ( 27);

• LZ-CT-TSH: N traj coupled nuclear trajectories are evolved according to adiabatic forces and electronic coefficients are evolved according to Eq. ( 12); the hopping probability is given by Eq. ( 28).

All algorithms are implemented in G-CTMQC 69 In what follows, the dynamics in the three models is monitored by looking at the time trace of the electronic excited-state population, at the time trace of the electronic indicator of coherence, and at nuclear distributions. Appendix B provides information about how the observables shown below are computed, as well as computational details for the dynamics simulations.

Let us comment in detail on the dynamics in Model I. The evolution starts with the nuclear wavepacket in the excited state S 1 , and after less than 30 fs (1250 h/E h ) it experiences a nonadiabatic population transfer to S 0 as it goes through the avoided crossing at R ac . Model I is weakly nonadiabatic, thus the main portion of the nuclear wavepacket remains in S 1 when moving in the region R > R ac . Owing to the fact that the S 1 PEC is very steep, the S 1 wavepacket is rapidly reflected back towards R ac . The second nonadiabatic event occurring at about 75 fs (3125 h/E h ) involves the reflected S 1 wavepacket which transfers once again population onto S 0 . This process creates an additional portion of the S 0 wavepacket moving in the region R < R ac , while the portion previously created -during the first nonadiabatic event -is being reflected back to R ac . At about 130 fs (5417 h/E h ), the S 1 wavepacket, once again reflected, from the left wall of the S 1 PEC and moving towards R ac , encounters the (first) S 0 wavepacket that has been reflected by the right wall of the S 0 PEC: they interfere when they meet at the avoided crossing, and we thus observe interstate interferences. Afterwards, at about 150 fs (6250 h/E h ), two wavepackets in S 0 meet in the region R < R ac and interfere, but such intrastate interferences clearly do not alter the electronic population evolution.

Just before the end of the simulated dynamics, at about 175 fs (7300 h/E h ) the S 0 and the S 1 wavepackets interfere once again when going through the avoided crossing.

As mentioned in the Introduction, the revivals of coherence and the interferences just discussed, which are observed as well in Model II and Model III, are not unrealistic quantum effects, and thus artefacts of the one-dimensional models studied here. In fact, recently, similar features have been reported in the relaxation dynamics of photo-excited molecules, such as DMABN and fulvene, which have been studied with AIMS and FS-TSH(-EDC). [START_REF] Ibele | A molecular perspective on Tully models for nonadiabatic dynamics[END_REF] The nonadiabatic population transfers described above for Model I will be described in detail in Section 4.1 (clearly identified in Figure 2). The frequent intrastate and interstate interferences might make it challenging for trajectory-based methods to correctly capture electronic dynamics. On the one hand, after every passage through the avoided crossing, the system experiences decoherence, as different portions of the nuclear wavepackets, in S 0 and in S 1 , evolve nearly independently. However, revivals of coherence occur every time population is transferred from one state to the other. On the other hand, if the S 0 and S 1 wavepackets meet at the avoided crossing, they interfere while exchanging population.

The simulated dynamics for Model II is very similar to the dynamics in Model I described above, but it is more classical, in the sense that interstate interferences occur only during the third nonadiabatic event at about 120 fs (5000 h/E h ). However, at around 175 fs (7300 h/E h ), quantum dynamics shows intrastate interferences between two portions of the S 1 wavepacket occurring at R ac while interstate interferences with S 0 occur as well.

In Model III, dynamics starts in S 0 , and since it is weakly nonadiabatic (similarly to Model I) only a small portion of electronic population is transferred to S 1 at about 45 fs (1875 h/E h ). Later on, the S 1 population slightly decreases at 105 fs (4375 h/E h ) when the excited wavepacket is reflected back towards R ac by the right wall of the S 1 PEC. At 150 fs the dynamics shows a new increase in S 1 population, when the S 0 wavepacket passes through R ac after being reflected by the right wall of the S 0 PEC. At this time, the S 1 wavepacket that is being created by S 0 interferes with another S 1 wavepacket which is reaching R ac .

Detailed analysis of electronic and nuclear dynamics for the models illustrated above is reported in Sections 4.1 and in Section 4.2, respectively, allowing to assess the performance on the new coupled-trajectory algorithms proposed here.

Electronic observables

In Figure 2 We show Model I in the top row, Model II in the middle row, and Model III in the bottom row. CT-MQC results with N traj = 1000 coupled trajectories in the left panels (blue lines) are compared to CT-MQC results obtained by producing n = 10 independent bundles of n traj = 100 coupled trajectories. Initial conditions are ordered according to their total energy (E-ordered, in red), according to their kinetic energy (T-ordered, in orange), or randomly organized (random, in green). FS-TSH results with N traj = 1000 trajectories are shown as light-blue and fuchsia lines, and FS-TSH-EDC results with the same number of trajectories are shown as purple and light-green lines. The label "cl." means that the population is estimated from the number of trajectories evolving on S 1 at time t, whereas "q." means that the average over the trajectories of the squared modulus of the electronic coefficient for S 1 at time t is shown. FS/LZ-CT-TSH results are obtained by running N traj = 1000 coupled trajectories using the fewest-switches or the Landau-Zener scheme. Populations (labelled "q.") are shown as fuchsia and as light-green lines, for FS-CT-TSH and LZ-CT-TSH, respectively. For completeness, the ratio of trajectories evolving in S 1 at each time step (labelled "cl.") is shown as light-blue lines and purple lines, for FS-CT-TSH and LZ-CT-TSH, respectively. that, if the initial sampling for the n traj trajectories in each bundle is biased via E-or Tordering, populations results improve even if compared to the fully couple CT-MQC (with The population of S 1 computed based on FS-TSH and FS-TSH-EDC is shown in the central panels of Figure 2 and is compared to exact results. In surface-hopping schemes, the population of the electronic states at time t has to be estimated as the ratio of trajectories evolving in each state at that time, and the total number of trajectories. If internal consistency holds, then such "classical" estimate agrees with the "quantum" estimate, which is given by the average over the trajectories of the squared moduli of the electronic coefficients propagated according to the surface-hopping electronic equation ( 13). In the central panels of Figure 2, we notice that for all models FS-TSH results for classical (cl.) and quantum (q.) populations are different (light-blue and fuchsia lines), whereas they agree when the energy-decoherence correction in FS-TSH-EDC is used (purple and light-green lines). Such energy-decoherence correction is, in fact, designed to restore surface-hopping internal consistency.

When the coupled-trajectory approach is combined with surface hopping, the procedures FS-CT-TSH and LZ-CT-TSH are derived, as discussed above. The populations of S 1 thus predicted are shown in the right panels of Figure 2 for the three models, and are compared to exact results. As discussed in Section 3, electronic population of a given state has to be estimated via the squared modulus of the corresponding expansion coefficient from Eq. ( 11), averaged over the trajectories. In fact, the distribution of classical trajectories "on" the PECs is only used to sample the energy landscape and, thus, the nuclear configuration space, since the trajectories are not propagated according to the "exact force". Therefore, by contrast with the standard formulation of surface hopping, it is the "quantum" estimate of the population which should be the reference (indicated as fuchsia and light-green lines in the right panels of Figure 2), not the "classical" estimate (indicated as light-blue and purple lines in the right panels of Figure 2). We find anyway interesting to compare the two populations, both in the fewest-switches procedure and in the Landau-Zener procedure, and to observe satisfactory agreement. Note that, as clarified in Appendix B, LZ-CT-TSH trajectories are propagated for shorter times than FS-CT-TSH. The comparison of LZ/FS-CT-TSH with exact results shows an improvement in predicting electronic populations over independent-trajectory surface hopping as well as over CT-MQC.

In summary, comparing quantum-classical results with exact calculations in Figure 2 we observe the following trends:

• In general, CT-MQC results agree very well with the reference, and in some cases, as for instance in Model III, IBA results slightly improve on the fully coupled approach.

At long times, CT-MQC results lose accuracy, which might be explained by the fast oscillations that develop in the potentials driving nuclear dynamics (see Figures 5,6 and 7 below), as consequence of interferences. Those oscillations are not correctly captured by the trajectories, thus their evolution deviates from the expected behavior as the simulation proceeds.

• For all models, decoherence corrections are essential to restore internal consistency in FS-TSH-EDC. However, only in Model II, after the decoherence correction is applied, surface-hooping populations converge to the correct results. It is interesting to notice, in fact, that when the model is weakly nonadiabatic, i.e., Model I and Model III, the energy-decoherence corrections act too strongly on the system, and FS-TSH-EDC populations deviate from the reference. The situation does not improve when using other values of the parameter that tunes the effect of the correction or when the results are averaged over different hop histories (results not shown).

• FS/LZ-CT-TSH populations remain at all times very close to the exact results, improving even on CT-MQC. Perhaps, as discussed above, this happens because the error in the calculation of the full force due to interference effects, from Eq. ( 17), is circumvented by propagating the trajectories adiabatically. In addition, and as expected, the hopping scheme does not seem to have an effect on the electronic populations (at least for times at which the comparison between fewest switches and Landau-Zener is shown).

In Figure 3 we show an indicator of coherence between the electronic states S 0 and S 1 , which, according to its definition given in Appendix B, is related to the overlap between nuclear wavepackets associated to the two electronic states. Therefore, every time a wavepacket passes through the nonadiabatic region and transfers amplitude to the other state, the indicator of coherence increases; coherence, instead, decays when the wavepackets separate in space, as consequence of the different slopes of the adiabatic PECs that drive their dynamics.

Analysis of the indicator of coherence shown in Figure 3 suggests that:

• Even though in CT-MQC with IBA, n bundles of coupled trajectories are propagated independetly from each other, decoherence is well-described by all sampling schemes, with small differences depending on the sampling. In general, decoherence effects applied via the quantum momentum are not as strong as they should be, and in fact the indicator of decoherence does not decay to zero as in exact calculations.

• As expected, without applying decoherence corrections, FS-TSH yields over-coherent results, and after the first crossing through the nonadiabatic coupling region, coherences are created and do not decay to zero when the trajectories leave the coupling region. Energy-decoherence corrections definitely cure such over-coherence, even though the corrections appear to be too strong.

• The effect of the coupling of trajectories and of the quantum momentum in FS/LZ-CT-TSH improves the results over CT-MQC and FS-TSH-EDC. Probably, the more stable nuclear dynamics with respect to CT-MQC (discussed above) yields a more stable quantum momentum acting on the electronic populations. Nevertheless, we can find similar features as in CT-MQC results, namely a slower decay of the coherence indicator if compared to exact results, even though it is not as severe as in CT-MQC.

With the aim to analyze the effect of decoherence induced by the quantum momentum and, thus, by the coupling among the trajectories, in (fully coupled) CT-MQC and in FS-CT-TSH, we show in Figure 4 results for three selected trajectories, for Model I in the left panels, for Model II in the central panels, and for Model III in the right panels. For each model, as functions of time, we plot: the active (or force) state in FS-CT-TSH along the selected trajectory (light-blue circles), which is either zero for S 0 or one for S 1 ; the population of S 1 along that trajectory, evolved according to FS-CT-TSH (fuchsia lines), which should be equal to one when the active state is S 1 if internal consistency holds valid; the population of S 1 along the trajectory evolved according to CT-MQC (blue lines) which has the same initial conditions, i.e., position and momentum, as the chosen FS-CT-TSH trajectory.

In FS-CT-TSH, and it is the case in LZ-CT-TSH as well, the quantum momentum acting in the electronic evolution equation ( 12) correctly accounts for decoherence effects by inducing a decay of the electronic coefficients/populations that are not associated to the active the active (or force) state in FS-CT-TSH along the selected trajectory (light-blue circles), which is either zero for S 0 or one for S 1 ; the population of S 1 along the selected trajectory, evolved according to FS-CT-TSH, which should be equal to one when the active state is S 1 (fuchsia lines); the population of S 1 along the trajectory evolved according to CT-MQC which has the same initial conditions (position and momentum) as the chosen FS-CT-TSH trajectory (blue lines). state. A clear example is shown for Model II in the central panel of the top row of Figure 4: just before 30 fs (1250 h/E h ), the active state switches from S 1 to S 0 , and, similarly, the electronic population of state S 1 switches from 1 to 0; later on, just after 120 fs (5000 h/E h ), the active state is S 1 again, and the electronic population of state S 1 switches from 0 to 1.

For a trajectory with the same initial conditions and propagated according to CT-MQC, the electronic population of state S 1 follows very closely the FS-CT-TSH population.

In some cases, the three quantities shown in Figure 4 are in reasonable good agreement with each other (see for instance the central panel in the top row just described). However, sometimes, CT-MQC results deviate from FS-CT-TSH, as is the case shown in the left panel of the top row. There, the quantum momentum probably acts differently towards the end of the simulated dynamics, since the whole ensemble of trajectories behaves differently in CT-MQC and in FS-CT-TSH. Furthermore, even though the two trajectories have the same initial conditions, it is very likely that their dynamics deviates during the simulation, as they are propagated according to different forces. In particular, the stochastic "ingredient" in FS-CT-TSH might be a major source of deviation. Note, instead, that the active state often agrees with the electronic population, suggesting some degree of internal consistency of the surface-hopping procedure along a single trajectory. However, it is worth recalling that, in FS-CT-TSH, the effect of the quantum momentum in the electronic evolution equation (12) is not derived by imposing internal consistency. Also, internal consistency is not required along a single trajectory, as a single trajectory completely loses physical meaning within such a coupled-trajectory scheme.

Nuclear observables

In order to finalize the assessment of the coupled-trajectory procedures, let us now turn towards the analysis of nuclear dynamics. To this end, we compare the quantum nuclear distribution, i.e., the quantum nuclear density, and the TDPES, both obtained from exact calculations, with the spatial distribution of the trajectories. In CT-MQC, we can

compute α GI (t) = m |C α m (t)| 2 E α m
, which is the CT-MQC-approximate expression of the gauge-invariant part of the TDPES given in Eq. ( 6). This same quantity can be computed in FS-CT-TSH. In addition, in FS-CT-TSH as well as in FS-TSH, we can look at the active state as function of the trajectories positions at some given time step. It is interesting, then, to compare the active states to the approximation of the TDPES and to the exact TDPES.

All those comparisons are presented in Figure 5 for Model I, in Figure 6 for Model II, and in Figure 7 for Model III, at three selected time steps along the simulated dynamics (the selected time steps are all the same for the three models and are indicated in the panels of the figures). Note that for these comparisons we analyze FS-TSH rather than FS-TSH-EDC because, in general, the used energy-decoherence correction seems to produce population dynamics with larger deviation from exact results than the uncorrected dynamics.

For Model I, shown in Figure 5, we observe that trajectories obtained based on IBA follow closely the fully-coupled CT-MQC trajectories (left panels). This is indeed a desired result, as it shows that using independent bunches provides a reasonable and accurate strategy to accelerate the simulations. Such an observation remains valid for Model II and Model III, as it is clearly shown in Figures 6 and7 (left panels).

Independent-trajectory surface-hopping dynamics, i.e., FS-TSH, reproduces reasonably well the quantum nuclear distribution. However, in the most challenging situation represented by Model I, already at intermediate times, i.e., at t = 136.8 fs (5700 h/E h ) shown in the central panel of the middle row of Figure 5, the trajectories miss the portion of the nuclear density localized in the vicinity of R = 2 a 0 . There, the exact TDPES shows that the nuclear density is mainly associated to S 0 , since the TDPES lies on the S 0 PEC. Instead, the trajectories are exclusively found on S 1 . The disagreement is even more severe at later times, i.e., at t = 172.8 fs (7200 h/E h ) in Figure 5, as the portion of the density in the vicinity of R = -2 a 0 is completely missed by FS-TSH trajectories.

The deviations between the TDPES and the active state in FS-TSH are completely eliminated when using FS-CT-TSH, since often the active state agrees with the shape of the exact and of the approximated TDPES. Nonetheless, some disagreement can still be observed at early times, at t = 74.4 fs (3100 h/E h ) in Figure 5 (right panel of the top row): this happens because when crossing the nonadiabatic region, surface-hopping trajectories are forced to hop (or not), but since in Model I the system is weakly nonadiabatic, the typical steps in the exact TDPES that bridge between adiabatic shapes appear slowly, [START_REF] Abedi | Dynamical steps that bridge piecewise adiabatic shapes in the exact time-dependent potential energy surface[END_REF][START_REF] Agostini | The exact forces on classical nuclei in non-adiabatic charge transfer[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] only some time after the nuclear density has left the nonadiabatic region.

Model II manifests strong nonadiabaticity, and all trajectory-based results agree pretty well with each other and with exact calculations, as shown in Figure 6. It is interesting to notice that, at the intermediate time t = 136.8 fs (5700 h/E h ), the approximate TDPES of CT-MQC (left panel of the middle row) is particularly noisy in the region around the avoided crossing. This behavior is due to the fact that the trajectories just left the coupling region after having "interfered" with each other (see discussion above on the dynamics in Model II).

The active state in FS-TSH and FS-CT-TSH at the same times, shown in the central and right panels of the middle row, respectively, capture much better the shape of the TDPES.

The approximate TDPES in FS-CT-TSH (right panel of the middle row) is definitely more neat than the CT-MQC approximation.

At the selected time steps in Model III, shown in Figure 7, we observe good agreement among all trajectory-based simulation methods. The main deviation appears at t = 136.8 fs (5700 h/E h ) for FS-TSH (central panel of the middle row in Figure 7). There, in fact, a larger amount of trajectories are found in S 0 around R = -2 a 0 than in S 1 around R = 2 a 0 .

Looking at the nuclear density, the inverse behavior is expected, and in fact, at that time, the S 1 population predicted by FS-TSH, shown in Figure 2, largely deviates from exact results.

The correct distribution is, instead, obtained at the same time for FS-CT-TSH (right panel of the middle row of Figure 7). Another deviation from the exact TDPES is observed at time t = 172.8 fs (7200 h/E h ) for FS-CT-TSH results (right panel in the bottom row), as the step developing approximately at R = 1 a 0 is not correctly captured. This issue is probably related to the fact that the quantum momentum induces decoherence effects in a smoother/softer way than what is expected from exact calculations. Such an observation can be confirmed by looking at the indicator of coherence in Figure 3, which shows that while for the quantum systems the indicator is zero, for the FS-CT-TSH system it has a (small) finite value.

Conclusions

In this work we proposed and tested different coupled-trajectory strategies to simulate excited-state dynamics based on the exact factorization of the electron-nuclear wavefunction.

In particular, we focused on capturing challenging quantum effects related to decoherence, revivals of coherence, and nonadiabatic (interstate) interferences. We proved, based on numerical case studies, the importance of the coupling among the trajectories to capture those effects.

In CT-MQC, i.e., the coupled-trajectory mixed quantum-classical algorithm derived from the exact factorization, since trajectories are propagated simultaneously, the computational cost rapidly increases with the number of trajectories and with the number of electronic states included in the simulation (especially due to the calculations of the nonadiabatic coupling vectors among all states for determining the force). Therefore, as first strategy, we introduced, tested and validated the suitability of the independent-bundle approximation, in order to propagate independently bundles of (a smaller number of) coupled trajectories, a procedure which can be easily and trivially parallelized. As a second strategy, we looked into combining the coupled-trajectory idea of CT-MQC with surface hopping, employing the fewest-switches and the Landau-Zener probabilities for the stochastic hops. Since in surface hopping nuclear trajectories are propagated adiabatically, these two new schemes allow one to circumvent the expensive calculation of the nonadiabatic coupling vectors of the CT-MQC force, while maintaining a coupled-trajectory perspective.

Combining ideas borrowed from surface hopping and from the exact factorization performs surprisingly well, thus we believe that the algorithms developed in this work will be able to open new avenues for stable, efficient and accurate nonadiabatic molecular-dynamics simulations of systems of growing complexity. In particular, the encouraging results obtained based on such a novel coupled-trajectory surface-hopping scheme motivate further investigation aiming, for instance, to include spin-orbit coupling or the effect of an external time-dependent perturbation and overcoming the limitations of current methods.

Gaussians and the third term from the spatial derivative of the classical action in the cosine function. In the numerators of the second and third line of the above expression, products of Gaussians centered at different positions, i.e., R β (t) and R γ (t), appear. In order for such products to be not negligible, the position of the trajectory R α (t), at which the quantum momentum is computed, shall lie close to the centroid position R c ν (t) = R β ν (t)+R γ ν (t) 2

(when the Gaussians have the same width, as is the case here). The centroid is the position where the Gaussian functions have the same value. Therefore, if R α ν (t) R c ν (t), the numerator of the second term goes to zero and, thus, can be neglected if compared to the first term.

In fact, the first term includes in the numerator only products of Gassians centered at the same position. The third term is negligible as well, if one considers the following relation: S β (t) = S[R β (t)] = S[R c (t) + ∆ β (t)] (and similarly for S γ (t) with ∆ γ (t) = -∆ β (t)). Taylorexpanding the action around the centroid position, one gets that the third term of Eq. ( 33) is small if R α (t) is close to R c (t), which is the point where the product of two Gaussianscentered at different positions -is reasonably different from zero. Therefore, we showed here that the first term in Eq. ( 33) is the leading one, based on the representation of the nuclear wavefunction as a combination of complex frozen Gaussians. Our simulations confirmed these analytical findings.

For algorithmic simplicity, since the diagonal numerator is the leading term in the expression of the quantum momentum, we neglect as well all terms β = γ in the denominator of Eq. (33). We have, thus, justified that, when using a semiclassical-like expression of the nuclear density, the quantum momentum can be written as In surface-hopping-based calculations, after a hop to a new state occurs, the velocity of that trajectory is rescaled along the direction of the corresponding nonadiabatic coupling to enforce energy conservation; in the case of frustrated hops, no action is taken to modify the velocity of the trajectory. In FS-TSH-EDC calculations, the coefficient that is not associated to the active state is exponentially damped over a typical decoherence time as prescribed in

P α ν (t) = Γ α ν (t)R α ν (t) -R α ν (t) (34)
Ref. [68] (the value of the used parameter for the decoherence correction is C = 0.1 E h ).

In LZ-CT-TSH the hopping probability is computed (only once) as soon as a trajectory R α (t) satisfies simultaneously a distance criterion, i.e., |R α (t) -R ac | < 0.1 a 0 , and an energy criterion, i.e., |E α S 0 -E α S 1 | < 0.04 E h . For each passage of the trajectory ensemble through the coupling region, centered at R ac , we enforced that each trajectory "decides" only once if a hop occurs or not. However, during the simulated dynamics, the trajectories cross several times the coupling region and, as time goes on, this can happen for different trajectories simultaneously. Therefore, we limit the propagation of LZ-CT-TSH trajectories to a shorter time than FS-CT-TSH. Even if over a shorter time frame, the aim of the comparison between LZ-CT-TSH and FS-CT-TSH is achieved, that is, we proved numerically that different hopping schemes yield similar results when coupled trajectories are employed within a surface-hopping procedure.

Electronic populations and the indicator of coherence along the dynamics are computed

The indicator of coherence η S 0 S 1 (t) between the states S 0 and S 1 relates to the overlap between the wavepackets χ S 0 (R, t) and χ S 1 (R, t), which explains why every time a nuclear wavepacket goes through the nonadiabatic region centered at R ac , a revival of coherence is observed. In order to define an indicator that can be evaluated in the same way using quantum properties and trajectory-based properties, we compute

η S 0 S 1 (t) = dR | χS 0 (R, t)χ S 1 (R, t)| 2 |χ(R, t)| 2 (40)
in the quantum calculations, and

η CT-MQC S 0 S 1 (t) = 1 N traj N traj α=1 Cα S 0 (t)C α S 1 (t) 2 (41)
in CT-MQC. The same expression is used in FS/LZ-CT-TSH. Even though in FS-TSH(-EDC) the electronic coefficients do not provide physical information about the dynamics, we use the same expression to estimate η FS-TSH S 0 S 1 (t) and η FS-TSH-EDC S 0 S 1 (t).

  ) which shows that the diabatic potential energy curves (PECs) are parabolas displaced in position and in energy. The values of the parameters in atomic units are: k = 0.020 E 2 h m e /h 2 , a = 3.0 a -2 0 , R 1 = 6.0 a 0 , R 2 = 2.0 a 0 , R 3 = 3.875 a 0 . The remaining parameters, i.e., b and ∆, are used to define three different situations: b = 0.01 E h and ∆ = 0.01 E h in Model I (weak nonadiabatic coupling in an asymmetric double-well potential), b = 0.005 E h and ∆ = 0.01 E h in Model II (strong nonadiabatic coupling in an asymmetric double-well potential), and b = 0.01 E h and ∆ = 0 in Model III (weak nonadiabatic coupling in a symmetric double-well potential). Figure 1 represents schematically Model I along with some of the processes analyzed below.

Figure 1 :

 1 Figure 1: Schematic representation of the model defined in Eq. (29). Dark-orange and darkgreen lines represent the adiabatic PECs of the ground state (S 0 ) and of the excited state (S 1 ). The parameters ∆ and b are tuned to produce different nonadiabatic situations; in the figure, the parameters and initial condition defining Model I are shown. Nuclear wavepackets are shown: before the photo-excitation (t < 0), just after the photo-excitation (t = 0), and at a later time (t = 48 fs). Coupled trajectories (indicated by filled circles with yellow shaded contours) follow the nuclear wavepackets, and at intermediate times (t = 25 -30 fs) they hop between the adiabatic PECs. They remain associated to one of the two states (dark-orange or dark-green contours, depending on the state) but the corresponding values of the electronic populations have intermediate values between 0 and 1 (filled circles with orange-to-green transition). The arrows indicate various processes: upward gray arrow for the instantaneous photo-excitation (from t < 0 to t = 0); downward light-green arrow (from t = 0 to t = 25-30 fs) for the wavepacket/trajectories sliding down S 1 ; rightward light-green and light-orange arrows (from t = 25 -30 fs to t = 48 fs) for the wavepacket/trajectories moving on S 1 and S 0 , respectively; light-green and light-orange arrows (for t > 50 fs) for wavepacket/trajectories turning back on S 1 and S 0 , respectively.

1 (

 1 and computational details are given in Appendix B. The adiabatic PECs and nonadiabatic coupling used in the calculations are obtained by diagonalizing the Hamiltonian of Eq. (29). The PECs are labeled S 0 (ground state) and S excited state), and present an avoided crossing at R ac (R ac = R 3 in Model I and Model II; R ac = 4.0 a 0 in Model III). When the nuclear wavepacket and the trajectories cross the nonadiabatic region, electronic population transfer between S 0 and S 1 occurs. Owing to the simplicity of the models, in the present study we can challenge the performance of the coupled-trajectory and independent-trajectory schemes by looking at long-time dynamics, which manifest coherence revivals and interferences as the result of multiple passages of the nuclear wavepacket/trajectories through the nonadiabatic region. If we estimate a characteristic time from the oscillation frequency in the harmonic wells of Eq. (29), which is approximately 25 -30 fs, then "long-time dynamics" means ∼ 200 fs.

  we show the population of the excited state S 1 calculated based on CT-MQC (left panels), FS-TSH(-EDC) (central panels) and FS/LZ-CT-TSH (right panels) for the three model systems studied here, namely Model I (top row), Model II (middle row) and Model III (bottom row).For CT-MQC, we propagate N traj = 1000 coupled trajectories (blue lines in the left panels of Figure2), and we compare those results to IBA, where n = 10 bundles of n traj = 100 trajectories each are propagated independently. The outcome of each independent run is then averaged over n. IBA results slightly change depending on the sampling procedures to select initial positions and momenta, which are described in Appendix B; in the left panels of Figure2, the different strategies are indicated as: red lines, when the initial conditions are ordered depending on their total energy (E-ordered); orange lines, when the initial conditions are ordered depending on their kinetic energy (T-ordered); green lines, when the initial conditions are randomly initialized in each bundle (random). It is interesting to observe -TSH cl.FS-CT-TSH q. LZ-CT-TSH cl. LZ-CT-TSH q.Exact vs FS/LZ-CT-

Figure 2 :

 2 Figure 2: S 1 population along the simulated dynamics: exact results (black lines) CT-MQC (left panels), FS-TSH(-EDC) (central panels), and FS/LZ-CT-TSH (right panels).We show Model I in the top row, Model II in the middle row, and Model III in the bottom row. CT-MQC results with N traj = 1000 coupled trajectories in the left panels (blue lines) are compared to CT-MQC results obtained by producing n = 10 independent bundles of n traj = 100 coupled trajectories. Initial conditions are ordered according to their total energy (E-ordered, in red), according to their kinetic energy (T-ordered, in orange), or randomly organized (random, in green). FS-TSH results with N traj = 1000 trajectories are shown as light-blue and fuchsia lines, and FS-TSH-EDC results with the same number of trajectories are shown as purple and light-green lines. The label "cl." means that the population is estimated from the number of trajectories evolving on S 1 at time t, whereas "q." means that the average over the trajectories of the squared modulus of the electronic coefficient for S 1 at time t is shown. FS/LZ-CT-TSH results are obtained by running N traj = 1000 coupled trajectories using the fewest-switches or the Landau-Zener scheme. Populations (labelled "q.") are shown as fuchsia and as light-green lines, for FS-CT-TSH and LZ-CT-TSH, respectively. For completeness, the ratio of trajectories evolving in S 1 at each time step (labelled "cl.") is shown as light-blue lines and purple lines, for FS-CT-TSH and LZ-CT-TSH, respectively.

Figure 3 :

 3 Figure 3: S 0 /S 1 indicator of coherence. The color code and the organization of the panels is the same as in Figure 2.

Figure 4 :

 4 Figure4: Analysis of decoherence for three selected trajectories for Model I (left panels), Model II (central panels), and Model III (right panels). The plots show, as function of time: the active (or force) state in FS-CT-TSH along the selected trajectory (light-blue circles), which is either zero for S 0 or one for S 1 ; the population of S 1 along the selected trajectory, evolved according to FS-CT-TSH, which should be equal to one when the active state is S 1 (fuchsia lines); the population of S 1 along the trajectory evolved according to CT-MQC which has the same initial conditions (position and momentum) as the chosen FS-CT-TSH trajectory (blue lines).

Figure 5 :

 5 Figure5: Comparison of the exact TDPES (only the gauge-invariant part is shown as black lines) with: the TDPES estimated with CT-MQC using fully coupled trajectories (blue circles in the left panels) and using n bunches of E-ordered trajectories (red dots in the left panels); the active state in FS-TSH dynamics (light-blue circles in the central panels); the active state in FS-CT-TSH dynamics (light-blue circles in the right panels) and the TDPES estimated with FS-CT-TSH (fuchsia dots in the right panels). The plots are shown at the times indicated in the panels for Model I. For reference, the adiabatic PECs are shown, S 0 in dark-orange and S 1 in dark-green, along with the exact nuclear density at the corresponding times (think black lines).

Figure 6 :

 6 Figure 6: Comparison analogous to Figure 5 but for Model II. The color code is the same as in Figure 5.

Figure 7 :

 7 Figure 7: Comparison analogous to Figure 5 but for Model III. The color code is the same as in Figure 5.

  Wigner transform of the initial nuclear probability density used in quantum calculations, i.e., the squared modulus of Eq.(30). The same initial conditions are used in all simulations. In FS-TSH(-EDC) and in FS-CT-TSH we run as well 100 trajectories 10 times with different seeds for the random number generator so as to sample different histories of jumps. The results thus obtained do not deviate significantly from the previous ones. The same initial conditions used in fully-coupled CT-MQC are organized in n = 10 independent bundles of n traj = 100 coupled trajectories, and are ordered according to their total energy (E-ordered) or according to their kinetic energy (T-ordered). In addition, n = 10 bundles of randomly organized n traj = 100 coupled trajectories are propagated as well (random).

We recall that when using a trajectory-based scheme like surface hopping or CT-MQC, the effect of decoherence is reflected on the electronic wavefunction via its "collapse" along a trajectory α to a single adiabatic state, say m, whose population, i.e., |C α m (t)| 2 , becomes equal to one. Norm conservation along that trajectory imposes that all other populations, for l = m, become zero.

A way to circumvent point (i) would be to use the "child" trajectories only as auxiliary trajectories. However, this procedure would use computational resources without extracting actual information.

Acknowledgement

This work is supported by a public grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program, through the "ADI 2019" project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02.

A. Semiclassical analysis of the nuclear wavefunction

In this Appendix we describe how to derive the expression of the quantum momentum based on the semiclassical, frozen-Gaussian representation [START_REF] Heller | Frozen Gaussians: A very simple semiclassical approximation[END_REF][START_REF] Heller | Time-dependent approach to semiclassical dynamics[END_REF] of the nuclear wavefunction of Eq. (2).

At the position R α (t) of the trajectory α, the value of the nuclear wavefunction can be approximated as

where the symbol

stands for the value of a normalized Gaussian function (with constant width σ) centered at R β (t), and S β (t) is the classical action computed along the trajectory β. Usually, the frozen width is chosen as the width of the initial distribution. Starting from this expression the nuclear density |χ(R α (t), t)| 2 can be computed as

Applying the gradient -h∇ ν to the nuclear density of Eq. ( 32), and dividing by (twice) the nuclear density to get the quantum momentum, we find

where the first two terms on the right-hand side follow from the spatial derivative of the and

Equation ( 34) looks like a linear function where Γ α ν (t) is the slope and R α ν (t) is the yintercept. On the one hand, the slope given in Eq. ( 35) is clearly independent of time and of the trajectory index, and is fully determined by the values of the variances of the frozen Gaussians in Eq. ( 31). In the original derivation of CT-MQC, [START_REF] Min | Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes[END_REF][START_REF] Agostini | Quantum-classical non-adiabatic dynamics: Coupled-vs. independent-trajectory methods[END_REF][START_REF] Min | Ab Initio Nonadiabatic Dynamics with Coupled Trajectories: A Rigorous Approach to Quantum (De)Coherence[END_REF] the slope was chosen time dependent and adapted along the dynamics. However, to stabilize the dynamics and reduce the number of parameters to be chosen for a CT-MQC simulation, in Ref. [95] we introduced the frozen-Gaussian expression (35). The expression of the y-intercept R α ν (t), on the other hand, will be modified so as to satisfy the condition (23), stating that no population transfer between two electronic states should be observed, when averaged of the trajectories, if the corresponding nonadiabatic couplings are zero.

B. Computational details

Quantum calculations are carried out by solving the TDSE in the diabatic basis using the split-operator technique. 104 The used time step is dt = 0.0024 fs (0.1 h/E h ) and the spatial grid is defined in the range R ∈ [-5.0, 15] bohr with 2000 grid points for all models.

Trajectory-based simulations are performed by solving nuclear equations with the velocity-Verlet algorithm and the electronic equations with the Runge-Kutta-Gill algorithm, as implemented in G-CTMQC. 69 The time step is dt = 0.0024 fs (0.1 h/E h ) as in quantum calculations, for both nuclear and electronic dynamics. For the assessment of the stability of CT-MQC, we tested as well dt = 0.012, 0.024, 0.048 fs (0.5, 1.0, 2.0 h/E h ) and we did not observe signifiant deviations.

We employ N traj = 1000 in fully-coupled CT-MQC, FS-TSH(-EDC) and LZ/FS-CT-TSH. Initial positions and momentum are sampled from the (harmonic) quantum distribution in different ways depending on the simulation method. Using the index m = S 0 , S 1 to label the electronic states, we have that exact populations are obtained as

where χ m (R, t) is the nuclear wavepacket associated to the electronic state m; it is computed by transforming the diabatic wavepackets, which are evolved when solving the TDSE, into adiabatic wavepackets along the quantum dynamics. In CT-MQC and FS/LZ-CT-TSH this same observable is determined from the average over the trajectories of the squared moduli of the electronic coefficients

whereas in FS-TSH(-EDC) we monitor the number of trajectories N m (t) associated to each electronic state m as a function of time

The quantity defined in Eq. ( 38) has been referred to as "quantum population" in the main text, whereas the quantity defined in Eq. ( 39) is the "classical population". Clearly, in CT-MQC, only the quantum population, thus Eq. ( 38), is accessible, because in this procedure one cannot identify the active state. Similarly, as we stressed in the main text, electronic populations in FS/LZ-CT-TSH have to be estimated according to the quantum expression (38).

By contrast, in a pure surface-hopping scheme, where the electronic equation is not derived from the exact equation ( 4), the classical expression is to be employed, thus Eq. (39).

In order to assess the internal consistency of the surface-hopping scheme, usually achieved when decoherence corrections are used, as in FS-TSH-EDC, the classical and the quantum electronic populations can be compared.