# Ultrafast dynamics with the exact factorization

Federica Agostini, E. Gross

► To cite this version:

Federica Agostini, E. Gross. Ultrafast dynamics with the exact factorization. The European Physical Journal B: Condensed Matter and Complex Systems, 2021, 94 (9), pp.179. 10.1140/epjb/s10051-021-00171-2. hal-04294574

HAL Id: hal-04294574

https://universite-paris-saclay.hal.science/hal-04294574v1

Submitted on 20 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Ultrafast dynamics with the exact factorization

Federica Agostini<fnref n="1" />,ᵃ and E. K. U. Gross<fnref n="2" />

Received 1 April 2021 / Accepted 19 July 2021 / Published online 9 September 2021

© The Author(s) 2021

**Abstract.** 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.

## 1 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–9], (de)coherence [10–15], interferences [16,17], energy relaxation [18–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 non-perturbative 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 light-induced ultrafast phenomena in isolated molecular systems, and employing the formalism of the exact factorization of the time-dependent electron-nuclear wavefunction [23–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 using classical-like trajectories to mimic nuclear dynamics while maintaining a quantum viewpoint for electronic dynamics [27–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 time-dependent 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–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–44]; on the other hand, nuclear quantum effects, such as interferences and branching in configuration space, can

<footnote n="1">¹ Institut de Chimie Physique UMR8000, Université Paris-Saclay, CNRS, 91405 Orsay, France</footnote>

<footnote n="2">² Fritz Haber Center for Molecular Dynamics, Institute of Chemistry, The Hebrew University of Jerusalem, Jerusalem 91904, Israel</footnote>

<footnote>ᵃ e-mail: federica.agostini@universite-paris-saclay.fr (corresponding author)</footnote>
manifest themselves as consequence of the coupling to electronic dynamics, and compete with other similar effects, such as tunneling [45–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 spin-forbidden 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 continuous-wave (cw) laser will be discussed). The possibility of including quantum nuclear effects by adopting a quantum trajectory approach is also mentioned.

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 Sect. 2, we briefly recall the exact factorization, but we mainly devote this section to the derivation of the trajectory-based procedure, dubbed coupled-trajectory mixed quantum-classical (CT-MQC), ultimately leading to an actual algorithm for quantum-classical molecular dynamics simulations. Section 3 focuses on a simple molecular application by presenting the performance of CT-MQC in describing the photo-dissociation of IBr, including spin-orbit coupling. A general assessment of the work done so far on the exact factorization and on future directions is presented in Sect. 4. Our conclusions are stated in Sect. 5.

## 2 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 external time-dependent field to the molecular Hamiltonian. Therefore, the system is described by the Hamiltonian $\hat{H}(\mathbf{x}, \mathbf{R}, t) = \hat{T}_n(\mathbf{R}) + \hat{H}_{BO}(\mathbf{x}, \mathbf{R}) + \hat{H}_{SO}(\mathbf{x}, \mathbf{R}) + \hat{V}(\mathbf{x}, \mathbf{R}, t)$, where $\hat{T}_n$ is the nuclear kinetic energy, $\hat{H}_{BO}$, i.e., the Born–Oppenheimer Hamiltonian, contains the electronic kinetic energy, together with the electronic and nuclear potentials, $\hat{H}_{SO}$ is the SOC, $\hat{V}$ is the external time-dependent potential. Electronic position-spin variables are labeled as $\mathbf{x} = [\mathbf{r}_1, \sigma_1], [\mathbf{r}_2, \sigma_2], \dots, [\mathbf{r}_{N_{el}}, \sigma_{N_{el}}]$, and nuclear positions are labeled as $\mathbf{R} = \mathbf{R}_1, \mathbf{R}_2, \dots, \mathbf{R}_{N_n}$. Note that, the electronic operators $\hat{H}_{BO}$ and $\hat{V}$ are block-diagonal in spin space, and the particular form [51–53] chosen for $\hat{H}_{SO}$ does not affect any of the results presented below.

The solution of the time-dependent Schrödinger equation (TDSE) $i\hbar\partial_t\Psi(\mathbf{x}, \mathbf{R}, t) = \hat{H}(\mathbf{x}, \mathbf{R}, t)\Psi(\mathbf{x}, \mathbf{R}, t)$ is factored as [24]

$$ \Psi(\mathbf{x}, \mathbf{R}, t) = \chi(\mathbf{R}, t)\Phi_{\mathbf{R}}(\mathbf{x}, t) \quad (1) $$

with $\chi(\mathbf{R}, t)$ the nuclear wavefunction, and $\Phi_{\mathbf{R}}(\mathbf{x}, t)$ the electronic conditional factor that parametrically depends on $\mathbf{R}$. The partial normalization condition $\int d\mathbf{x} |\Phi_{\mathbf{R}}(\mathbf{x}, t)|^2 = 1 \quad \forall \mathbf{R}, t$ [54,55] guarantees that $|\chi(\mathbf{R}, t)|^2$ reproduces at all times the (exact) nuclear density obtained from $\Psi(\mathbf{x}, \mathbf{R}, t)$. Here, the integral over $\mathbf{x}$ stands for an integral over the $3N_{el}$-dimensional electronic configuration space and $N_{el}$ sums over electronic spins. The evolution of $\chi(\mathbf{R}, t)$ and $\Phi_{\mathbf{R}}(\mathbf{x}, t)$ is given by

$$ i\hbar\partial_t\chi(\mathbf{R}, t) = \left[ \sum_{\nu=1}^{N_n} \frac{[-i\hbar\nabla_\nu + \mathbf{A}_\nu(\mathbf{R}, t)]^2}{2M_\nu} + \epsilon(\mathbf{R}, t) + v(\mathbf{R}, t) \right] \chi(\mathbf{R}, t) \quad (2) $$

$$ i\hbar\partial_t\Phi_{\mathbf{R}}(\mathbf{x}, t) = \left[ \hat{H}_{BO}(\mathbf{x}, \mathbf{R}) + \hat{H}_{SO}(\mathbf{x}, \mathbf{R}) + \hat{V}(\mathbf{x}, \mathbf{R}, t) + \hat{U}_{en}[\Phi_{\mathbf{R}}, \chi] - \epsilon(\mathbf{R}, t) - v(\mathbf{R}, t) \right] \Phi_{\mathbf{R}}(\mathbf{x}, t) \quad (3) $$

The coupled nuclear and electronic evolution equations are derived by inserting the product form (1) of the molecular wavefunction into the TDSE and using the partial normalization condition, as described in detail in Refs. [25, 46, 54, 55]. In Eq. (2), the symbol $\nabla_\nu$ indicates the gradient taken with respect to the positions of the $\nu$-th nucleus and $M_\nu$ are the nuclear masses. The time-dependent vector potential is a three-dimensional vector for each value of the index $\nu$, and it is defined as [31, 32, 56, 57]

$$ \mathbf{A}_{\nu}(\mathbf{R}, t) = \langle \Phi_{\mathbf{R}}(t) | -i\hbar\nabla_{\nu} \Phi_{\mathbf{R}}(t) \rangle_{\mathbf{x}} \quad (4) $$

The time-dependent scalar potentials are

$$ \epsilon(\mathbf{R}, t) = \langle \Phi_{\mathbf{R}}(t) | \hat{H}_{BO}(\mathbf{x}, \mathbf{R}) + \hat{H}_{SO}(\mathbf{x}, \mathbf{R}) + \hat{U}_{en}[\Phi_{\mathbf{R}}, \chi] - i\hbar\partial_t | \Phi_{\mathbf{R}}(t) \rangle_{\mathbf{x}} \quad (5) $$

