# Exact Factorization Adventures: A Promising Approach for Non-Bound States

Evaristo Villaseco Arribas, Federica Agostini, Neepa T. Maitra

► To cite this version:

Evaristo Villaseco Arribas, Federica Agostini, Neepa T. Maitra. Exact Factorization Adven-
tures: A Promising Approach for Non-Bound States. Molecules, 2022, 27 (13), pp.4002.
10.3390/molecules27134002 . hal-04295721

HAL Id: hal-04295721

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

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
Article

# Exact Factorization Adventures: A Promising Approach for Non-Bound States

Evaristo Villaseco Arribas<fnref n="1" />, Federica Agostini<fnref n="2" />* and Neepa T. Maitra <fnref n="1" />,*

**Abstract:** Modeling the dynamics of non-bound states in molecules requires an accurate description of how electronic motion affects nuclear motion and vice-versa. The exact factorization (XF) approach offers a unique perspective, in that it provides potentials that act on the nuclear subsystem or electronic subsystem, which contain the effects of the coupling to the other subsystem in an exact way. We briefly review the various applications of the XF idea in different realms, and how features of these potentials aid in the interpretation of two different laser-driven dissociation mechanisms. We present a detailed study of the different ways the coupling terms in recently-developed XF-based mixed quantum-classical approximations are evaluated, where either truly coupled trajectories, or auxiliary trajectories that mimic the coupling are used, and discuss their effect in both a surface-hopping framework as well as the rigorously-derived coupled-trajectory mixed quantum-classical approach.

**Keywords:** non-adiabatic dynamics; light-matter interactions; molecular dynamics algorithms

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

check for
updates

**Citation:** Villaseco Arribas, E.; Agostini, F.; Maitra, N.T. Exact Factorization Adventures: A Promising Approach for Non-Bound States. *Molecules* **2022**, *27*, 4002.
https://doi.org/10.3390/molecules27134002

Academic Editor: Sergey I. Bokarev

Received: 13 May 2022

Accepted: 9 June 2022

Published: 22 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2022 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

A challenge in the theoretical simulation of molecules driven out of their ground-state by a strong field is how to account for, and interpret, the coupling between the electrons and ions. If a molecule was only gently perturbed at frequencies below electronic excitations, the Born-Oppenheimer (BO) picture can still be useful: the disparity in the time-scales of nuclear and electronic motion allows for the nuclei to evolve on a single BO surface for much of the time, and regions where this breaks down tend to be somewhat localized in space and involving the interaction of only a small number of BO states. In those regions, a fully quantum treatment of the molecule based on the BO picture would involve propagating nuclear wavefunctions projected on the different surfaces, including coupling terms between them. If the nuclear degrees of freedom were to be treated classically, the effective force driving their motion would deviate from that of a single surface, or even a mean surface: different parts of the underlying nuclear wavepacket would be associated with different electronic surfaces, even within the same spatial region. The problems with capturing this electron-nuclear correlation in such mixed quantum-classical methods are well-known [1,2]. But with strong fields, or with high-frequency fields, the problem is even more challenging: the nuclear wavepacket straddles many BO states, a continuum beyond the ionization threshold. Field-dressed surfaces such as quasi-static [3–5] or Floquet [6–8] can be instructive in this regard; although both serve equally well as a basis for the nuclear dynamics, as does the BO basis, approximations inherent in mixed quantum-classical approaches may perform better in one or the other, for physical reasons or computational reasons, but for general field parameters, the jury is out. Further, the observables of interest depend on the application: for example, we may want to correlate electron ionization and nuclear dissociation channels to get a full understanding of the dynamics [9], or, we may be primarily interested in the behavior of the ionizing electron itself, with a time-dependent potential driving the electron helping to guide our understanding. As intuitive

<footnote n="1">¹ Department of Physics, Rutgers University, Newark, NJ 07102, USA; ev259@newark.rutgers.edu</footnote>

<footnote n="2">² Institut de Chimie Physique UMR8000, Université Paris-Saclay, CNRS, 91405 Orsay, France; federica.agostini@universite-paris-saclay.fr</footnote>

<footnote>* Correspondence: neepa.maitra@rutgers.edu</footnote>
and instrumental the BO picture has been, it does not provide such a potential, nor is it convenient in strong fields when many surfaces get involved.

An alternative, but equally fundamental, picture is provided by the exact factorization (XF) approach [10–18]. XF re-casts the dynamics of interacting quantum subsystems in a way that provides unique potentials that drive the motion of each subsystem and account for complete coupling to others. Applied to the electron-nuclear case, the exact molecular wavefunction $\Psi(\underline{r}, \underline{R}, t)$ is factored into a single correlated product of a marginal factor and a conditional factor. Choosing the marginal to be a function of $\underline{R}$, one obtains a time-dependent Schrödinger equation (TDSE) for the nuclear wavefunction, in which the potentials contain the complete effect of coupling to the electrons [16]. Instead, choosing the marginal to be a function of $\underline{r}$ only, a TDSE is obtained for the electronic wavefunction, in which the potentials contain the complete effect of coupling to the nuclei [19]. These potentials drive the nuclear and electronic subsystem respectively, and have been instructive in interpreting correlated electron-ion dynamics, including for dissociation [16] and ionization processes [20] in strong fields. To be a practical approach, approximations must be made, and recently two mixed quantum-classical (MQC) approaches have been derived, in which a key ingredient in the electron-nuclear coupling terms is the so-called nuclear quantum momentum. The quantum momentum measures the spatial variation of the nuclear density during non-adiabatic dynamics, and couples strongly to the electronic evolution giving large corrections to mean-field approximations. Different schemes have proposed different ways to calculate this term, using either the support of coupled trajectories or of auxiliary trajectories to track the delocalization of the nuclear density, but until now a detailed comparison has not been made. Here we compare the use of auxiliary trajectories within an XF-based surface-hopping framework (SHXF) [21,22] to using coupled trajectories within the surface-hopping procedure (CTSH) [23] as well as within the original coupled-trajectory mixed quantum-classical approach (CTMQC) [24–26] that was derived directly from the XF. We further study the impact of a modified definition of the quantum momentum in the coupled-trajectories calculations CTSH and CTMQC, that ensures the physical condition of zero population transfer in regions of no non-adiabatic coupling.

The paper is organized as follows. In Section 2, we review the formalism of the XF approach, and briefly meander down various avenues in which the XF idea has been extended and explored. Section 3 demonstrates its usefulness for the dynamics of non-bound systems, briefly reviewing the results of refs. [16,17] for dissociation and refs. [20,27] for ionization, in a one-dimensional model of $H_2^+$. The most practical impact of the XF approach is the development of rigorous MQC methods, and in Section 4 we outline the CTMQC, CTSH, and SHXF schemes, before delving into the different ways the quantum momentum is computed in different implementations of these schemes in Section 5. We offer some conclusions in Section 6.

## 2. The XF Approach in a Nutshell

Consider the time evolution of a system of coupled quantum subsystems, possibly subject to some external fields:

$$ \hat{H}\Psi(q, \underline{Q}, \dots, t) = i\partial_t\Psi(q, \underline{Q}, \dots, t) \qquad (1) $$

(Atomic units, $e^2 = \hbar = m_e = 1$, are used throughout this article unless otherwise stated). These could be different types of particles, for example with $\underline{q}$ representing electronic coordinates and $\underline{Q}$ nuclear coordinates, and perhaps another set of coordinates for photons. Or, for example, we could have a system of the same type of particle, say electrons, but with $\underline{q}$ being occupation numbers of some set of orbitals, e.g., weakly-correlated and $\underline{Q}$ being occupation numbers of strongly correlated orbitals. The problem is quite general,
and in the XF the full wavefunction is factored into a marginal term and a conditional term [10,15–17]:

$$ \Psi(q, \underline{Q}, t) = \chi(\underline{Q}, t)\Phi_{\underline{Q}}(q, t), \text{ where } \langle\Phi_{\underline{Q}}(t)|\Phi_{\underline{Q}}(t)\rangle_{\underline{q}} = 1 \text{ for all } t, \underline{Q} \quad (2) $$

Here $\langle...|...\rangle_q$ represents an inner product over all $q$ variables only. The factorization Equation (2) is exact and unique up to a $\underline{Q}$- and $t$-dependent phase, where $e^{i\theta(\underline{Q},t)}$ multiplies $\chi$ while $e^{-i\theta(\underline{Q},t)}$ multiplies $\Phi_{\underline{Q}}$, yielding the same product. For Hamiltonians that have the form

$$ \hat{H} = -\frac{1}{2} \sum_{I} \frac{1}{M_I} \nabla_{\underline{Q}_I}^2 - \frac{1}{2} \sum_{i} \frac{1}{m_i} \nabla_{q_i}^2 + V(q, \underline{Q}, t) \quad (3) $$

where $V(q, \underline{Q}, t)$ is a purely multiplicative operator in $q, \underline{Q}$, applying the Dirac-Frenkel variational principle shows that the marginal factor $\chi(\underline{Q}, t)$ satisfies a TDSE with a scalar and vector potential that depend on $\Phi_{\underline{Q}}(q, t)$ while the conditional factor satisfies a more complicated equation, with coupling terms dependent on $\chi(\underline{Q}, t)$ [16,17,28]. In Equation (3) we used the symbols $M_I$ and $m_i$ to indicate the masses—in atomic units—of the two sets of particles with positions $\underline{Q}$ and $q$, respectively. Specifically, for the electron-nuclear problem, with $\underline{Q} = \underline{R}$ the nuclear coordinates, and $q = \underline{r}$ the electronic coordinates, and the partial normalization condition $\int d\underline{r}|\Phi_{\underline{R}}(\underline{r}, t)|^2 = 1$, the wavefunctions $\Phi_{\underline{R}}(\underline{r}, t)$ and $\chi(\underline{R}, t)$ satisfy:

$$ (\hat{H}_{el}(\underline{r}, \underline{R}, t) - \epsilon(\underline{R}, t))\Phi_{\underline{R}}(\underline{r}, t) = i\partial_t \Phi_{\underline{R}}(\underline{r}, t) \quad (4) $$

$$ \left( \sum_{\nu=1}^{N_n} \frac{1}{2M_\nu} (-i\nabla_\nu + \mathbf{A}_\nu(\underline{R}, t))^2 + \hat{V}_n^{\text{ext}}(\underline{R}, t) + \epsilon(\underline{R}, t) \right) \chi(\underline{R}, t) = i\partial_t \chi(\underline{R}, t) \quad (5) $$

with the electronic Hamiltonian $\hat{H}_{el}$ defined as

$$ \hat{H}_{el}(\underline{r}, \underline{R}, t) = \hat{H}_{BO}(\underline{r}, \underline{R}) + \hat{U}_{en}[\Phi_{\underline{R}}, \chi] + \hat{V}_{e}^{\text{ext}}(\underline{r}, t) \quad (6) $$

and with electron-nuclear coupling term

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

Nuclear masses in atomic units are indicated with the symbol $M_\nu$, with $\nu$ labelling the nuclei. Here $\hat{H}_{BO}$ is the usual BO Hamiltonian (the sum of the electron kinetic energy, electron-electron, nuclear-nuclear, and electron-nuclear interaction operators). The $\hat{V}_n^{\text{ext}}(\underline{R}, t)$ or $\hat{V}_e^{\text{ext}}(\underline{r}, t)$ potential is any externally applied potential, such as a laser field, acting on the nuclei or the electrons, respectively. We denote the time-dependent scalar potential appearing in Equations (4) and (5)

$$ \epsilon(\underline{R}, t) = \langle \Phi_{\underline{R}}(t) | \hat{H}_{el}(\underline{R}, t) - i\partial_t | \Phi_{\underline{R}}(t) \rangle_{\underline{r}} \quad (8) $$

as the time-dependent potential energy surface (TDPES). Together with the time-dependent vector potential appearing in Equations (5) and (7)

$$ \mathbf{A}_\nu(\underline{R}, t) = \langle \Phi_{\underline{R}}(t) | -i\nabla_\nu \Phi_{\underline{R}}(t) \rangle_{\underline{r}}, \quad (9) $$

and the electron-nuclear coupling term $\hat{U}_{en}$, these three terms embody the exact electron-nuclear coupling. Notice that the potentials appearing in the nuclear equation, $\epsilon(\underline{R}, t)$ and $\mathbf{A}_\nu(\underline{R}, t)$, involve the electronic wavefunction $\Phi_{\underline{R}}(\underline{r}, t)$, while the electron-nuclear coupling
term appearing in the electronic equation, $\hat{U}_{en}$, involves the nuclear wavefunction $\chi$: thus the evolution of each factor is intimately correlated with the other. Refs. [16,17] showed Equations (4)–(9) are form-invariant under a gauge-like transformation that multiplies $\Phi_R$ by an $\underline{R}, t$-dependent phase $e^{i\theta(\underline{R},t)}$, and $\chi$ by its inverse; the potentials transform correspondingly as $\underline{A} \rightarrow \underline{A} + \nabla\theta(\underline{R},t)$ and $\epsilon(\underline{R},t) \rightarrow \epsilon(\underline{R},t) + \partial_t\theta(\underline{R},t)$. Further, it was shown that we can identify $\chi(\underline{R},t)$ as the nuclear wavefunction, in the sense that its modulus-square gives the $N_n$-body nuclear density, and the exact $N_n$-body nuclear current-density can be extracted in the usual way from its phase together with the vector potential [16].

These equations emphasize the distinction from the far simpler BO approximation, that also involves a single correlated product, but with a time-independent electronic wavefunction and much simpler potential terms. The TDSE-form of the nuclear equation means that the potentials appearing in it can be directly interpreted as driving the nuclear motion, fully incorporating the effect of non-adiabatic and time-dependent coupling to the electronic system.

## 2.1. *XF Adventures in Different Realms*

Although originally presented for coupled electron-nuclear dynamics, the XF approach applies quite generally to interacting quantum subsystems. It has been explored and extended along myriad different paths over the past decade. We summarize some of these here.

*Exact potentials driving the nuclear motion* The first applications [16,17] showed the usefulness of the XF approach in providing potentials that directly correlate with nuclear motion in laser-driven dynamics, aiding in the interpretation of the physical processes. We defer further discussion of this to Section 3. The exact surfaces have been instructive for analysis of basis-dependence of MQC surface-hopping schemes for laser-driven bond-softening and bond-hardening processes [29], and control of laser-induced electron localization [30]. For field-free dynamics following a photo-excitation, the exact potential was found to track the BO surface initially occupied, develop a diabatic-like character as the wavepacket passed through a non-adiabatic coupling region and then form steps that bridge the different BO surfaces associated to the local nuclear wavepacket [31–34]. More precisely, the steps in the gauge-independent term bridge the BO surfaces, while there is a partially compensating step in the gauge-dependent term. The latter, when gauge-transformed to a vector potential, could be interpreted as a momentum-rescaling. Interestingly, in situations where nuclear wavepackets associated with different electronic states meet, the exact TDPES develops features such that independent classical nuclear trajectories evolving on it mimic the true nuclear distribution of interfering nuclear wavepackets [35]. In all these studies, the exact potentials were obtained from an inversion of the numerical solution of the full molecular TDSE.

