# Nonadiabatic dynamics with classical trajectories: The problem of an initial coherent superposition of electronic states

Evaristo Villaseco Arribas, Neepa T Maitra, Federica Agostini

► To cite this version:

Evaristo Villaseco Arribas, Neepa T Maitra, Federica Agostini. Nonadiabatic dynamics with classical trajectories: The problem of an initial coherent superposition of electronic states. The Journal of Chemical Physics, 2024, 160 (5), 10.1063/5.0186984. hal-05287016

HAL Id: hal-05287016

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

Submitted on 27 Sep 2025

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 Attribution 4.0 International License
# Nonadiabatic dynamics with classical trajectories: The problem of an initial coherent superposition of electronic states

Evaristo Villaseco Arribas

*Department of Physics, Rutgers University, Newark 07102, New Jersey, USA and*
*Université Paris-Saclay, CNRS, Institut de Chimie Physique UMR8000, 91405 Orsay, France*

Neepa T. Maitra

*Department of Physics, Rutgers University, Newark 07102, New Jersey, USA*

Federica Agostini

*Université Paris-Saclay, CNRS, Institut de Chimie Physique UMR8000, 91405 Orsay, France*

(*Electronic address: federica.agostini@universite-paris-saclay.fr)

(*Electronic address: neepa.maitra@rutgers.edu)

(*Electronic address: evaristo.villaseco@rutgers.edu)

(Dated: January 10, 2024)

Advances in coherent light sources and development of pump-probe techniques in recent decades have opened the way to study electronic motion in its natural time-scale. When an ultrashort laser pulse interacts with a molecular target a coherent superposition of electronic states is created and the triggered electron dynamics is coupled to the nuclear motion. A natural and computationally efficient choice to simulate this correlated dynamics is a trajectory-based method where the quantum-mechanical electronic evolution is coupled to a classical-like nuclear dynamics. These methods must approximate the initial correlated electron-nuclear state by associating an initial electronic wavefunction to each classical trajectory in the ensemble. Different possibilities exist that reproduce the initial populations of the exact molecular wavefunction when represented in a basis. We show that different choices yield different dynamics, and explore the effect of this choice in Ehrenfest, surface-hopping, and exact-factorization-based coupled-trajectory schemes in a one-dimensional two-electronic-state model system that can be solved numerically exactly. This work aims to clarify the problems that standard trajectory-based techniques might have when a coherent superposition of electronic states is created to initialize the dynamics, to discuss what properties and observables are affected by different choices of electronic initial conditions, and to point out the importance of the quantum-momentum-induced electronic transitions in the coupled-trajectory schemes.

## I. INTRODUCTION

While at the end of the 20th century femtosecond laser pulses made it possible to probe and image molecular motion [1, 2], the development of extreme ultra-violet (XUV) attosecond laser pulses from high harmonic generation in the beginning of the 21st century [3, 4] enabled real-time observation of electron dynamics [5, 6] sparking the field of attochemistry [7–10]
and culminating with the 2023 Nobel prize in physics awarded to Pierre Agostini, Ferenc Krausz and Anne L'Huillier for the development of experimental methods that generate attosecond pulses of light for the study of electron dynamics in matter. When a molecular system gets ionized by an XUV pulse, a coherent superposition of cationic states is created (electronic wavepacket) which will then evolve in time in the ultrafast time regime (10⁻¹⁵ – 10⁻¹⁸ s). While in some cases the electronic charge redistribution occurs purely driven by electron correlation [11–13] before the nuclear dynamics takes place (charge migration), many are strongly affected by the nuclear motion [14–19] (charge transfer). Thus, accurate simulation of these processes involve describing the coupled motion of electrons and nuclei.

The theoretical description of these ultrafast processes remains very challenging. Ideally, the simulation would include the laser pulse as part of the dynamics, starting from the ground state of the molecule. However, if a correlated wavefunction method is used, this becomes numerically costly; if instead time-dependent density functional theory is used, approximate functionals perform poorly when the physical state evolves away from being dominated by a single Slater determinant [20–22]. Even if one were to use an efficient electronic structure method, difficulties in dealing with electronic continuum states [23–25] tend to discourage simulations of the actual ionization/excitation step. Therefore, usually an instantaneous excitation is assumed to initiate the evolution, invoking the so called "sudden approximation" [26–29]. The dynamics is launched from a chosen superposition of electronic eigenstates selected from some physical arguments, e.g., in ionization, the wavepacket corresponds to the sudden removal of an electron from an orbital in the correlated ground state wavefunction of the neutral molecule [30]. More generally, in many cases one could obtain this superposition from first-order perturbation theory, where in the length gauge and dipole approximation we write $\hat{H}(\underline{\mathbf{r}}, \underline{\mathbf{R}}, t) = \hat{H}_0(\underline{\mathbf{r}}, \underline{\mathbf{R}}) + \mathcal{E}(t)\mathbf{n} \cdot \hat{\mathbf{d}}(\underline{\mathbf{r}}, \underline{\mathbf{R}})$, where $\hat{H}_0(\underline{\mathbf{r}}, \underline{\mathbf{R}})$ is the molecular Hamiltonian with eigenstates $\Psi_k(\underline{\mathbf{r}}, \underline{\mathbf{R}})$ satisfying $\hat{H}_0(\underline{\mathbf{r}}, \underline{\mathbf{R}})\Psi_k(\underline{\mathbf{r}}, \underline{\mathbf{R}}) = E_k\Psi_k(\underline{\mathbf{r}}, \underline{\mathbf{R}})$. The perturbation is mediated by the dipole operator $\hat{\mathbf{d}} = \sum_i^{N_e} \hat{\mathbf{r}}_i + \sum_v^{N_v} \hat{\mathbf{R}}_v$ and by the field parameters, $\mathbf{n}$ being the polarization direction of the incident field, and $\mathcal{E}(t)$ a temporal function describing the laser pulse. Then the distribution of excited states after the pulse is directly related to the frequency distribution of the pulse and the dipole couplings between the field-free molecular eigenstates. In this way, the contribution of the k-th state to the excitation reads [31]

$$c_k(t) = -\frac{i}{\hbar} \langle \Psi_k(\underline{\mathbf{r}}, \underline{\mathbf{R}}) | \mathbf{n} \cdot \hat{\mathbf{d}} | \Psi_0(\underline{\mathbf{r}}, \underline{\mathbf{R}}) \rangle_{\underline{\mathbf{r}}=\underline{\mathbf{R}}} \int_{t_0}^{t} \mathcal{E}(t') e^{i\omega_{k0}t'} dt' \quad (1)$$

where $\omega_{k0} = (E_k - E_0)/\hbar$. However even within this modeling of the initial state, we encounter a challenge when trajectory-based methods are used, as we will describe shortly.

After the excitation process, the ultrafast electron-nuclear non-adiabatic dynamics beyond the Born-Oppenheimer approximation [32] begins. Since its quantum-mechanical treatment is limited to a few degrees of freedom, a large number of ultrafast simulations of charge dynamics are done with fixed nuclei, but nuclear motion may in reality play a relevant role even in the early stages. The coupling to the nuclear motion is a defining factor in attochemistry, where the goal is to control chemical reactions by steering the electron dynamics using attosecond laser pulses. While the large bandwidth of the attosecond pump pulse creates a coherent superposition of electronic eigenstates, the control of the evolution of such an electronic wavepacket can be achieved potentially using a second laser pulse, i.e., a probe pulse [33], or via a coupling with nuclei that recouples the states, leading
to coherence-dependent dynamics, with coherences lasting long enough to be able to direct the reactivity. Thus, an accurate description of the evolution of the coherence, including accounting for decoherence, is essential for understanding and controlling the resulting reaction products. The electron-nuclear coupling causes decoherence affecting the nature and time-scale of the oscillations in the electron density [34–36]. Approximate methods that treat the nuclei in terms of classical-like trajectories are the natural choice for simulations of molecular systems due to their affordable computational cost. In these schemes, an ensemble of classical-like nuclear trajectories is propagated, each associated with an instantaneous time-dependent electronic wavefunction typically expressed in a Born-Huang (BH) expansion.

Returning now to the challenge mentioned of the initial state, we note that while the choice of initial nuclear positions and momenta is often determined from sampling the Wigner distribution of the ground-state nuclear wavepacket [37] under the assumption of a vertical excitation (which neglects the dynamics occurring during the pulse and the position dependence of the transition dipoles [28, 29]) suitably constructing the initial electronic state in such an approximate picture of the process is a critical step to realistically simulate the coupled electron-nuclear dynamics after the photo-excitation/ionization step. In most non-adiabatic simulations, all the trajectories are initialized in a given excited state. However, in a realistic photo-excitation process a laser pulse of certain width interacts with the molecular target and can create a coherent superposition of electronic eigenstates. Even with knowledge of the exact wavefunction from an evaluation of all the excitation and ionization probability amplitudes after the pulse, various ways to initialize the electronic state(s) associated to the trajectories ensemble can be proposed, and, as we will discuss in the following sections, this choice is strongly dependent of the non-adiabatic dynamics method used to simulate the coupled electron-nuclear dynamics.

Aiming to address the issue of initial coherent superpositions of electronic states in a trajectory-based description of ultrafast processes, we study the impact that different initialization schemes have on the dynamics of a two-electronic-state one-dimensional model. We compare the performance of Ehrenfest (Eh) dynamics, a widely-used approach in these situations [38–44], and the coupled-trajectory mixed quantum-classical algorithm (CTMQC) [45–47], with the following two choices of initial electronic state. In one, which we refer to as “pure”, all trajectories are initially associated to the same coherent superposition of the two electronic states. In the other, which we refer to as a “mixed”, each trajectory is associated with a single adiabatic electronic state, such that the distribution of trajectories reproduces the same adiabatic populations of the two electronic states as in the pure case. Note that the Ehrenfest approach is often used when employing real-time time-dependent density functional theory (TDDFT) for the electronic dynamics. Specifically, the key advantage of Ehrenfest dynamics is that it can be formulated in terms of the electron density avoiding the explicit calculation of adiabatic states, and, thus, can be naturally combined with TDDFT. This approach (TDDFT-Eh) has been widely used to study coupled electron-nuclear dynamics initiated by a photo-ionization process that creates a coherent superposition of cationic states [40, 42–44]. We also test the performance of the Tully surface hopping (TSH) scheme with and without energy-decoherence correction [48], and its coupled-trajectory version [49, 50], CTTSH, derived from the exact factorization [51, 52].

