Ultrafast dynamics with the exact factorization

The exact factorization of the time-dependent electron-nuclear wavefunction has been employed successfully in the field of quantum molecular dynamics simulations for interpreting and simulating light-induced ultrafast processes. In this work, we summarize the major developments leading to the formulation of a trajectory-based approach, derived from the exact factorization equations, capable of dealing with nonadiabatic electronic processes, and including spin-orbit coupling and the non-perturbative effect of an external time-dependent field. This trajectory-based quantum-classical approach has been dubbed coupled-trajectory mixed quantum-classical (CT-MQC) algorithm, whose performance is tested here to study the photo-dissociation dynamics of IBr.


I. INTRODUCTION
The theoretical description of nonequilibrium processes at the microscopic scale poses continuous challenges in many fields, such as molecular and chemical physics, condensed matter physics, and theoretical chemistry. Theory needs to be able to describe purely quantum effects such as electronic transitions [1][2][3][4][5][6][7][8][9], (de)coherence [10][11][12][13][14][15], interferences [16,17], energy relaxation [18][19][20], or phase transitions [21,22], that are often the result of complex interactions between electrons and nuclei on ultrafast time scales and of the nonperturbative effect of time-dependent external fields. While quantum mechanics is the key to unravel these processes, actual applications and studies relying on computational methods require the introduction of approximations. In the domain of quantum molecular dynamics, approximations can be intended in various way, (i) to reduce the complexity of the original problem via simplified models, (ii) to make the underlying equations of motion computationally tractable based on mathematical or physical observations, or (iii) to neglect some effects in favor of others depending on the situations.
In the present work, we will provide examples on these three strategies, limiting our study to lightinduced ultrafast phenomena in isolated molecular systems, and employing the formalism of the exact factorization of the time-dependent electron-nuclear wavefunction [23][24][25]. Using this theoretical framework, particular attention is devoted to present the procedures yielding simplified equations of motion, i.e., point (ii) above, that have made in recent years the exact factorization a suitable "tool" to perform quantum molecular dynamics simulations [4,26].
Introducing the factored form of a molecular -electronic and nuclear -wavefunction allows one to rewrite the time-dependent Schrödinger equation as coupled evolution equations for the electronic and the nuclear components of the full wavefunction. Once the problem is decomposed, the idea that has motivated and driven our efforts in the last years is to circumvent the exponential scaling of quantum mechanical equations of motion by using classical-like trajectories to mimic nuclear dynamics while maintaining a quantum viewpoint for electronic dynamics [27][28][29][30].
The exact factorization naturally lends itself for such a quantum-classical scheme because the original problem, i.e., the molecular time-dependent Schrödinger equation, is decomposed into a (single) electronic and nuclear part without invoking any approximations. Furthermore, nuclear dynamics appears to be fairly standard, in the sense that the nuclei evolve according to a (new) time-dependent Schrödinger equation where the effect of the electrons is accounted for via timedependent vector [31,32] and scalar potentials [33,34]. Being explicitly dependent on time, these potentials are able to describe effects due to electronic excited states and induced or driven by external fields [24,[35][36][37][38]. Aside from being a suitable framework for developing approximate numerical schemes, the exact factorization offers a framework to interpret, and perhaps even disentangle, complex effects: on one hand, nonadiabaticity and quantum decoherence are often strongly related to each other and affected by nuclear motion [30,[39][40][41][42][43][44]; on the other hand, nuclear quantum effects, such as interferences and branching in configuration space, can manifest themselves as consequence of the coupling to electronic dynamics, and compete with other similar effects, such as tunneling [45][46][47].
While a comprehensive discussion on these points is beyond the scope of this work, we focus here on the possibility of capturing nonadiabatic effects and quantum decoherence during photo-induced ultrafast processes, by means of the quantum-classical scheme derived from the exact factorization. For completeness, we treat at the same level spin-allowed electronic nonadiabatic transitions, induced by nuclear motion, and spinforbidden electronic transitions mediated by spin-orbit coupling [48,49]. Following recent developments [50], we present as well how to explicitly account for external time-dependent fields (only the case of a continuouswave (cw) laser will be discussed). The possibility of including quantum nuclear effects by adopting a quantum trajectory approach is also mentioned.
In order to present the exact factorization, together with its quantum-classical formulation, and to apply the theory to the study of a molecular process, we organize the paper as follows. In Section II we briefly recall the exact factorization, but we mainly devote this section to the derivation of the trajectory-based procedure, dubbed coupled-trajectory mixed quantumclassical (CT-MQC), ultimately leading to an actual algorithm for quantum-classical molecular dynamics simulations. Section III focuses on a simple molecular application by presenting the performance of CT-MQC in describing the photo-dissociation of IBr, including spinorbit coupling. A general assessment of the work done so far on the exact factorization and on future directions is presented in Section IV. Our conclusions are stated in Section V.

II. EXCITED-STATE DYNAMICS WITH THE EXACT FACTORIZATION FORMALISM
This section is devoted to the presentation of the exact factorization formalism to set the basis for the derivation of the CT-MQC algorithm.
We consider a system of interacting electrons and nuclei, including spin-orbit coupling (SOC) and an ex-ternal time-dependent field to the molecular Hamiltonian. Therefore, the system is described by the HamiltonianĤ(x, R, t) =T n (R) +Ĥ BO (x, R) + H SO (x, R) +V (x, R, t), whereT n is the nuclear kinetic energy,Ĥ BO , i.e., the Born-Oppenheimer Hamiltonian, contains the electronic kinetic energy, together with the electronic and nuclear potentials,Ĥ SO is the SOC,V is the external time-dependent potential. Electronic position-spin variables are labelled as x = [r 1 , σ 1 ], [r 2 , σ 2 ], . . . , [r N el , σ N el ], and nuclear positions are labelled as R = R 1 , R 2 , . . . , R Nn . Note that, the electronic operatorsĤ BO andV are block-diagonal in spin space, and the particular form [51][52][53] chosen forĤ SO does not affect any of the results presented below.
The solution of the time-dependent Schrödinger with χ(R, t) the nuclear wavefunction, and Φ R (x, t) the electronic conditional factor that parametrically depends on R. The partial normalization condition [54,55] guarantees that |χ(R, t)| 2 reproduces at all times the (exact) nuclear density obtained from Ψ(x, R, t). Here, the integral over x stands for an integral over the 3N el -dimensional electronic configuration space and N el sums over electronic spins. The evolution of χ(R, t) and Φ R (x, t) is given by The coupled nuclear and electronic evolution equations are derived by inserting the product form (1) of the molecular wavefunction into the TDSE and by using the partial normalization condition, as described in detail in Refs. [25,46,54,55]. In Eq. (2) the symbol ∇ ν indicates the gradient taken with respect to the positions of the ν-th nucleus and M ν are the nuclear masses. The time-dependent vector potential is a three-dimensional vector for each value of the index ν, and it is defined as [31,32,56,57] The time-dependent scalar potentials are and mediate the coupling between electrons and nu-clei, beyond the adiabatic regime. We distinguish two scalar potential contributions in order to separate the effect of the external time-dependent field. The first contribution, given by Eq. (5), will be referred to as timedependent potential energy surface (TDPES) [16,27,33,34,58,59]; the second contribution, Eq. (6), is associated to the external field. The integration over x is indicated as · x . Like the vector potential and the TDPES, the external contribution, Eq. (6), is, in general, an N-body operator, even though the external potentialV (x, R, t) (typically representing a laser pulse) is a one-body operator. The electron-nuclear coupling operatorÛ en [Φ R , χ] is [28,29,60,61] and depends explicitly on the nuclear wavefunction χ(R, t), and implicitly on the electronic factor Φ R (x, t), via its dependence on the vector potential, given by Eq. (4). The product form (1) of the molecular wavefunction is invariant under a gauge-like R, t-dependent phase transformation of the electronic and nuclear components. 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. (2) and (3) with any given choice of gauge leads to a unique solution of the TDSE, as expected.
In the next sections we describe the procedure and approximations ultimately leading to the CT-MQC equations. In Section II A we rewrite the nuclear TDSE (2) employing the so-called quantum hydrodynamic formalism [62] and devise a trajectory-based scheme to solve it by invoking a classical limit (on the nuclear degrees of freedom). In Section II B, we introduce the ex-pansion of the electronic wavefunction on a basis and we show how to solve the electronic evolution equation (3) along the classical nuclear trajectories. In Section II C we summarize the algorithm and briefly discuss the computational procedure for its implementation. Note that, recently, an attempt to solve numerically exactly Eqs. (2) and (3) has been made [63], but numerical instabilities seem to develop at the early stages of the propagation even for a simple one-dimensional model. This study confirms the need for the development of approximate schemes, as the one described below.

A. Solution based on characteristics of the nuclear equation
The polar representation of the (complexvalued) nuclear wavefunction, χ(R, t) = |χ(R, t)| exp [(i/ )S(R, t)], when inserted into the nuclear TDSE (2), yields the two coupled (real) equations [46] The evolution equation for the phase S(R, t), Eq. (8), has a Hamilton-Jacobi form, where the Hamiltonian on the right-hand side contains a kinetic term, with ∇ ν S(R, t) + A ν (R, t) being the momentum field, two "classical" time-dependent potential terms, (R, t) and v(R, t), and a quantum potential term [62,64], Since the quantum potential depends on the modulus of the nuclear wavefunction, it couples the Hamilton-Jacobi equation to the continuity equation, Eq. (9), which describes the evolution of the nuclear probability density |χ(R, t)| 2 . A quantum trajectory scheme to solve the partial differential equation (PDE) Eq. (8) has been derived in Ref. [46] and relies on the method of characteristics. The PDE is replaced by a set of ordinary differential equations (ODEs), i.e., the characteristic equations, for the "variables" R(s), t(s), S(s),P(s) = ∇S(s), ∂ t S(s) that always satisfy the original PDE, as functions of the parameter s. We indicate here the 3N n -dimensional gradient of S as ∇S(s) = {∇ ν S(s)} ν=1,...,Nn to simplify the notation. We introduce the function F (R, t, S,P, S t ) = S t + H n (P, R, t) = 0, whose total differential is dF = F R · dR + F t dt + F S dS + FP · dP + F St dS t . The quantities F R , F t , F S , FP, and F St indicate the partial derivatives of the function F with respect to its variables. If a set of "initial conditions" is chosen such that F = 0, the (nontrivial) characteristic equations representing the evolution of the variables (R, t, S,P, S t ) arė with ν = 1, . . . , N n , where the Hamiltonian H n (P, R, t) is the right-hand side of Eq. (8). The symbolsṘ ν anḋ P ν are intended as total derivatives with respect to the parameter s for each value of the nuclear index ν. The characteristic equationṫ = 1 shows that the parameter s can be identified as the physical time t. The characteristic equations for S and S t , not reported here, can be easily derived fromṖ ν , sinceP ν = ∇ ν S. Equations (11) and (12) guarantee that the vector (dR, dt, dS, dP, dS t ) is orthogonal to the gradient (F R , F t , F S , FP, F St ) of F, and thus along the characteristics dF = 0. If at any time t the characteristic ODEs are solved for any initial condition, the field S(R, t) (and ∇ ν S(R, t) as well) can be reconstructed.
We proposed different procedures to solve Eqs. (8) and (9): (i) The classical procedure [26,30,41,43,44] relies on neglecting the quantum potential Q(R, t), such that the equations decouple, in the sense that the equation for S(R, t) can be solved independently from the continuity equation. The nuclear density is reconstructed simply from the distribution of trajectories. Quantum effects such as tunnelling cannot be captured in this case.
(ii) The pseudo-quantum procedure [46,60] accounts for the quantum potential Q(R, t) in the evolution for S(R, t). However, the nuclear density is only approximated as a sum of Gaussians centered at the position of the trajectories, without solving the continuity equation (9). Tunneling can be captured, even though the fine details of the nuclear density, and thus of the quantum potential, cannot be correctly reproduced, posing issues to capture effects such as interferences.
(iii) The quantum procedure [45,65,66] is devised to solve the coupled equations for the phase S(R, t) and for the density |χ(R, t)| 2 according to Eqs. (8) and (9). This procedure is fully equivalent to the solution of the nuclear TDSE. However, so far, only a proof-of-principle study has been conducted in this direction [45], due to the evident numerical instabilities caused by the quantum potential.
Henceforth, all quantities depending on R become functions of R(t) as we solve the dynamics along the flow of characteristics. This is independent of the procedure chosen among the three possibilities just presented. Note that the term "trajectory" is used to indicate the collection of 3N n nuclear coordinates that evolve in time.
The derivation of CT-MQC follows the classical procedure, thus Q(R(t), t) = 0. In addition, after having calculated the gradient of the Hamiltonian on the right-hand side of Eq. (12), we impose the condition where we introduced the new symbol P ν =P ν + A ν for the classical momentum. Equation (14) shows that the force to be used in CT-MQC to propagate classical nuclear trajectories is derived from the time-dependent vector potential. Note that the SOC contribution to the TDPES, i.e., the second term in the definition given in Eq. (5), has been gauged away and does not appear in Eq. (14).
In Section II C we provide a more explicit expression of the classical force, using the dependence of A ν (R(t), t) on the electronic wavefunction. To this end, we describe in Section II B how to derive the evolution equation for Φ R(t) (x, t) along the nuclear characteristics.

B. Solution of the electronic equation along the characteristics
Solving the electronic PDE (3) along the flow of nuclear trajectories -the characteristics -requires to switch from the Eulerian frame to the Lagrangian frame. In the Lagrangian frame, only total time derivatives can be evaluated along the flow, that is why the symbol Φ R(t) (x, t) has to be used, rather than ∂ t Φ R(t) (x, t). This is done by applying the chain rule anywhere the partial time derivative is found.
The derivation of CT-MQC relies on the expansion of the electronic wavefunction Φ R(t) (x, t) on a basis, namely The electronic representation can be chosen in different ways, as illustrated here: (i) In the absence of an external time-dependent field and of SOC, the adiabatic representation is usually preferred [44]. The adiabatic basis is formed by the eigenstates of the Born-Oppenheimer Hamil-tonianĤ BO . Transitions among electronic adiabatic states during the dynamics, that will emerge as a result of the coupled electron-nuclear motion, are induced by the nonadiabatic coupling (NAC). NAC is referred to as kinetic coupling, consequence of the action of the nuclear kinetic energy operator on the parametric dependence of the electronic states on the nuclear positions. Quantum chemistry softwares usually provide electronic structure information in the adiabatic basis, a feature that is crucial when combining CT-MQC with different approaches to electronic structure theory.
(ii) In the presence of SOC, the spin-diabatic representation is often the chosen one [49,52,67]. Spindiabatic states can be labeled depending on their spins, e.g., singlets, or triplets, and are, thus, easily employed for interpreting spectroscopic results. The spin-diabatic basis is formed by the eigenstates of the Born-Oppenheimer Hamilto-nianĤ BO . Transitions among electronic spindiabatic states of different spin-multiplicity are induced by the SOC, whereas NAC can mediate transitions within the same spin multiplet. Alternatively, the spin-adiabatic representation can be used, formed by the eigenstates ofĤ BO +Ĥ SO . Spin-adiabatic states are combinations of different spin multiplicities and transition among them is purely of kinetic nature (NAC).
(iii) In the presence of an external time-dependent field, in particular of a laser pulse, a possible choice is the representation based on the eigenstates of the HamiltonianĤ BO (orĤ BO +Ĥ SO if SOC is present), which has been called field-diabatic representation in Ref. [68]. Transitions among electronic states are of NAC or SOC character, depending on the physical problem at hand. In addition, transitions can be induced by the field itself, respecting the spin selection rules, and are mediated by the transition dipole moment when only an external electric field is considered.
(iv) In the case of a cw laser, the Floquet diabatic representation [50,[69][70][71][72][73][74][75] appears to be a viable option especially for interpreting absorption/emission processes [35,[76][77][78][79][80][81][82] in terms of exchanges of photons between the system and the field. Floquet representation exploits the periodicity induced by the laser in the Hamiltonian [83], which is encoded in the mixed character of field-electronic states, defined as the eigenstates ofĤ BO − i ∂ t . Transitions among Floquet diabatic states are either of NAC nature, conserving the number of photons, or are induced by photon exchanges. For completeness we mention as well the Floquet adiabatic representation [50], formed by the eigenstates of H BO +V − i ∂ t , even though its applications in the field of quantum molecular dynamics simulations have not been thoroughly investigated so far.
In order to present a derivation of CT-MQC that is as general as possible we use here the spin-diabatic representation. Once the final equations are obtained, it is easy to see how they can be modified to accommodate alternative representations. In the following derivation, we do not consider the external field because including the effect of a cw laser via the Floquet formalism requires the use of a slightly more involved representation [50], whose implementation in the CT-MQC algorithm is still ongoing work. The electronic CT-MQC equation describes the time evolution of the expansion coefficients introduced in Eq. (15). The purpose of the following discussion is, thus, to present a derivation of the quantityĊ m R(t), t . However, we only describe the critical points of such derivation, since the details can be found in Refs. [19,43,44,49].
In the definition (7) ofÛ en we identify two terms: the first term has been shown [61,84,85] to be smaller if compared to the second term, and it is thus neglected. The corresponding term in the TDPES will be neglected as well to maintain gauge invariance of the exact factorization equations in their trajectory-based formulation [30,44]. The electron-nuclear coupling operator depends on the nuclear wavefunction. When its polar representation is used, the expression ofÛ en reduces tô The first term in parenthesis of Eq. (16) is the velocity of the trajectory from the characteristic equation (13), whereas the second term contains the quantum momentum [30], P ν (R(t), t) = − |χ(R(t), t)| −1 ∇ ν |χ(R(t), t)|, that has been defined in previous work on CT-MQC. The quantum momentum [86,87] induces quantum decoherence effects by tracking the spatial delocalization over time of the nuclear density (or, equivalently, of its modulus) [13,40]. Since the nuclear density at each time has to be reconstructed to evaluate its spatial derivative, and thus the quantum momentum, CT-MQC trajectories are "coupled", in the sense that they cannot be propagated independently from each other. CT-MQC equations encode some non-local information about the dynamics, which is the key to correctly capture the quantum mechanical effect of decoherence. The expression of the electronic wavefunction along the trajectories in the chosen representation is inserted in Eq. (15). Note that when the total time derivative acts on the electronic states, only the gradient term survives, because they depend on time implicitly, via their dependence on the trajectory.
In order to isolate the evolution of one of the coefficients, say C m R(t), t , one has to project the electronic equation ontoφ R α (t) x , respectively. Note that we used the superscript α to label the dependence on the trajectory R α (t).
The action of the Born-Oppenheimer Hamiltonian on its eigenstates yields the eigenvalues E α m , clearly depending on the position of the trajectory.
As consequence of the action of the gradient operator ∇ ν in Eq. (16) on the electronic wavefunction, the equation contains a term like ∇ ν C α m (t), which is approximated by neglecting the spatial dependence of the modulus of C α m (t) in favor of the spatial dependence of its phase. Therefore, ∇ ν C α m (t) (i/ )f α ν,m C α m (t) and we use a simple expression for the gradient of the phase, ; this quantity is a spindiabatic force accumulated over time along the trajectory α.
With these definitions and approximations in mind, we can rewrite the electronic evolution equation given in Eq. (3) in the spin-diabatic basis along the trajectory Some quantities in this equation depend explicitly on time, and such dependence has been indicated together with the dependence on R α (t); some other quantities only depend on time via their dependence on the trajectory, which is indicated with the superscript α.
The first diagonal term in Eq. (17) contains in square brackets a purely imaginary part, depending on the spin-diabatic energy, and a purely real part, depending on the quantum momentum. The former is responsible for the oscillating phase of C α m (t), the latter affects the modulus of C α m (t) and is, thus, source of decoherence effects. It is worth mentioning here that the norm of the electronic wavefunction is preserved along the dynamics by Eq. (17) The time-dependent vector potential appears in Eq. (17). According to its definition given in Eq. (4), it depends on the electronic wavefunction, therefore, it can be expressed in the chosen electronic basis along a trajectory R α (t) as A α Recalling the approximation used above to express the gradient of the electronic coefficients, and observing that the first term of this ex-pression accumulates over time, while the second is localized in space due to the dependence on the NAC vectors, the CT-MQC expression of the (real-valued) timedependent vector potential is The classical force of Eq. (14) used to propagate CT-MQC trajectories contains the vector potential, and its total time derivative has to be computed to give an explicit expression of the force in terms of electronic coefficients and the other electronic properties introduced in Eq. (17). The final expression in given in Section II C.

C. CT-MQC algorithm
Sections II A and II B described the procedure and the approximations used to derive CT-MQC electronic and nuclear evolution equations. Essentially, those equations are a set of ODEs, thus referred to as quantumclassical, representing a way to reformulate the quantum mechanical equations of the exact factorization, i.e., Eqs. (2) and (3). In particular, the electronic equation (17) is solved along a trajectory, thus it lends itself for an onthe-fly approach, where electronic structure properties are only computed at the instantaneous nuclear posi-tions. In turn, the nuclei evolve according to a classicallike force determined by the instantaneous state of the electrons.
Nuclear forces of Eq. (14) are determined by computing the total time derivative of the vector potential. Writing explicitly this derivative yields the following expression of the force for the trajectory R α (t) The Ehrenfest-like term (Eh) is a standard mean-field force The quantum momentum term (qm) depends on P α ν (t), namely and couples the trajectories through the presence of the quantum momentum, as described in Section II B. The SOC term contains the matrix elements of the spin-orbit Hamiltonian It is clear from this expression that it is crucial to correctly capture decoherence effects in the presence of SOC. In fact, usually, SOC is very delocalized in (nuclear) space and, in some situations, it is even considered spatially constant [88]. This feature is clearly different from the behavior of NAC, which is typically localized in the regions of avoided crossings and conical intersections. If the SOC contribution to the force, Eq. (22), goes to zero, it does so as effect of decoherence, and not because the SOC itself goes to zero. Decoherence manifests itself as the "collapse" of the electronic time-dependent wavefunction along a given trajectory to a single electronic state, with the corresponding coefficient becoming one. Norm conservation along that trajectory imposes that all other coefficients become zero, finally yielding no contribution from the SOC force as consequence of the decoherence process. It should be noted that adapting standard trajectorybased algorithms for nonadiabatic dynamics to the delocalized nature of SOC requires to revisit those algorithms if calculations are performed in the spin-diabatic basis. In fact, methods such as ab initio multiple spawning [67,89,90] and surface hopping [52,[91][92][93] rely on the fact that NACs are spatially localized. The choice between spin-diabatic or spin-adibatic basis when dealing with SOC depends on the quantum-chemistry package chosen for the electronic-structure calculations; quantum-chemistry codes typically provide electronicstructure information in the spin-diabatic basis. However, if calculations are possible in the spin-adiabatic basis, all couplings become NACs, thus spatially localized, and the algorithms can be employed in their original form. The interested reader is referred to Ref. [49] for a detailed discussion on this topic with connections to various trajectory-based approaches considering SOC.
Similarly to the expression of the force given in Eq. (19), we can rewrite the electronic evolution equation (17) aṡ where we identify the following termṡ In the absence of SOC, the term [H α SO ] ml in Eq. (24c) is not present, and the electronic equation reduces to the one in the original derivation of CT-MQC in the adiabatic representation reported in Refs. [30]. In the presence of NAC and SOC, the algorithm is dubbed generalized CT-MQC (G-CT-MQC) in the spin-diabatic basis [48]. As discussed earlier, in the presence of spinorbit effects, one can transform the coupling mediated by the SOC into NAC in the spin-adiabatic representation. In this case, [H α SO ] ml of Eq. (24c) does not appear and the electronic equation is identical with the original CT-MQC equation.
In the presence of a cw laser field, if the Floquet diabatic representation is used, the periodic time dependence induced by the external field is expressed in Fourier space, via the harmonics of the driving frequency. An additional term appears in the expression of the force of Eq. (19) , and in the expression of the evolution of the electronic coefficients of Eq.(23), namelyĊ α m (t) ext = −i −1 l V α ml C α l (t). Furthermore, the energies E α m of Eq. (24b) are shifted by a fixed amount depending on the photon energy of the corresponding harmonic. Note that the indices m, l in the Floquet picture label the electronic physical states as well as the harmonic of the driving frequency. Thus, V α ml induces transition between electronic states and between harmonics. We recall that the matrix elements of the external field in the Floquet diabatic basis are expressed in Fourier space (that is why V α ml does not depend on time). Using the Floquet diabatic representation, the algorithm is dubbed F-CT-MQC [50].
In the next section we apply the CT-MQC algorithm to study the photo-dissociation reaction of IBr. The molecule IBr is known to have strong spin-orbit effects, and thus, it is selected here to show the performance of (G-)CT-MQC in the spin-adiabatic and spin-diabatic flavors presented above.

III. PHOTO-DISSOCIATION OF IBR
In this section we simulate the photo-induced dissociation of IBr based on CT-MQC and we focus on the calculation of the branching ratio of the products.
As reported in the literature [94], a strong transition dipole moment couples the ground-state of the molecule X( 1 Σ 0 + ) to its excited electronic state B( 3 Π 0 + ). In the Franck-Condon region, in particular, this transition dipole moment is the dominant one, thus we will consider that photo-dissociation is initiated by an excitation X( 1 Σ 0 + ) → B( 3 Π 0 + ). The electronic states considered here are denoted using typical spin-diabatic labels [91,95] even though we will work in the spinadiabatic representation as well, based on the model potentials of Ref. [96]. In particular, the spin-diabatic characters B( 3 Π 0 + ) and Y(0 + ), discussed below, describe the first-excited and second-excited electronic states, respectively, in the Franck-Condon region R < 5 bohr.
The molecule IBr manifests strong SOC, which is responsible for the appearance of an avoided crossing between the first-excited and second-excited electronic states, ultimately leading to the dissociation of the molecule via two channels: within 300 fs from photoexcitation, the products of such an ultrafast photoreaction are I + Br( 2 P 3/2 ), if the molecule dissociates via the first-excited state, and I + Br * ( 2 P 1/2 ), if the dissociation takes place via the second-excited state. In the dissociation limit far away from the Franck-Condon region, the first-excited (spin-adiabatic) state acquires a Y(0 + ) (spin-diabatic) character, whereas the secondexcited (spin-adiabatic) state acquires a B( 3 Π 0 + ) (spindiabatic) character. Therefore, determining the populations of the electronic states at the end of the process allows us to determine the branching ratio of the dissociation products In our simulations we use the model Hamiltonian of Ref. [96], expressed in a spin-diabatic basis aŝ where the nuclear coordinate R is the IBr internuclear distance. The diagonal elements of the electronic Hamiltonian are the spin-diabatic potential energy curves (PECs), whereas the off-diagonal elements are the SOC, which are chosen to be constant. Following Ref. [96] we have The spin-diabatic PECs are shown in Fig. 1 as yellow dashed lines. The curves corresponding to H 1 (R) and H 2 (R) cross at R c = 6.2 bohr, and in fact, after diagonalization of the Hamiltonian in Eq. (26), an avoided crossing appears at R c . The ground state is not coupled to the excited states, as clearly shown in the Hamiltonian of Eq. (26) . However, as the transition dipole moment between X( 1 Σ 0 + ) and B( 3 Π 0 + ) is strong, an optical transition can be induced via an external field.
In this studied case, we suppose an instantaneous initial excitation. The molecule is prepared in the ground state, and the initial density is shown in Fig. 1. Since the ground state PEC is nearly harmonic around the equilibrium position at R 0 = 4.666 bohr, the initial nuclear with variance σ = 0.096 bohr. The nuclear mass used in the calculations is the IBr reduced mass M = 90023 amu. The instantaneous excitation promotes the ground state wavepacket to the first excited state without geometric rearrangements, thus the first excited state is fully populated at time t = 0, and all along the dynamics, the ground state remains non-populated because it is not coupled to the excited states. CT-MQC dynamics is performed in the spin-adiabatic and in the spin-diabatic basis considering explicitly the 3 electronic states. Initial conditions are randomly sampled in position-momentum space from the (Gaussian) Wigner distribution associated to the initial nuclear state of Eq. (30). For all calculations, we run N tr = 1000 coupled trajectories. For all trajectoires, the electronic coefficients corresponding to the initially populated state are set equal to one.
When the molecule reaches R c , population is partially transferred from the first excited state to the second excited state, and transfer is concluded just after 100 fs. In Fig. 2, the populations of the spin-adiabatic states are shown as continuous lines as functions of time, whereas dashed lines are used to indicate spin-diabatic populations. Comparison with exact quantum dynamics (thin black lines) shows that CT-MQC results reproduce extremely well the spin-adiabatic populations. As shown in Fig. 1, the spin-diabatic PECs are identical with the spin-adiabatic curves in the asymptotic regions, but they differ in the vicinity of R c . Therefore, between 50 fs and 100 fs, when the population transfer process takes place in Fig. 2, the behavior of the populations differs. However, at long time, the spin-diabatic and spin-adiabatic populations agree: the population of the first excited state (red line in Fig. 2) is the same as the population of the state Y(0 + ) in the dissociating limit (dashed cyan line in Fig. 2); the population of the second excited state (green line in Fig. 2) is the same as the population of the state B( 3 Π 0 + ) in the dissociating limit (dashed darkyellow line in Fig. 2). Therefore, the branching ratio of the dissociation products can be determined by using either set of results; numerical values are reported in Table I. exact spin-adiabatic spin-diabatic   (25). The branching ratio is also calculated based on CT-MQC simulations in the spin-adiabatic and spin-diabatic basis.
In Table I  Finally, we show the comparison between nuclear quantum wavepacket dynamics and the trajectories. In Fig. 3, the distributions of CT-MQC trajectories along the PECs that guide their evolution is superimposed to the nuclear wavepackets corresponding to the first excited state (red) and to the second excited state (green). Three snapshots at the time steps indicated in the figure are shown. CT-MQC equations are solved in the spinadiabatic (black dots) and in the spin-diabatic (colored dots) basis, thus two sets of classical results are shown at the selected time steps. A very good agreement between the two sets of CT-MQC results, as well as with exact results, is observed.
In general, based on the above observations and on previous work, we conclude that CT-MQC is a reliable and flexible algorithm for trajectory-based calculations of light-induced ultrafast phenomena, involving nonadiabatic and spin-orbit coupling. Additional developments can be envisaged, especially aiming to refine the way the coupling of the trajectories is treated in the calculation of the quantum momentum, and to implement the effect of an external time-dependent field beyond the Floquet formalism.

IV. STATE OF THE ART AND PERSPECTIVES
In the past few years, the exact factorization has been developed in different flavors to describe various kinds of systems and processes.
In its original time-dependent electron-nuclear formulation proposed by Gross and co-workers [24], studies addressed the question as to what forces drive nuclear dynamics in nonadiabatic regimes [31,33,34], ultimately leading to the derivation of CT-MQC [30]. Those studies focused on the characterization of the timedependent potential energy surface and of the timedependent vector potential of the theory in key situations manifesting strong nonadiabaticity with decoherence [33], quantum interferences [16], and conical intersections [31]. In a similar spirit, and inverting the role of electronic and nuclear coordinates, the inverse factorization has been proposed to analyze the time-dependent Schrödinger equation for electronic dynamics with nonclassical nuclei [37]. In situations of weak nonadiabaticity, instead, for instance when the nuclei move slower than the electrons, the electronic equation can be solved perturbatively [61], the Born-Oppenheimer regime being the unperturbed state of the system. This idea has been successfully applied to compute electronic flux densities "within" the Born-Oppenheimer approximation [97], the response of chiral molecules to infrared left and right circularly polarized light, known as vibrational circular dichroism [84], and to estimate nonadiabatic corrections to infrared spectra of hydrogen-based molecules [85].
As discussed previously, the exact factorization naturally lends itself to the inclusion of (classical) time-dependent external fields, for instance, laser pulses or cw lasers, to describe phenomena such as photo-induced ultrafast dynamics [43], strong field ionization [36,98], or periodically driven processes [35]. This is because the formalism already accounts for time-dependent potentials to describe the coupling between electrons and nuclei. In fact, the time-dependent potentials are modified by the presence of the external field, and have been analyzed to unravel dynamical details of dissociation [25], electron localization [37,99], charge-resonance enhanced ionization [36] in H + 2 , or to challenge singleelectron pictures [100] to describe molecules in strong lasers [38,101]. Recently, the exact factorization and CT-MQC (F-CT-MQC algorithm) have been combined with the Floquet formalism to interpret laser-driven dynamics in terms of single-or multi-photon absorption and emission processes [50].
Interesting work has been conducted in the framework of the stationary Schrödinger equation [102][103][104], by proposing the static, time-independent version of the exact factorization of the electron-nuclear wavefunction. In this context it is worth mentioning three major contributions. One is the study of the behavior of the scalar potential and of the vector potential in the presence of conical intersections and in relation to geometric and topological phases [56][57][58]. This work is oriented towards finding an answer to the question: Is the appearance of the molecular geometric phase effects within the Born-Oppenheimer approximation an artefact? The traditional molecular Berry phase is intimately connected to the requirement of infinitely slow nuclear motion. While this notion of adiabatic transport is a beautiful mathematical concept, real-world nuclei move the way move (i.e., not adiabatically) and this makes it difficult to relate the traditional molecular Berry phase to actual physical observables. By contrast, the geometric phase associated with the exact vector potential does not require any adiabaticity and, hence, it might be measurable. Even though conclusions of general validity are perhaps difficult to draw, enlightening case studies have been reported. Another topic that has been approached recently, and has been already applied to LiF, is the combination of the exact factorization formalism with density functional theory to solve the electronic equation, which is coupled to the nuclear Schrödinger equation [105,106]. Finally, a quantum electronic embedding method has been derived from the exact factorization in order to calculate static properties of a many-electron system [107,108].
Motivated by experimental and theoretical advances of the past few years in the domain of physics and chemistry in cavities, the exact factorization of the electronphoton wavefunction and the exact factorization of the electron-nuclear-photon wavefunction have been proposed. Currently, studies are conducted focusing on the analysis of the properties of the scalar potential acting on the electrons as effect of the photons [109] (reminiscent of the inverse factorization), as well as on the dynamic effect of the cavity on the electron-nuclear timedependent potential energy surface [59,110], in analogy with polaritonic and cavity-Born-Oppenheimer surfaces. Application of the exact factorization in this do-main is still preliminary, but the success of the original theory in many diverse fields strongly suggest that insightful novel developments are to be expected.
In the field of quantum molecular dynamics simulations, so far CT-MQC has proven its strength and flexibility. CT-MQC is a trajectory-based method, thus it is suitable for on-the-fly molecular dynamics calculations where electronic structure properties are determined along the dynamics only at the visited nuclear geometries [26]. The necessary electronic structure information to solve CT-MQC equations are energies, gradients and derivative couplings, along with spin-orbit coupling, when working in the spin-diabatic basis, and with the transition dipole moment, if the external time-dependent electric field is explicitly considered (beyond the Condon approximation). This information is available in standard quantum chemistry packages, thus CT-MQC can be combined with different approaches to electronic structure theory. It has been shown, in fact, that this is a viable route, employing time-dependent density functional theory and the CPMD software [111]. The coupled nature of CT-MQC trajectories, however, remains the bottleneck even when parallellization strategies are employed, as was done in the CPMD implementation. To circumvent this issue, recently [42,[112][113][114] the surface-hopping algorithm has been combined with the exact factorization including decoherence effects while maintaining (as much as possible) an independent-trajectory perspective. Further studies are ongoing to propose an algorithm/implementation able to access systems of "exper-imental" complexity. Interesting new developments are clearly envisaged, focusing on the theoretical formulation of CT-MQC, for instance aiming to capture nuclear quantum effects based on the quantum trajectory formulation of nuclear dynamics, on its applications, addressing systems of growing complexity, and on the algorithmic aspect, with the purpose to improve the computational performance of such a coupled-trajectorybased scheme.

V. CONCLUSIONS
We reported an overview of the state of the art and perspectives on the application of the exact factorization of the electron-nuclear wavefunction in the field of quantum molecular dynamics.
We showed that the CT-MQC scheme is a reliable and flexible numerical procedure to solve the exact factorization equations with the support of classical-like trajectories. CT-MQC lends itself to on-the-fly ab initio molecular dynamics calculations to simulate and interpret highly nonequilibrium processes governed by various effects. Here, we presented a general method to treat nonadiabatic ultrafast dynamics, with spin-orbit effects, all within the same formalism.
Owing to the diverse achievements obtained so far and documented here, we envisage interesting new results of the ongoing developments.