*Mixed quantum-classical and quantum trajectory methods* The most practical impact that the XF has had so far has been in the development of rigorous MQC methods [24–26,36,37]. We defer a discussion of this to Section 4. While a purely quantum-classical scheme is not able to capture correctly nuclear quantum effects, such as tunnelling through a potential barrier of light particles [38–45] or scattering processes at low temperatures [46], the use of *quantum trajectories* within an XF framework [47–49] holds the promise to overcome this limitation while maintaining a general non-adiabatic viewpoint. Proof-of-principle studies in this direction have been proposed employing a Bohmian perspective [50–53] on nuclear dynamics with time-dependent XF potentials, and developments are currently ongoing especially aiming to address well-known numerical instabilities [54–59] inherent a quantum-trajectory method. Similarly, using a hydrodynamic description of the XF equations, a trajectory-based representation for the electrons had been introduced in ref. [60] coupled to quantum nuclei.

*Inclusion of spin-orbit coupling* Non-adiabatic processes mediated by spin-orbit coupling (SOC), i.e., intersystem crossings, play a key role in unraveling many excited-state processes such as photo-induced dynamics in organometallic complexes [61–66] and even collision reactions involving systems of light elements [67–69]. Refs. [70,71] extended the XF approach to derive a MQC approach that addresses both spin-allowed electronic transitions
mediated by non-adiabatic coupling (NAC) as well as spin-forbidden transitions mediated by the SOC in the same theoretical construction.

**Geometric phase** The concept of a geometric phase arises naturally in XF and can be defined as the circulation integral of the vector potential along a closed path, both in the stationary formulation of XF [72–74] and in its time-dependent version [34,75–77]. It is worth underlining, though, that a non-vanishing geometric phase (not topological, as it is not quantized) in the framework of XF does not require invoking the adiabatic separation of electronic and nuclear motion in a molecule, thus it is a more general concept that does not rely on the particular representation (adiabatic vs. diabatic, for instance) used to describe the electronic subsystem.

**Density functionalization** The XF has been used to extend density functional theory (DFT) to non-adiabatic situations, where the conditional electron density $n_R(\underline{r}) = N_e \int |\Phi_R(\underline{r}, \underline{r}_2...\underline{r}_N)|^2 d^3r_2...d^3r_N$ replaces the $N_e$-body wavefunction as basic variable [78] of the theory. The full nuclear wavefunction $\chi(\underline{R})$ is retained, and the exchange-correlation functional of the usual adiabatic DFT is generalized to include nonadiabatic contributions arising from the $\hat{U}_{en}$ coupling term; it has functional dependence on the conditional electronic density and paramagnetic current-density, nuclear wavefunction, vector potential, and the quantum geometric tensor. Ref. [79] proposed a local conditional density approximation, and applied the XF-based DFT to study dissociation in a Hubbard-model of the LiF molecule, showing that non-adiabatic electron-nuclear coupling results in a 0.5 Bohr elongation of the bond-length at which electron-transfer is onset.

**Exact electron factorization** The XF formalism can be applied within a purely electronic system, where the marginal factor is a function of one electronic coordinate $r_1$, with the conditional factor a function of the rest, conditionally dependent on $r_1$ [80]. From the partial normalization condition and the antisymmetry of the full electronic wavefunction, it follows that the marginal density gives the exact one-body electronic density and current-density. The exact XF potentials for an electron ionized by a field have been studied [80,81], and approximated, and they provide a new perspective for one-electron approaches commonly used in strong-field physics, such as the single active electron approximation.

**Electronic embedding: strong correlation** The XF idea can also be applied to an electronic wavefunction in Fock space [82–84], where some orbitals (e.g., the more strongly correlated ones) are selected to be those of the marginal (the "fragment") while all others are conditionally dependent on the occupation numbers of the marginal orbitals. XF then yields an embedding Hamiltonian on the fragment, which is approximated by solving the full system at a low-level; the embedded fragment is then solved with a high-level method. This has been demonstrated to be accurate for ground-state energies over a full range of weak to strong correlation in Hubbard model systems [83]. A generalized Kohn-Sham approach has also been used to find the conditional factor [84], using an orbital-dependent functional approximation, which accurately captured the topological phase diagram of a multiband Hubbard model.

**Reverse factorization** There is nothing in the factorization presented above that requires us to choose the marginal as the nuclear coordinate, and the conditional as the electronic. It would be formally equally possible to choose the marginal factor to be the wavefunction for the electronic coordinate, and then the nuclear wavefunction is conditionally dependent on the electronic one. Although less intuitive, this is particularly useful for problems where we are primarily interested in electronic dynamics, since it yields a TDSE for the electrons, in which the potentials contain the effect of coupling to the nuclei. This was first studied in ref. [19] with the example of laser-control of electron localization: the time-delay between the excitation of a nuclear wavepacket from the ground to the first excited state of $H_2^+$ and an applied infrared pulse controls the localization asymmetry, i.e., the probability the electron ends up on the "left" or "right" nucleus. We will return to this briefly in Section 3 [19,20,27]. For static problems, the idea of a marginal-conditional decomposition was applied to relative coordinates in two-electron atoms, to define an exact central field
model and analyze the Coulomb hole [85]. Recently, the idea was extended to three-particle Coulomb systems to characterize the atom-molecule transition [86].

*Mathematical analysis* Formally, the XF expression of the time-dependent molecular wavefunction closely resembles the BO approximation. However, in the latter case, the electronic conditional term is constrained to be an eigenstate of the BO Hamiltonian $\hat{H}_{BO}(\underline{r}, \underline{R})$ all along the dynamics and, thus, it does not depend on time. Ref. [87] showed that the BO approximation can be recovered from XF as a limiting case for the electron-nuclear mass ratio $\mu$ tending to zero. Not only is this procedure an elegant scheme to connect XF to the well-known BO approximation, but it has been employed as well to develop algorithms that solve the coupled electronic (4) and nuclear (5) equations of XF. More specifically, the expansion of those equations in powers of $\mu$ helped in identifying leading non-adiabatic contributions in the electronic Equation (4) to be retained in the approximate CTMQC algorithm (see Section 4), as well as in suggesting strategies for a perturbative treatment of non-adiabatic effects [88–91] (discussed below).

It is of fundamental interest whether the set of equations, Equations (4)–(9), are numerically stable to propagate, but it is also of practical interest, since in making approximations, we would like to know which terms might need care in developing methods that converge stably. In fact the mathematical form of the equations is unprecedented, and the usual numerical methods for TDSE’s fail when applied in a straightforward way to the exact equations [92]; a deeper mathematical analysis can be found in ref. [93].

*Polaritonic chemistry* The application of XF to systems of coupled electrons, nuclei, and photons [94], has provided new perspectives in the burgeoning field of polaritonic chemistry [95]. Confining a molecule to an optical or plasmonic cavity enhances the light-matter coupling strength that scales as the square-root of the volume, and distorts potential energy surfaces according to how close the energy gaps at a given nuclear configuration are to the resonant frequencies of the cavity, thus altering reaction pathways. Choosing the marginal factor as the nuclear coordinate provides the exact TDPES that drives the nuclear motion, directly correlating with the nuclear dynamics in contrast to the polaritonic surfaces [96] which are generalizations of the BO surfaces; ref. [97] showed how its features result in the suppression of proton-coupled electron-transfer in a model system, unlike the polaritonic surfaces. To properly model the complex systems used in the experiments, approximations such as MQC would be required, and running classical dynamics on the exact TDPES [98] gives a big improvement than when using simply the weighted polaritonic surface. Work is in progress to extend the self-consistent XF-based MQC methods to the polaritonic case for practical applications. To account for many cavity modes, a classical Ehrenfest treatment for the dynamics of the photonic system has been explored [99], and using the XF with the marginal as the photonic displacement coordinate could explain the general underestimation of photon numbers in this approach [100] through comparing the exact potential driving the photons with that underlying the Ehrenfest approach. A further useful factorization could be $\Psi(\underline{r}, \underline{R}, \underline{q}, t) = \chi(\underline{R}, \underline{q}, t)\Phi_{\underline{R}, \underline{q}}(\underline{r}, t)$ since this would yield a TDSE for both the nuclear and photonic systems in which the potential “exactifies” the concept of cavity-BO surfaces introduced in ref. [101]. Also, in general, when several identical, or non-identical particles compose the subsystems, nested commutators along the lines of the formalism of ref. [102] could be considered.

*Perturbative limit* The XF has been employed to address diverse classes of static and dynamic problems in molecules [88–90] and solids [91] manifesting essential, though not strong, non-adiabatic effects. For instance, it has been observed that, in determining vibrational circular dichroism of chiral molecules [90,103,104] or in computing electronic current densities [88,105,106], the BO approximation fails to fully account for electronic effects during the dynamics. Even though non-adiabaticity is not strong in these examples because the system does not evolve close to avoided crossings or conical intersections, including excited-state effects as a perturbation to the BO (unperturbed) state is essential—and sufficient—to recover the correct behavior. In the context of XF, the electron-nuclear coupling term of Equation (7) has been treated as a perturbation to the BO Hamiltonian [90],
where the nuclear velocity is the small parameter. This results in a vector potential that can be seen as a position-dependent mass correction to the bare nuclear masses and that can be used to correct vibrational harmonic frequencies of hydrogen-based molecules [89] or phonon frequencies in solids [91].

*Quantum-mechanical clock* The separation of degrees of freedom at the basis of the stationary version of the XF was used in refs. [107,108] to introduce the concept of time in quantum mechanics via the classical limit of the quantum-mechanical clock's degrees of freedom $\underline{R}$ represented by the marginal wavefunction. A clock-dependent continuity equation for the conditional wavefunction was then derived.

### 3. Exact Potentials Driving Dissociation and Ionization

Returning to the theme of dynamics of non-bound states, the exact TDPES of the XF has proved to be a useful interpretative tool. Consider a laser field causing dissociation and ionization in a diatomic molecule. If we are primarily interested in the dissociation aspect, the original factorization of Equations (4)–(9) might be most helpful for analysis since it provides a rigorous definition of the potentials that drive the nuclear motion. If we are interested in the ionization aspect, we may instead want to utilize the reverse factorization, to study the exact potential driving the electronic system. In this section, we discuss examples of these two potentials, previously studied in refs. [16,17] and refs. [20,27] respectively, for the case of a one-dimensional (1D) model of $H_2^+$, subject to a laser field.

The model used in refs. [16,17] is from Refs. [109–112], and involves one electronic coordinate and one internuclear (relative) coordinate with soft-Coulomb-interactions. We note that the laser does not couple directly to the relative nuclear coordinate: the dissociation is driven solely by the electronic motion. In the XF picture, from the nuclear viewpoint this is represented entirely by the TDPES and the time-dependent vector potential. For one-dimensional problems, we can always pick a gauge in which the vector potential vanishes everywhere, and then the TDPES represents the complete potential driving the nuclear system. This is the gauge chosen in this study [16,17], and the TDPES is obtained by inversion of Equation (5). The two-dimensional TDSE for the electron-nuclear wavefunction is first solved, starting in the ground state of the molecule, and then $\chi$ is constructed: its magnitude is the nuclear density obtained by integrating the full wavefunction over electronic coordinates and its phase is settled by the condition that $A(R,t) = 0$ [17].

The laser field has a wavelength of 228 nm, which corresponds to a frequency of about twice the dissociation energy. A sketch is shown in the top middle panel of Figure 1. We see from the black curves in the upper left panel (b), that for an intensity of $I_1 = 10^{14}$ W/cm$^2$, the system dissociates, accompanied by significant ionization. This suggests a Coulomb explosion mechanism. On the other hand, for the weaker intensity of $I_2 = 2.5 \times 10^{13}$ W/cm$^2$, over the same time period there is dissociation with $\langle R \rangle$ reaching about half the distance as for the stronger field, but very little ionization (0.008 of an electron by the end of the time shown, as opposed to a little more than 0.5 for the stronger field). The other curves in these plots show the predictions from approximate methods: Ehrenfest, time-dependent Hartree, and exact-Ehrenfest. In the Ehrenfest approach [113,114], a single nuclear classical trajectory, initially located at the $\langle R \rangle$ of the initial nuclear probability density and with zero momentum, is propagated by the classical Hamilton's equations with a force given by the average of the soft-Coulomb nuclear-nuclear interaction and electron-density-weighted electron-nuclear interaction. In exact-Ehrenfest, on the other hand, a classical nuclear trajectory is propagated under the gradient of the exact TDPES [16,17]. These methods both treat the nuclear dynamics classically, while time-dependent Hartree (time-dependent self-consistent field) involves a quantum propagation of a nuclear wavefunction based on an uncorrelated factorization of the molecular wavefunction [114] resulting in the potential driving the nuclei determined self-consistently through the electronic density. We observe that while all three of the approximate methods predict dissociation in the stronger field $I_1$, with exact-Ehrenfest being closest to the exact result, none of them is able to predict the dissociation in the weaker field $I_2$ case. The exact TDPES can help us understand why.
The TDPES, shown in the lower set of panels in Figure 1b,c, shows a strong deviation from the BO surface, with large oscillations in time in the region from near the equilibrium distance on out. The tail of the TDPES alternately falls sharply and returns in correspondence with the field, which softens the bond, releasing the nuclear density. In this way the TDPES mediates the transfer of energy from the field-accelerated electronic motion to the nuclei. The features are similar for both the stronger and weaker field cases, and why a classical particle evolving on the TDPES does not therefore dissociate for the case of $I_2$ is perhaps not immediately obvious. The answer is revealed when we take a closer look at the structure of the surfaces close to the equilibrium distance and plot the classical energy of the particle evolving in the TDPES (the exact-Ehrenfest calculation) as a function of its position. These are black dots shown in the time snapshots. While for the stronger field $I_1$, the particle quickly becomes loose of the binding potential, being kicked up by the TDPES oscillations, in the case of the weaker field it gets trapped in a narrow potential well not so far from the equilibrium bond length. Without tunneling, it cannot escape, and therefore dissociation does not occur in the classical calculation. This shows that the mechanism for dissociation of the quantum system in the weaker field is tunneling: the TDPES transfers field energy to the nuclei via tunneling, and although the classical exact-Ehrenfest calculation shows a larger amplitude of oscillation of $\langle R \rangle$ than the other approximation methods, it ultimately cannot tunnel through the barrier. Despite treating the nuclei quantum mechanically, so in principle being able to capture effects like tunneling, the time-dependent Hartree method still fails. This is because of its uncorrelated nature: the electronic density determining the potential for the nuclear motion is the same for every $R$, and does not deviate far from the exact electronic density near the equilibrium bond-length, due to the weak field strength. The reader is referred to ref. [17] where the Hartree surfaces are explicitly shown. This example emphasizes that the mechanism of dissociation via tunneling requires not only an quantum mechanical description of the nuclei but also an adequate description of electron-nuclear correlation.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** 1D model of laser-driven $H_2^+$ dissociation: (a) Shows a sketch of the model molecule, and the electric field of the $\lambda = 228$ nm laser applied to the system, with the shaded part indicating the optical cycle for which snapshots for two different intensities are shown in the other panels. (b) The top panel shows $\langle R \rangle(t)$ as a measure of dissociation, for the stronger field intensity of $I_1$, and the lower panel shows $I_p = 1 - \int_{\text{box}_e} dz \int dR |\Psi(z, R, t)|^2$ which measures ionization through the number of electrons outside a box chosen of size $|z| \le 10$ a.u; the exact is shown in black, while predictions of the traditional classical Ehrenfest method is in red, the quantum time-dependent Hartree in blue, and a classical evolution on the exact TDPES (labelled as exact-Ehrenfest) as green. The lowest panel shows time-snapshots of the nuclear density and the TDPES, and the black dot shows the position and energy of a classical particle evolving under the TDPES. (c) This shows the same quantities as in (b) but for the weaker field intensity $I_2$. Reproduced from ref. [17] with the permission of AIP Publishing.</figcaption></figure>
We now briefly turn to the electron’s viewpoint, and ask what is the potential that is
driving its motion? In this case, we choose the reverse factorization where the marginal
has the electronic coordinate: $\Psi(\underline{r}, \underline{R}, t) = \Phi_{\underline{r}}(\underline{R}, t)\chi(\underline{r}, t)$, so that the equations are exactly
the same as Equations (4)–(9) but with $\underline{r} \leftrightarrow \underline{R}$. The electronic wavefunction satisfies a
TDSE, and, again picking a gauge where $\underline{A} = 0$, the TDPES appearing in the equation for
$\Phi(\underline{r}, t)$ (e-TDPES) together with the external laser field account completely for the electronic
system’s motion.