The paper is organized as follows. In Section II we begin by reviewing the density matrix formalism and by defining the reduced electronic density matrix of an ensemble of trajectories, distinguishing between a pure and a mixed case. In Section III
we provide a brief overview of the exact factorization approach introducing the coupled-trajectory methods derived from it. Section IV is devoted to the numerical simulations comparing the effect that the choice of the initial state has on the dynamics of an exactly solvable model system. Finally, section V gathers the main conclusions of our study.

## II. DEFINITION OF PURE ENSEMBLE AND MIXED ENSEMBLE VIA THE INITIAL CONDITIONS

In this section, we discuss the notion of a pure ensemble and a mixed ensemble of electronic initial conditions in the context of a trajectory-based description of the nuclear dynamics. To this end, we will employ the density-matrix formalism. Therefore, let us first define the density operator of the full molecular system, which is described by the molecular wavefunction $\Psi(\underline{\mathbf{R}}, \underline{\mathbf{R}}, t) = \langle \underline{\mathbf{R}}, \underline{\mathbf{R}} | \Psi(t) \rangle$. Using the exact factorization [51–60], we express this as $\Psi(\underline{\mathbf{R}}, \underline{\mathbf{R}}, t) = \chi(\underline{\mathbf{R}}, t) \phi_{\underline{\mathbf{R}}}(\underline{\mathbf{R}}, t)$ with the partial normalization condition $\langle \phi_{\underline{\mathbf{R}}}(t) | \phi_{\underline{\mathbf{R}}}(t) \rangle_{\underline{\mathbf{R}}} = 1$ for all $\underline{\mathbf{R}}, t$. The density operator reads

$$ \hat{\Gamma}(t) = |\Psi(t)\rangle\langle\Psi(t)| \qquad (2) $$

whose matrix elements in the nuclear configuration space are