$$ v(\mathbf{R}, t) = \langle \Phi_{\mathbf{R}}(t) | \hat{V}(\mathbf{x}, \mathbf{R}, t) | \Phi_{\mathbf{R}}(t) \rangle_{\mathbf{x}} \quad (6) $$

and mediate the coupling between electrons and nuclei, beyond the adiabatic regime. We distinguish two scalar potential contributions to separate the effect of the external time-dependent field. The first contribution, given by Eq. (5), will be referred to as time-dependent potential energy surface (TDPES) [16, 27, 33, 34, 58, 59]; the second contribution, Eq. (6), is associated to the external field. The integration over $\mathbf{x}$ is indicated as $\langle \cdot | \cdot \rangle_{\mathbf{x}}$. Like the vector potential and the TDPES, the external contribution, Eq. (6), is, in general, an N-body operator, even though the external potential
$\hat{V}(\mathbf{x}, \mathbf{R}, t)$ (typically representing a laser pulse) is a one-body operator. The electron-nuclear coupling operator $\hat{U}_{en}[\Phi_{\mathbf{R}}, \chi]$ is [28, 29, 60, 61]

$$ \hat{U}_{en}[\Phi_{\mathbf{R}}, \chi] = \sum_{\nu=1}^{N_n} \frac{1}{M_{\nu}} \left( \frac{[-i\hbar\nabla_{\nu} - \mathbf{A}_{\nu}(\mathbf{R}, t)]^2}{2} + \left( \frac{-i\hbar\nabla_{\nu}\chi(\mathbf{R}, t)}{\chi(\mathbf{R}, t)} + \mathbf{A}_{\nu}(\mathbf{R}, t) \right) \left( -i\hbar\nabla_{\nu} - \mathbf{A}_{\nu}(\mathbf{R}, t) \right) \right) \quad (7) $$