Although formally the same as the direct factorization equations Equations (4)–(9), the structure of the TDPES is significantly different, with different terms in Equation (8) contributing in different ways. In particular, the term in the e-TDPES that comes from the electron-nuclear coupling operator $\hat{U}_{en}$, $\frac{1}{2m_e}\langle\Phi_{\underline{r}}|\nabla_j^2\Phi_{\underline{r}}\rangle_{\underline{R}}$ (where $\nabla_j$ is the gradient with respect to the jth electronic coordinate) is significantly larger in the reverse factorization than in the direct factorization, due to the division by the much smaller electron mass. This term tends to give peak structures in the e-TDPES, which can have a significant impact on the electron dynamics. For example, in the case of laser-controlled electron localization [115–117], mentioned in Section 2.1, the peak led to an interatomic barrier that enhanced the one in the traditional pictures that comes from Coulomb interaction, resulting in an earlier localization time and smaller localization asymmetry than predicted by traditional methods [19]. The traditional methods use a “quasistatic” picture where the electron-nuclear electrostatic potential is added to the laser field to determine the electronic wavefunction’s evolution. This work showed that the quasistatic approach misses important effects coming from dynamical electron-nuclear correlation, that can have a notable influence on the predicted observables.

Dynamical electron-nuclear correlation was also shown to be important in under-
standing the process of charge-resonance enhanced ionization (CREI). CREI refers to the
phenomenon that the ionization rate from a molecule can be several orders of magnitude
higher than the rate from the constituent atoms at a critical range of internuclear sepa-
rations [110,118–122]. This was originally explained by a quasistatic argument, where
nuclei are instantaneously fixed point particles. The combined electrostatic electron-nuclear
attraction, together with the field, leads to Stark-shifted atomic levels, and the potential
acting on the electrons involves both an outer barrier on the down-field side, and an in-
ternuclear barrier that grows as internuclear separation increases. Provided that the field
is turned on fast enough such that population in the up-field level remains and does not
tunnel back to the down-field atomic level, the molecule can rapidly ionize over both
the inner and outer barriers, at an optimal internuclear separation. This gives rise to the
enhanced ionization rate. By requiring that the Stark-shifted up-field atomic level exceeds
the top of both the inner and outer field modified Coulombic barriers, one finds the critical
internuclear separation for CREI as 4.07 divided by the ionization potential. The analysis
can be generalized to the case of a laser field whose period is shorter than the tunneling
time [118,123–125].

This physical picture is very instructive, but in modeling the experiment, the premise
of a quasi-static picture is questionable. A time-resolved study showed that even with static
nuclei, the ionization tends to occur in multiple sub-cycle bursts [124,125], not at the peak
of the field cycles, which was assumed in the prior analysis. Furthermore, to observe CREI,
the molecule needs to be stretched to the CREI region rapidly enough that appreciable
electron density remains un-ionized. In fact, in some experiments, CREI is not observed
because there is too much ionization before the critical distance is reached [126–128].
Typically, a fraction of the nuclear density remains near equilibrium separations, while a
fraction dissociates, and it is the electron density associated with the dissociating fragment
that is responsible for CREI. Thus not only the point particle picture for the nuclei is not
accurate, but also the nuclear dynamics means the quasistatic picture is not appropriate,
especially as in some dissociating channels the fragment velocities can be comparable to the
electronic velocities. Refs. [20,27] showed that the dynamical electron-nuclear correlation
terms in the exact e-TDPES are needed to accurately predict the ionization characteristics
when modeling a real CREI experiment. Going beyond the quasistatic treatment by only accounting for the width and splitting of the nuclear wave packet is generally not enough to get the correct dynamics of CREI.

These examples show that the exact TDPES for nuclei or e-TDPES for electrons can provide an illuminating picture on the nature of dynamical electron-nuclear correlation in laser-driven dissociation and ionization, but such studies can only be done when the full molecular TDSE can be somehow solved exactly or highly accurately. In the next section, we discuss approximations based on the XF, which move the XF to be a useful approach for real systems.

## 4. XF-Based Mixed Quantum-Classical Methods

Being an exact reformulation of the molecular TDSE, the XF equations are no easier to solve than the original molecular TDSE; in fact, due to numerical stability issues, they appear to be harder [92,93]! Although the form of the exact coupling terms can give us insight and help us interpret phenomena related to electron-nuclear correlation, to have a practical impact, the XF requires approximations. MQC approximations are a natural direction for this, justified by the large mass of the nuclei [87].

By taking a classical limit for the nuclear motion, and discarding terms in the equations justified by the exact studies to have a small effect [33], refs. [24,25] derived the coupled-trajectory MQC (CTMQC) approach. In CTMQC the electrons are propagated quantum-mechanically whereas the nuclear dynamics consists of classical trajectories coupled through the nuclear quantum momentum

$$Q_v(\underline{\mathbf{R}}, t) = - \frac{\nabla_v |\chi(\underline{\mathbf{R}}, t)|}{|\chi(\underline{\mathbf{R}}, t)|} = - \frac{\nabla_v |\chi(\underline{\mathbf{R}}, t)|^2}{2|\chi(\underline{\mathbf{R}}, t)|^2} \quad (10)$$

When writing the nuclear wavefunction in terms of its modulus $|\chi(\underline{\mathbf{R}}, t)|$ and phase $S(\underline{\mathbf{R}}, t)$, the term in parenthesis in Equation (7) that depends on $\chi(\underline{\mathbf{R}}, t)$ can be rewritten as

$$\frac{-i\nabla_v \chi(\underline{\mathbf{R}}, t)}{\chi(\underline{\mathbf{R}}, t)} + \mathbf{A}_v(\underline{\mathbf{R}}, t) = \nabla_v S(\underline{\mathbf{R}}, t) + \mathbf{A}_v(\underline{\mathbf{R}}, t) + i\mathbf{Q}_v(\underline{\mathbf{R}}, t) \quad (11)$$

The first two terms on the right-hand side sum up to the nuclear momentum field, which has a clear classical interpretation. The last term, i.e., the quantum momentum, appears as an imaginary correction to the classical momentum, with no classical counterpart.

The conditional electronic wavefunction associated with each trajectory $\alpha$ is expanded in the BO basis $\{\Phi_{\underline{\mathbf{R}}}^l(\underline{\mathbf{r}})\}$ such that $\Phi_{\underline{\mathbf{R}}}^l(\underline{\mathbf{r}}, t) \to \Phi_{\underline{\mathbf{R}}^\alpha(t)}(\underline{\mathbf{r}}, t) = \sum_l C_l^\alpha(t) \Phi_{\underline{\mathbf{R}}^\alpha(t)}^l(\underline{\mathbf{r}})$. Then the CTMQC equations for the electronic and nuclear evolution are

$$\dot{C}_l^{(\alpha)}(t) = \dot{C}_{l,Eh}^{(\alpha)}(t) + \dot{C}_{l,XF}^{(\alpha)}(t) \quad (12)$$

$$\mathbf{F}_v^{(\alpha)}(t) = \mathbf{F}_{v,Eh}^{(\alpha)}(t) + \mathbf{F}_{v,XF}^{(\alpha)}(t) \quad (13)$$

where we note the shorthand $f(\underline{\mathbf{R}}^{(\alpha)}(t), t) = f^{(\alpha)}(t)$. The first terms in the electronic and nuclear equations are Ehrenfest-like terms

$$\dot{C}_{l,Eh}^{(\alpha)}(t) = -iE_l^{(\alpha)} C_l^{(\alpha)}(t) - \sum_k \sum_{\nu}^{N_n} \dot{\mathbf{R}}_v^{(\alpha)}(t) \cdot \mathbf{d}_{v,lk}^{(\alpha)} C_k^{(\alpha)}(t) \quad (14)$$

$$\mathbf{F}_{v,Eh}^{(\alpha)}(t) = \sum_l |C_l^{(\alpha)}(t)|^2 (-\nabla_v E_l^{(\alpha)}) + \sum_{l,k} C_l^{(\alpha)*}(t) C_k^{(\alpha)}(t) (E_l^{(\alpha)} - E_k^{(\alpha)}) \mathbf{d}_{v,lk}^{(\alpha)} \quad (15)$$

with $\mathbf{d}_{v,lk}^{(\alpha)}$ being the non-adiabatic coupling vector (NAC) along the $v$'th nuclear coordinate between BO states $l$ and $k$ evaluated at the coordinate $\underline{\mathbf{R}}^\alpha(t)$, i.e., $\langle \Phi_{\underline{\mathbf{R}}}^l(\underline{\mathbf{r}}) | \nabla_v \Phi_{\underline{\mathbf{R}}}^k(\underline{\mathbf{r}}) \rangle |_{\underline{\mathbf{R}}^\alpha}$ and
$E_k^{(\alpha)}$ the k-th BO electronic surface. The second terms in Equations (12) and (13) are the corrections coming from XF:

$$ \dot{C}_{l,XF}^{(\alpha)}(t) = \sum_{\nu}^{N_{\eta}} \frac{\mathbf{Q}_{\nu}^{(\alpha)}(t)}{M_{\nu}} \cdot \left( \mathbf{f}_{\nu,l}^{(\alpha)} - \sum_{k} |C_{k}^{(\alpha)}|^{2} \mathbf{f}_{\nu,k}^{(\alpha)}(t) \right) C_{l}^{(\alpha)}(t) \quad (16) $$

$$ \mathbf{F}_{\nu,XF}^{(\alpha)}(t) = \sum_{\mu} \frac{2\mathbf{Q}_{\mu}^{(\alpha)}(t)}{M_{\mu}} \cdot \left( \sum_{l} |C_{l}^{(\alpha)}(t)|^2 \mathbf{f}_{\mu,l}^{(\alpha)} \right) \left( \mathbf{f}_{\nu,l}^{(\alpha)} - \sum_{k} |C_{k}^{(\alpha)}|^{2} \mathbf{f}_{\nu,k}^{(\alpha)}(t) \right) \quad (17) $$

The term $\mathbf{f}_{\nu,l}^{(\alpha)}$ is a momentum that equals the time-integrated adiabatic force accumulated on the lth surface along the trajectory $\alpha$:

$$ \mathbf{f}_{\nu,l}^{(\alpha)} = - \int_{0}^{t} \nabla_{\nu} E_{l}^{(\alpha)} d\tau \quad (18) $$

and $\mathbf{Q}_{\nu}^{(\alpha)}(t)$ is the quantum momentum evaluated at the position of the trajectory $\mathbf{R}^{\alpha}(t)$. In a trajectory-based scheme, at any given time, information of the position of all trajectories is required to reconstruct the nuclear spatial distribution and compute the quantum momentum. Hence, this is a coupled-trajectory scheme, not an independent-trajectory one. CTMQC has been successfully employed to simulate several excited-state processes, such as the ring-opening process of oxirane [26,129] and the photo-isomerization of 2-cis-penta-2,4-dienimium cation, a retinal photoreceptor [130]. Key to its success in these problems is its ability to capture decoherence from first principles, unlike traditional MQC methods as we will shortly discuss.

Another approach is to use the electronic equation derived from XF, i.e., Equation (12), in a surface-hopping scheme [21,22] where the nuclei evolve on a single BO surface at any time, making hops between them according to a stochastic algorithm. This XF-based surface-hopping (SHXF) has been applied to a range of complex systems, from organic molecules to molecular motors [131–136]. Part of what has made SHXF able to treat such large systems, is that the quantum momentum term is computed via auxiliary trajectories thus enabling an independent trajectory scheme; these auxiliary trajectories approximately mimic the local coupling of a trajectory to nearby ones. Instead, computing the quantum momentum using coupled trajectories in the same way as in CTMQC, but within a surface-hopping scheme, yields the coupled-trajectory surface-hopping (CTSH) method [23]. More details on the different ways to compute the quantum momentum are given in Section 4.1, and an investigation of their effect on dynamics in Section 5.

In traditional MQC schemes such as Ehrenfest and surface-hopping [1,2,114,137], the electronic system suffers from “overcoherence”: after passing through a region of non-adiabatic interactions, the nuclear wavepacket splits with different parts of the wavepacket being correlated with different electronic surfaces. However, the Ehrenfest nuclear distribution follows a mean-field surface and does not split (even if it approximates averaged observables adequately), while, although the surface-hopping nuclear distribution is able to split with different branches evolving on different electronic surfaces, the electronic coefficients remain coherent. In a sense, surface-hopping is more unsettling than Ehrenfest, because of the disconnect between electronic and nuclear systems. Although it is clear how the Ehrenfest approach can be derived from a time-dependent self-consistent field ansatz [114], a solid derivation of surface-hopping has proved to be challenging with recent progress in ref. [138]. Various ad hoc decoherence corrections have been imposed on the surface-hopping algorithm [1,139,140], most often involving an empirical parameter; they are often adequate to reproduce experimental results, but not always, and not reliably.