$$ \hat{\Gamma}(\underline{\mathbf{R}}, \underline{\mathbf{R}}, t) = \langle \underline{\mathbf{R}}' | \hat{\Gamma}(t) | \underline{\mathbf{R}} \rangle = \chi^*(\underline{\mathbf{R}}, t) \chi(\underline{\mathbf{R}}', t) |\phi_{\underline{\mathbf{R}}'}(t)\rangle \langle \phi_{\underline{\mathbf{R}}}(t)| \qquad (3) $$

Considering only the diagonal part in $\underline{\mathbf{R}}$-space and expressing it in the adiabatic basis, i.e., the set of electronic eigenstates of the Born-Oppenheimer Hamiltonian that parametrically depend on the nuclear positions, $\{|\phi_{\underline{\mathbf{R}}}(t)\rangle\}, i = 1, \dots, n_{\text{states}}$, we find

$$ \Gamma_{ij}(\underline{\mathbf{R}}, t) = \langle \phi_{\underline{\mathbf{R}}}^i | \hat{\Gamma}(\underline{\mathbf{R}}, \underline{\mathbf{R}}, t) | \phi_{\underline{\mathbf{R}}}^j \rangle = |\chi(\underline{\mathbf{R}}, t)|^2 C_i(\underline{\mathbf{R}}, t) C_j^*(\underline{\mathbf{R}}, t) = |\chi(\underline{\mathbf{R}}, t)|^2 \rho_{ji}(\underline{\mathbf{R}}, t) \qquad (4) $$

where we introduced a BH-like expansion of the electronic time-dependent state of the exact factorization, i.e., $|\phi_{\underline{\mathbf{R}}}(t)\rangle = \sum_i C_i(\underline{\mathbf{R}}, t) |\phi_{\underline{\mathbf{R}}}^i\rangle$, and we defined the conditional electronic density matrix elements $\rho_{ij}(\underline{\mathbf{R}}, t) = C_i^*(\underline{\mathbf{R}}, t) C_j(\underline{\mathbf{R}}, t)$, usually referred to as “coherences”. In the adiabatic basis, the elements of the reduced density matrix for the electronic subsystem $\hat{\gamma}(t) = \text{Tr}_{\underline{\mathbf{R}}} (\hat{\Gamma}(t))$ are

$$ \gamma_{ij}(t) = \text{Tr}_{\underline{\mathbf{R}}} (\Gamma_{ij}(\underline{\mathbf{R}}, t)) = \int d\underline{\mathbf{R}} |\chi(\underline{\mathbf{R}}, t)|^2 \rho_{ji}(\underline{\mathbf{R}}, t) \qquad (5) $$

In a trajectory-based mixed quantum-classical (MQC) description, the nuclear density is represented by an ensemble of trajectories, such that

$$ |\chi(\underline{\mathbf{R}}, t)|^2 \stackrel{\text{MQC}}{=} \frac{1}{N_{\text{tr}}} \sum_{\alpha=1}^{N_{\text{tr}}} \delta(\underline{\mathbf{R}} - \underline{\mathbf{R}}^{(\alpha)}(t)) \qquad (6) $$
with $N_{tr}$ the total number of trajectories and $\underline{\mathbf{R}}^{(\alpha)}(t)$ indicating their positions. This yields the MQC electronic reduced density matrix as

$$ \gamma_{ij}^{MQC}(t) = \frac{1}{N_{tr}} \sum_{\alpha=1}^{N_{tr}} \rho_{ji}^{(\alpha)}(t) \qquad (7) $$

Henceforth, we will use the shorthand notation $f^{\alpha}(t)$ and $f^{\alpha}$, to indicate the values of a function $f(\underline{\mathbf{R}}, t)$ at the position of the trajectory $\alpha$ at time $t$ and for the value of a function $f(\underline{\mathbf{R}})$ at the position of the trajectory $\alpha$, respectively.

Let us now consider the purity of the electronic subsystem through

$$ \mathrm{Tr}(\hat{\gamma}^2(t)) = \sum_{i,k} \gamma_{ik}(t)\gamma_{ki}(t) = \int d\underline{\mathbf{R}} d\underline{\mathbf{R}}' |\chi(\underline{\mathbf{R}}, t)|^2 |\chi(\underline{\mathbf{R}}', t)|^2 \sum_{i,k} \rho_{ki}(\underline{\mathbf{R}}, t)\rho_{ik}(\underline{\mathbf{R}}', t) \qquad (8) $$

$$ \stackrel{MQC}{=} \frac{1}{N_{tr}^2} \sum_{i,k} \sum_{\alpha=1}^{N_{tr}} \rho_{ki}^{(\alpha)}(t) \sum_{\beta=1}^{N_{tr}} \rho_{ik}^{(\beta)}(t) \qquad (9) $$

Equation (8), i.e., the quantum-mechanical expression of the purity, can yield a value smaller than one, as a result of the partial-trace operation in the definition of Eq. (5), attesting to a mixed-state character of the electronic subsystem, even when the underlying molecular state is pure, (i.e., $\mathrm{Tr}_{\underline{\mathbf{R}},\underline{\mathbf{R}}'}(\hat{\Gamma}^2) = \mathrm{Tr}_{\underline{\mathbf{R}},\underline{\mathbf{R}}'}(\hat{\Gamma}) = 1$). It may be instructive to note that if the conditional electronic coefficients are somehow independent of $\underline{\mathbf{R}}$, then $\mathrm{Tr}\hat{\gamma}^2 = 1$, as if the electronic subsystem were a pure state.

Equation (9), i.e., the MQC expression of the purity, yields a purity equal to one at time $t$ if all the trajectories have the same electronic coefficients, i.e., the coherences $\rho_{ij}^{(\alpha)}(t) = \rho_{ij}(t)$ are independent of $\alpha$, which is in a sense analogous to the special case of the fully quantum case just considered. Then, in such a case

$$ \sum_{\alpha=1}^{N_{tr}} \rho_{ki}^{(\alpha)}(t) \sum_{\beta=1}^{N_{tr}} \rho_{ik}^{(\beta)}(t) = N_{tr}^2 \rho_{ki}(t) \rho_{ik}(t) \qquad (10) $$

and this situation allows us to define a pure ensemble since

$$ \mathrm{Tr}(\hat{\gamma}_{\mathrm{pure}}^2(t)) = \frac{1}{N_{tr}^2} \sum_{i,k} N_{tr}^2 \rho_{ki}(t) \rho_{ik}(t) = \sum_i |C_i(t)|^2 \sum_k |C_k(t)|^2 = 1. \qquad (11) $$

Note that $\sum_i |C_i^{\alpha}(t)|^2 = 1 \forall \alpha, t$ follows from the partial normalization condition of the electronic conditional wavefunction of the exact factorization. Such a pure-ensemble situation occurs at $t = 0$, for instance, when all the trajectories are initialized in a single electronic eigenstate, or when they are all initialized in the same coherent superposition of adiabatic eigenstates. On the other hand, if at a certain time $t$ different trajectories are associated with different electronic coefficients, then the concept of a mixed ensemble arises, as follows. Consider first a two-component ensemble, and let $N_I$ and $N_J$ be the numbers of trajectories in two different (superpositions of) electronic states, say $I$ and $J$, such that $N_I + N_J = N_{tr}$. For the superposition $I(J)$, we indicate
the coherences as $\rho_{ki}^I(t)$ ($\rho_{ki}^J(t)$), with all coherences the same within the ensemble I (J). For this situation, we have

$$ \sum_{\alpha=1}^{N_{tr}} \rho_{ki}^{(\alpha)}(t) \sum_{\beta=1}^{N_{tr}} \rho_{ik}^{(\beta)}(t) = (N_I \rho_{ki}^I + N_J \rho_{ki}^J)(N_I \rho_{ik}^I + N_J \rho_{ik}^J) \quad (12) $$

by explicitly calculating the sums over the trajectories. The purity in such a situation is

$$ \mathrm{Tr} (\hat{\gamma}_{\mathrm{mixed}}^2(t)) = \frac{1}{N_{\mathrm{tr}}^2} \sum_{i,k} [N_I^2 \rho_{ki}^I(t) \rho_{ik}^I(t) + N_J^2 \rho_{ki}^J(t) \rho_{ik}^J(t) + N_I N_J (\rho_{ki}^I \rho_{ik}^J + \rho_{ki}^J \rho_{ik}^I)] \quad (13) $$

In the first two terms in the right-hand side, the sums over i and k can be performed using the partial normalization condition, i.e., $\sum_{i,k} \rho_{ki}^I(t) \rho_{ik}^I(t) = \sum_i |C_i^I(t)|^2 \sum_k |C_k^I(t)|^2 = 1$ as in Eq. (11), thus we get

$$ \mathrm{Tr} (\hat{\gamma}_{\mathrm{mixed}}^2(t)) = \frac{N_I^2 + N_J^2}{N_{\mathrm{tr}}^2} + \frac{N_I N_J}{N_{\mathrm{tr}}^2} \sum_{i,k} (\rho_{ki}^{(I)} \rho_{ik}^{(J)} + \rho_{ki}^{(J)} \rho_{ik}^{(I)}) \neq 1 \quad (14) $$

which is in general different from unity, thus allowing us to define a mixed ensemble. In general (more than two components in the ensemble),

$$ \mathrm{Tr} (\hat{\gamma}_{\mathrm{mixed}}^2(t)) = \frac{1}{N_{\mathrm{tr}}^2} \left( \sum_I N_I^2(t) + 2 \sum_{I,J \neq I} N_I(t) N_J(t) \sum_{ik} \mathrm{Re} (\rho_{ki}^{(I)} \rho_{ik}^{(J)}) \right) \neq 1 \quad (15) $$

Initializing the trajectories in a pure ensemble or in a mixed ensemble may yield the same (trajectory-averaged) electronic populations with, however, different electronic coherences and different value of the purity.

When the MQC dynamics is simulated using Ehrenfest-based methods, the most natural choice for the initialization of the dynamics is the pure ensemble, such that all trajectories are associated to the same coherent superposition of electronic states with coefficients independent of the trajectories positions. In surface-hopping-based methods, in contrast, one would typically choose the mixed state option, where each trajectory is associated to a single electronic state (not a coherent superposition of states), consistent with the adiabatic force driving the nuclear trajectory's evolution. In this latter choice, the off-diagonal elements of the electronic density matrix, i.e., the coherences, are zero.

For a concrete example, compare (1) the case of a two-level system in which the ensemble of trajectories is associated to the same initial superposition of electronic states, e.g., $|\phi^{(\alpha)}(0)\rangle = \frac{1}{\sqrt{2}} (|\phi^{0,(\alpha)}\rangle + e^{i\phi}|\phi^{1,(\alpha)}\rangle)$, to (2) the case where 50% of the trajectories in the ensemble are associated to one state, e.g., $|\phi^{(\alpha)}(0)\rangle = |\phi^{0,(\alpha)}\rangle$, and the other 50% to the other state, $|\phi^{(\alpha)}(0)\rangle = |\phi^{1,(\alpha)}\rangle$. The electronic density matrices of these two systems at time $t=0$ are

$$ \hat{\gamma}_{\mathrm{pure}}^{\mathrm{MQC}}(0) = \frac{1}{2} \begin{pmatrix} 1 & e^{i\phi} \\ e^{-i\phi} & 1 \end{pmatrix} \qquad \hat{\gamma}_{\mathrm{mixed}}^{\mathrm{MQC}}(0) = \frac{1}{2} \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \quad (16) $$

The diagonal elements, i.e., the populations, show that both systems have the same trajectory-averaged populations ($\frac{1}{2}$ in each state) but different coherences (off-diagonal elements). Furthermore, one can easily show that $\mathrm{Tr}(\hat{\gamma}_{\mathrm{pure}}^2) = 1$ while $\mathrm{Tr}(\hat{\gamma}_{\mathrm{mixed}}^2) =$
½. Thus, one can prepare different initial states of an ensemble of trajectories that have the same populations but differing in the coherences. As we will see in the next sections, this freedom in choosing the form of the initial state will determine the fate of the dynamics, since the time evolution of the electronic populations and the force that the trajectories experience depend on the coherences.

It is important to note that, even at the initial time, molecular properties will differ if they depend on the coherences. An example is the one-body electronic current-density

$$ 
\begin{aligned}
\mathbf{j}^e(\mathbf{r}, t) &= \frac{\hbar N_e}{m} \int d\underline{\mathbf{R}} d\mathbf{r}_2 \dots d\mathbf{r}_{N_e} \mathfrak{S} \{ \Psi^*(\underline{\mathbf{r}}, \underline{\mathbf{R}}, t) \nabla_{\mathbf{r}} \Psi(\underline{\mathbf{r}}, \underline{\mathbf{R}}, t) \} \\
&= \int d\underline{\mathbf{R}} |\chi(\underline{\mathbf{R}}, t)|^2 \mathfrak{S} \left( \sum_{l,k} C_l^*(\underline{\mathbf{R}}, t) C_k(\underline{\mathbf{R}}, t) \mathbf{j}_{lk, \underline{\mathbf{R}}}^e(\mathbf{r}) \right)
\end{aligned}
\quad (17) $$

where $N_e$ is the number of electrons, and $\mathbf{j}_{lk, \underline{\mathbf{R}}}^e(\mathbf{r}) = \frac{\hbar N_e}{m} \int d\mathbf{r}_2 \dots d\mathbf{r}_{N_e} \phi_{\underline{\mathbf{R}}}^{l*}(\mathbf{r}, \mathbf{r}_2 \dots \mathbf{r}_{N_e}) \nabla_{\mathbf{r}} \phi_{\underline{\mathbf{R}}}^k(\mathbf{r}, \mathbf{r}_2 \dots \mathbf{r}_{N_e})$. The corresponding MQC expression is

$$ \mathbf{j}^{e,MQC}(\mathbf{r}, t) = \frac{1}{N_{tr}} \sum_{\alpha=1}^{N_{tr}} \sum_{l,k} \mathfrak{S} \{ \rho_{lk}^{(\alpha)}(t) \mathbf{j}_{lk}^{e,(\alpha)}(\mathbf{r}) \} \quad (18) $$

Considering the two situations discussed above, i.e., the pure and mixed cases in Eq. (16), with equally populated states, the initial current-densities are

$$ \mathbf{j}_{\text{pure}}^{e,MQC}(\mathbf{r}, 0) = \frac{1}{2} \frac{\rho_{lk}(0)}{N_{tr}} \sum_{\alpha=1}^{N_{tr}} \mathfrak{S} \{ \mathbf{j}_{12}^{e,(\alpha)}(\mathbf{r}) \} \neq 0 \quad \text{and} \quad \mathbf{j}_{\text{mixed}}^{e,MQC}(\mathbf{r}, 0) = 0 \quad (19) $$

In the expression on the left, which is, in general, different from zero, the coherence for each trajectory is the same, i.e., independent of the trajectories index/position, $\rho_{lk}^{(\alpha)}(0) = \rho_{lk}(0) \forall \alpha$. Instead, the expression on the right is identically zero, since in the mixed-ensemble case the coherences are initially zero. The choice of an initial mixed ensemble, for which the electronic density matrix has zero off-diagonal elements in the adiabatic basis will always produce a strictly zero electronic current-density, even though it may capture the correct electronic populations.

The above discussion shows that when selecting the initial conditions for the electronic state of the system in trajectory-based simulations of non-adiabatic processes, it is important not only to consider the populations, but also how these populations are “distributed” over the trajectories, as described above, and further also the phases of the coefficients of the exact initial state that the MQC dynamics is attempting to simulate (which is a topic that emerged previously in the literature [61] and was tackled in the context of the multistate semiclassical Liouville equation [62, 63]). These affect not only the coherences and the equations of motion, as we will see in the next section, but also directly the observables, which may be wrong from the very beginning if not properly considered.
### III. EXACT-FACTORIZATION-BASED COUPLED-TRAJECTORY METHODS

In this section, we briefly review the coupled-trajectory methods derived from the exact factorization of the molecular wave-function, namely coupled-trajectory mixed quantum-classical (CTMQC) algorithm and coupled-trajectory Tully surface hopping (CTTSH). Since they are variations of the trajectory-independent Eh method and TSH method respectively, we first discuss these two schemes.

In both Eh and TSH, the quantum-mechanical electronic evolution along a trajectory $\alpha$ is given by the evolution of the expansion coefficients of the electronic wavefunction in a BH (adiabatic) basis

$$ \dot{C}_{l, \text{Eh}/\text{TSH}}^{(\alpha)}(t) = -i\epsilon_l^{(\alpha)} C_l^{(\alpha)}(t) - \sum_k \sum_{v=1}^{N_n} \dot{\mathbf{R}}_v^{(\alpha)}(t) \cdot \mathbf{d}_{v,lk}^{(\alpha)} C_k^{(\alpha)}(t) \quad (20) $$

with $\epsilon_l^{(\alpha)}$ the l-th adiabatic or Born-Oppenheimer (BO) potential energy surface (PES) evaluated at $\mathbf{R}^{(\alpha)}$ and $\mathbf{d}_{v,lk}^{(\alpha)}$ is the non-adiabatic coupling (NAC) vector along the v nuclear coordinate between BO states l and k, i.e., $\langle \phi_{\underline{R}}^l | \nabla_v \phi_{\underline{R}}^k \rangle |_{\underline{R}^{(\alpha)}}$. In TSH, the trajectories evolve on a given adiabatic surface (active state), i.e., $\mathbf{F}_{v, \text{TSH}}^{(\alpha)} = -\nabla_v \epsilon_{\text{active}(t)}^{(\alpha)}$, where the active state stochastically changes in general according to the fewest-switches probability [64]. Since each trajectory is “collapsed” onto a given surface, in order to avoid internal inconsistency, the dynamics is typically started in a mixed ensemble with zero coherences with respect to the adiabatic basis (using the language introduced in Section II) [65, 66]. Note that while TSH is generally performed in the adiabatic representation, there exist alternative representations that allow one to account explicitly for a laser pulse [67–71], though the dynamics is strongly dependent on the chosen electronic basis and the used TSH parameters [72]. On the other hand, in Eh the trajectories evolve on a mean-field surface. Defining $\hat{H}_{BO}$ as the electronic Hamiltonian, namely $\hat{H}_{BO} = \hat{H} - \hat{T}_n$ with $\hat{H}$ the molecular Hamiltonian and $\hat{T}_n$ the nuclear kinetic energy, the mean-field Eh force is

$$ \mathbf{F}_{v, \text{Eh}}^{(\alpha)}(t) = -\langle \Phi^{(\alpha)}(t) | \nabla_v \hat{H}_{BO} | \Phi^{(\alpha)}(t) \rangle = -\sum_l \rho_{ll}^{(\alpha)}(t) \nabla_v \epsilon_l^{(\alpha)} + \sum_{l,k} \rho_{lk}^{(\alpha)}(t) \Delta \epsilon_{lk}^{(\alpha)} \mathbf{d}_{v,lk}^{(\alpha)} \quad (21) $$

where $|\Phi^{(0)}(t)\rangle$ is the instantaneous electronic state along the trajectory $\alpha$. Such a force reflects a coherent evolution, thus it makes more sense, at least formally, that all trajectories are initialized in the same coherent superposition state (pure ensemble, using the language introduced in Section II), if this is the case in the “real” dynamics or, anyway, in the dynamics that one aims to reproduce.

The exact factorization (EF) of the molecular wavefunction [51–60], i.e., $\Psi(\underline{\mathbf{R}}, \underline{\mathbf{r}}, t) = \chi(\underline{\mathbf{R}}, t) \phi_{\underline{\mathbf{R}}}(\underline{\mathbf{r}}, t)$ with $\int d\underline{\mathbf{r}} |\phi_{\underline{\mathbf{R}}} (\underline{\mathbf{r}}, t)|^2 = 1 \forall \underline{\mathbf{R}}, t$, provides a different view of photochemical processes than the standard BH description [73–75]. In EF the time evolution of the nuclear wavefunction is given by a time-dependent Schrödinger equation in the presence of time-dependent scalar and vector potentials which capture the exact electron-nuclear correlation and are functionals of the time-dependent conditional
electronic wavefunction

$$i\partial_t \chi(\underline{\mathbf{R}}, t) = \left( \sum_{\mathbf{v}=1}^{N_v} \frac{(-i\nabla_{\mathbf{v}} + \mathbf{A}_{\mathbf{v}}[\phi_{\underline{\mathbf{R}}}(t)])^2}{2M_{\mathbf{v}}} + \varepsilon[\phi_{\underline{\mathbf{R}}}(t)] \right) \chi(\underline{\mathbf{R}}, t) \quad (22)$$

while the time evolution of the electronic subsystem follows a much more complicated equation [51, 59, 76]. Hence, the BH description of photochemistry in which different wavepacket components evolve on various static PESs, is replaced by one in which a time-dependent nuclear wavefunction evolves on a single time-dependent potential energy surface (TDPES) $\varepsilon(\underline{\mathbf{R}}, t)$ under the influence of a vector potential $\mathbf{A}_{\mathbf{v}}(\underline{\mathbf{R}}, t)$. Note that this factorization is unique up to a $(\underline{\mathbf{R}}, t)$-dependent gauge transformation $\chi(\underline{\mathbf{R}}, t) \rightarrow e^{-i\theta(\underline{\mathbf{R}}, t)}\chi(\underline{\mathbf{R}}, t)\Phi_{\underline{\mathbf{R}}}(\underline{\mathbf{r}}, t) \rightarrow e^{i\theta(\underline{\mathbf{R}}, t)}\Phi_{\underline{\mathbf{R}}}(\underline{\mathbf{r}}, t)$ with the potentials transforming under the gauge according to $\varepsilon(\underline{\mathbf{R}}, t) \rightarrow \varepsilon(\underline{\mathbf{R}}, t) + \partial_t \theta(\underline{\mathbf{R}}, t)$ and $\mathbf{A}_{\mathbf{v}}(\underline{\mathbf{R}}, t) = \mathbf{A}_{\mathbf{v}}(\underline{\mathbf{R}}, t) + \nabla_{\mathbf{v}}\theta(\underline{\mathbf{R}}, t)$.

The presence of these potentials, that drive the nuclear motion incorporating the coupling to the electrons, makes EF an ideal framework to develop MQC schemes for non-adiabatic dynamics. The method derived first was CTMQC [45–47, 52], which relies in expanding the time-dependent electronic conditional wavefunction in the BH basis

$$\phi_{\underline{\mathbf{R}}}(\underline{\mathbf{r}}, t) = \sum_{l} C_l(\underline{\mathbf{R}}, t) \phi_{\underline{\mathbf{R}}}^l(\underline{\mathbf{r}}) \quad (23)$$

and taking the classical limit of the nuclear motion, neglecting some terms in the equations justified by exact studies [77, 78]. In CTMQC the electronic evolution and nuclear force acting on the trajectories have an Eh component, plus a coupled-trajectory (CT) term that arises within EF:

$$\dot{C}_{l}^{(\alpha)}(t) = \dot{C}_{l,Eh/TSH}^{(\alpha)}(t) + \dot{C}_{l,CT}^{(\alpha)}(t) \quad (24)$$

$$\mathbf{F}_{\mathbf{v}}^{(\alpha)}(t) = \mathbf{F}_{\mathbf{v},\text{Eh}}^{(\alpha)}(t) + \mathbf{F}_{\mathbf{v},\text{CT}}^{(\alpha)}(t) \quad (25)$$

The coupled-trajectory terms read

$$\dot{C}_{l,CT}^{(\alpha)}(t) = \sum_{\mathbf{v}=1}^{N_v} \sum_{k} \frac{\mathbf{Q}_{\mathbf{v}}^{(\alpha)}(t)}{M_{\mathbf{v}}} \cdot \Delta\mathbf{f}_{\mathbf{v},lk}^{(\alpha)}(t) \rho_{kk}^{(\alpha)}(t) C_{l}^{(\alpha)}(t) \quad (26)$$

$$\mathbf{F}_{\mathbf{v},\text{CT}}^{(\alpha)}(t) = \sum_{\mu=1}^{N_v} \sum_{l,k} \frac{2\mathbf{Q}_{\mu}^{(\alpha)}(t)}{M_{\mu}} \cdot \mathbf{f}_{\mu,l}^{(\alpha)}(t) \rho_{ll}^{(\alpha)}(t) \rho_{kk}^{(\alpha)}(t) \Delta\mathbf{f}_{\mathbf{v},lk}^{(\alpha)}(t) \quad (27)$$

with $\Delta\mathbf{f}_{\mathbf{v},lk}^{(\alpha)}(t) = \mathbf{f}_{\mathbf{v},l}^{(\alpha)}(t) - \mathbf{f}_{\mathbf{v},k}^{(\alpha)}(t)$, where $\mathbf{f}_{\mathbf{v},l}^{(\alpha)}(t)$ is the time-integrated adiabatic force on nucleus v accumulated on the l-th BO
PES along the trajectory α (i.e., an adiabatic momentum)

$$ \mathbf{f}_{v,l}^{(\alpha)}(t) = - \int_{0}^{t} \nabla_{v} \epsilon_{l}^{(\alpha)} d\tau, \qquad (28) $$

and $\mathbf{Q}_{v}^{(\alpha)}(t)$ is the nuclear quantum momentum evaluated at the position of the trajectory $\mathbf{R}^{(\alpha)}(t)$

$$ \mathbf{Q}_{v}^{(\alpha)}(t) = - \frac{\nabla_{v} |\chi^{(\alpha)}(t)|^2}{2 |\chi^{(\alpha)}(t)|^2} \qquad (29) $$

The quantum momentum measures the spatial (de)localization of the nuclear density over time. This quantity is normally computed as a linear function by imposing the net contribution of Eq. (26) to be zero over the whole ensemble, which avoids net population transfer away from NAC regions [46, 47, 60]. For a given trajectory, the nuclear quantum momentum requires non-local information on the position and electronic coefficients of the rest of trajectories in the ensemble, and therefore the trajectories are coupled and cannot be propagated independently. Note that although the trajectory coupling adds computational cost with respect to independent-trajectory methods, efficient parallelization schemes could mitigate this. Further note that the coupling among the trajectories does not require additional electronic structure calculations, which is generally the computational bottleneck in non-adiabatic simulations of molecules. The extra term in the electronic evolution Eq. (26), when compared to TSH or Eh, is not only able to rigorously recover decoherence effects lacking in those methods, but also yield quantum-momentum-induced population transfers which are key to capture multi-state dynamics [79, 80]. The accumulated force can be redefined to satisfy energy conservation over the ensemble of trajectories as done in the CTMQC-E algorithm [81, 82].

The CTMQC electronic evolution equation has been employed within a TSH framework for the nuclei, resulting in the CTTSH algorithm [49, 50, 83] and surface hopping based on the exact factorization (SHXF) [84–86] methods

$$ \dot{\mathbf{C}}_{l}^{(\alpha)}(t) = \dot{\mathbf{C}}_{l,Eh/TSH}^{(\alpha)}(t) + \dot{\mathbf{C}}_{l,CT}^{(\alpha)} \qquad (30) $$

$$ \mathbf{F}_{v}^{(\alpha)}(t) = -\nabla_{v} \epsilon_{\text{active}(t)}^{(\alpha)} \qquad (31) $$

In SHXF, auxiliary trajectories are used to approximate the quantum momentum in order to produce an independent-trajectory scheme. On the other hand, in CTTSH the quantum momentum is evaluated through coupled trajectories in the same way as in CTMQC. CTTSH scales linearly with system size and quadratically with the number of trajectories [87]. Studies on the properties of the exact TDPES [77, 88, 89] revealed it develops adiabatic shapes after passing through a coupling region, with steps connecting regions of different adiabatic slopes and while acquiring a more diabatic shape in coupling regions. Therefore the main motivation behind these EF-based TSH schemes is to approximate the nuclear propagation on the TDPES via purely adiabatic dynamics, where the surface hops are set to approximately mimic the effect of the steps. Since in CTTSH the electronic dynamics is given by an evolution equation that is derived from EF, one may argue that the electronic coefficients
should be taken to be the measure of the electronic populations instead of the more usual distribution of the active states for the trajectories. One may expect that, unlike in TSH, in CTTSH trajectories should be initialized in a pure ensemble, such that each trajectory in the ensemble carries a coherent superposition, which would eliminate the problem of initially incorrect observables such as vanishing electronic current-densities, that beginning in a mixed ensemble suffers from. The active states can, then, be stochastically selected to match the quantum populations, which not only achieves internal consistency at the initial time, but also yields better sampling of the energy landscape, as we will demonstrate with our numerical study in the next section. It is important to note that whether a pure or mixed ensemble is chosen has a significant impact on the dynamics from the very start, since the time-evolution of the electronic populations and the forces that the trajectories experience (Eq. (21) and Eq. (27)) depend on the electronic coherences.

## IV. NUMERICAL RESULTS

Our model system consists on two coupled one-dimensional harmonic oscillators shifted in position and displaced in energy whose BO Hamiltonian is given in the diabatic basis as

$$H_{BO}(R) = \begin{pmatrix} \frac{1}{2}k(R-R_1)^2 & be^{-a(R-R_3)^2} \\ be^{-a(R-R_3)^2} & \frac{1}{2}k(R-R_2)^2 + \Delta \end{pmatrix} \qquad (32)$$

The values of the parameters chosen for this system are in atomic units: $k = 0.02 \text{ Ha}^2\text{m}_e/\hbar^2$, $a = 3.0 \text{ bohr}^{-2}$, $b = 0.01 \text{ Ha}$,

<figure><img src="image_2.png" /><figcaption>FIG. 1: Diabatic (dashed) and adiabatic (solid) potential energy curves, scaled NAC between ground and excited states (blue dashed) and ground (green) and excited (magenta) nuclear densities at the initial time.</figcaption></figure>

$R_1 = 6.0 \text{ bohr}$, $R_2 = 2.0 \text{ bohr}$ (identifying the Franck-Condon region where the instantaneous excitation takes place), $R_3 = 3.875 \text{ bohr}$, $\Delta = 0.01 \text{ Ha}$ and the nuclear mass $M = 20000 \text{ m}_e$. Figure 1 shows the diabatic and adiabatic potential energy curves for this system together with the NAC between ground and excited state. This model was used to test different coupled-trajectory exact-factorization-based schemes, when beginning purely in the excited state, showing the importance of the trajectory coupling to capture decoherence behavior and non-adiabatic interferences [49]. Propagating the coupled trajectories with a purely adiabatic
force (CTTSH) improved the long-time dynamics behavior of CTMQC due to noisy scalar and vector potentials [90] resulting from accumulated errors in the phase of the electronic coefficients. We will study the evolution of this system in a situation where we supposed that after interaction with a short laser pulse 20% of the probability density initially located in the ground state, gets vertically promoted to the first excited state creating a coherent superposition (see below), and simulate the dynamics starting from this state, with the field turned off.

For the quantum dynamics (QD) simulations, the time-dependent Schrödinger equation is solved on a grid in the diabatic basis using the split-operator method [91, 92]. The initial state (just after the pulse) is taken to be:

$$|\Psi(R, 0)\rangle = e^{-\frac{(R-R_0)}{2\sigma^2}} \left( \sqrt{0.8}|\phi_R^0\rangle + \sqrt{0.2}|\phi_R^1\rangle \right) \quad (33)$$

where two Gaussian wavepackets with centers at $R_0 = 2$ bohr and variances $\sigma = 0.2236$ bohr are located in the two states and integrate to 0.8 in the ground state and to 0.2 in the excited state (they are green and magenta densities in Fig. 1). In our analysis, we will show the exact TDPES as determined by our QD simulations. This is obtained by computing the time-dependent conditional electronic wavefunction via inversion from the full molecular wavefunction [77] in a gauge in which $A_V(\underline{R}, t) = 0$ [88, 89], i.e., the time-evolution of $\chi(\underline{R}, t)$ is determined entirely by the TDPES $\epsilon(\underline{R}, t)^1$.

In our trajectory-based description we have a choice of whether we use a pure ensemble or a mixed ensemble to initialize the simulations (Section II), the only constraint being that the trajectory-averaged populations at the initial time are

$$\frac{1}{N_{tr}} \sum_{\alpha}^{N_{tr}} \rho_{00}^{(\alpha)}(0) = 0.8 \quad \text{and} \quad \frac{1}{N_{tr}} \sum_{\alpha}^{N_{tr}} \rho_{11}^{(\alpha)}(0) = 0.2. \qquad (34)$$

We will first study the performance of Eh-based methods, i.e., of CTMQC and pure Eh, in the two situations identified in Section II as pure ensemble and mixed ensemble: In the former case, all trajectories in the ensemble are associated to the same coherent superposition of electronic states, thus to the same electronic coefficients at the initial time as given in Eq. (33). In the latter case, we have a two-component ensemble, with 80% of the trajectories in the ensemble with electronic coefficients being 1 in the electronic ground state and 0 otherwise and the remaining 20% with electronic coefficient being 1 in the excited state and 0 otherwise; thus the two sets of trajectories are associated to different electronic wavefunctions. Further, we also test the performance of TSH-based methods, namely CTTSH, TSH and TSH with energy decoherence correction (TSHEDC) with decoherence parameter of $C = 0.1$ Ha [48], using again a pure ensemble and a mixed ensemble to initialize the dynamics. In the case of a pure ensemble, the active states are stochastically selected to match the net populations at the initial time. For Eh and CTMQC, $N_{tr} = 2000$ trajectories were used for the mixed ensemble calculations while convergence was achieved with $N_{tr} = 1000$ trajectories for the pure ensemble. $N_{tr} = 4000$ trajectories were used for the TSH-based mixed ensemble calculations, while $N_{tr} = 2000$ trajectories were run for the pure ensemble. The initial positions and momenta of the trajectories are sampled

<footnote n="1">¹ The TDPES has a gauge-dependent and a gauge-independent contribution $\epsilon(\underline{R}, t) = \epsilon_{GD}(\underline{R}, t) + \epsilon_{GI}(\underline{R}, t)$. In our analysis we just show the gauge-invariant part of the TDPES, due to the different gauge choices used in the exact calculations and in CTMQC.</footnote>
from a Wigner distribution of the ground state nuclear wavepacket. The time-step used in the simulations was $dt = 0.1 \hbar/Ha$, and all the trajectory simulations were performed in the G-CTMQC package [93].

The excited state populations and coherences in the QD simulations are obtained via numerical integration

$$ 
\begin{aligned}
P_1^{\mathrm{QD}}(t) &= \int dR |\chi_1(R,t)|^2 \\
\sigma_{01}^{\mathrm{QD}}(t) &= \int dR \frac{|\chi_0(R,t)|^2 |\chi_1(R,t)|^2}{|\chi(R,t)|^2}
\end{aligned} 
\qquad (35) $$

where $\{\chi_l(R,t)\}$ are the BH expansion coefficients of the full state of the system and $|\chi(R,t)|^2 = |\chi_0(R,t)|^2 + |\chi_1(R,t)|^2$ is the total nuclear density. The MQC corresponding quantities are trajectory averages

$$ 
\begin{aligned}
P_1^{\mathrm{MQC}}(t) &= \frac{1}{N_{tr}} \sum_{\alpha}^{N_{tr}} \rho_{11}^{(\alpha)}(t) \\
\sigma_{01}^{\mathrm{MQC}}(t) &= \frac{1}{N_{tr}} \sum_{\alpha}^{N_{tr}} |\rho_{01}^{(\alpha)}(t)|^2.
\end{aligned} 
\qquad (36) $$

For the TSH-based calculations the fraction of trajectories on the excited state is also computed

$$ \Pi_1(t) = \frac{N_1(t)}{N_{tr}} \qquad (37) $$

Let us first analyze the QD results. The wavepacket component in the ground state remains trapped in the potential well in the Franck-Condon region, while the excited state wavepacket reaches the NAC region at around 25 fs and splits. This model is weakly nonadiabatic and the main portion of the nuclear wavepacket remains in excited state but there is some transfer to the ground state, creating a small coherence that then decays. After the splitting event at the NAC, the part of the wavepacket that remains in the excited state reflects and re-encounters the NAC region at around 75 fs creating a wavepacket in the ground state that interferes with the trapped wavepacket at 90 fs. The increase in the coherence observed at 90 fs occurs as the excited-state wavepacket traveling back to the left after reflection overlaps with the ground-state wavepacket still localized in the Franck-Condon region. Since the indicator of coherence of Eq. (35) is related to the spatial overlap of the ground-state wavepacket with the excited-state wavepacket, the peak in the coherence at around 100 fs is observed when the overlap of the two wavepackets is maximum. At 125 fs there is an additional peak in the indicator of coherence when the wavepacket component that was created in the ground state in the first splitting event at the NAC, meets the excited state wavepacket travelling to the right after a second reflection at the avoided crossing. There are further interference events later on as wavepacket components on the same surface with different histories meet in nuclear space.
<figure><img src="image_3.png" /><figcaption>FIG. 2: Electronic excited state populations (lower panel) and coherences (upper panel) as a function of time obtained via QD simulations (Exact, in black), Eh starting in a pure ensemble (green), Eh starting in a mixed ensemble (blue), CTMQC starting in a pure ensemble (red) and CTMQC starting in a mixed ensemble (magenta).</figcaption></figure>

A. Discussion on Ehrenfest-based dynamics

Figure 2 shows the excited state populations (lower panel) and coherences (upper panel) as a function of time. QD results are used to benchmark Eh and CTMQC results starting in a pure ensemble or in a mixed ensemble. We observe that starting from a pure ensemble with Eh dynamics yields negligible net population transfer and the electronic wavefunction retains its coherence throughout the dynamics. On the other hand, CTMQC captures the transfer events and decoherence behavior, although all methods miss the recoherence at 90 fs, resulting from the interferences.

In the trajectory-based methods, the coherence is defined as a trajectory average of the individual coherences for each tra-
jectory. In the case of Eh dynamics, this recoherence cannot be captured due to the lack of a quantum-momentum-induced
mechanism that yields an *R*-dependent effect in the electronic populations and coherences (see Eq. (36)). Specifically, in Eh,
this can only be achieved when the trajectories encounter a region of non-zero NAC. In the case of CTMQC, instead, a change
in populations or coherences in *R* does exist through the quantum-momentum term, but it only acts for the trajectories associ-
ated to an instantaneous time-dependent electronic wavefunction that is a superposition of electronic eigenstates (as evidenced
from Eq. (38) below). For the trajectories associated to the trapped ground-state wavepacket in the Franck-Condon region, the
<figure><img src="image_4.png" /><figcaption>FIG. 3: Upper panel: time snapshots of the exact density (orange), BO PESs (red lines), NAC (magenta), TDPES (black) and distribution of CTMQC trajectories starting in a pure ensemble (blue dots) and Eh trajectories starting in a pure ensemble (green dots). Lower panel: Excited state populations as function of R or $R^{(\alpha)}(t)$ for QD (black line), CTMQC starting in a pure ensemble (blue dots) and Eh starting in a pure ensemble (green dots).</figcaption></figure>

electronic coefficients are all “collapsed” to the ground electronic state, i.e., the electronic coefficients of the ground state are one while the coefficients of the excited states are zero, making this quantum-momentum term inactive; for the majorities of the trajectories describing the excited-state wavepacket travelling back to the left, the electronic coefficients are collapsed as well, but to the excited state, making the quantum-momentum term inactive for a large number of trajectories. Thus, very few
<figure><img src="image_5.png" /><figcaption>FIG. 4: Time snapshots of the exact density (orange), density reconstructed from the distribution of trajectories starting in a pure ensemble with CTMQC (blue), and Eh (green), TDPES (black line), and distribution of trajectories of CTMQC (blue dots), and Eh (green dots).</figcaption></figure>

trajectories contribute to the net rise in coherence starting at around 90 fs, yielding an insignificant overall rise. Note that this peak of coherence at 100 fs is different from the other peaks, which are captured by CTMQC (at 20, 70, 120 and 180 fs), because these events occur as a consequence of the population transfer in a NAC region.

In particular, CTMQC starting in the pure ensemble gives the closest agreement to the exact populations overall but captures the dynamics in the transition regions a bit less accurately than the mixed ensemble. The mixed ensembles initially incorrectly begin with zero coherence, as expected from the earlier discussions; they would consequently fail to capture observables such as the electronic current-density, regardless of the propagation method.

We now take a closer look at the pure ensemble case. Figure 3 shows time snapshots of the quantum nuclear density and distribution of positions and excited state populations of the CTMQC and Eh classical-like trajectories at early times starting in a pure ensemble together with the BO PESs, NAC and TDPES. Let us first analyze the Eh results. The trajectories propagate on a mean-field surface and their velocity is determined by its slope and their initial momenta. Since the ensemble of trajectories are away from the interaction region each trajectory experiences zero population transfer in the first 14 fs and thus Eh cannot capture the first decoherence event and splitting of the nuclear density as the excited state wavepacket moves away from the trapped ground wavepacket. Instead, CTMQC corrects both of these problems, thanks to the key role played by the quantum momentum. For this two-electronic-state one-dimensional system, Eq. (26) and Eq. (27) reduce to
$$ \dot{\rho}_{CT,11}(t) = \frac{2}{N_{tr}} \sum_{\alpha} \frac{Q^{(\alpha)}(t)}{\hbar M} \rho_{00}^{(\alpha)}(t) \rho_{11}^{(\alpha)}(t) \Delta f_{10}^{(\alpha)}(t) \quad (38) $$

$$ F_{CT}(t) = \frac{2}{N_{tr}} \sum_{\alpha} \frac{Q^{(\alpha)}(t)}{\hbar M} \rho_{00}^{(\alpha)}(t) \rho_{11}^{(\alpha)}(t) \Delta f_{10}^{(\alpha)2}(t) \quad (39) $$

Although, as evident from Figure 2, there is no net population transfer in the first 25 fs where the trajectories are in a region of zero NAC, each individual trajectory experiences quantum-momentum-induced population transfer responsible for capturing decoherence: the lower panel of Figure 3 shows the position-resolved populations of the excited state at each time curving away from the initial uniform value, with the leading edge of the distribution transferring to the upper surface and the trailing edge transferring to the lower, while the average of the populations over the distribution remains constant. This can be explained by the following. The quantum momentum has a linear shape (by construction), positive on the leading edge of the trajectory distribution and and negative on the trailing edge. Since the trajectories are close enough to each other in $R$ so that the accumulated force difference has the same sign, the coupled-trajectory term (Eq. (38)) changes the populations of the trajectories on the leading and trailing edge in opposite ways leading to branching of the electronic coefficients achieving decoherence. The net population change averaged over all trajectories yet remains zero by design (see comment below Eq. (26)). On the other hand, the sign of the force acting on the trajectories (Eq. (39)) depends solely on the quantum momentum, and as we can observe the force acts on opposite ways on the leading and trailing trajectories splitting the trajectory distribution. This early-time behavior dramatically exemplifies the crucial role of quantum-momentum-driven transitions: in regions of negligible NAC, this mechanism is essential for wavepacket splitting to be even qualitatively captured when the system is in a pure ensemble that involves more than one surface/state.

Figure 4 shows additional time snapshots over the duration of the dynamics, now showing also the density reconstructed from the distribution of classical-like trajectories. We observe at 19 fs that the energy of the CTMQC trajectories follows the shape of the exact TDPES with some trajectories corresponding to the trapped wavepacket in the ground state remaining at 2 bohr, while the others, corresponding to the vertically excited initial wavepacket, evolve to the right, but the splitting of the density is not quite as clean as in the QD case. For Eh, on the other hand, all the trajectories evolve together and slowly as a consequence of the propagation in the mean-field surface. After the first splitting event, we observe at 53 fs the CTMQC trajectories have three parts corresponding to the three parts of the true nuclear density; the left part trapped on the lower surface, the middle part transmitted past the NAC region on the upper surface (largely), and a very small part to the right (around 9 bohr) which CTMQC incorrectly associates more closely to the upper surface instead of the lower surface that the exact density lives on in this region.

In general, we observe that the distributions of CTMQC trajectories at long time are particularly broad, more than the quantum density. This behavior can be very likely attributed to numerical instabilities produced in the classical nuclear force that accumulate over time as a consequence of “noisy” time-dependent scalar and vector potentials in the presence of quantum interferences
<figure><img src="image_6.png" /><figcaption>FIG. 5: Time snapshots of the exact density (orange), density reconstructed from the distribution of trajectories starting in a mixed ensemble with CTMQC (blue), and Eh (green), TDPES (black line), and distribution of trajectories of CTMQC (blue dots), and Eh (green dots).</figcaption></figure>

and revivals of coherence. The CTTSH method of Ref. [49] has been mainly introduced to circumvent these problems since the forces used for the classical propagation are purely adiabatic in this case.

We now turn to the results obtained with Eh and CTMQC starting from the mixed ensemble. As noted earlier (Fig. 2), both Eh and CTMQC capture the transfer events, with Eh deviating from the QD results at around 120 fs. In terms of coherences we note that while the initial mixed ensemble cannot capture the initial coherence behavior as noted earlier, CTMQC dynamics redeems it after the initial decoherence event in the exact system, where it does almost as well as the pure ensemble, describing the periodic recoherence/decoherence well but again missing the large recoherence resulting from the interference event. Eh predicts the small recoherence occurring at 25 fs, but remains overcoherent in the remaining non-adiabatic events. Figure 5 shows time snapshots of the quantum nuclear density and distribution of classical-like trajectories together with the reconstructed density form the trajectories, the BO PESs and the TDPES. Since 20% of the trajectories in both methods start in the excited state surface and 80% on the ground surface, experiencing the different respective forces, we observe at 19 fs that the first splitting of the nuclear density is well captured in both methods. The snapshot at 53 fs indicates that both methods miss the part of the density at 8 bohr corresponding to the wavepacket created on the ground surface in the first interaction region, although few CTMQC trajectories get pushed to the ground surface. In the case of CTMQC this might be a consequence of the imposition that the quantum momentum have the same slope for all trajectories. At 145 fs we observe that CTMQC indeed predicts some density at −1 bohr while Eh completely misses it.
## B. Discussion on surface-hopping-based dynamics

Finally, we turn to the TSH-based methods. In this case, we will analyze only the excited-state populations and the coherences. Figure 6 shows the electronic populations and fraction of trajectories on the excited state (middle and lower panels) and coherences (upper panel) as a function of time, once again benchmarking the trajectory-based methods TSH, TSHEDC and CTTSH against QD results. The dynamics is initialized in the pure ensemble and in a mixed ensemble.

<figure><img src="image_7.png" /><figcaption>FIG. 6: Fraction of trajectories on the excited state (lower panel), excited state electronic populations (middle panel) and coherences (upper panel) as a function of time obtained via QD simulations (Exact, in black), TSH starting in a pure ensemble (orange) and in a mixed ensemble (violet), TSHEDC starting in a pure ensemble (green) and in a mixed ensemble (cyan), and CTTSH starting in a pure ensemble (blue) and in a mixed ensemble (red). Note that the green and cyan lines are superimposed in the lower panel.</figcaption></figure>

Let us first focus on the pure ensemble results. We observe that while TSH and TSHEDC perform very poorly (as expected),
CTTSH shows very good agreement with the QD reference. In the case of TSH the overcoherence of the trajectories yields wrong dynamics when the first interaction region is encountered, and thereon; in the case of CTTSH, due to the quantum-momentum-induced population transfers for each trajectory (as we saw for CTMQC), decoherence is recovered (see upper panel of Fig. 6) and capturing this early time behavior is key to yield accurate dynamics as the trajectories propagate through the successive interaction regions. On the other hand, the poor and notably unusual performance of TSHEDC starting from the pure ensemble can be explained by the following. In TSHEDC at every time-step the coefficients of the non-active states are exponentially dampened with the active-state coefficient being renormalized as

$$C_i^{(\alpha)\prime}(t) = e^{-\frac{\Delta\Delta\epsilon^{(\alpha)}}{\hbar(1+C/\epsilon_{kin}^{(\alpha)})}} C_i^{(\alpha)}(t) \quad i \neq \text{active} \qquad (40)$$

$$C_{active}^{(\alpha)\prime}(t) = C_{active}^{(\alpha)}(t) \sqrt{\frac{\sum_{j \neq active} |C_j^{(\alpha)\prime}(t)|^2}{|C_{active}^{(\alpha)}|^2}} \qquad (41)$$

in order to impose internal consistency, namely to enforce that the quantum population of Eq. (36) equals at all times the fraction of trajectories from Eq. (37). Here, $\Delta t$ is the electronic time-step used in the simulation, $E_{kin}^{(\alpha)}$ the kinetic energy of the $\alpha$-th trajectory, $\Delta\epsilon^{(\alpha)} = \epsilon_i^{(\alpha)} - \epsilon_{active}^{(\alpha)}$ is the energy difference between the active state PES and the $i$-th PES at $\underline{R}(t) = \underline{R}^{(\alpha)}(t)$ and $C$ is a constant that we set to 0.1 Ha [48, 49, 94].

In the pure ensemble, all the trajectories start in the same superposition state, with a distribution of active states that matches the trajectory-averaged electronic populations at the initial time, i.e., 80% of the trajectories (denoted $\alpha$-trajectories in the following) will have the ground state as the active state $F^{(\alpha)}(0) = -\nabla\epsilon_0^{(\alpha)}$ while the remaining 20% (denoted $\beta$-trajectories in the following) will start in the excited adiabatic state $F^{(\beta)}(0) = -\nabla\epsilon_1^{(\beta)}$. For the $\alpha$-trajectories the excited state is the inactive state and so the corresponding electronic coefficient will be exponentially dampened, whereas for the $\beta$-trajectories it is the inactive ground-state coefficient that will be dampened. However, the rate of decay will be faster for trajectories with higher kinetic energies (see Eq. (40)), i.e., for trajectories propagating on the excited state ($\beta$-trajectories) due to the larger initial force they experience. Hence, since for 20% of the trajectories in the ensemble the excited-state population increases at a higher exponential rate than the excited state population decay of the remaining 80%, we see a rapid increase in the excited state population between 0 and 20 fs. It is worth stressing that, in general, in standard TSH methods the choice of initiating the dynamics in a pure ensemble is not preferred. Specifically, with this choice of initial conditions, internal consistency is violated at the single trajectory level from the very start: the motion of each trajectory is determined by a single active surface, and yet it is associated to more than one non-zero state coefficient/population and with independent trajectories it is difficult to enforce the internal consistency over the average, as this example dramatically demonstrates.

Now let us turn to the mixed ensemble results. In terms of coherences, all methods miss the early coherence behavior since the ensemble already starts decohered; again, they will suffer from poor initial observables as discussed earlier. Also, they miss the strong recoherence at 90 fs result of the first quantum interference event. TSH remains overcoherent after the first crossing through the NAC region, while CTTSH and TSHEDC predict the decoherence and recoherence events with TSHEDC slightly underestimating the recoherences. In terms of populations, we observe that the fraction of trajectories in the excited
state obtained with TSH underestimates the net population transfer to the ground state after the second interaction region. This behavior is similar for the electronic populations of CTTSH. TSHEDC on the other hand starts deviating in the first NAC region already significantly underestimating the transfer to the ground state.

V. CONCLUSIONS AND OUTLOOK

In this work we have discussed the problem of the initialization of the electronic state in non-adiabatic dynamics when simulations are performed using trajectory-based methods, such as in Ehrenfest, surface hopping and coupled-trajectory schemes derived from the exact factorization. Such methods describe the electron-nuclear dynamics beyond the Born-Oppenheimer approximation using classical-like nuclear trajectories coupled to quantum-mechanical electrons. Often in the literature the effect of different sampling schemes to select the initial conditions for the trajectories have been discussed, but here we focused on the other aspect of the problem, namely the initialization of the electronic dynamics.

We have introduced the concepts of pure ensemble and mixed ensemble for the initialization of the electronic wavefunc-
tion in the context of a trajectory-based approach. These concepts naturally arise when the “real” (exact, in our study) initial
state of the system is a coherent superposition of electronic states, a situation that can be created experimentally upon photo-
excitation/ionization of a molecular system with an attosecond pulse.

Based on these concepts, we have performed calculations for a model system where the exact molecular wavefunction begins
in a coherent superposition of adiabatic electronic states. We compare both standard and exact-factorization-based trajectory
methods starting the dynamics in a pure or in a mixed ensemble, against the exact quantum dynamics. While the pure ensemble
choice for Ehrenfest dynamics is the natural one, it performs very poorly from the start even if one is just interested in the
electronic populations, due to the mean-field nature of the force on the nuclei and the resulting erroneous feedback on the
electronic equations. As expected, standard surface-hopping methods perform better with the mixed ensemble, since internal
consistency is valid since the beginning of the dynamics at the single trajectory level, however, it is pointed out that with this
choice some observables that are directly related to coherences (e.g., electronic current-density) start off on the wrong foot. The
coupled-trajectory methods seem to perform quite well with both choices of initial conditions, although when starting in the
mixed ensemble they suffer from the same problem just mentioned. The quantum-momentum-induced transitions are crucial in
the initial stages of dynamics starting in a pure state to yield the correct forces for the trajectories correlated with each of the
electronic states, and throughout the dynamics both the full CTMQC as well as the CTTSH scheme do the best job for electronic
and nuclear observables overall out of the methods studied.

The present work studied only a one-dimensional two-electronic-state system, while the differences between the pure and the
mixed ensembles of initial electronic state for more degrees of freedom are discussed only theoretically. It would be interesting
to study the performance of the various methods for these different choices of the ensembles as a function of both the number of
electronic states and the number of vibrational degrees of freedom. For example, in the limit of a manifold of dense electronic
states, the choice of initiating pure or mixed ensembles of electronic states is expected to become less important, thus, the
difference between the different methods likely decreases. Likewise, in the limit of a large number of vibrational degrees
of freedom, we might expect that any initial coherence in the electronic state decays rapidly. Also in this case, it would be
interesting to analyze the role of the quantum-momentum-induced transitions in bringing about the decoherence and how the
different methods perform when starting the dynamics with the pure and mixed ensembles. Performing simulations in the
condensed phase is, however, very challenging for all the trajectory methods discussed here, except for Ehrenfest because they
require resolving electronic coefficients on the near-continuum of electronic states. In this regard, we note that a recent recasting
of the exact-factorization-based methods that operates with real-time TDDFT may be particularly useful [95].

VI. ACKNOWLEDGEMENTS

This research was supported by the ANR Q-DeLight project, Grant No. ANR-20-CE29-0014 of the French Agence Nationale de la Recherche, the National Science Foundation Award CHE-2154829, the Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences under Award No. DE-SC0020044, the Computational Chemistry Center: Chemistry in Solution and at Interfaces funded by the U.S. Department of Energy, Office of Science Basic Energy Sciences, under Award DE-SC0019394 and the Chateaubriand Fellowship from the Mission Scientifique et Technologique of the Embassy of France in the United States.

VII. DATA AVAILABILITY

The data presented in this study are available on reasonable request from the corresponding author.

[1] A. H. Zewail, Angew. Chem., Int. Ed. **39**, 2586 (2000).

[2] A. H. Zewail, J. Phys. Chem. A **104**, 5660 (2000).

[3] P. M. Paul, E. S. Toma, P. Breger, G. Mullot, F. Augé, P. Balcou, H. G. Muller, and P. Agostini, Science **292**, 1689 (2001).

[4] M. Hentschel, R. Kienberger, C. Spielmann, G. A. Reider, N. Milosevic, T. Brabec, P. Corkum, U. Heinzmann, M. Drescher, and F. Krausz, Nature **414**, 509 (2001).

[5] G. Sansone, F. Kelkensberg, J. F. Pérez-Torres, F. Morales, M. F. Kling, W. Siu, O. Ghafur, P. Johnsson, M. Swoboda, E. Benedetti, et al., Nature **465**, 763 (2010).

[6] F. Calegari, D. Ayuso, A. Trabattoni, L. Belshaw, S. D. Camillis, S. Anumula, F. Frassetto, L. Poletto, A. Palacios, P. Decleva, et al., Science **346**, 336 (2014).

[7] F. Krausz and M. Ivanov, Rev. Mod. Phys. **81**, 163 (2009).

[8] F. Lépine, M. Y. Ivanov, and M. J. J. Vrakking, Nat. Photonics **8**, 195 (2014).

[9] M. Nisoli, P. Decleva, F. Calegari, A. Palacios, and F. Martín, Chem. Rev. **117**, 10760 (2017).

[10] A. Palacios and F. Martín, Wiley Interdiscip. Rev. Comput. Mol. Sci. **10**, 1430 (2020).

[11] L. Cederbaum and J. Zobeley, Chem. Phys. Lett. **307**, 205 (1999).
[12] F. Calegari, D. Ayuso, A. Trabattoni, L. Belshaw, S. D. Camillis, S. Anumula, F. Frassetto, L. Poletto, A. Palacios, P. Decleva, et al., *Science* **346**, 336 (2014).

[13] A. I. Kuleff¹ and L. Cederbaum, *J. Phys. B: At. Mol. Opt. Phys.* **47**, 124002 (2014).

[14] M. Vacher, D. Mendive-Tapia, M. J. Bearpark, and M. A. Robb, *J. Chem. Phys.* **142**, 094105 (2015).

[15] M. Vacher, L. Steinberg, A. J. Jenkins, M. J. Bearpark, and M. A. Robb, *Phys. Rev. A* **92**, 040502 (2015).

[16] M. Vacher, D. Mendive-Tapia, M. J. Bearpark, and M. A. Robb, *J. Chem. Phys.* **142**, 094105 (2015).

[17] J. Ajay, J. Šmydke, F. Remacle, and R. D. Levine, *J. Phys. Chem. A* **120**, 3335 (2016).

[18] C. Arnold, O. Vendrell, and R. Santra, *Phys. Rev. A* **95**, 033425 (2017).

[19] A. Scheidegger, J. Vaníček, and N. V. Golubev, *J. Chem. Phys.* **156**, 034104 (2022).

[20] A. Bruner, S. Hernandez, F. Mauger, P. M. Abanador, D. J. LaMaster, M. B. Gaarde, K. J. Schafer, and K. Lopata, *J. Phys. Chem. Lett.* **8**, 3991 (2017).

[21] L. Lacombe and N. T. Maitra, *J. Phys. Chem. Lett.* **12**, 8554 (2021).

[22] L. Lacombe and N. T. Maitra, *Npj Comput. Mater.* **9**, 124 (2023).

[23] I. Sanchez and F. Martin, *J. Phys. B: At. Mol. Opt. Phys.* **30**, 679 (1997).

[24] C. Marante, M. Klinker, I. Corral, J. González-Vázquez, L. Argenti, and F. Martín, *J. Chem. Theory Comput.* **13**, 499 (2017).

[25] M. Ruberti, P. Decleva, and V. Averbukh, *J. Chem. Theory Comput.* **14**, 4991 (2018).

[26] B. T. Pickup, *Chem. Phys.* **19**, 193 (1977).

[27] L. S. Cederbaum, W. Domcke, J. Schirmer, and W. V. Niessen, *Adv. Chem. Phys.* pp. 115–159 (1986).

[28] J. Suchan, D. Hollas, B. F. E. Curchod, and P. Slavíček, *Faraday Discuss.* **212**, 307 (2018).

[29] D. V. Makhov and D. V. Shalashilin, *Chem. Phys.* **515**, 46 (2018).

[30] J. Breidbach and L. S. Cederbaum, *Phys. Rev. Lett.* **94**, 033901 (2005).

[31] M. Persico and G. Granucci, *Photochemistry. A Modern Theoretical Perspective* (Springer, 2018).

[32] M. Born and R. Oppenheimer, *Ann. Phys.* **389**, 457 (1927).

[33] I. C. D. M. D. Jacquermin and M. Vacher, *J. Phys. Chem. Lett.* **12**, 8404 (2021).

[34] M. Vacher, D. Mendive-Tapia, M. J. Bearpark, and M. A. Robb, *J. Chem. Phys.* **142**, 094105 (2015).

[35] A. Jenkins, M. Vacher, M. J. Bearpark, and M. A. Robb, *J. Chem. Phys.* **144**, 104110 (2016).

[36] M. Vacher, M. J. Bearpark, M. A. Robb, and J. P. Malhado, *Phys. Rev. Lett.* **118**, 083001 (2017).

[37] M. Barbatti and K. Sen, *Int. J. Quantum Chem.* **116**, 762 (2016).

[38] M. Vacher, D. Mendive-Tapia, M. J. Bearpark, and M. A. Robb, *Theor. Chem. Acc.* **133**, 1505 (2014).

[39] M. Vacher, F. E. A. Albertani, A. J. Jenkins, I. Polyak, M. J. Bearpark, and M. A. Robb, *Faraday Discuss.* **194**, 95 (2016).

[40] M. Lara-Astiaso, A. Palacios, P. Decleva, I. Tavernelli, and F. Martín, *Chem. Phys. Lett.* **683**, 357 (2017).

[41] D. Danilov, A. J. Jenkins, M. J. Bearpark, G. A. Worth, and M. A. Robb, *J. Phys. Chem. Lett.* **14**, 6127 (2023).

[42] M. Lara-Astiaso, D. Ayuso, I. Tavernelli, P. Decleva, A. Palacios, and F. Martín, *Faraday Discuss.* **194**, 41 (2016).

[43] D. Ayuso, A. Palacios, P. Decleva, and F. Martín, *Phys. Chem. Chem. Phys.* **19**, 19767 (2017).

[44] P. López-Tarifa, D. Grzegorz, Piekarski, E. Rossich, M.-A. H. du Penhoat, R. Vuilleumier, M.-P. Gaigeot, I. Tavernelli, M.-F. Politis, Y. Wang, et al., *J. Phys. Conf. Ser.* **488**, 012037 (2014).

[45] S. K. Min, F. Agostini, and E. K. U. Gross, *Phys. Rev. Lett.* **115**, 073001 (2015).

[46] F. Agostini, S. K. Min, A. Abedi, and E. K. U. Gross, *J. Chem. Theory Comput.* **12**, 2127 (2016).
[47] S. K. Min, F. Agostini, I. Tavernelli, and E. K. U. Gross, J. Phys. Chem. Lett. **8**, 3048 (2017).

[48] G. Granucci and M. Persico, J. Chem. Phys. **126**, 134114 (2007).

[49] C. Pieroni and F. Agostini, J. Chem. Theory Comput. **17**, 5969 (2021).

[50] F. Talotta, D. Lauvergnat, and F. Agostini, J. Chem. Phys. **156**, 184104 (2022).

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

[52] F. Agostini and E. K. U. Gross, Eur. Phys. J. B **94**, 179 (2021).

[53] G. Hunter, Int. J. Quantum Chem. **9**, 237 (1975).

[54] G. Hunter, Int. J. Quantum Chem. **9**, 311 (1975).

[55] G. Hunter, Int. J. Quantum Chem. **9**, 133 (1980).

[56] G. Hunter, Int. J. Quantum Chem. **19**, 755 (1981).

[57] G. Hunter and C. C. Tai, Int. J. Quantum Chem. **21**, 1041 (1982).

[58] N. I. Gidopoulos and E. K. U. Gross, Philos. Trans. Royal Soc. A. **372** (2014), ISSN 1364-503X.

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

[60] E. Villaseco Arribas, F. Agostini, and N. T. Maitra, Molecules **27**, 13 (2022).

[61] A. Donoso, D. Kohen, and C. C. Martens, J. Chem. Phys. **112**, 7345 (2000).

[62] C. C. Martens and J.-Y. Fang, J. Chem. Phys. **106**, 4918 (1997).

[63] A. Donoso and C. C. Martens, J. Phys. Chem. A **102**, 4291 (1998).

[64] J. C. Tully, J. Chem. Phys. **93**, 1061 (1990).

[65] J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi, Ann. Rev. Phys. Chem. **67**, 387 (2016).

[66] M. Barbatti, G. Granucci, M. Persico, M. Ruckenbauer, M. Vazdar, M. Eckert-Maksić, and H. Lischka, J. Photochem. Photobiol. A **190**, 228 (2007).

[67] R. Mitrić, J. Petersen, and V. Bonačić Koutecký, Phys. Rev. A **79**, 053416 (2009).

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

[69] J. J. Bajo, J. González-Vázquez, I. Sola, J. Santamaría, M. Richer, P. Marquetand, and L. González, J. Phys. Chem. A **116**, 2800 (2012).

[70] T. Fiedlschuster, J. Handt, and R. Schmidt, Phys. Rev. A **93**, 053409 (2016).

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

[72] M. Heindl and L. Gonzalez, J. Chem. Phys. **154**, 144102 (2021).

[73] A. D. McLachlan, Mol. Phys. **8**, 39 (1964).

[74] L. M. Ibele, B. F. E. Curchod, and F. Agostini, J. Phys. Chem. A **126**, 1263 (2022).

[75] F. Agostini and B. F. E. Curchod, Wiley Interdiscip. Rev. Comput. Mol. Sci. **0**, e1417 (2019).

[76] G. H. Gossel, L. Lacombe, and N. T. Maitra, J. Chem. Phys. **150**, 154112 (2019).

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

[78] F. G. Eich and F. Agostini, J. Chem. Phys. **145**, 054110 (2016).

[79] P. Vindel-Zandbergen, S. Matsika, and N. T. Maitra, J. Phys. Chem. Lett. **13**, 1785 (2022).

[80] E. Villaseco Arribas, P. Vindel-Zandbergen, S. Roy, and N. T. Maitra, Phys. Chem. Chem. Phys. **25**, 26380 (2023).

[81] E. Villaseco Arribas and N. T. Maitra, J. Chem. Phys. **158** (2023).

[82] E. Villaseco Arribas, L. M. Ibele, D. Lauvergnat, N. T. Maitra, and F. Agostini, J. Chem. Theory Comput. (2023).

[83] C. Pieroni, E. Sangiogo Gil, L. M. Ibele, M. Persico, G. Granucci, and F. Agostini, ChemRxiv (2023).
[84] J.-K. Ha, I. S. Lee, and S. K. Min, J. Phys. Chem. Lett. **9**, 1097 (2018).

[85] P. Vindel-Zandbergen, L. M. Ibele, J.-K. Ha, S. K. Min, B. F. E. Curchod, and N. T. Maitra, J. Chem. Theory Comput. **17**, 3852 (2021).

[86] J.-K. Ha and S. K. Min, J. Chem. Phys. **156**, 174109 (2022).

[87] A. Pereira, J. Knapik, A. Chen, D. Lauvergnat, and F. Agostini, Eur. Phys. J. Spec. Top. **232**, 1917 (2023).

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

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

[90] B. F. E. Curchod, F. Agostini, and E. K. U. Gross, J. Chem. Phys. **145**, 034103 (2016).

[91] M. Feit, J. Fleck, and A. Steiger, J. Comp. Phys. **47**, 412 (1982).

[92] A. D. Bandrauk and H. Shen, J. Chem. Phys. **99**, 1185 (1993).

[93] F. Agostini, E. Marsili, F. Talotta, and E. Villaseco Arribas, G-CTMQC (Last accessed Oct 2023), https://gitlab.com/agostini.work/g-ctmqc.

[94] C. Zhu, S. Nangia, A. W. Jasper, and D. G. Truhlar, J. Chem. Phys. **121**, 7658 (2004).

[95] D. Han, J. K. Ha, and S. Min, J. Chem. Theory Comput. **19**, 2186 (2023).