and depends explicitly on the nuclear wavefunction $\chi(\mathbf{R}, t)$, and implicitly on the electronic factor $\Phi_{\mathbf{R}}(\mathbf{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 $\mathbf{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 Sect. 2.1, 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 Sect. 2.2, we introduce the expansion of the electronic wavefunction on a basis and we show how to solve the electronic evolution equation (3) along the classical nuclear trajectories. In Sect. 2.3, 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. In addition, some (successful) studies are currently ongoing [64] aiming to solve the exact equations employing their quantum hydrodynamic formulation. Since the original work on the exact factorization [24], the main effort has been devoted on exploiting Eqs. (2) and (3) to develop approximations, as done in CT-MQC, rather than studying their fundamental mathematical properties. Nonetheless, we can comment here on some of these properties, like the existence and uniqueness of their solutions. On the one hand, since the solution of the TDSE exists, and Ref. [24] proved the existence of the exact factorization product form, then the solutions of equations (2) and (3) exist. Uniqueness, on the other hand, has not been thoroughly investigated, and non-physical solutions might emerge, which are not related to the physical molecular wavefunction. Future theoretical devel-

opments might, indeed, focus on circumventing the numerical instabilities encountered when solving the coupled partial differential equations (2) and (3), as those observed in Ref. [63] (even though this might require to transform back the exact factorization equations into the molecular TDSE), or on curing eventual instabilities emerging in the trajectory-based CT-MQC scheme (see Sect. 2.2). In more general terms, we can say that the numerically-exact solution of the molecular TDSE is prohibitive already for small molecules like H₂. Nonetheless, techniques able to deal with the quantum mechanical evolution of nuclear wavepackets [19,20] including the effect of multiple electronic states can access up to few tens of nuclear degrees of freedom, but the potential energy surfaces and their couplings have to be pre-computed (or even expressed in analytical forms). In approaches such as the multi-configurational time-dependent Hartree (MCTDH) method [65–69], up to hundreds, or even thousands, nuclear degrees of freedom can be handled, however, the potential energy surfaces need to be pre-computed and expressed as sums of single-particle functions (and often, the majority of the degrees of freedom are harmonic oscillators that mimic the effect of a bath).

## 2.1 Solution based on characteristics of the nuclear equation

The polar representation of the (complex-valued) nuclear wavefunction, $\chi(\mathbf{R}, t) = |\chi(\mathbf{R}, t)|\exp[(i/\hbar)S(\mathbf{R}, t)]$, when inserted into the nuclear TDSE (2), yields the two coupled (real) equations [46]:

$$ -\partial_t S(\mathbf{R}, t) = \sum_{\nu=1}^{N_n} \frac{|\nabla_{\nu} S(\mathbf{R}, t) + \mathbf{A}_{\nu}(\mathbf{R}, t)|^2}{2M_{\nu}} + \epsilon(\mathbf{R}, t) + v(\mathbf{R}, t) + Q(\mathbf{R}, t) \quad (8) $$

$$ -\partial_t |\chi(\mathbf{R}, t)|^2 = \sum_{\nu=1}^{N_n} \left[ \nabla_{\nu} \cdot \frac{\nabla_{\nu} S(\mathbf{R}, t) + \mathbf{A}_{\nu}(\mathbf{R}, t)}{M_{\nu}} |\chi(\mathbf{R}, t)|^2 \right] \quad (9) $$

The evolution equation for the phase $S(\mathbf{R}, t)$, Eq. (8), has a Hamilton-Jacobi form, where the Hamiltonian on the right-hand side contains a kinetic term, with $\nabla_{\nu}S(\mathbf{R}, t) + \mathbf{A}_{\nu}(\mathbf{R}, t)$ being the momentum field, two “classical” time-dependent potential terms, $\epsilon(\mathbf{R}, t)$ and $v(\mathbf{R}, t)$, and a quantum potential term [62,70],

$$ Q(\mathbf{R}, t) = \sum_{\nu=1}^{N_n} \frac{-\hbar^2 \nabla_{\nu}^2 |\chi(\mathbf{R}, t)|}{2M_{\nu} |\chi(\mathbf{R}, t)|} \quad (10) $$

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 $|\chi(\mathbf{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” $\mathbf{R}(s)$, $t(s)$, $S(s)$, $\tilde{\mathbf{P}}(s) = \nabla S(s)$, $\partial_t S(s) = S_t(s)$ that always satisfy the original PDE, as functions of the parameter $s$. We indicate here the $3N_n$-dimensional gradient of $S$ as $\nabla S(s) = \{\nabla_\nu S(s)\}_{\nu=1,...,N_n}$ to simplify the notation. We introduce the function $F(\mathbf{R}, t, S, \tilde{\mathbf{P}}, S_t) = S_t + H_n(\tilde{\mathbf{P}}, \mathbf{R}, t) = 0$, whose total differential is $dF = F_{\mathbf{R}} \cdot d\mathbf{R} + F_t dt + F_S dS + F_{\tilde{\mathbf{P}}} \cdot d\tilde{\mathbf{P}} + F_{S_t} dS_t$. The quantities $F_{\mathbf{R}}$, $F_t$, $F_S$, $F_{\tilde{\mathbf{P}}}$, and $F_{S_t}$ 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 (non-trivial) characteristic equations representing the evolution of the variables $(\mathbf{R}, t, S, \tilde{\mathbf{P}}, S_t)$ are

$$ \dot{\mathbf{R}}_{\nu} = \frac{\tilde{\mathbf{P}}_{\nu} + \mathbf{A}_{\nu}}{M_{\nu}} \qquad (11) $$

$$ \dot{\tilde{\mathbf{P}}}_{\nu} = -\nabla_{\nu} H_n(\tilde{\mathbf{P}}, \mathbf{R}, t) \qquad (12) $$

with $\nu = 1, ..., N_n$, where the Hamiltonian $H_n(\tilde{\mathbf{P}}, \mathbf{R}, t)$ is the right-hand side of Eq. (8). The symbols $\dot{\mathbf{R}}_{\nu}$ and $\dot{\tilde{\mathbf{P}}}_{\nu}$ are intended as total derivatives with respect to the parameter $s$ for each value of the nuclear index $\nu$. The characteristic equation $\dot{t} = 1$ shows that the parameter $s$ can be identified as the physical time $t$: the map $t(s) = s + cnst$ is, thus, invertible for a given choice of cnst (we define the map such that $t(s = 0) = 0$). The characteristic equations for $S$ and $S_t$, not reported here, can be easily derived from $\dot{\tilde{\mathbf{P}}}_{\nu}$, since $\tilde{\mathbf{P}}_{\nu} = \nabla_{\nu} S$. Equations (11) and (12) guarantee that the vector $(d\mathbf{R}, dt, dS, d\tilde{\mathbf{P}}, dS_t)$ is orthogonal to the gradient $(F_{\mathbf{R}}, F_t, F_S, F_{\tilde{\mathbf{P}}}, F_{S_t})$ 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(\mathbf{R}, t)$ (and $\nabla_{\nu} S(\mathbf{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(\mathbf{R}, t)$, such that the equations decouple, in the sense that the equation for $S(\mathbf{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 tunneling cannot be captured in this case.

(ii) The pseudo-quantum procedure [46, 60] accounts for the quantum potential $Q(\mathbf{R}, t)$ in the evolution for $S(\mathbf{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, 64, 71] is devised to solve the coupled equations for the phase $S(\mathbf{R}, t)$ and for the density $|\chi(\mathbf{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 $\mathbf{R}$ become functions of $\mathbf{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(\mathbf{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 $\epsilon(\mathbf{R}(t), t) + v(\mathbf{R}(t), t) + \sum_{\nu} \dot{\mathbf{R}}_{\nu}(t) \cdot \mathbf{A}_{\nu}(\mathbf{R}(t), t) = 0$ to fix the gauge freedom. The characteristic equations thus reduce to

$$ \dot{\mathbf{R}}_{\nu}(t) = \frac{\mathbf{P}_{\nu}(t)}{M_{\nu}} \qquad (13) $$

$$ \dot{\mathbf{P}}_{\nu}(t) = \dot{\mathbf{A}}_{\nu}(\mathbf{R}(t), t) \qquad (14) $$

where we introduced the new symbol $\mathbf{P}_{\nu} = \tilde{\mathbf{P}}_{\nu} + \mathbf{A}_{\nu}$ 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 Sect. 2.3, we provide a more explicit expression of the classical force, using the dependence of $\mathbf{A}_{\nu}(\mathbf{R}(t), t)$ on the electronic wavefunction. To this end, we describe in Sect. 2.2 how to derive the evolution equation for $\Phi_{\mathbf{R}(t)}(\mathbf{x}, t)$ along the nuclear characteristics.

## 2.2 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 $\dot{\Phi}_{\mathbf{R}(t)}(\mathbf{x}, t)$ has to be used, rather than $\partial_t \Phi_{\mathbf{R}(t)}(\mathbf{x}, t)$. This is done by applying the chain rule $\partial_t \Phi_{\mathbf{R}(t)}(\mathbf{x}, t) = \dot{\Phi}_{\mathbf{R}(t)}(\mathbf{x}, t) - \sum_{\nu} \dot{\mathbf{R}}_{\nu}(t) \cdot \nabla_{\nu} \Phi_{\mathbf{R}(t)}(\mathbf{x}, t)$ anywhere the partial time derivative is found.

The derivation of CT-MQC relies on the expansion of the electronic wavefunction $\Phi_{\mathbf{R}(t)}(\mathbf{x}, t)$ on a basis, namely
$$ \Phi_{\mathbf{R}(t)}(\mathbf{x}, t) = \sum_{m} C_m(\mathbf{R}(t), t) \varphi_{\mathbf{R}(t)}^{(m)}(\mathbf{x}) \quad (15) $$

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 Hamiltonian $\hat{H}_{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,72]. Spin-diabatic 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 Hamiltonian $\hat{H}_{BO}$. Transitions among electronic spin-diabatic 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 $\hat{H}_{BO} + \hat{H}_{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 $\hat{H}_{BO}$ (or $\hat{H}_{BO} + \hat{H}_{SO}$ if SOC is present), which has been called field-diabatic representation in Ref. [73]. 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,74–80] appears to be a viable option especially for interpreting absorption/emission processes [35,81–87] in terms of exchanges of photons between the system and the field. Floquet representation exploits the periodicity induced by the laser in the Hamiltonian [88], which is encoded in the mixed character of field-electronic states, defined as the eigenstates

of $\hat{H}_{BO} - i\hbar\partial_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 $\hat{H}_{BO} + \hat{V} - i\hbar\partial_t$, even though its applications in the field of quantum molecular dynamics simulations have not been thoroughly investigated so far.

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 $\dot{C}_m(\mathbf{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 $\hat{U}_{en}$ we identify two terms: the first term has been shown [61,89,90] 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 $\hat{U}_{en}$ reduces to

$$ \hat{U}_{en} \approx \sum_{\nu=1}^{N_n} \left( \dot{\mathbf{R}}_{\nu}(t) + i \frac{\mathbf{P}_{\nu}(\mathbf{R}(t), t)}{M_{\nu}} \right) (-i\hbar\nabla_{\nu} - \mathbf{A}_{\nu}(\mathbf{R}(t), t)) \quad (16) $$

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], $\mathbf{P}_{\nu}(\mathbf{R}(t), t) = -\hbar|\chi(\mathbf{R}(t), t)|^{-1}\nabla_{\nu}|\chi(\mathbf{R}(t), t)|$, that has been defined in previous work on CT-MQC. The quantum momentum [91,92] 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. If we use an index $\alpha$ to label the dependence on the tra-
trajectory $\mathbf{R}^{\alpha}(t)$, the quantum momentum at the position of the trajectory is

$$ \mathbf{P}_{\nu}^{\alpha}(t) = \frac{-\hbar \nabla_{\nu} |\chi(\mathbf{R}^{\alpha}(t), t)|^2}{2 |\chi(\mathbf{R}^{\alpha}(t), t)|^2} \quad (17) $$

where $|\chi(\mathbf{R}^{\alpha}(t), t)|^2 = [|\chi(\mathbf{R}, t)|^2]_{\mathbf{R}=\mathbf{R}^{\alpha}(t)}$. The nuclear density at time $t$ is constructed as the sum of (normalized real) Gaussians centered at the positions of the traj trajectories, thus

$$ |\chi(\mathbf{R}^{\alpha}(t), t)|^2 = \frac{1}{N_{\text{traj}}} \sum_{\beta=1}^{N_{\text{traj}}} G_{\sigma}(\mathbf{R}^{\alpha}(t) - \mathbf{R}^{\beta}(t)) \quad (18) $$

$G_{\sigma}$ is a multi-dimensional Gaussian, which can be expressed as the product of $3N_n$ uni-dimensional Gaussians. Each of the $3N_n$ nuclear degrees of freedom can be associated to a different width, which is kept constant along the dynamics, and each of the $N_{\text{traj}}$ replicas has the same width. The widths of the Gaussians are indicated in Eq. (18) with the symbol $\sigma$. Equation (18) clearly shows that, to determine the nuclear density, and thus, the quantum momentum at time $t$, information about the positions of all trajectories $\mathbf{R}^{\beta}(t)$ at that time is required. It is in this sense that CT-MQC trajectories are coupled.

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.

To isolate the evolution of one of the coefficients, say $C_m(\mathbf{R}(t), t)$, one has to project the electronic equation onto $\varphi_{\mathbf{R}(t)}^{(m)}(\mathbf{x})$, by integrating over $\mathbf{x}$. The NAC vectors and the matrix elements of the SOC appear in the final equation, and are defined as $\mathbf{d}_{\nu,ml}^{\alpha} = \langle \varphi_{\mathbf{R}^{\alpha}(t)}^{(m)} | \nabla_{\nu} \varphi_{\mathbf{R}^{\alpha}(t)}^{(l)} \rangle_{\mathbf{x}}$ and $[H_{SO}^{\alpha}]_{ml} = \langle \varphi_{\mathbf{R}^{\alpha}(t)}^{(m)} | \hat{H}_{SO}(\mathbf{R}^{\alpha}(t)) | \varphi_{\mathbf{R}^{\alpha}(t)}^{(l)} \rangle_{\mathbf{x}}$, respectively.

The action of the Born-Oppenheimer Hamiltonian on its eigenstates, i.e., $\hat{H}_{BO}(\mathbf{x}, \mathbf{R}^{\alpha}(t))\varphi_{\mathbf{R}^{\alpha}(t)}^{(m)}(\mathbf{x}) = E_m(\mathbf{R}^{\alpha}(t)) \varphi_{\mathbf{R}^{\alpha}(t)}^{(m)}(\mathbf{x})$, yields the eigenvalues $E_m^{\alpha} = E_m(\mathbf{R}^{\alpha}(t))$, clearly depending on the position of the trajectory.

As consequence of the action of the gradient operator $\nabla_{\nu}$ in Eq. (16) on the electronic wavefunction, the equation contains a term like $\nabla_{\nu}C_m^{\alpha}(t)$, which is approximated by neglecting the spatial dependence of the modulus of $C_m^{\alpha}(t)$ in favor of the spatial dependence of its phase. Therefore, $\nabla_{\nu}C_m^{\alpha}(t) \simeq (i/\hbar)\mathbf{f}_{\nu,m}^{\alpha}C_m^{\alpha}(t)$ and we use a simple expression for the gradient of the phase, namely $\mathbf{f}_{\nu,m}^{\alpha} = \int_{0}^{t} d\tau [-\nabla_{\nu}E_m^{\alpha}]$; this quantity is a spin-diabatic force accumulated over time along the trajectory $\alpha$.

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 $\mathbf{R}^{\alpha}(t)$ as

$$ \dot{C}_m^\alpha(t) = \left[ -\frac{i}{\hbar} E_m^\alpha + \sum_{\nu=1}^{N_n} \frac{\mathbf{P}_\nu^\alpha(t)}{\hbar M_\nu} \cdot (\mathbf{f}_\nu^\alpha - \mathbf{A}_\nu^\alpha(t)) \right] C_m^\alpha(t) \quad (19a) $$

$$ -\sum_l \left( \frac{i}{\hbar} [H_{SO}^{\alpha}]_{ml} + \sum_{\nu} \dot{\mathbf{R}}_{\nu}^{\alpha}(t) \cdot \mathbf{d}_{\nu,ml}^{\alpha} \right) C_{m}^{\alpha}(t) \quad (19b) $$

Some quantities in this equation depend explicitly on time, and such dependence has been indicated together with the dependence on $\mathbf{R}^{\alpha}(t)$; some other quantities only depend on time via their dependence on the trajectory, which is indicated with the superscript $\alpha$.

In Eqs. (19) electronic information is encoded in the energies, i.e., $E_m^{\alpha}$, in the NACs, i.e., $\mathbf{d}_{\nu,ml}^{\alpha}$, and in SOC, i.e., $[H_{SO}^{\alpha}]_{ml}$. Those quantities are usually computed on-the-fly at each time by solving the electronic (many-body) problem for a given nuclear configuration, i.e., for the configurations visited by the trajectories during the dynamics. In general, the (finite) number of electronic states considered in a simulation is chosen a priori, so as to truncate the expansion in Eq. (15). Selecting the electronic states requires to assess electronic structure properties of the system by scanning the potential energy landscape along different nuclear paths, to identify minima, transitions states, and degeneracies. In particular, potential energies can become degenerate, for instance at conical intersections, causing (i) discontinuous changes in the electronic wavefunction characters and (ii) nonadiabatic singularities. However, in trajectory-based schemes, issue (i) is circumvented since only local electronic structure quantities are required, and issue (ii) seldom appears because the trajectories do not often encounter the real singularity.

The first diagonal term, Eq. (19a), 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^{\alpha}(t)$, the latter affects the modulus of $C_m^{\alpha}(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 Eqs. (19), i.e., $\sum_m |C_m^{\alpha}(t)|^2 = 1 \forall t, \alpha$. The off-diagonal contributions in Eq. (19b) drive electronic transitions from state $m$ to state $l$, mediated by the SOC, if states with different spin multiplicity are concerned, and by the NAC vectors, if states with the same spin multiplicity are involved.

The time-dependent vector potential appears in Eq. (19a). 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 $\mathbf{R}^{\alpha}(t)$ as $\mathbf{A}_{\nu}^{\alpha}(t) = -i\hbar \sum_m \dot{C}_m^{\alpha}(t) \nabla_{\nu} C_m^{\alpha}(t)$ -
$i\hbar \sum_{m,l} \bar{C}_m^\alpha(t) C_l^\alpha(t) \mathbf{d}_{\nu,ml}^\alpha$. Recalling the approximation used above to express the gradient of the electronic coefficients, and observing that the first term of this expression 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) time-dependent vector potential is

$$ \mathbf{A}_{\nu}^{\alpha}(t) = \sum_{m} |C_{m}^{\alpha}(t)|^{2} \mathbf{f}_{\nu,m}^{\alpha} \quad (20) $$

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 Eqs. (19). The final expression in given in Sect. 2.3.

## 2.3 CT-MQC algorithm

Sections 2.1 and 2.2 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 quantum-classical, representing a way to reformulate the quantum mechanical equations of the exact factorization, i.e., Eqs. (2) and (3). In particular, the electronic equation (19) is solved along a trajectory, thus it lends itself for an on-the-fly approach, where electronic structure properties are only computed at the instantaneous nuclear positions. In turn, the nuclei evolve according to a classical-like 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 $\mathbf{R}^\alpha(t)$

$$ \mathbf{F}_{\nu}^{\alpha}(t) = \mathbf{F}_{\nu,\text{Eh}}^{\alpha}(t) + \mathbf{F}_{\nu,\text{qm}}^{\alpha}(t) + \mathbf{F}_{\nu,\text{SOC}}^{\alpha}(t) \quad (21) $$

The Ehrenfest-like term (Eh) is a standard mean-field force

$$ \mathbf{F}_{\nu,\text{Eh}}^{\alpha}(t) = \sum_{m} |C_{m}^{\alpha}(t)|^{2} (-\nabla_{\nu} E_{m}^{\alpha}) + \sum_{m,l} \bar{C}_{m}^{\alpha}(t) C_{l}^{\alpha}(t) (E_{m}^{\alpha} - E_{l}^{\alpha}) \mathbf{d}_{\nu,ml} \quad (22) $$

The quantum momentum term (qm) depends on $\mathbf{P}_{\nu}^{\alpha}(t)$, namely

$$ \mathbf{F}_{\nu,\text{qm}}^{\alpha}(t) = \frac{2}{\hbar} \sum_{m} |C_{m}^{\alpha}(t)|^{2} \left[ \sum_{\mu=1}^{N_{n}} \mathbf{P}_{\mu}^{\alpha}(t) \cdot \mathbf{f}_{\mu,m}^{\alpha} \right] (\mathbf{f}_{\nu,m}^{\alpha} - \mathbf{A}_{\nu}^{\alpha}(t)) \quad (23) $$

and couples the trajectories through the presence of the quantum momentum, as described in Sect. 2.2. The

SOC term contains the matrix elements of the spin-orbit Hamiltonian

$$ \mathbf{F}_{\nu,\text{SOC}}^{\alpha}(t) = \frac{1}{\hbar} \sum_{m,l} \operatorname{Im} (\bar{C}_{m}^{\alpha}(t) C_{l}^{\alpha}(t) [H_{\text{SO}}^{\alpha}]_{ml}) (\mathbf{f}_{\nu,l}^{\alpha}(t) - \mathbf{f}_{\nu,m}^{\alpha}(t)) \quad (24) $$

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 [93]. 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. (24), 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 trajectory-based 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 [72,94,95] and surface hopping [52,96–98] rely on the fact that NACs are spatially localized. The choice between spin-diabatic or spin-adiabatic basis when dealing with SOC depends on the quantum-chemistry package chosen for the electronic structure calculations; quantum-chemistry codes typically provide electronic structure 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. (21), we can rewrite the electronic evolution equation (19) as

$$ \dot{\mathbf{C}}_{m}^{\alpha}(t) = \dot{\mathbf{C}}_{m}^{\alpha}(t)|_{\text{Eh}} + \dot{\mathbf{C}}_{m}^{\alpha}(t)|_{\text{qm}} + \dot{\mathbf{C}}_{m}^{\alpha}(t)|_{\text{SOC}} \quad (25) $$

where we identify the following terms:

$$ \dot{C}_{m}^{\alpha}(t)|_{\text{Eh}} = -\frac{i}{\hbar} E_{m}^{\alpha} C_{m}^{\alpha}(t) - \sum_{l} \sum_{\nu=1}^{N_{n}} \dot{\mathbf{R}}_{\nu}^{\alpha}(t) \cdot \mathbf{d}_{\nu,ml} C_{l}^{\alpha}(t) \quad (26a) $$

$$ \dot{C}_{m}^{\alpha}(t)|_{\text{qm}} = \sum_{\nu=1}^{N_{n}} \frac{\mathbf{P}_{\nu}^{\alpha}(t)}{\hbar M_{\nu}} \cdot (\mathbf{f}_{\nu,m}^{\alpha} - \mathbf{A}_{\nu}^{\alpha}(t)) C_{m}^{\alpha}(t) \quad (26b) $$
$$ \dot{C}_m^\alpha(t)|_{\text{SOC}} = -\frac{i}{\hbar} \sum_l [H_{SO}^\alpha]_{ml} C_l^\alpha(t) \quad (26c) $$

In the absence of SOC, the term $[H_{SO}^\alpha]_{ml}$ in Eq. (26c) is not present, and the electronic equation reduces to the one in the original derivation of CT-MQC in the adiabatic representation reported in Ref. [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 spin-orbit effects, one can transform the coupling mediated by the SOC into NAC in the spin-adiabatic representation. In this case, $[H_{SO}^\alpha]_{ml}$ of Eq. (26c) 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. (21), namely, $F_{\nu,\text{ext}}^\alpha(t) = \hbar^{-1} \sum_{m,l} \text{Im} [C_m^\alpha(t) C_l^\alpha(t)] V_{ml}^\alpha(f_{\nu,l}^\alpha - f_{\nu,m}^\alpha)$, and in the expression of the evolution of the electronic coefficients of Eq. (25), namely $\dot{C}_m^\alpha(t)|_{\text{ext}} = -i\hbar^{-1} \sum_l V_{ml}^\alpha C_l^\alpha(t)$. Furthermore, the energies $E_m^\alpha$ of Eq. (26b) 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}^\alpha$ 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}^\alpha$ 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.

## 3 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 [99], a strong transition dipole moment couples the ground state of the molecule $X(^1\Sigma_0^+)$ to its excited electronic state $B(^3\Pi_{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\Sigma_{0+}) \rightarrow B(^3\Pi_{0+})$. The electronic states considered here are denoted using typical spin-diabatic labels [96, 100] even though we will work in the spin-adiabatic representation as well, based on the model potentials of Ref. [101]. In particular, the spin-diabatic characters

$B(^3\Pi_{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 photo-excitation, the products of such an ultrafast photo-reaction are $I + \text{Br}(^2\text{P}_{3/2})$, if the molecule dissociates via the first-excited state, and $I + \text{Br}^*(^2\text{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 second-excited (spin-adiabatic) state acquires a $B(^3\Pi_{0+})$ (spin-diabatic) 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

$$ Q = \frac{[I + \text{Br}^*]}{[I + \text{Br}] + [I + \text{Br}^*]} \quad (27) $$

In our simulations, we use the model Hamiltonian of Ref. [101], expressed in a spin-diabatic basis as

$$ \hat{H}_{BO}(R) = \begin{pmatrix} H_0(R) & 0 & 0 \\ 0 & H_1(R) & V_{12} \\ 0 & V_{12} & H_2(R) \end{pmatrix} \quad (28) $$

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. [101], we have

$$ H_0(R) = A_0 \left[ \left( 1 - e^{-\alpha_0 (R - R_0)} \right)^2 - 1 \right] \quad (29) $$

$$ H_1(R) = A_1 \left[ \left( 1 - e^{-\alpha_1 (R - R_1)} \right)^2 - 1 \right] + \Delta \quad (30) $$

$$ H_2(R) = A_2 e^{-\alpha_2 R} + B_2 e^{-\beta_2 R} \quad (31) $$

where $A_0 = 0.067$ Ha, $\alpha_0 = 0.996$ bohr$^{-1}$, $R_0 = 4.666$ bohr, $A_1 = 0.01019$ Ha, $\alpha_1 = 1.271$ bohr$^{-1}$, $R_1 = 5.3479$ bohr, $\Delta = 0.01679$ Ha, $A_2 = 2.82$ Ha, $\alpha_2 = 0.9186$ bohr$^{-1}$, $B_2 = 3.0 \times 10^7$ Ha, $\beta_2 = 4.3$ bohr$^{-1}$, and $V_{12} = 0.0006834$ Ha.

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. (28), 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. (28). However, as the transition dipole moment between $X(^1\Sigma_{0+})$ and $B(^3\Pi_{0+})$ is strong, an optical transition can be induced via an external field.
<figure><img src="image_2.png" /><figcaption>Fig. 1 Potential energy curves (PECs) corresponding to the spin-(a)diabatic states considered to describe the photo-dissociation process of IBr. Continuous colored lines are the spin-adiabatic PECs: ground state in blue, first-excited state in red, second-excited state in green. They are labeled using the notation of Ref. [99] in the Franck–Condon region. Dashed yellow lines are spin-diabatic PECs. The initial nuclear density centered at the equilibrium position of the ground state PEC is shown as a thin black line. The arrows indicate the dissociation channels along the two electronic states</figcaption></figure>

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 wavefunction $\chi(R, t = 0)$ is a (real) Gaussian function

$$ \chi(R, t = 0) = \sqrt[4]{\frac{1}{\pi\sigma^2}} e^{-\frac{(R-R_0)^2}{2\sigma^2}} \quad (32) $$

with variance $\sigma = 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. Equation (32) is the vibrational ground state of lowest electronic PEC around the equilibrium position. Therefore, Eq. (32) is a suitable initial nuclear wavefunction before the instantaneous excitation promotes it to the first-excited state. In addition, electronic population of that state can safely be chosen equal to unity since the Franck–Condon region is far away from the avoided crossing located at $R_c$.

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. (32). The Wigner transform of a Gaussian function in $R$ (see Eq. (32)) yields a position-

<figure><img src="image_3.png" /><figcaption>Fig. 2 Populations of the spin-(a)diabatic states as functions of time. Continuous lines are used to indicate spin-adiabatic populations: ground state in blue, first-excited state in red, second-excited state in green. Dashed lines are used to indicate spin-diabatic populations: $X(^1\Sigma_{0+})$ in gray, $B(^3\Pi_{0+})$ in dark-yellow, $Y(0^+)$ in cyan. Exact spin-adiabatic results are shown as thin black lines</figcaption></figure>

momentum distribution in the form of an uncorrelated product of a Gaussian in $R$ and a Gaussian in $P$; both distributions are strictly positive and standard sampling techniques [102] can be applied to select the initial conditions. For this reason, in general situations the Wigner sampling is usually employed by invoking the harmonic approximation of the potential energy in the region where the sampling is performed. In the presence of strong anharmonicities, a classical Boltzmann sampling is, instead, preferred [19], to avoid dealing with negative portions of the Wigner distribution. For all calculations, we run $N_{traj} = 1000$ coupled trajectories. For all trajectories, 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.

Populations of the electronic states, denoted here $\rho_m(t)$, are computed based on the CT-MQC procedure as averages over the trajectories, namely

$$ \rho_m(t) = \frac{1}{N_{traj}} \sum_{\alpha=1}^{N_{traj}} |C_m^{\alpha}(t)|^2 \quad (33) $$

with $m = \text{ground state}$, 1st excited state, 2nd excited state in the spin-adiabatic representation, or $m = X(^1\Sigma_{0+})$, $B(^3\Pi_{0+})$, $Y(0^+)$ in the spin-diabatic representation. Note that, when performing an average operation to compute observable quantities, as it is done in Eq. (33), the trajectories have all the same weight, because the initial Wigner sampling is an importance
<table><caption>**Table 1** Branching ratio of the products of the photo-dissociation reaction of IBr. The reference values (exact) are determined based on quantum mechanical simulations using Eq. (27). The branching ratio is also calculated based on CT-MQC simulations in the spin-adiabatic and spin-diabatic basis</caption><thead><tr><th></th><th>Exact</th><th>Spin-adiabatic</th><th>Spin-diabatic</th></tr></thead><tbody><tr><td>Q</td><td>75%</td><td>76(1)%</td><td>76(1)%</td></tr></tbody></table>

sampling which initializes the trajectories where the initial density is large.

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\Pi_{0^+})$ in the dissociating limit (dashed dark-yellow line in Fig. 2). Therefore, the branching ratio of the dissociation products can be determined using either set of results; numerical values are reported in Table 1.

In Table 1 the branching ratios calculated based on CT-MQC results are compared to the exact values calculated from the results shown in Fig. 2. For CT-MQC, the errors on the branching ratios (in parentheses) are determined from the standard deviations on the predictions of the final populations, and, within the error, CT-MQC values are in perfect agreement with the reference.

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 spin-adiabatic (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 cal-

<figure><img src="image_4.png" /><figcaption>**Fig. 3** PECs corresponding to the spin-adiabatic states are shown as thin colored lines. The nuclear densities corresponding to the first-excited state (thick red lines) and to the second-excited state (thick green lines) are shown at three time steps: $t = 12, 60, 180$ fs, as indicated in the figure. At the same times, the distributions of CT-MQC trajectories evolved in the spin-adiabatic (black dots) and in the spin-diabatic (colored dots) basis are shown. In particular, we use the value of the TDPES from Eq. (5) at the position of each trajectory</figcaption></figure>

culation of the quantum momentum, and to implement the effect of an external time-dependent field beyond the Floquet formalism.

## 4 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 time-dependent potential energy surface and of the time-dependent 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 non-classical 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 [103], the response of chiral molecules to infrared left and right circularly polarized light, known as vibrational circular dichroism [89],
and to estimate nonadiabatic corrections to infrared spectra of hydrogen-based molecules [90].

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 photoinduced ultrafast dynamics [43], strong field ionization [36,104], 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,105], charge-resonance enhanced ionization [36] in $H_2^+$, or to challenge single-electron pictures [106] to describe molecules in strong lasers [38,107]. 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 [108–110], 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–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 artifact? 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 [111,112]. Finally, a quantum electronic embedding method has been derived from the exact factorization to calculate static properties of a many-electron system [113,114].

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 electron–photon 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 [115] (reminiscent of the inverse factorization), as well as on the dynamic effect of the cavity on the electron–nuclear time-dependent potential energy surface [59,116], in analogy with polaritonic and cavity-Born–Oppenheimer surfaces. Application of the exact factorization in this domain is still preliminary, but the success of the original theory in many diverse fields strongly suggests 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 [117]. The coupled nature of CT-MQC trajectories, however, remains the bottleneck even when parallelization strategies are employed, as was done in the CPMD implementation. To circumvent this issue, recently [42,118–120] 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 “experimental” 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-trajectory-based scheme.

## 5 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 var-
ious 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.

**Data Availability Statement** This manuscript has no associated data or the data will not be deposited. [Authors' comment: Data is available upon request from the corresponding author.]

**Acknowledges** EKUG acknowledges financial support from the European Research Council Advanced Grant FACT (ERC-2017-AdG-788890).

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

1. A.H. Zewail, Science **242**, 1645 (1988)

2. B.F.E. Curchod, T.J. Martínez, Chem. Rev. **118**, 3305 (2018)

3. R. Crespo-Otero, M. Barbatti, Chem. Rev. **118**, 7026 (2018)

4. F. Agostini, B.F.E. Curchod, WIREs Comput. Mol. Sci. **9**, e1417 (2019)

5. J. Fregoni, G. Granucci, M. Persico, S. Corni, Chem. **6**, 250 (2019)

6. J.J. Bajo, G. Granucci, M. Persico, J. Chem. Phys. **140**, 044113 (2014)

7. G.A. Worth, L.S. Cederbaum, Annu. Rev. Phys. Chem. **55**, 127 (2004)

8. J.E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, N. Bellonzi, Annu. Rev. Phys. Chem. **67**, 387 (2016)

9. L.M. Ibele, B.F.E. Curchod, Phys. Chem. Chem. Phys. **22**, 15062 (2020)

10. G.S. Engel, T.R. Calhoun, E.L. Read, T.K. Ahn, Y.C.C.T. Mančal, R.E. Blankenship, G.R. Fleming, Nature **446**, 782 (2007)

11. G.D. Scholes, G.R. Fleming, L.X. Chen, A. Aspuru-Guzik, A. Buchleitner, D.F. Coker, G.S. Engel, R. van Grondelle, A. Ishizaki, D.M. Jonas et al., Nature **543**, 647 (2017)

12. C.A. Rozzi, S.M. Falke, N. Spallanzani, A. Rubio, E. Molinari, D. Brida, M. Maiuri, G. Cerullo, H. Schramm, J. Christoffers et al., Nat. Commun. **4**, 1602 (2013)

13. F. Agostini, Euro. Phys. J. B **91**, 143 (2018)

14. V. Cantatore, G. Granucci, M. Persico, Comput. Theo. Chem. **1040**, 126 (2014)

15. S. Giannini, A. Carof, M. Ellis, H. Yang, O.G. Zioshos, S. Ghosh, J. Blumberger, Nat. Commun. **10**, 3843 (2019)

16. B.F.E. Curchod, F. Agostini, E.K.U. Gross, J. Chem. Phys. **145**, 034103 (2016)

17. G. Granucci, M. Persico, Chem. Phys. Lett. **246**, 228 (1995)

18. S. Gozem, H.L. Luk, I. Schapiro, M. Olivucci, Chem. Rev. **117**, 13502 (2017)

19. C. Pieroni, E. Marsili, D. Lauvergnat, F. Agostini, J. Chem. Phys. **154**, 034104 (2021)

20. E. Marsili, M. Olivucci, D. Lauvergnat, F. Agostini, J. Chem Theory Comput. **16**, 6032 (2020)

21. D. Fausti, R.I. Tobey, N. Dean, S. Kaiser, A. Dienst, M.C. Hoffmann, S. Pyon, H.T.T. Takayama, A. Cavalleri, Science **189**, 331 (2011)

22. D. Nicoletti, A. Cavalleri, Adv. Opt. Photon. **8**, 401 (2016)

23. F. Agostini, E.K.U. Gross, in *Quantum chemistry and dynamics of excited states: Methods and applications*, ed. by L. González, R. Lindh (Wiley, 2021), pp. 531–562

24. A. Abedi, N.T. Maitra, E.K.U. Gross, Phys. Rev. Lett. **105**, 123002 (2010)

25. A. Abedi, N.T. Maitra, E.K.U. Gross, J. Chem. Phys. **137**, 22A530 (2012)

26. F. Agostini, B.F.E. Curchod, R. Vuilleumier, I. Taver-nelli, E.K.U. Gross, in *Handbook of Materials Modeling*, ed. by W. Andreoni, S. Yip (Springer, 2018), (pp. 1–47)

27. F. Agostini, A. Abedi, Y. Suzuki, E.K.U. Gross, Mol. Phys. **111**, 3625 (2013)

28. A. Abedi, F. Agostini, E.K.U. Gross, Europhys. Lett. **106**, 33001 (2014)

29. F. Agostini, A. Abedi, E.K.U. Gross, J. Chem. Phys. **141**, 214101 (2014)

30. S.K. Min, F. Agostini, E.K.U. Gross, Phys. Rev. Lett. **115**, 073001 (2015)

31. B.F.E. Curchod, F. Agostini, J. Phys. Chem. Lett. **8**, 831 (2017)

32. F. Agostini, B.F.E. Curchod, Euro. Phys. J. B **91**, 141 (2018)

33. A. Abedi, F. Agostini, Y. Suzuki, E.K.U. Gross, Phys. Rev. Lett **110**, 263001 (2013)

34. F. Agostini, A. Abedi, Y. Suzuki, S.K. Min, N.T. Maitra, E.K.U. Gross, J. Chem. Phys. **142**, 084303 (2015)

35. T. Fiedlschuster, J. Handt, E.K.U. Gross, R. Schmidt, Phys. Rev. A **95**, 063424 (2017)

36. E. Khosravi, A. Abedi, N.T. Maitra, Phys. Rev. Lett. **115**, 263002 (2015)

37. Y. Suzuki, A. Abedi, N.T. Maitra, K. Yamashita, E.K.U. Gross, Phys. Rev. A **89**, 040501(R) (2014)

38. A. Schild, E.K.U. Gross, Phys. Rev. Lett. **118**, 163202 (2017)

39. F. Agostini, E.K.U. Gross, B.F.E. Curchod, Comput. Theo. Chem. **1151**, 99 (2019)

40. G. Gossel, F. Agostini, N.T. Maitra, J. Chem. Theory Comput. **14**, 4513 (2018)
41. B.F.E. Curchod, F. Agostini, I. Tavernelli, Euro. Phys. J. B **91**, 168 (2018)

42. J.K. Ha, I.S. Lee, S.K. Min, J. Phys. Chem. Lett. **9**, 1097 (2018)

43. S.K. Min, F. Agostini, I. Tavernelli, E.K.U. Gross, J. Phys. Chem. Lett. **8**, 3048 (2017)

44. F. Agostini, S.K. Min, A. Abedi, E.K.U. Gross, J. Chem. Theory Comput. **12**, 2127 (2016)

45. F. Talotta, F. Agostini, G. Ciccotti, J. Phys. Chem. A **124**, 6764 (2020)

46. F. Agostini, I. Tavernelli, G. Ciccotti, Euro. Phys. J. B **91**, 139 (2018)

47. Y. Suzuki, K. Watanabe, Phys. Rev. A **94**, 032517 (2016)

48. F. Talotta, S. Morisset, N. Rougeau, D. Lauvergnat, F. Agostini, Phys. Rev. Lett. **124**, 033001 (2020)

49. F. Talotta, S. Morisset, N. Rougeau, D. Lauvergnat, F. Agostini, J. Chem Theory Comput. **16**, 4833 (2020)

50. M. Schirò, F.G. Eich, F. Agostini, J. Chem. Phys. **154**, 114101 (2021)

51. E. van Lenthe, E.J. Baerends, J.G. Snijders, J. Chem. Phys. **99**, 4597 (1993)

52. F.F. de Carvalho, I. Tavernelli, J. Chem. Phys. **143**, 224105 (2015)

53. C.M. Marian, WIREs Comput. Mol. Sci. **2**, 187 (2012)

54. J.L. Alonso, J. Clemente-Gallardo, P. Echeniche-Robba, J.A. Jover-Galtier, J. Chem. Phys. **139**, 087101 (2013)

55. A. Abedi, N.T. Maitra, E.K.U. Gross, J. Chem. Phys. **139**, 087102 (2013)

56. R. Requist, F. Tandetzky, E.K.U. Gross, Phys. Rev. A **93**, 042108 (2016)

57. R. Requist, C.R. Proetto, E.K.U. Gross, Phys. Rev. A **96**, 062503 (2017)

58. S.K. Min, A. Abedi, K.S. Kim, E.K.U. Gross, Phys. Rev. Lett. **113**, 263004 (2014)

59. L. Lacombe, N.M. Hoffmann, N.T. Maitra, Phys. Rev. Lett. **123**, 083201 (2019)

60. F. Agostini, S.K. Min, E.K.U. Gross, Ann. Phys. **527**, 546 (2015)

61. F.G. Eich, F. Agostini, J. Chem. Phys. **145**, 054110 (2016)

62. R.E. Wyatt, in *Quantum Dynamics with Trajectories: Introduction to Quantum Hydrodynamics* (Interdisciplinary Applied Mathematics, Springer, 2005)

63. G.H. Gossel, L. Lacombe, N.T. Maitra, J. Chem. Phys. **150**, 154112 (2019)

64. D.D. Holm, J.I. Rawlinson, C. Tronci. arXiv:2012.03569 [physics.chem-ph] (2021)

65. H.D. Meyer, U.M.L.S. Cederbaum, Chem. Phys. Lett. **120**, 73 (1990)

66. R. Binder, D. Lauvergnat, I. Burghardt, Phys. Rev. Lett. **120**, 227401 (2018)

67. H. Wang, J. Phys. Chem. A **119**, 7951 (2015)

68. H.D. Meyer, G.A. Worth, Theor. Chem. Acc. **109**, 251 (2003)

69. P.J. Castro, A. Perveaux, D. Lauvergnat, M. Reguero, B. Lasorne, Chem. Phys. **509**, 30 (2018)

70. C.L. Lopreore, R.E. Wyatt, Phys. Rev. Lett. **82**, 5190 (1999)

71. B. Gu, I. Franco, J. Chem. Phys. **146**, 194104 (2017)

72. B.F.E. Curchod, C. Rauer, P. Marquetand, L. González, T.J. Martínez, J. Chem. Phys. **144**, 101102 (2016)

73. B. Mignolet, B.F.E. Curchod, J. Phys. Chem. A **123**, 3582 (2019)

74. V. Makhova, D.V. Shalashilin, Chem. Phys. **515**, 46 (2018)

75. I. Horenko, B. Schmidt, C. Schütte, J. Chem. Phys. **115**, 5733 (2001)

76. Z. Zhou, H.T. Chen, A. Nitzan, J.E. Subotnik, J. Chem. Theory Comput. **16**, 821 (2020)

77. T. Fiedlschuster, J. Handt, R. Schmidt, Phys. Rev. A **93**, 053409 (2016)

78. J.J. Bajo, J. González-Vázquez, I.R. Sola, J. Santa-maria, M. Richter, P. Marquetand, L. González, J. Phys. Chem. A **116**, 2800 (2012)

79. M.K.G. Subramanian, R. Santra, R. Welsch, Phys. Rev. A **102**, 013107 (2020)

80. J. Kim, H. Tao, T.J. Martinez, P. Bucksbaum, J. Phys. B: At. Mol. Opt. Phys. **48**, 164003 (2015)

81. A. Giusti-Suzo, X. He, O. Atabek, F.H. Mies, Phys. Rev. Lett. **64**, 515 (1990)

82. J.H. Posthumus, J. Plumridge, L.J. Frasinski, K. Codling, E.J. Divall, A.J. Langley, P.F. Taday, J. Phys. B: At. Mol. Opt. Phys. **33**, L563 (2000)

83. X. He, O. Atabek, A. Giusti-Suzor, Phys. Rev. A **38**, 5596 (1988)

84. G. Jolicard, O. Atabek, Phys. Rev. A **46**, 5845 (1992)

85. P.H. Bucksbaum, A. Zavriyev, H.G. Muller, D.W. Schumacher, Phys. Rev. Lett. **64**, 1883 (1990)

86. L.J. Frasinski, J. Plumridge, J.H. Posthumus, K. Codling, P.F. Taday, E.J. Divall, A.J. Langley, Phys. Rev. Lett. **86**, 2541 (2001)

87. M. Uhlmann, T. Kunert, R. Schmidt, Phys. Rev. A **72**, 045402 (2005)

88. H. Sambe, Phys. Rev. A **7**, 2203 (1973)

89. A. Scherrer, F. Agostini, D. Sebastiani, E.K.U. Gross, R. Vuilleumier, J. Chem. Phys. **143**, 074106 (2015)

90. A. Scherrer, F. Agostini, D. Sebastiani, E.K.U. Gross, R. Vuilleumier, Phys. Rev. X **7**, 031035 (2017)

91. S. Garashchuk, V.A. Rassolov, J. Chem. Phys. **120**, 1181 (2004)

92. S. Garashchuk, V.A. Rassolov, Chem. Phys. Lett. **376**, 358 (2003)

93. B. Fu, Y.C. Han, J.M. Bowman, L. Angelucci, N. Balucani, F. Leonori, P. Casavecchia, Proc. Nat. Acad. Sci. USA **109**, 9733 (2012)

94. D.A. Fedorov, A.O. Lykhin, S.A. Varganov, J. Phys. Chem. A **122**, 3480 (2018)

95. D.A. Fedorov, S.R. Pruitt, K. Keipert, M.S. Gordon, S.A. Varganov, J. Phys. Chem. A **120**, 2911 (2016)

96. M. Richter, P. Marquetand, J. González-Vázquez, I. Sola, L. González, J. Chem. Theory Comput. **7**, 1253 (2011)

97. S. Mai, P. Marquetand, L. González, Int. J. Quant. Chem. **115**, 1215 (2015)

98. G. Cui, W. Thiel, J. Chem. Phys. **141**, 124101 (2014)

99. Y. Kobayashi, K.F. Chang, T. Zeng, D.M. Neumark, S.R. Leone, Science **365**, 79 (2019)

100. C. Sanz-Sanz, G.W. Richings, G.A. Worth, Faraday Discuss. **153**, 275 (2011)

101. H. Guo, J. Chem. Phys. **99**, 1685 (1993)
102. G.E.P. Box, M.E. Muller, Ann. Math. Stat. **29**, 610 (1958)

103. A. Schild, F. Agostini, E.K.U. Gross, J. Phys. Chem. A **120**, 3316 (2016)

104. E. Khosravi, A. Abedi, A. Rubio, N.T. Maitra, Phys. Chem. Chem. Phys. **19**, 8269 (2017)

105. Y. Suzuki, A. Abedi, N.T. Maitra, E.K.U. Gross, Phys. Chem. Chem. Phys. **17**, 29271 (2015)

106. J. Kocák, E. Kraisler, A. Schild. J. Phys. Chem. Lett. **12**, 3204 (2021)

107. J. Kocák, A. Schild, Phys. Rev. Res. **2**, 043365 (2020)

108. N.I. Gidopoulos, E.K.U. Gross, Philos. Trans. R. Soc. A **372**, 20130059 (2014)

109. G. Hunter, Int. J. Quant. Chem. **9**, 237 (1975)

110. L.S. Cederbaum, J. Chem. Phys. **138**, 224110 (2013)

111. R. Requist, E.K.U. Gross, Phys. Rev. Lett. **117**, 193001 (2019)

112. C. Li, R. Requist, E.K.U. Gross, J. Chem. Phys. **148**, 084110 (2018)

113. E.K.U.G. Ryan Requist. arXiv:1909.07933 [cond-mat.str-el] (2019)

114. L. Lacombe, N.T. Maitra, Phys. Rev. Lett. **124**, 206401 (2020)

115. A. Abedi, E. Khosravi, I. Tokatly, Euro. Phys. J. B **91**, 194 (2018)

116. N.M. Hoffmann, H. Appel, A. Rubio, N. Maitra, Euro. Phys. J. B **91**, 180 (2018)

117. CPMD, copyright 2000-2019, jointly by IBM Corp. and by Max Planck Institute, Stuttgart. http://www.cpmd.org/

118. M. Filatov, S.K. Min, K.S. Kim, J. Chem. Theory Com-put. **14**, 4499 (2018)

119. M. Filatov, S.K. Min, C.H. Choi, Phys. Chem. Chem. Phys. **21**, 2489 (2019)

120. M. Filatov, S.K. Min, K.S. Kim, Mol. Phys. **117**, 1128 (2019)