On the other hand, due to its origin in the rigorous XF approach, the quantum momentum term in Equation (12) has been shown to give a first-principles description of decoherence in a number of complex systems as well as in detailed studies on model systems [21,25,26,130,135,141,142]; in some cases, the XF-based correction to surface-hopping
performs similarly, or better, than the traditional decoherence-corrected methods [135]. A significant improvement over traditional MQC methods was recently found in ref. [136]: here, the quantum-momentum term has a role that goes beyond just decoherence. It was shown to be crucial to describe processes where multiple electronic states are simultaneously coupled in some regions of the nuclear configuration space via non-adiabatic coupling.

## 4.1. Computation of the Quantum Momentum

Currently, available codes to perform CTMQC simulations are G-CTMQC [23,143] and PyUNIxMD [22]. G-CTMQC includes the coupled-trajectory schemes CTMQC and CTSH, with hopping calculated either via fewest-switches hopping, or Landau-Zener, and it is interfaced with the potential library ModelLib [144] Furthermore, it has been extended to treat spin-forbidden electronic transitions mediated by spin-orbit coupling, both in the spin-diabatic and spin-adiabatic representations [70,71]. PyUNIxMD has CTMQC and SHXF (with fewest switches hopping probability) capabilities, and is interfaced with a number of electronic structure codes. Recently, the simulation of photo-isomerization dynamics of a protonated Schiff base was calculated using both SHXF and CTMQC, achieving results that are in good agreement with each other [145].

In the CTMQC algorithm there are two main ways to compute the quantum momentum. The first is based on the idea of reconstructing the nuclear density using a sum of gaussians centered at the position of each classical trajectory

$$|\chi^{(\alpha)}(t)|^2 = \frac{1}{N_{tr}} \sum_{\beta}^{N_{tr}} \prod_{\nu=1}^{N_n} g_{\sigma_{\nu}^{(\beta)}}(\mathbf{R}_{\nu}^{(\alpha)}(t) - \mathbf{R}_{\nu}^{(\beta)}(t)) \quad (19)$$

Note that we have indicated with the symbol $g_{\sigma_{\nu}^{(\beta)}}(\mathbf{R}_{\nu}^{(\alpha)}(t) - \mathbf{R}_{\nu}^{(\beta)}(t))$ a normalized three-dimensional gaussian (for each nucleus $\nu$) centered at $\mathbf{R}_{\nu}^{(\beta)}(t)$ and with variance $\sigma_{\nu}^{(\beta)}$. Then, the quantum momentum is computed analytically applying Equation (10). In the original CTMQC implementation [24–26] the variance of the gaussians was time-dependent and determined from the distribution of classical trajectories along the dynamics. However, to stabilize the dynamics and for practical reasons, a frozen gaussian approach with time-independent widths is currently used. This results in the following expression for the quantum momentum:

$$\mathbf{Q}_{\nu}^{(\alpha)}(t) = \Gamma_{\nu}^{(\alpha)}(\mathbf{R}_{\nu}^{(\alpha)}(t) - \mathbf{R}_{\nu}^{0(\alpha)}(t)) \quad (20)$$

where the slope of the quantum momentum is computed as $\Gamma_{\nu}^{(\alpha)} = \sum_{\beta} W_{\nu}^{\alpha\beta}$ where

$$W_{\nu}^{\alpha\beta} = \frac{\prod_{\nu}^{N_n} g_{\sigma_{\nu}^{(\beta)}}(\mathbf{R}_{\nu}^{(\alpha)}(t) - \mathbf{R}_{\nu}^{(\beta)}(t))}{2\sigma_{\nu}^{(\beta)} \sum_{\beta'}^{N_{tr}} \prod_{\nu'}^{N_n} g_{\sigma_{\nu'}^{(\beta')}}(\mathbf{R}_{\nu'}^{(\alpha)}(t) - \mathbf{R}_{\nu'}^{(\beta')}(t))}, \quad (21)$$

and the center of the quantum momentum is defined as

$$\mathcal{R}_{\nu}^{0(\alpha)}(t) = \frac{\sum_{\beta'}^{N_{tr}} W_{\nu}^{\alpha\beta'} \mathbf{R}_{\nu}^{(\beta')}(t)}{\sum_{\beta} W_{\nu}^{\alpha\beta}} \quad (22)$$

Note that if the gaussians for all the trajectories carry the same nucleus-dependent width $\sigma_{\nu}$, given e.g., by the widths of the initial nuclear wavepacket, then the slope becomes independent of the trajectory index:

$$\Gamma_{\nu}^{(\alpha)} \rightarrow \Gamma_{\nu} = \frac{1}{2\sigma_{\nu}^{2}} \quad (23)$$
Because of the approximations introduced to derive CTMQC, the expression for the quantum momentum Equation (20) does not satisfy the physical condition that no population transfer among electronic states should be observed in regions where the NAC vectors are zero, once averaged over all trajectories. The CTMQC equation for the time-variation of the trajectory-averaged population of the lth electronic state, $P_l(t) = N_{tr}^{-1} \sum_{\alpha} |C_l^{(\alpha)}(t)|^2$ is:

$$ \frac{dP_l(t)}{dt} = \frac{1}{N_{tr}} \sum_{\nu\alpha} \left( -\sum_k \mathrm{Re}(\rho_{lk}^{(\alpha)}) \mathbf{d}_{\nu,lk}^{(\alpha)} \cdot \dot{\mathbf{R}}_{\nu}^{(\alpha)} + \frac{2\mathbf{Q}_{\nu}^{(\alpha)}}{M_{\nu}} \cdot \left( \mathbf{f}_{\nu,l}^{(\alpha)} - \sum_k \rho_{kk}^{(\alpha)} \mathbf{f}_{\nu,k}^{(\alpha)} \right) \rho_{ll}^{(\alpha)} \right) \quad (24) $$

where $\rho_{lk}^{(\alpha)} = C_l^{(\alpha)*} C_k^{(\alpha)}$ are electronic density-matrix elements. To enforce the condition of zero net population transfer when $\mathbf{d}_{\nu,lk}^{(\alpha)} = 0$, the second term on the right-hand-side of Equation (24) is required to be zero, and this gives an expression for the quantum momentum [26]. Note that this means the overall population transfer of the whole swarm of trajectories comes solely from the NAC term in Equation (24) and the quantum momentum affects the electronic population dynamics only indirectly through the evolution of the coefficients appearing in the first term. The condition is imposed pairwise in the electronic states and for each degree of freedom, which results in a different quantum momentum center for each pair of electronic states which is required to be independent of the trajectory index, $\mathcal{R}_{\nu}^{0(\alpha)}(t) \to \mathcal{R}_{lk,\nu}^{0}$:

$$ \sum_{\nu\alpha} \frac{2\mathbf{Q}_{lk,\nu}^{(\alpha)}(t)}{M_{\nu}} \cdot (\mathbf{f}_{\nu,l}^{(\alpha)} - \mathbf{f}_{\nu,k}^{(\alpha)}(t)) \rho_{ll}^{(\alpha)}(t) \rho_{kk}^{(\alpha)}(t) = 0 \quad (25) $$

Imposing the condition, the numerical expression satisfies

$$ \mathbf{Q}_{lk,\nu}^{(\alpha)}(t) = \Gamma_{\nu}^{(\alpha)}(\mathbf{R}_{\nu}^{(\alpha)}(t) - \mathcal{R}_{lk,\nu}^{0}(t)) \quad (26) $$

where now the center of the quantum momentum is computed as

$$ \mathcal{R}_{lk,\nu}^{0} = \frac{\sum_{\alpha} \Gamma_{\nu}^{(\alpha)} R_{\nu}^{(\alpha)} (f_{k,\nu}^{(\alpha)}(t) - f_{l,\nu}^{(\alpha)}(t)) |C_{l}^{(\alpha)}(t)|^2 |C_{k}^{(\alpha)}(t)|^2}{\sum_{\alpha'} \Gamma_{\nu}^{(\alpha')} (f_{k,\nu}^{(\alpha')}(t) - f_{l,\nu}^{(\alpha')}(t)) |C_{l}^{(\alpha')}(t)|^2 |C_{k}^{(\alpha')}(t)|^2} \quad (27) $$

A alternative approach is to choose the y-intercept of Equation (20), $\Gamma_{\nu}^{(\alpha)} \mathcal{R}_{\nu}^{0(\alpha)}(t)$, as independent of the trajectory index instead [26].

In the implementation of the CTMQC algorithm the expression of the quantum momentum center to be used, either the original Equation (22) or the modified Equation (27), is selected according to a given threshold (referred to as the M-parameter in the code). If the population decoheres and the wavefunction collapses to a single BO state the quantum momentum is set to zero. If the denominator of Equation (27) is very small for other reasons different than decoherence, then the analytical expression (22) for the center of the quantum momentum is used.

The CTSH algorithm also requires the quantum momentum in the electronic equation, and it is implemented with the same two options as possibilities.

Instead, in the independent-trajectory surface-hopping scheme SHXF, auxiliary trajectories are used to compute the quantum momentum. In the SHXF algorithm, when the propagated nuclear trajectory enters a region of non-zero non-adiabatic coupling an auxiliary trajectory is created on the BO state in which the population appears. This virtual trajectory propagates with a momentum given by isotropic energy conservation until the BO population in that state becomes negligible and the auxiliary trajectory is killed. At each
time step, the quantum momentum is computed placing a gaussian centered at the position of each trajectory and applying the analytical definition Equation (10)

$$ \mathbf{Q}_{v}^{(\alpha)}(t) = -\frac{1}{2\sigma_{v}^{2}} \left( \mathbf{R}_{v}^{(\alpha)}(t) - \sum_{l} |C_{l}^{(\alpha)}(t)|^{2} \mathbf{R}_{l,v}^{(\alpha)}(t) \right) \quad (28) $$

The variance of the gaussian centered on the auxiliary trajectory $\sigma_v$ in each degree of freedom $v$ is an input parameter, and, to eliminate empiricism is chosen to be that of the initial distribution. This has given good results for the real molecules studied. Recently, a new implementation of SHXF [146] was proposed in which a time-dependent width $\sigma_v(t)$ is computed by imposing maximum overlap between nuclear wavefunctions on different BO states in phase space.

Both SHXF and CTSH follow the same equations, except for the computation of the quantum momentum, and, in particular, Equation (28) of SHXF approximately mimics the computation with coupled trajectories of CTSH using the original definition Equation (22) in Equation (20). Although similar in spirit, we see that there is a significant difference in centers of the quantum momentum (second term on the right of Equation (28) c.f. Equation (20)). For two-state problems, while all trajectories of the ensemble contribute in CTSH, only two trajectories enter in SHXF: the running state, and one auxiliary trajectory. In fact, in the two-state case, Equation (28) reduces to

$$ \mathbf{Q}_{v}^{(\alpha)} = -\frac{1}{2\sigma_{v}^{2}} \left( \mathbf{R}_{v,d}^{(\alpha)}(t) - \mathbf{R}_{v,in}^{(\alpha)}(t) \right) (1 - |C_{d}^{(\alpha)}(t)|^{2}) \quad (29) $$

where $\mathbf{R}_{v,d}^{(\alpha)}(t)$ denotes the position of the active trajectory (running state), with electronic coefficient $C_{d}^{(\alpha)}(t)$, and $\mathbf{R}_{v,in}^{(\alpha)}(t)$ the auxiliary trajectory. In contrast, the center of the quantum momentum in CTSH is vastly different in Equation (22), including contributions from all trajectories on all surfaces, even though the originating equation is the same.

The main advantages of SHXF are computational efficiency and numerical stability. However, both SHXF and CTSH suffer from deficiencies inherited from the ad hoc nature of the surface-hopping procedure itself: the ambiguity of velocity-adjustments after a hop [135,147–149], and the incorrect forces in the vicinity of a NAC region where the exact force tends to have a diabatic or averaged character rather than that from one surface. Further, although the disconnect between the single-surface evolution of the nuclei and the coherent electronic evolution is partially ameliorated with the XF term, reducing the difference between the electronic population at a given time and the fraction of trajectories evolving on a given surface, this so-called "internal consistency" is not guaranteed. Finally, from the XF side, we note that the modified computation of the quantum momentum Equation (27) is not possible within SHXF because the trajectories are not coupled to each other.

## 5. Computing the Quantum Momentum

We now investigate the effect on the dynamics of a one-dimensional system of the three ways of calculating the quantum momentum term: with coupled-trajectories (CTSH and CTMQC) using the original definition Equation (22), and using the modified implementation Equation (27), and with the auxiliary trajectories method (as done in SHXF). Further we compare the effect of the XF terms in the nuclear equation of CTMQC versus the surface-hopping schemes SHXF and CTSH.

We choose our model system as the Tully model with an extended coupling region with reflection (ECR) [137], where an incoming nuclear wavepacket enters a region of extended non-adiabatic coupling between two BO surfaces that are asymptotically parallel. Along with other Tully models, the performance of CTMQC was studied in detail in ref. [25]. In that work, a time-dependent variance was computed for the gaussians. Nevertheless the
results of ref. [25] are in good agreement with our work in which a frozen gaussian approach Equation (23) was used instead.

The analytical form of the electronic Hamiltonian matrix elements, the BO surfaces, and the NAC are shown in Figure 2.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Analytical form of the electronic Hamiltonian on the **left**, with the BO potential energy surfaces (red and black) and NAC (dashed green) shown on the **right** for ECR model.</figcaption></figure>

The parameters were chosen to be $A = 6 \times 10^{-4}$ a.u., $B = 0.1$ a.u. and $C = 0.90$ a.u. and the atomic mass 2000 a.u. We will study two cases, a Gaussian nuclear wavepacket incoming on the lower surface from the left with a higher ($k_0 = 30$ a.u.) and a lower ($k_0 = 10$ a.u.) initial momentum. In the NAC region, some population is transferred to the upper surface and in both cases, the off-diagonal elements of the density matrix in the BO basis (coherences) rise and then eventually vanish after passing the non-adiabatic region: For the higher initial momentum, the nuclear wavepacket associated with the lower surface moves faster, separating in nuclear space from the part that has transferred to the upper surface, while for the low initial momentum case one branch of the nuclear wavepacket is reflected and the other is transmitted. The decay of the coherence cannot be accounted for by uncorrected surface hopping since the electronic equations are propagated fully coherently along each independent trajectory [137], and the method becomes “internally inconsistent” in that the BO population determined from the modulus-square of the associated electronic coefficient differs from that determined by the fraction of nuclear trajectories evolving on that surface [1,2,114,137,139,140].

We ran 1000 Wigner-distributed trajectories starting on the lower BO surface for CTSH and SHXF, while CTMQC calculations converge already with 200 trajectories. The initial position was $R_0 = -15$ a.u. and the variance $\sigma_0$ of the gaussian wavepacket at time zero was chosen to be 20 times the inverse of the initial momentum. The time-step used in the calculations is 0.1 a.u.

5.1. ECR Model $k_0 = 30$ a.u.

As we can see from panel A) in Figure 3, both CTMQC and the CTSH electronic populations with the modified definition of the quantum momentum (Equation (27)) are on top of each other and closely reproduce the exact result. For CTSH, measuring the population instead by the fraction of trajectories evolving on the surface yields small oscillations around the exact results. On the other hand, panel A) in Figure 3 shows that both CTMQC and CTSH with the original definition of the quantum momentum predict a small population transfer back to the ground state after the first interaction region when considering the electronic populations. Interestingly, the fraction of trajectories measure, however, does not show this, and instead has a similar trend as the exact but with an underestimate of CTSH. Comparing with panel B) of Figure 3, we see that the trend in the electronic populations is the same in SHXF. That is, the electronic populations yield a back-transfer when the original definition of the quantum momentum is used, whether it
is computed via auxiliary trajectories as in SHXF or with coupled-trajectories as in CTSH. The two ways yield similar results for the populations, although it should be noted that the SHXF results show a dependence on the width-parameter of Equation (28). Taking $\sigma = \sigma_0$ arguably gives an overall closer result to the parameter-free CTSH population, although SHXF with $\sigma = \sigma_0/10$ matches closely the CTSH populations at earlier times. The fraction of trajectories in SHXF does not show the back-transfer, just as that in CTSH did not; this measure in SHXF with $\sigma = \sigma_0/10$ gives a larger underestimate of the population transfer. But the fraction of trajectories measure for either $\sigma$-parameter in SHXF as well as either definition of the quantum momentum in CTSH are quite close to the exact, and close to that of CTMQC with the modified definition. At the same time we observe that uncorrected SH electronic populations are close to that of the exact while the fraction of trajectories closely resembles that of CTSH. However the uncorrected SH coherences do not decay in time and would lead to errors in later evolution if the system encounters another interaction region.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Upper panels: electronic populations (A,B), and of the excited state $\rho_{22}(t)$ (solid) and fraction of trajectories on the excited state $\Pi_2(t) = N_2(t)/N_{\text{traj}}$ (dotted), with $N_2(t)$ being the number of trajectories running on the upper surface at a given time and $N_{\text{traj}}$ the total number of trajectories, are plotted for ECR model with $k_0 = 30$ a.u., as functions of time starting at 400 a.u. which is when the trajectories and the quantum wavepacket approach the NAC region. The label $Q_m$ indicates the use of Equation (27), i.e., the modified definition of the quantum momentum, while $Q_0$ indicates the use of Equation (22), i.e., the original definition. Lower panels: electronic coherences (C,D), $|\rho_{12}|^2 = \sum_{\alpha}^{N_{\text{traj}}} |C_1^{\alpha*}(t)C_2^{\alpha}(t)|^2/N_{\text{traj}}$. The color code is the same as in the upper panels. The coupled-trajectory results are compared with exact results in the left panels, whereas auxiliary-trajectory results are shown in the right panels. For reference, surface-hopping (SH) results with no decoherence corrections are shown in the right panels.</figcaption></figure>

The incorrect back-transfer of population seen in the electronic populations when using the original definition of quantum momentum can be related directly to the violation of the condition of zero net population transfer in regions of zero NAC. To demonstrate this, we plot the spatial distribution of trajectories at different times in Figure 4 for CTSH and SHXF along with the NAC. The top panels show that the nuclear distributions for CTSH with the two definitions are similar at these times, and that between 1000 and 1200 a.u., only a small fraction of trajectories are in a region of non-negligible non-adiabatic coupling; so the population transfer seen with the original definition Equation (22) is a violation of the condition. The lower panels show that a similar violation occurs for SHXF $\sigma = \sigma_0$ and $\sigma_0/10$, consistent with the unphysical population transfer seen there. However, using the fraction
of trajectories measure of populations appears to remedy this problem: the fewest-switches hopping probability is proportional to the NAC so if this is too small, then the system is very unlikely to hop. Thus, in a sense, the fewest-switches hopping probability somewhat takes care of the zero net population transfer in regions of zero NAC, since the probability is zero when the NAC is. We trade one violation for another: the violation of the internal inconsistency enables us to obtain physically reasonable populations through the fraction of trajectories measure of populations instead of through the electronic coefficients, which avoids us violating the condition of zero net population transfer in regions of zero NAC.

Looking at the coherences plotted in panel C) of Figure 3, both CTSH and CTMQC with the modified definition of the quantum momentum are in agreement with each other and show faster decoherence rates than the exact calculations, while using the original definition of the quantum momentum predicts slower decoherence times. On the other hand, as we see in panel D) the coherences for SHXF $\sigma = \sigma_0$ are similar to those of CTSH with the original quantum momentum definition, and the decoherence rate increases with decreasing $\sigma$ as expected from Equation (28).

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Spatial distribution of trajectories for the ECR model with $k_0 = 30$ a.u.: CTSH with original definition of the quantum momentum using Equation (22) (panel A), CTSH with the modified definition of the quantum momentum using Equation (27) (panel B), SHXF with $\sigma = \sigma_0$ (panel C) and SHXF $\sigma = \sigma_0/10$ (panel D).</figcaption></figure>

## 5.2. ECR Model $k_0 = 10$ a.u.

As we can see in Figure 5, in contrast to the higher momentum case, there is substantial disagreement between the methods for the low incoming momentum case.

First let us consider the coupled-trajectory methods with the modified quantum momentum expression, Equation (27). CTMQC with Equation (27) closely reproduces the exact results, slightly overestimating the population transfer to the excited state and underestimating the back-transfer from the part of the nuclear wavepacket that reflects. Using Equation (27) the electronic populations of CTSH however do not reproduce the transfer back to the ground state that would be associated with reflecting trajectories. On the other hand, the CTSH fraction of trajectories reproduces the populations better than the electronic populations, even though, in a sense the latter are more fundamentally derived [23]; we had seen this also for the higher momentum case in Section 5.1 when the original definition of the quantum momentum was used, but now here we see it also for the modified definition. Unlike the higher momentum case, this cannot be attributed to violation of the zero net population transfer condition since the modified definition respects this condition.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Upper panels: electronic populations (A,B) of the excited state (solid) and fraction of trajectories on the excited state (dotted) are plotted for ECR model with *k*₀ = 10 a.u. Lower panels: electronic coherences (C,D).</figcaption></figure>

With the modified definition of the quantum momentum, although the CTSH trajectory-
averaged electronic populations are similar to CTMQC until about 4000 a.u., the coefficients
are different from about 3000 a.u. onwards, as is evident from the coherences shown in
the lower panel. CTMQC coherences decay after the first passage on a similar but shorter
time-scale as the exact before rising later when the reflected trajectories pass the NAC.
But CTSH coherences rise after a much shorter time: the CTSH trajectories spend more
time in the NAC region, reflecting earlier, and the coefficients of the individual trajectories
evolve, even if the population averaged over all trajectories remains constant. The quantum
momentum term (which depends on |C₁⁽α⁾(t)|²|C₂⁽α⁾(t)|²) in the equation of motion for the
coefficients remains activated in CTSH from 3200 a.u. onwards, unlike in CTMQC where
it is essentially zero from 3200 a.u. to about 4000 a.u. These observations are illustrated
in Figure 6. The figure also shows the electronic energies associated with each trajectory,
which, for the CTMQC maps out the effective approximate TDPES. While the CTSH trajecto-
ries dwell in the NAC region, there is a significant hopping probability, yielding changes in
which surface a given trajectory evolves on. Even at large times, the internal inconsistency
is severe in both the reflected and transmitted region. In particular, in the first non-adiabatic
event, some trajectories do not hop but the associated electronic coefficient changes to the
upper surface leading to a transmitted part of the nuclear distribution with poor internal
consistency as we observe on the trailing edge of the transmitted trajectories (Figure 6).
Arguably the trend of the fraction of trajectories measure of populations in Figure 5 is closer
to the exact than that from the electronic populations, just as the running state population
is closer to that of CTMQC in Figure 6. The reason why the CTSH trajectories stay longer in
the NAC region and reflect earlier than the CTMQC trajectories is related to the effective
potential energies they experience: the ones which have hopped to the upper surface evolve
on the more sharply rising upper BO surface, with a velocity penalty to conserve energy
from making the hop, while the CTMQC ones evolve on an approximate TDPES which
tends to be bounded by the upper BO surface but can be below it. The CTMQC trajectories
(black dots moving on the thin red lines representing the BO energy curves in the left
panels) leave the NAC region earlier than the CTSH ones (black dots in the right panels),
and move further to the right, before the reflection occurs. While at the beginning of the
non-adiabatic event (shown at time *t* = 2875 a.u. in Figure 6) the distributions of CTMQC
and CTSH trajectories are similar, the issues just described can be identified at later times,
for instance at *t* = 3215 a.u. and *t* = 3395 a.u. in Figure 6. The distribution generated by the CTSH evolution at *t* = 3795 a.u. and *t* = 4675 a.u. has a splitting that agrees better with CTMQC than at previous times. Movies describing the full dynamics as documented in Figure 6 are provided as Supplementary Material.

In Figure 6, comparing the CTSH electronic populations (green dots) in the top right-hand plots at each time with the running state (red dots, or, compare with the bottom plots) attests to the breakdown of the internal consistency of the algorithm, discussed above. Specifically, we note that in some regions of space and already at early times, the running state for each trajectory (the red dots are at zero for the ground electronic state and one for the excited state), differs from the value of the excited-state population for those same trajectories. When a trajectory is on the ground-state potential energy curve, the corresponding population associated to the excited state ($\rho_{22}^{\alpha}(t)$) should become zero as the trajectory leaves the interaction region. However, this is not always the case, as evident from Figure 6. For CTMQC, this issue of internal consistency does not exist. We show for completeness, in the same figure but in the middle plots of each panel, the contribution of the quantum-momentum XF term in the electronic evolution equation (in particular, we show as blue dots the second term in the right-hand side of Equation (24) but without the trajectory-sum). We note that at times $t = 3215, 3395, 3795$ a.u. the effect of such term is non-zero, in contrast to what is predicted by CTMQC.

We turn now to the CTMQC and CTSH results using the original definition of the quantum momentum Equation (22). With this definition we see from the lower left panel of Figure 5 that neither of these decohere during this time period. The CTMQC populations continue to evolve throughout, beginning to differ from the exact result already around about 2400 a.u. Soon after this time, many of the trajectories have moved passed the NAC region, as evident in Figure 7 (black dots in the bottom panels of each plot on the left shown at different times along the dynamics—those times are different than Figure 6 since the overall evolution is different if compared to the previous case where the modified quantum momentum was used), and there is violation of the zero net population transfer condition in this part of the nuclear wavepacket. In fact the contribution of the quantum momentum term to the equation of motion for the coefficient (blue dots on the left panels in Figure 7) can be larger in the transmitted part than it is in the NAC region, which is a big distinction with the modified quantum momentum result where it was zero in the transmitted region (compare Figures 6 and 7). Driven by the quantum momentum (which is predominantly negative in this leading edge [141]) populations in this region start transferring back to the ground state after previously having transferred to the excited state, earlier than when the modified quantum momentum definition is used. The effective TDPES shown by the electronic energies associated with each CTMQC trajectory shows a significant difference from that when using the modified definition of the quantum momentum in Figure 6, and does not reduce to the BO surfaces away from the interaction region as they should. The distribution of nuclear trajectories do not show a clear splitting. Finally, note that the population of the electronic excited state at the positions of the trajectories (green dots on the left in Figure 7) clearly attests to the lack of spatial splitting and efficient decoherence observed for the dynamics of the trajectories.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Top panels in each plot: Time snapshots of the electronic populations of the excited state for CTMQC (**left**) and CTSH (**right**), with the modified definition of the quantum momentum, at the positions of the trajectories (green dots) and running state (red dots). Middle panels in each plot: XF contribution to the time-variation of the electronic populations (blue dots), i.e., the second-term on the right-hand-side of Equation (24) but without the trajectory-sum. Lower panels in each plot: BO surfaces (thin red lines), NAC (thin dashed lines), and spatial distribution of the nuclear trajectories on the BO surfaces (black dots).</figcaption></figure>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Same as in Figure 6 but using the original definition of the quantum momentum.</figcaption></figure>

The electronic populations of CTSH with the original quantum momentum definition Equation (22) show a strikingly different behavior than that of CTMQC, with many continuing to rise after the first transfer event, before leveling off (see the right panels of Figure 7 for CTSH results to be compared to the CTMQC results just described). As in CTMQC with this definition, there is violation of the zero net population transfer condition in the transmitted wavepacket, but the populations behave very differently: each of the transmitted and reflected wavepackets have populations spanning both upper and lower surfaces in CTSH, and the unphysical growth of the trajectory-averaged populations in Figure 5 comes primarily with trajectories associated with the transmitted part of the wavepacket. In contrast, in CTMQC, there is a linear shape to the population going from more on the upper surface on the reflected side to more on the lower surface on the transmitted packet. Interestingly, the CTSH fraction of trajectories is similar to that with the modified definition
of the quantum momentum, and has a trend closer to the exact. Although the electronic populations differ, because the nuclear trajectories follow adiabatic forces and are in similar regions of space, the hopping probabilities are similar, since hopping can only occur in the region of the NAC, so, as observed in the higher-momentum case, in a sense ameliorates the violation of the condition of zero net transfer in regions of zero NAC. Movies describing the full dynamics as documented in Figure 7 are provided as Supplementary Material.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Spatial distribution of trajectories for the ECR $k = 10$ a.u. model. CTSH with original definition of the QM (panel A), CTSH with the modified definition of the QM (panel B), SHXF with $\sigma = \sigma_0$ (panel C) and SHXF $\sigma = \sigma_0/10$ (panel D).</figcaption></figure>

Turning now to the independent trajectory surface-hopping calculations shown in the upper right panel of Figure 5. We see that uncorrected surface hopping in fact looks very similar to CTSH with the modified definition of the quantum momentum, in both the populations as well as the fraction of trajectories but gives poor coherences (lower right panel). Including the XF correction with $\sigma = \sigma_0$, leads to a large back-transfer of electronic populations, but again with a similar fraction of trajectories as the uncorrected surface hopping and the CTSH cases. Taking instead $\sigma = \sigma_0/10$ further causes the system to decohere faster initially and improves the internal consistency but at the expense of worse agreement with the exact. A comparison of the histogram of the nuclear trajectories at different time-snapshots is shown in Figure 8 for CTSH with the modified and original definitions of quantum momentum and SHXF with the two different values of $\sigma$, showing they yield similar nuclear distributions as a function of time, consistent with the fraction of trajectories being similar. At the later times when the transmitted part of the wavepacket has left the NAC region, the associated populations of CTSH $Q_0$ and SHXF however continue to evolve (compare with Figure 7 and the movies), consistent with the violation of the condition of zero net population transfer.

Interestingly, the trend in the populations of SHXF is more similar to that in CTMQC with Equation (22) and opposite to that of CTSH (with Equation (22)), after the first interaction region. Despite being similar in the underlying equations, this may not be surprising given the very different implementation of Equation (22) in the two methods (see discussion at the end of Section 4.1). In the case of the higher-momentum $k = 30$ a.u., the population trends in CTSH and SHXF were in fact similar: likely because there is essentially one localized wavepacket on the upper surface and one on the lower, so the coupled-trajectory computation of the quantum momentum with Equation (22) is similar to the auxiliary-trajectory way which takes simply each independent active trajectory plus one auxiliary trajectory on the other surface. For the lower momentum case $k = 10$ a.u., two wavepackets develop on the lower surface going in opposite directions, as well as a
wavepacket on the upper surface: while CTSH computes the quantum momentum from
contributions from all of these, again the SHXF just takes the independent active surface
trajectory plus one auxiliary trajectory on the other surface, trajectory-by-trajectory.

6. Discussion

The XF approach offers new insights into the interaction of coupled quantum subsystems, through its definition of exact coupling terms between the wavefunctions. For the degrees of freedom chosen as the marginal in the factorization, one obtains an effective TDSE in which the potentials contain the full effects of the correlated motion of the other subsystems. This has proven instructive in revealing molecular dissociation mechanisms under different laser parameters, as well as highlighting the effects of electron-nuclear correlation in ionization processes going beyond the usual quasi-static picture of the nuclear potentials. The generality of the XF idea has led to adventures in a wide range of applications outside coupled electron-nuclear dynamics, including electron-electron factorization to define a new single-active electron picture, quantum embedding for strongly correlated electronic systems, and inclusion of photons for molecular polaritonics.

The most practical impact of the XF approach so far has been in the development
of rigorous MQC methods, and CTMQC and SHXF have both successfully predicted
dynamics in a number of complex situations [26,129–135], giving improved results from
first-principles compared to traditional MQC methods [136]. These methods involve
computing a term that depends on the nuclear quantum momentum, for which three
different implementations exist. In particular, with coupled trajectories, there is the original
definition coming directly from $\nabla_v |\chi(\mathbf{R}, t)| / |\chi(\mathbf{R}, t)|$, as well as a modified definition
that ensures the physical condition that there should be no net population transfer in
regions of zero NAC. Further these two definitions can be used either in the full CTMQC
scheme where there are XF corrections in both electronic and nuclear equations of motion,
or in a surface-hopping scheme CTSH where they appear only in the electronic equation.
The quantum momentum may also be computed via auxiliary trajectories rather than
coupled trajectories, which is done in the surface-hopping scheme SHXF.

The effects of these different methods were studied here on two cases of dynamics
through Tully's model of an extended coupling region between two BO surfaces. In one
case, the incoming momentum of the nuclear wavepacket is high enough that there is a
single pass through the NAC region with significant transfer of population to the other BO
surface, and we showed that while CTMQC with the modified definition yields accurate
population transfer and coherence properties, CTMQC with the original definition of the
quantum momentum suffers from the violation of the condition. The original definition
however yields accurate population transfer and coherences when used within surface-
hopping framework, provided the fraction of trajectories measure is used, but at the
expense of violating internal consistency. This is because the hopping probability is only
non-zero in regions of the NAC, so in this sense, takes care of the condition. This is true
both for when the quantum momentum is computed using coupled trajectories, or from
auxiliary trajectories.

In the second case, the incoming momentum is low enough that part of the wavepacket
that transfers to the upper surface, reflects and recrosses the NAC region a second time.
Here CTMQC with the modified definition does a good job, although underestimates
the second rise in coherence, while the other methods are less accurate and also less in
agreement with each other. Using the modified definition in a surface-hopping approach
however yields poor electronic populations likely due to the incorrect forces on the nuclei
from the BO surfaces in the region of interaction, and incorrect imposition of energy
conservation on an individual trajectory level, that lead the trajectories to spend longer in
these regions. Using the original definition of the quantum momentum leads to unphysical
population transfer. Again the fraction of trajectories measure reduces the error and
temporarily, gives a reasonable trend.
The method that is most well-grounded from a theoretical point of view is CTMQC, and, with the modified definition of the quantum momentum, gives reliable results for both the populations and coherences. From a practical standpoint, it has the advantage of needing fewer trajectories for convergence, and no adjustable parameters. However, the method can be more numerically challenging to deal with than methods like surface hopping. In CTSH the forces that propagate the trajectories are purely adiabatic, and the NAC vectors $d_{v,lk}^{(\alpha)}$, which are computational expensive for molecular systems, are not needed along the dynamics unless the velocity rescaling after a surface hop has occurred is performed along the direction of the corresponding NAC vector. Note that in the electronic Equation (14) only the time-derivative coupling, i.e. $d_{v,lk}^{(\alpha)} \cdot \dot{R}_{v}^{(\alpha)}(t) = \langle \Phi_{\underline{R}}^{l}(\underline{r}) | \Phi_{\underline{R}}^{k}(\underline{r}) \rangle_{\underline{R}^{a}}$ is needed, and this term can be efficiently computed for instance as the overlap between the wavefunctions at subsequent time steps [150]. However with coupled trajectories, patience is required to wait for all the trajectories at each time-step in order to compute the coupled-trajectory term. The computational advantage of SHXF is evident due to its use of independent trajectories supplemented with auxiliary trajectories. So, in addition to the first-principles nature of the electronic equation, SHXF has an attraction due to its independent trajectories, but surface-hopping in itself is unsatisfactory due to the ad hoc aspects such as the very hopping procedure, choice of velocity adjustment, energy conservation on the independent trajectory level, and the choice of how to deal with forbidden hops. These problems contribute to the internal consistency, which is generally partly reduced by the XF-based quantum momentum term as seen in previous studies on decoherence, but we have seen here that it can also exacerbate it. Still, even then, we can sometimes *temporarily* take advantage of this internal inconsistency by using the fraction of trajectories measure, since it takes care of the condition that there should be zero net population transfer when there is zero NAC; this is only temporary since the evolving populations will continue to be wrong, and will give poor results when the system then meets another interaction region.

In conclusion, the XF offers a new playground for the development of non-adiabatic dynamics of electrons and ions, and we hope this work will be helpful to inform future adventures in this area.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/molecules27134002/s1, Video S1: Results for the ECR model for the lower initial momentum case. Top panel shows the time evolution of the electronic populations of the excited state for CTMQC with the modified definition of the quantum momentum at the positions of the trajectories (green dots), middle panel shows the XF contribution to the time-variation of the electronic populations (blue dots), and lower panel shows the BO surfaces (thin red lines), NAC (thin dashed lines), and evolution of the spatial distribution of the nuclear trajectories on the BO surfaces (black dots). Video S2: Middle and lower panels are the same as in video S1 but for CTSH with the modified definition of the quantum momentum. Middle panel also includes the running state in CTSH (red dots). Video S3: Same as in video S1 but for CTMQC with the original definition of the quantum momentum. Video S4: Same as in video S2 but for CTSH with the original definition of the quantum momentum.

**Author Contributions:** Conceptualization, N.T.M., E.V.A. and F.A.; methodology, F.A. and N.T.M.; software, F.A.; validation, N.T.M., E.V.A. and F.A.; formal analysis, N.T.M., E.V.A. and F.A.; investigation, E.V.A.; resources, N.T.M.; data curation, E.V.A.; writing—original draft preparation, E.V.A., N.T.M. and F.A.; writing—review and editing, N.T.M. and F.A.; visualization, E.V.A.; supervision, N.T.M. and F.A.; project administration, N.T.M.; funding acquisition, N.T.M. and F.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** Financial support from the National Science Foundation Award CHE-1940333 (EVA), and by 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 (N.T.M.) as part of the Computational Chemical Sciences Program, are gratefully acknowledged. FA acknowledges financial support by the ANR Q-DeLight project, Grant No. ANR-20-CE29-0014.
**Data Availability Statement:** The data presented in this study are available on reasonable request from the corresponding author.

**Acknowledgments:** We warmly thank Patricia Vindel Zandbergen for helpful discussions.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

## References

1. Crespo-Otero, R.; Barbatti, M. Recent Advances and Perspectives on Nonadiabatic Mixed Quantum-Classical Dynamics. *Chem. Rev.* **2018**, *118*, 7026–7068. [CrossRef] [PubMed]

2. Agostini, F.; Curchod, B.F.E. Different flavors of nonadiabatic molecular dynamics. *Wiley Interdiscip. Rev. Comput. Mol. Sci.* **2019**, *9*, e1417. [CrossRef]

3. Thachuk, M.; Ivanov, M.Y.; Wardlaw, D.M. A semiclassical approach to intense-field above-threshold dissociation in the long wavelength limit. *J. Chem. Phys.* **1996**, *105*, 4094–4104. [CrossRef]

4. Thachuk, M.; Ivanov, M.Y.; Wardlaw, D.M. A semiclassical approach to intense-field above-threshold dissociation in the long wavelength limit. II. Conservation principles and coherence in surface hopping. *J. Chem. Phys.* **1998**, *109*, 5747–5760. [CrossRef]

5. Bajo, J.J.; González-Vázquez, J.; Sola, I.; Santamaria, J.; Richer, M.; Marquetand, P.; González, L. Mixed quantum-classical dynamics in the adiabatic representation to simulate molecules driven by strong laser pulses. *J. Phys. Chem. A* **2012**, *116*, 2800. [CrossRef]

6. Horenko, I.; Schmidt, B.; Schütte, C. A theoretical model for molecules interacting with intense laser pulses: The Floquet-based quantum-classical Liouville equation. *J. Chem. Phys.* **2001**, *115*, 5733–5743. [CrossRef]

7. Fiedlschuster, T.; Handt, J.; Schmidt, R. Floquet surface hopping: Laser-driven dissociation and ionization dynamics of $H_2^+$. *Phys. Rev. A* **2016**, *93*, 053409. [CrossRef]

8. Schirò, M.; Eich, F.G.; Agostini, F. Quantum-classical nonadiabatic dynamics of Floquet driven systems. *J. Chem. Phys.* **2021**, *154*, 114101. [CrossRef]

9. Lépine, F.; Ivanov, M.Y.; Vrakking, M.J.J. Attosecond molecular dynamics: Fact or fiction? *Nat. Photonics* **2014**, *8*, 195. [CrossRef]

10. Hunter, G. Conditional probability amplitudes in wave mechanics. *Int. J. Quantum Chem.* **1975**, *9*, 237. [CrossRef]

11. Hunter, G. Ionization potentials and conditional amplitudes. *Int. J. Quantum Chem.* **1975**, *9*, 311. [CrossRef]

12. Hunter, G. Nodeless wave function quantum theory. *Int. J. Quantum Chem.* **1980**, *9*, 133. [CrossRef]

13. Hunter, G. Nodeless wave functions and spiky potentials. *Int. J. Quantum Chem.* **1981**, *19*, 755. [CrossRef]

14. Hunter, G.; Tai, C.C. Variational marginal amplitudes. *Int. J. Quantum Chem.* **1982**, *21*, 1041. [CrossRef]

15. Gidopoulos, N.I.; Gross, E.K.U. Electronic non-adiabatic states: Towards a density functional theory beyond the Born–Oppenheimer approximation. *Philos. Trans. R. Soc. Lond. A: Math. Phys. Eng. Sci.* **2014**, *372*, 20130059. [CrossRef]

16. Abedi, A.; Maitra, N.T.; Gross, E.K.U. Exact Factorization of the Time-Dependent Electron-Nuclear Wave Function. *Phys. Rev. Lett.* **2010**, *105*, 123002. [CrossRef]

17. Abedi, A.; Maitra, N.T.; Gross, E.K.U. Correlated electron-nuclear dynamics: Exact factorization of the molecular wavefunction. *J. Chem. Phys.* **120**, 137, 22A530. [CrossRef]

18. Agostini, F.; Gross, E.K.U. Ultrafast dynamics with the exact factorization. *Eur. Phys. J. B* **2021**, *94*, 179. [CrossRef]

19. Suzuki, Y.; Abedi, A.; Maitra, N.T.; Yamashita, K.; Gross, E.K.U. Electronic Schrödinger equation with nonclassical nuclei. *Phys. Rev. A* **2014**, *89*, 040501(R). [CrossRef]

20. Khosravi, E.; Abedi, A.; Maitra, N.T. Exact Potential Driving the Electron Dynamics in Enhanced Ionization of $H_2^+$. *Phys. Rev. Lett.* **2015**, *115*, 263002. [CrossRef]

21. Ha, J.K.; Lee, I.S.; Min, S.K. Surface Hopping Dynamics beyond Nonadiabatic Couplings for Quantum Coherence. *J. Phys. Chem. Lett.* **2018**, *9*, 1097–1104. [CrossRef]

22. Lee, I.S.; Ha, J.K.; Han, D.; Kim, T.I.; Moon, S.W.; Min, S.K. PyUNIxMD: A Python-based excited state molecular dynamics package. *J. Comput. Chem.* **2021**, *42*, 1755–1766. [CrossRef] [PubMed]

23. Pieroni, C.; Agostini, F. Nonadiabatic Dynamics with Coupled Trajectories. *J. Chem. Theory Comput.* **2021**, *17*, 5969–5991. [CrossRef] [PubMed]

24. Min, S.K.; Agostini, F.; Gross, E.K.U. Coupled-Trajectory Quantum-Classical Approach to Electronic Decoherence in Nonadiabatic Processes. *Phys. Rev. Lett.* **2015**, *115*, 073001. [CrossRef] [PubMed]

25. Agostini, F.; Min, S.K.; Abedi, A.; Gross, E.K.U. Quantum-Classical Nonadiabatic Dynamics: Coupled- vs Independent-Trajectory Methods. *J. Chem. Theory Comput.* **2016**, *12*, 2127–2143. [CrossRef]

26. Min, S.K.; Agostini, F.; Tavernelli, I.; Gross, E.K.U. Ab Initio Nonadiabatic Dynamics with Coupled Trajectories: A Rigorous Approach to Quantum (De)Coherence. *J. Phys. Chem. Lett.* **2017**, *8*, 3048–3055. [CrossRef]

27. Khosravi, E.; Abedi, A.; Rubio, A.; Maitra, N.T. Electronic non-adiabatic dynamics in enhanced ionization of isotopologues of hydrogen molecular ions from the exact factorization perspective. *Phys. Chem. Chem. Phys.* **2017**, *19*, 8269–8281. [CrossRef] [PubMed]
28. Abedi, A.; Maitra, N.T.; Gross, E.K.U. Response to "Comment on 'Correlated electron-nuclear dynamics: Exact factorization of the molecular wavefunction'

29. Fiedlschuster, T.; Handt, J.; Gross, E.K.U.; Schmidt, R. Surface hopping in laser-driven molecular dynamics. Phys. Rev. A 2017, 95, 063424. [CrossRef]

30. Suzuki, Y.; Abedi, A.; Maitra, N.T.; Gross, E.K.U. Laser-induced electron localization in H₂: Mixed quantum-classical dynamics based on the exact time-dependent potential energy surface. Phys. Chem. Chem. Phys. 2015, 17, 29271-29280. [CrossRef]

31. Abedi, A.; Agostini, F.; Suzuki, Y.; Gross, E.K.U. Dynamical steps that bridge piecewise adiabatic shapes in the exact time-dependent potential energy surface. Phys. Rev. Lett. 2013, 110, 263001. [CrossRef] [PubMed]

32. Agostini, F.; Abedi, A.; Suzuki, Y.; Gross, E.K.U. Mixed quantum-classical dynamics on the exact time-dependent potential energy surfaces: A novel perspective on non-adiabatic processes. Mol. Phys. 2013, 111, 3625. [CrossRef]

33. Agostini, F.; Abedi, A.; Suzuki, Y.; Min, S.K.; Maitra, N.T.; Gross, E.K.U. The exact forces on classical nuclei in non-adiabatic charge transfer. J. Chem. Phys. 2015, 142, 084303. [CrossRef] [PubMed]

34. Curchod, B.F.E.; Agostini, F. On the Dynamics through a Conical Intersection. J. Phys. Chem. Lett. 2017, 8, 831-837. [CrossRef]

35. Curchod, B.F.E.; Agostini, F.; Gross, E.K.U. An exact factorization perspective on quantum interferences in nonadiabatic dynamics. J. Chem. Phys. 2016, 145, 034103. [CrossRef] [PubMed]

36. Agostini, F.; Abedi, A.; Gross, E.K.U. Classical nuclear motion coupled to electronic non-adiabatic transitions. J. Chem. Phys. 2014, 141, 214101. [CrossRef] [PubMed]

37. Abedi, A.; Agostini, F.; Gross, E.K.U. Mixed quantum-classical dynamics from the exact decomposition of electron-nuclear motion. Europhys. Lett. 2014, 106, 33001. [CrossRef]

38. Davis, M.J.; Heller, E.J. Quantum dynamical tunneling in bound states. J. Chem. Phys. 1981, 75, 246. [CrossRef]

39. Hughes, K.H.; Parry, S.M.; Parlant, G.; Burghardt, I. A hybrid hydrodynamic-liouvillian approach to mixed quantum-classical dynamics: Application to tunneling in a double well. J. Phys. Chem. A 2007, 111, 10269-10283. [CrossRef]

40. Basire, M.; Borgis, D.; Vuilleumier, R. Computing Wigner distributions and time correlation functions using the quantum thermal bath method: Application to proton transfer spectroscopy. Phys. Chem. Chem. Phys. 2013, 15, 12591-12601. [CrossRef]

41. Litman, Y.; Behler, J.; Rossi, M. Temperature dependence of the vibrational spectrum of porphycene: A qualitative failure of classical-nuclei molecular dynamics. Faraday Discuss. 2020, 221, 526-546. [CrossRef] [PubMed]

42. Lawrence, J.E.; Manolopoulos, D.E. Path integral methods for reaction rates in complex systems. Faraday Discuss. 2020, 221, 9-29. [CrossRef] [PubMed]

43. Ghosh, S.; Giannini, S.; Lively, K.; Blumberger, J. Nonadiabatic dynamics with quantum nuclei: Simulating charge transfer with ring polymer surface hopping. Faraday Discuss. 2020, 221, 501-525. [CrossRef] [PubMed]

44. Gu, B.; Franco, I. Partial hydrodynamic representation of quantum molecular dynamics. J. Chem. Phys. 2017, 146, 194104. [CrossRef]

45. Shushkov, P.; Li, R.; Tully, J.C. Ring polymer molecular dynamics with surface hopping. J. Chem. Phys. 2012, 137, 22A549. [CrossRef]

46. Dupuy, L.; Lauvergnat, D.; Scribano, Y. Smolyak representations with absorbing boundary conditions for reaction path Hamiltonian model of reactive scattering. Chem. Phys. Lett. 2022, 787, 139241.: 10.1016/j.cplett.2021.139241. [CrossRef]

47. Suzuki, Y.; Watanabe, K. Bohmian mechanics in the exact factorization of electron-nuclear wave functions. Phys. Rev. A 2016, 94, 032517. [CrossRef]

48. Talotta, F.; Agostini, F.; Ciccotti, G. Quantum trajectories for the dynamics in the exact factorization framework: A proof-of-principle test. J. Phys. Chem. A 2020, 124, 6764-6777. [CrossRef]

49. Agostini, F.; Tavernelli, I.; Ciccotti, G. Nuclear Quantum Effects in Electronic (Non)Adiabatic Dynamics. Eur. Phys. J. B 2018, 91, 139. [CrossRef]

50. Lopreore, C.L.; Wyatt, R.E. Quantum Wave Packet Dynamics with Trajectories. Phys. Rev. Lett. 1999, 82, 5190. [CrossRef]

51. Wyatt, R.E. Quantum wavepacket dynamics with trajectories: Wavefunction synthesis along quantum paths. Chem. Phys. Lett. 1999, 313, 189-197. [CrossRef]

52. Wyatt, R.E.; Na, K. Quantum trajectory analysis of multimode subsystem-bath dynamics. Phys. Rev. E 2001, 65, 016702. [CrossRef] [PubMed]

53. Wyatt, R.E. Quantum Dynamics with Trajectories: Introduction to Quantum Hydrodynamics; Interdisciplinary Applied Mathematics; Springer: New York, NY, USA, 2005.

54. Garashchuk, S.; Rassolov, V. Quantum Trajectory Dynamics Based on Local Approximations to the Quantum Potential and Force. J. Chem. Theory Comput. 2019, 15, 3906-3916. [CrossRef] [PubMed]

55. Garashchuk, S.; Vazhappilly, T. Multidimensional Quantum Trajectory Dynamics in Imaginary Time with Approximate Quantum Potential. J. Phys. Chem. C 2010, 114, 20595-20602. [CrossRef]

56. Wyatt, R.E.; Bittner, E.R. Quantum wave packet dynamics with trajectories: Implementation with adaptive Lagrangian grids. J. Chem. Phys. 2000, 119, 8898. [CrossRef]

57. Hughes, K.H.; Wyatt, R.E. Wavepacket dynamics on dynamically adapting grids: Application of the equidistribution principle. Chem. Phys. Lett. 2002, 366, 336-342. [CrossRef]

58. Kendrick, B.K. A new method for solving the quantum hydrodynamic equations of motion. J. Chem. Phys. 2003, 119, 5805. [CrossRef]
59. Trahan, C.J.; Wyatt, R.E. An arbitrary Lagrangian-Eulerian approach to solving the quantum hydrodynamic equations of motion: Equidistribution with "smart" springs. *J. Chem. Phys.* **2003**, *4784*, 336–342. [CrossRef]

60. Schild, A. Electronic quantum trajectories with quantum nuclei. *arXiv* **2021**, arXiv:2109.13632v1.

61. Tavernelli, I.; Curchod, B.F.E.; Rothlisberger, U. Nonadiabatic molecular dynamics with solvent effects: A LR-TDDFT QM/MM study of ruthenium (II) tris (bipyridine) in water. *Chem. Phys.* **2011**, *391*, 101–109. [CrossRef]

62. Talotta, F.; Boggio-Pasqua, M.; González, L. Early relaxation dynamics in the photoswitchable trans-[RuCl(NO)(py)₄]²⁺. *Chem.: Eur. J.* **2020**, *26*, 11522–11528. [CrossRef] [PubMed]

63. Atkins, A.J.; Talotta, F.; Freitag, L.; Boggio-Pasqua, M.; González, L. Assessing Excited State Energy Gaps with Time-Dependent Density Functional Theory on Ru(II) Complexes. *J. Chem. Theory Comput.* **2017**, *13*, 4123–4145. [CrossRef] [PubMed]

64. García, J.S.; Talotta, F.; Alary, F.; Dixon, I.M.; Heully, J.L.; Boggio-Pasqua, M. A Theoretical Study of the N to O Linkage Photoisomerization Efficiency in a Series of Ruthenium Mononitrosyl Complexes. *Molecules* **2017**, *22*, 1667.

65. Ando, H.; Iuchi, S.; Sato, H. Theoretical study on ultrafast intersystem crossing of chromium(III) acetylacetonate. *Chem. Phys. Lett.* **2012**, *535*, 177–181. [CrossRef]

66. Brahim, H.; Daniel, C. Structural and spectroscopic properties of Ir(III) complexes with phenylpyridine ligands: Absorption spectra without and with spin-orbit-coupling. *Comput. Theo. Chem.* **2014**, *1040–1041*, 219–229. [CrossRef]

67. Hu, W.; Lendvay, G.; Maiti, B.; Schatz, G.C. Trajectory Surface Hopping Study of the O(³P) + Ethylene Reaction Dynamics. *J. Phys. Chem. A* **2008**, *112*, 2093–2103. [CrossRef]

68. Fu, B.; Han, Y.C.; Bowman, J.M.; Angelucci, L.; Balucani, N.; Leonori, F.; Casavecchia, P. Intersystem crossing and dynamics in O(³P)+C₂H₄ multichannel reaction: Experiment validates theory. *Proc. Natl. Acad. Sci. USA* **2012**, *109*, 9733–9738. [CrossRef]

69. Hu, W.; Lendvay, G.; Maiti, B.; Schatz, G.C. Electronic Structure and Excited States of the Collision Reaction O(³P)+C₂H₄: A Multiconfigurational Perspective. *J. Phys. Chem. A* **2021**, *125*, 6075–6088.

70. Talotta, F.; Morisset, S.; Rougeau, N.; Lauvergnat, D.; Agostini, F. Spin-Orbit Interactions in Ultrafast Molecular Processes. *Phys. Rev. Lett.* **2020**, *124*, 033001. [CrossRef]

71. Talotta, F.; Morisset, S.; Rougeau, N.; Lauvergnat, D.; Agostini, F. Internal Conversion and Intersystem Crossing with the Exact Factorization. *J. Chem. Theory Comput.* **2020**, *16*, 4833–4848. [CrossRef]

72. Min, S.K.; Abedi, A.; Kim, K.S.; Gross, E.K.U. Is the molecular Berry phase an artefact of the Born-Oppenheimer approximation? *Phys. Rev. Lett.* **2014**, *113*, 263004. [CrossRef] [PubMed]

73. Requist, R.; Tandetzky, F.; Gross, E.K.U. Molecular geometric phase from the exact electron-nuclear factorization. *Phys. Rev. A* **2016**, *93*, 042108. [CrossRef]

74. Requist, R.; Proetto, C.R.; Gross, E.K.U. Asymptotic analysis of the Berry curvature in the $E \otimes e$ Jahn-Teller model. *Phys. Rev. A* **2017**, *96*, 062503. [CrossRef]

75. Agostini, F.; Curchod, B.F.E. When the exact factorization meets conical intersections. *Eur. Phys. J. B* **2018**, *91*, 141. [CrossRef]

76. Ibele, L.M.; Curchod, B.F.E.; Agostini, F. A photochemical reaction in different theoretical representations. *J. Phys. Chem. A* **2022**, *126*, 1263–1281. [CrossRef]

77. Hader, K.; Albert, J.; Gross, E.K.U.; Engel, V. Electron-nuclear wave-packet dynamics through a conical intersection. *J. Chem. Phys.* **2017**, *146*, 074304. [CrossRef]

78. Requist, R.; Gross, E.K.U. Exact Factorization-Based Density Functional Theory of Electrons and Nuclei. *Phys. Rev. Lett.* **2016**, *117*, 193001. [CrossRef]

79. Li, C.; Requist, R.; Gross, E.K.U. Density functional theory of electron transfer beyond the Born-Oppenheimer approximation: Case study of LiF. *J. Chem. Phys.* **2018**, *148*, 084110. [CrossRef]

80. Schild, A.; Gross, E.K.U. Exact Single-Electron Approach to the Dynamics of Molecules in Strong Laser Fields. *Phys. Rev. Lett.* **2017**, *118*, 163202. [CrossRef]

81. Kocák, J.; Schild, A. Many-electron effects of strong-field ionization described in an exact one-electron theory. *Phys. Rev. Res.* **2020**, *2*, 043365. [CrossRef]

82. Gonze, X.; Zhou, J.S.; Reining, L. Variations on the "exact factorization" theme. *Eur. Phys. J. B* **2018**, *91*, 224. [CrossRef]

83. Lacombe, L.; Maitra, N.T. Embedding via the Exact Factorization Approach. *Phys. Rev. Lett.* **2020**, *124*, 206401. [CrossRef] [PubMed]

84. Requist, R.; Gross, E.K.U. Fock-Space Embedding Theory: Application to Strongly Correlated Topological Phases. *Phys. Rev. Lett.* **2021**, *127*, 116401. [CrossRef] [PubMed]

85. Salas, L.D.; Arce, J.C. Potential energy surfaces in atomic structure: The role of Coulomb correlation in the ground state of helium. *Phys. Rev. A* **2017**, *95*, 022502. [CrossRef]

86. Salas, L.D.; Zamora-Yusti, B.; Arce, J.C. Characterization of the continuous transition from atomic to molecular shape in the three-body Coulomb system. *Phys. Rev. A* **2022**, *105*, 012808. [CrossRef]

87. Eich, F.G.; Agostini, F. The adiabatic limit of the exact factorization of the electron-nuclear wave function. *J. Chem. Phys.* **2016**, *145*, 054110. [CrossRef]

88. Schild, A.; Agostini, F.; Gross, E.K.U. Electronic Flux Density beyond the Born-Oppenheimer Approximation. *J. Phys. Chem. A* **2016**, *120*, 3316. [CrossRef]

89. Scherrer, A.; Agostini, F.; Sebastiani, D.; Gross, E.K.U.; Vuilleumier, R. On the mass of atoms in molecules: Beyond the Born-Oppenheimer approximation. *Phys. Rev. X* **2017**, *7*, 031035. [CrossRef]
90. Scherrer, A.; Agostini, F.; Sebastiani, D.; Gross, E.K.U.; Vuilleumier, R. Nuclear velocity perturbation theory for vibrational circular dichroism: An approach based on the exact factorization of the electron-nuclear wave function. *J. Chem. Phys.* **2015**, *143*, 074106. [CrossRef]

91. Requist, R.; Proetto, C.R.; Gross, E.K.U. Exact factorization-based density functional theory of electron-phonon systems. *Phys. Rev. B* **2019**, *99*, 165136. [CrossRef]

92. Gossel, G.H.; Lacombe, L.; Maitra, N.T. On the numerical solution of the exact factorization equations. *J. Chem. Phys.* **2019**, *150*, 154112. [CrossRef] [PubMed]

93. Lorin, E. Numerical analysis of the exact factorization of molecular time-dependent Schrödinger wavefunctions. *Commun. Nonlinear Sci. Numer. Simul.* **2021**, *95*, 105627. [CrossRef]

94. Hoffmann, N.M.; Appel, H.; Rubio, A.; Maitra, N.T. Light-matter interactions via the exact factorization approach. *Eur. Phys. J. B* **2018**, *91*, 180. [CrossRef]

95. Yuen-Zhou, J.; Xiong, W.; Shegai, T. Polariton chemistry: Molecules in cavities and plasmonic media. *J. Chem. Phys.* **2022**, *156*, 030401. [CrossRef] [PubMed]

96. Galego, J.; Garcia-Vidal, F.J.; Feist, J. Cavity-Induced Modifications of Molecular Structure in the Strong-Coupling Regime. *Phys. Rev. X* **2015**, *5*, 041022. [CrossRef]

97. Lacombe, L.; Hoffmann, N.M.; Maitra, N.T. Exact Potential Energy Surface for Molecules in Cavities. *Phys. Rev. Lett.* **2019**, *123*, 083201. [CrossRef]

98. Martinez, P.; Rosenzweig, B.; Hoffmann, N.M.; Lacombe, L.; Maitra, N.T. Case studies of the time-dependent potential energy surface for dynamics in cavities. *J. Chem. Phys.* **2021**, *154*, 014102. [CrossRef]

99. Hoffmann, N.M.; Schäfer, C.; Rubio, A.; Kelly, A.; Appel, H. Capturing vacuum fluctuations and photon correlations in cavity quantum electrodynamics with multitrajectory Ehrenfest dynamics. *Phys. Rev. A* **2019**, *99*, 063819. [CrossRef]

100. Rosenzweig, B.; Hoffmann, N.M.; Lacombe, L.; Maitra, N.T. Analysis of the classical trajectory treatment of photon dynamics for polaritonic phenomena. *J. Chem. Phys.* **2022**, *156*, 054101. [CrossRef]

101. Flick, J.; Ruggenthaler, M.; Appel, H.; Rubio, A. Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry. *Proc. Natl. Acad. Sci. USA* **2017**, *114*, 3026-3034. [CrossRef]

102. Cederbaum, L.S. The exact wavefunction of interacting N degrees of freedom as a product of N single-degree-of-freedom wavefunctions. *Chem. Phys.* **2015**, *457*, 129. [CrossRef]

103. Scherrer, A.; Vuilleumier, R.; Sebastiani, D. Vibrational circular dichroism from ab initio molecular dynamics and nuclear velocity perturbation theory in the liquid phase. *J. Chem. Phys.* **2016**, *145*, 084101. [CrossRef] [PubMed]

104. Scherrer, A.; Vuilleumier, R.; Sebastiani, D. Nuclear Velocity Perturbation Theory of Vibrational Circular Dichroism. *J. Chem. Theory Comput.* **2013**, *9*, 5305-5312. [CrossRef]

105. Diestler, D.J.; Kenfack, A.; Manz, J.; Paulus, B.; Pérez-Torres, J.F.; Pohl, V. Computation of the Electronic Flux Density in the Born-Oppenheimer Approximation. *J. Phys. Chem. A* **2013**, *117*, 8519-8527. [CrossRef]

106. Diestler, D.J. Beyond the Born-Oppenheimer Approximation: A Treatment of Electronic Flux Density in Electronically Adiabatic Molecular Processes. *J. Phys. Chem. A* **2013**, *117*, 4698-4708. [CrossRef] [PubMed]

107. Arce, J.C. Unification of the conditional probability and semiclassical interpretations for the problem of time in quantum theory. *Phys. Rev. A* **2012**, *85*, 042108. [CrossRef]

108. Schild, A. Time in quantum mechanics: A fresh look at the continuity equation. *Phys. Rev. A* **2018**, *98*, 052113. [CrossRef]

109. Kulander, K.C.; Mies, F.H.; Schafer, K.J. Model for studies of laser-induced nonlinear processes in molecules. *Phys. Rev. A* **1996**, *53*, 2562-2570. [CrossRef]

110. Chelkowski, S.; Foisy, C.; Bandrauk, A.D. Electron-nuclear dynamics of multiphoton $H_2^+$ dissociative ionization in intense laser fields. *Phys. Rev. A* **1998**, *57*, 1176-1185. [CrossRef]

111. Walsh, T.D.G.; Ilkov, F.A.; Chin, S.L.; Châteauneuf, F.; Nguyen-Dang, T.T.; Chelkowski, S.; Bandrauk, A.D.; Atabek, O. Laser-induced processes during the Coulomb explosion of $H_2$ in a Ti-sapphire laser pulse. *Phys. Rev. A* **1998**, *58*, 3922-3933. [CrossRef]

112. Lein, M.; Kreibich, T.; Gross, E.K.U.; Engel, V. Strong-field ionization dynamics of a model $H_2$ molecule. *Phys. Rev. A* **2002**, *65*, 033403. [CrossRef]

113. McLachlan, A.D. A variational solution of the time-dependent Schrödinger equation. *Mol. Phys.* **1964**, *8*, 39-44. [CrossRef]

114. Tully, J.C. Mixed quantum-classical dynamics. *Faraday Discuss.* **1998**, *110*, 407. [CrossRef]

115. Kelkensberg, F.; Sansone, G.; Ivanov, M.Y.; Vrakking, M. A semi-classical model of attosecond electron localization in dissociative ionization of hydrogen. *Phys. Chem. Chem. Phys.* **2011**, *13*, 8647-8652. [CrossRef]

116. Sansone, G.; Kelkensberg, F.; Pérez-Torres, J.F.; Morales, F.; Kling, M.F.; Siu, W.; Ghafur, O.; Johnsson, P.; Swoboda, M.; Benedetti, E.; et al. Electron localization following attosecond molecular photoionization. *Nature* **2010**, *465*, 763-766. [CrossRef]

117. He, F.; Ruiz, C.; Becker, A. Control of Electron Excitation and Localization in the Dissociation of $H_2^+$ and Its Isotopes Using Two Sequential Ultrashort Laser Pulses. *Phys. Rev. Lett.* **2007**, *99*, 083002. [CrossRef]

118. Zuo, T.; Bandrauk, A.D. Charge-resonance-enhanced ionization of diatomic molecular ions by intense lasers. *Phys. Rev. A* **1995**, *52*, R2511. [CrossRef]

119. Chelkowski, S.; Bandrauk, A.D. Two-step Coulomb explosions of diatoms in intense laser fields. *J. Phys. B: At. Mol. Opt. Phys.* **1995**, *28*, L723-L731. [CrossRef]
120. Chelkowski, S.; Zuo, T.; Atabek, O.; Bandrauk, A.D. Dissociation, ionization, and Coulomb explosion of H₂⁺ in an intense laser field by numerical integration of the time-dependent Schrödinger equation. *Phys. Rev. A* **1995**, *52*, 2977. [CrossRef]

121. Seideman, T.; Ivanov, M.Y.; Corkum, P.B. Role of Electron Localization in Intense-Field Molecular Ionization. *Phys. Rev. Lett.* **1995**, *75*, 2819–2822. [CrossRef]

122. Bandrauk, A.D.; Légaré, F. Enhanced Ionization of Molecules in Intense Laser Fields. In *Progress in Ultrafast Intense Laser Science VIII*; Yamanouchi, K., Nisoli, M., Hill, W.T., Eds.; Springer: Berlin/Heidelberg, Germany, 2012; pp. 29–46. [CrossRef]

123. Zuo, T.; Chelkowski, S.; Bandrauk, A.D. Harmonic generation by the H₂⁺ molecular ion in intense laser fields. *Phys. Rev. A* **1993**, *48*, 3837–3844. [CrossRef] [PubMed]

124. Takemoto, N.; Becker, A. Multiple Ionization Bursts in Laser-Driven Hydrogen Molecular Ion. *Phys. Rev. Lett.* **2010**, *105*, 203004. [CrossRef] [PubMed]

125. Takemoto, N.; Becker, A. Time-resolved view on charge-resonance-enhanced ionization. *Phys. Rev. A* **2011**, *84*, 023401. [CrossRef]

126. Beylerian, C.; Saugout, S.; Cornaggia, C. Non-sequential double ionization of H₂ using ultrashort 10 fs laser pulses. *J. Phys. B: At. Mol. Opt. Phys.* **2006**, *39*, L105–L112. [CrossRef]

127. Bocharova, I.; Karimi, R.; Penka, E.F.; Brichta, J.P.; Lassonde, P.; Fu, X.; Kieffer, J.C.; Bandrauk, A.D.; Litvinyuk, I.; Sanderson, J.; et al. Charge Resonance Enhanced Ionization of CO₂ Probed by Laser Coulomb Explosion Imaging. *Phys. Rev. Lett.* **2011**, *107*, 063201. [CrossRef]

128. Légaré, F.; Litvinyuk, I.V.; Dooley, P.W.; Quéré, F.; Bandrauk, A.D.; Villeneuve, D.M.; Corkum, P.B. Time-Resolved Double Ionization with Few Cycle Laser Pulses. *Phys. Rev. Lett.* **2003**, *91*, 093002. [CrossRef]

129. Curchod, B.F.E.; Agostini, F.; Tavernelli, I. CT-MQC—A Coupled-Trajectory Mixed Quantum/Classical method including nonadiabatic quantum coherence effects. *Eur. Phys. J. B* **2018**, *91*, 168. [CrossRef]

130. Marsili, E.; Olivucci, M.; Lauvergnat, D.; Agostini, F. Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model. *J. Chem. Theory Comput.* **2020**, *16*, 6032–6048. [CrossRef]

131. Filatov, M.; Min, S.K.; Kim, K.S. Non-adiabatic dynamics of ring opening in cyclohexa-1,3-diene described by an ensemble density-functional theory method. *Mol. Phys.* **2019**, *117*, 1128–1141. [CrossRef]

132. Filatov, M.; Paolino, M.; Min, S.K.; Kim, K.S. Fulgides as Light-Driven Molecular Rotary Motors: Computational Design of a Prototype Compound. *J. Phys. Chem. Lett.* **2018**, *9*, 4995–5001. [CrossRef]

133. Filatov, M.; Paolino, M.; Min, S.K.; Choi, C.H. Design and photoisomerization dynamics of a new family of synthetic 2-stroke light driven molecular rotary motors. *Chem. Commun.* **2019**, *55*, 5247–5250. [CrossRef] [PubMed]

134. Filatov, M.; Min, S.K.; Choi, C.H. Theoretical modelling of the dynamics of primary photoprocess of cyclopropanone. *Phys. Chem. Chem. Phys.* **2019**, *21*, 2489–2498. [CrossRef] [PubMed]

135. Vindel-Zandbergen, P.; Ibele, L.M.; Ha, J.K.; Min, S.K.; Curchod, B.F.E.; Maitra, N.T. Study of the Decoherence Correction Derived from the Exact Factorization Approach for Nonadiabatic Dynamics. *J. Chem. Theory Comput.* **2021**, *17*, 3852–3862. [CrossRef] [PubMed]

136. Vindel-Zandbergen, P.; Matsika, S.; Maitra, N.T. Exact-Factorization-Based Surface Hopping for Multistate Dynamics. *J. Phys. Chem. Lett.* **2022**, *13*, 1785–1790. [CrossRef]

137. Tully, J.C. Molecular dynamics with electronic transitions. *J. Chem. Phys.* **1990**, *93*, 1061. [CrossRef]

138. Lu, J.; Zhoud, Z. Frozen Gaussian Approximation with surface hopping for mixed quantum-classical dynamics: A mathematical justification of fewest switches surface hopping algorithms. *Math. Comp.* **2018**, *87*, 2189–2232. [CrossRef]

139. Wang, L.; Akimov, A.; Prezhdo, O.V. Recent Progress in Surface Hopping: 2011–2015. *J. Phys. Chem. Lett.* **2016**, *7*, 2100. [CrossRef]

140. Subotnik, J.E.; Jain, A.; Landry, B.; Petit, A.; Ouyang, W.; Bellonzi, N. Understanding the Surface Hopping View of Electronic Transitions and Decoherence. *Ann. Rev. Phys. Chem.* **2016**, *67*, 387–417. [CrossRef]

141. Gossel, G.H.; Agostini, F.; Maitra, N.T. Coupled-Trajectory Mixed Quantum-Classical Algorithm: A Deconstruction. *J. Chem. Theory Comput.* **2018**, *14*, 4513–4529. [CrossRef]

142. Agostini, F. An exact-factorization perspective on quantum-classical approaches to excited-state dynamics. *Eur. Phys. J. B* **2018**, *91*, 143. [CrossRef]

143. Agostini, F.; Marsili, E.; Talotta, F. G-CTMQC. 2021. Available online: https://gitlab.com/agostini.work/g-ctmqc (accessed on 13 May 2022).

144. Lauvergnat, D. ModelLib. 2018. Available online: https://github.com/lauvergn/QuantumModelLib/tree/OOP2022).

145. Kim, T.I.; Ha, J.K.; Min, S.K. Coupled- and Independent-Trajectory Approaches Based on the Exact Factorization Using the PyUNIxMD Package. *Top. Curr. Chem.* **2022**, *380*, 8. [CrossRef] [PubMed]

146. Ha, J.K.; Min, S.K. Independent Trajectory Mixed Quantum-Classical Approaches Based on the Exact Factorization. *J. Chem. Phys.* **2022**, *156*, 174109. [CrossRef] [PubMed]

147. Barbatti, M. Velocity Adjustment in Surface Hopping: Ethylene as a Case Study of the Maximum Error Caused by Direction Choice. *J. Chem. Theor. Comput.* **2021**, *17*, 3010–3018. [CrossRef] [PubMed]

148. Carof, A.; Giannini, S.; Blumberger, J. Detailed balance, internal consistency, and energy conservation in fragment orbital-based surface hopping. *J. Chem. Phys.* **2017**, *147*, 214113. [CrossRef]
149. Tang, D.; Shen, L.; Fang, W.h. Evaluation of Mixed Quantum-Classical Molecular Dynamics on cis-Azobenzene Photoisomerization. *Phys. Chem. Chem. Phys.* **2021**, *23*, 13951–13964. [CrossRef]

150. Hammes-Schiffer, S.; Tully, J.C. Proton transfer in solution: Molecular dynamics with quantum transitions. *J. Chem. Phys.* **1994**, *101*, 4657–4667. [CrossRef]