# A Photochemical Reaction in Different Theoretical Representations

Lea M. Ibele, Basile F. E. Curchod,* and Federica Agostini*

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

Cite This: *J. Phys. Chem. A* 2022, 126, 1263–1281

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

Read Online

## ACCESS

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

Metrics & More

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

Article Recommendations

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

Supporting Information

**ABSTRACT:** The Born–Oppenheimer picture has forged our representation and interpretation of photochemical processes, from photoexcitation down to the passage through a conical intersection, a funnel connecting different electronic states. In this work, we analyze a full in silico photochemical experiment, from the explicit electronic excitation by a laser pulse to the formation of photoproducts following a nonradiative decay through a conical intersection, by contrasting the picture offered by Born–Oppenheimer and that proposed by the exact factorization. The exact factorization offers an alternative understanding of photochemistry that does not rely on concepts such as electronic states, nonadiabatic couplings, and conical intersections. On the basis of nonadiabatic quantum dynamics performed for a two-state 2D model system, this work allows us to compare Born–Oppenheimer and exact factorization for (i) an explicit photoexcitation with and without the Condon approximation, (ii) the passage of a nuclear wavepacket through a conical intersection, (iii) the formation of excited stationary states in the Franck–Condon region, and (iv) the use of classical and quantum trajectories in the exact factorization picture to capture nonadiabatic processes triggered by a laser pulse.

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

## 1. INTRODUCTION

Our way of picturing molecules and chemical processes has been greatly shaped by the Born–Oppenheimer approximation, the assumption that the motion of electrons and nuclei can be treated separately in a molecule.¹,² The discussion of chemical structures, properties, and reactivity usually intrinsically assumes that the molecule remains in a given electronic eigenstate or, in other words, that electrons can adapt instantaneously to any nuclear motion, which is a direct consequence of the Born–Oppenheimer approximation. Photochemistry inherently goes beyond this picture because photon absorption by a molecule makes transitions to different electronic states possible.³⁻⁹ The subsequent out-of-equilibrium evolution of the excited molecule usually means that regions where one or more electronic states come close in energy can be visited. When the nuclear dynamics drives the molecule into such regions, the Born–Oppenheimer approximation breaks down. The *nonadiabatic* coupling between nuclear and electronic motions cannot be neglected anymore, and one needs to accept that nuclear motion can lead to a change of electronic states.¹⁰

A legitimate strategy to move beyond the Born–Oppenheimer approximation would be to try to supplement the Born–Oppenheimer picture with more electronic eigenstates and account for their mutual couplings mediated by nuclear motion. Importantly, this post-Born–Oppenheimer picture, often called the Born–Huang representation,¹¹ relies on the use of quantities originally defined within the framework of the Born–Oppenheimer approximation, such

as the potential energy surfaces obtained from the electronic Schrödinger equation for fixed nuclear degrees of freedom, in a context that is at the opposite of their initial raison d'être: strong electron–nuclear couplings.¹² This post-Born–Oppenheimer picture is at the heart of our way to regard photochemical processes, and the vocabulary used for such processes is intrinsically shaped by Born–Oppenheimer concepts: potential energy surfaces (eigenvalues from the time-independent Schrödinger equation),⁴ conical intersections (points of degeneracy between two adiabatic potential energy surfaces),¹³,¹⁴ the Berry phase (a phase picked up in the adiabatic representation when encircling a conical intersection),¹⁵,¹⁶ and the transition dipole moment (the electric dipole moment for the transition between two electronic states). Within this picture, the molecular wave function is expanded in the adiabatic basis, that is, using the set of time-independent eigenstates of the (electronic) Born–Oppenheimer Hamiltonian. An alternative choice for this basis expansion is to use the diabatic representation, which diagonalizes, when possible, the nuclear kinetic energy operator. In the diabatic representation, the electronic states preserve their character upon varying the nuclear config-

Received: November 6, 2021

Revised: January 28, 2022

Published: February 14, 2022

<figure><img src="image_7.png" /></figure>
uration. Conical intersections and the Berry phase are not present in this representation because the electron-nuclear coupling is mediated here via off-diagonal elements of the electronic Hamiltonian matrix. For systems with more than two electronic states, where diabatic states might not be defined, one needs to rely on quasi-diabatic states that correspond to a set of electronic states that minimize the nonadiabatic couplings.

An exact theoretical description of a photochemical process necessitates the solution of the full time-dependent Schrödinger equation, a task that is possible only for the smallest molecular systems. In the vast majority of cases, the in silico simulation of excited-state (or nonadiabatic) molecular dynamics is performed under the (post-) Born-Oppenheimer picture. In the simplest case, the molecule of interest is considered to be described by its ground vibrational and electronic states, an external time-dependent field (for example, a laser pulse) is applied to the molecule, and the light-matter interaction induces an electronic excitation leading to the formation of a nuclear wavepacket in an excited electronic state. This nuclear wavepacket subsequently evolves and can branch in regions of strong nonadiabaticity, leading to its (possibly partial) transfer to different electronic states, a nonradiative decay, until the ground state is potentially reached. The out-of-equilibrium dynamics experienced by the molecule during this overall process can make it explore regions of the nuclear configuration space that would not be thermally accessible from the ground electronic state, leading to the formation of so-called photoproducts.<sup>17-21</sup> The formation of these photoproducts can be triggered in the excited electronic states or during subsequent athermal dynamics in the ground electronic state following the nonradiative decay.<sup>22,23</sup> Accessing the photochemistry of molecules in their full dimensionality usually requires approximations for the dynamics, which are in many cases based on representing the nuclear wave functions through trajectory basis functions or substituting it with swarms of classical trajectories. (See refs 24-26 for examples of such strategies.) A consequence of the Born-Oppenheimer picture is that one needs to find strategies to account for nonadiabatic effects in these approximate methods and, more specifically, to describe the branching of nuclear wavepackets between different electronic states. A plethora of ideas and algorithms were proposed, where coupled or uncoupled trajectories would hop, evolve on mean-field surfaces, spawn, or clone other trajectories.<sup>27</sup>

In this work, we propose to explore the different steps of a typical photochemical experiment described above with an alternative representation of the coupled electron-nuclear dynamics of a molecule: the exact factorization.<sup>28,29</sup> The exact factorization introduces a framework exempt from the ideas emanating from the Born-Oppenheimer picture; there is no mention of electronic states or static potential energy surfaces. Instead, the exact factorization depicts the dynamics of a molecular system by a nuclear wave function whose dynamics is dictated by a single time-dependent vector<sup>30</sup> and scalar<sup>31</sup> potential. Earlier work<sup>30,32-34</sup> showed that this formalism changes our way of representing nonadiabatic molecular dynamics. In particular, the concepts of the transition between electronic states, conical intersections, and the (topological) Berry phase do not appear in the exact factorization. Hence, there is a real curiosity in unraveling how the exact factorization would describe a full photochemical experiment

for a 2D two-state molecular model, from photoexcitation with a laser pulse to the formation of photoproducts, and in comparing this picture to the more conventional Born-Oppenheimer representation. In addition, these simulations will allow us to shed light on other interesting aspects of an in silico photochemical experiment, such as (i) the effect of the Condon approximation, (ii) the analysis of the dynamics using representation-free quantities, and (iii) the use of classical and quantum trajectories to depict the entire nuclear dynamics during a photochemical process.

This article is organized as follows. We propose in section 2 a brief review of the Born-Oppenheimer (section 2.1) and exact-factorization (section 2.2) pictures to describe coupled electron-nuclear dynamics such as that observed in a photochemical process. We then define our two-state 2D model Hamiltonian and the light-matter interaction Hamiltonian with and without the Condon approximation in section 2.3. In the same section, we also highlight some important considerations on our model and provide the computational details. We then present the results of our in silico photochemical experiment in section 3, starting with the more conventional Born-Oppenheimer picture and then moving to the exact factorization. We also discuss the use of classical and quantum trajectories within the exact factorization formalism. Our conclusions are finally stated in section 4.

## 2. METHODS

The evolution of the time-dependent wave function describing the state of a molecule, $\Psi(\mathbf{r}, \mathbf{R}, t)$, follows the time-dependent Schrödinger equation

$$i\hbar \frac{\partial}{\partial t} \Psi(\mathbf{r}, \mathbf{R}, t) = \hat{H}(\mathbf{r}, \mathbf{R}, t) \Psi(\mathbf{r}, \mathbf{R}, t) \quad (1)$$

where the full Hamiltonian operator $\hat{H}(\mathbf{r}, \mathbf{R}, t)$ includes the nuclear kinetic energy $\hat{T}_n(\mathbf{R})$, a Born-Oppenheimer Hamiltonian $\hat{H}_{BO}(\mathbf{r}, \mathbf{R})$, and an external (time-dependent) potential $\hat{V}(\mathbf{r}, \mathbf{R}, t)$:

$$ \begin{aligned} \hat{H}(\mathbf{r}, \mathbf{R}, t) &= \hat{T}_n(\mathbf{R}) + \hat{H}_{BO}(\mathbf{r}, \mathbf{R}) + \hat{V}(\mathbf{r}, \mathbf{R}, t) \\ &= \sum_{\nu}^{N_n} \frac{-\hbar^2}{2M_{\nu}} \nabla_{\mathbf{R}_{\nu}}^2 + \hat{T}_e(\mathbf{r}) + \hat{V}_{ee}(\mathbf{r}) + \hat{V}_{nn}(\mathbf{R}) + \hat{V}_{en}(\mathbf{r}, \mathbf{R}) \\ &\quad - \hat{\mu}(\mathbf{r}, \mathbf{R}) \cdot \mathbf{E}(t) \end{aligned} \quad (2) $$

$\mathbf{R}$ and $\mathbf{r}$ are collective variables for the coordinates of the nuclei and electrons forming the molecule. The summation for the nuclear kinetic energy runs over all $N_n$ nuclei, each labeled by $\nu$ and with corresponding mass $M_{\nu}$. The Born-Oppenheimer Hamiltonian consists of the electronic kinetic energy $\hat{T}_e(\mathbf{r})$, the electron-electron $\hat{V}_{ee}(\mathbf{r})$, the nucleus-nucleus $\hat{V}_{nn}(\mathbf{R})$, and the electron-nucleus $\hat{V}_{en}(\mathbf{r}, \mathbf{R})$ interaction. In the dipole approximation, the effect of an external time-dependent electric field $\mathbf{E}(t)$ coupled to the electric dipole moment operator $\hat{\mu}(\mathbf{r}, \mathbf{R})$ is encoded in $\hat{V}(\mathbf{r}, \mathbf{R}, t)$.

### 2.1. Born-Oppenheimer Picture of Nonadiabatic Dynamics
Within the Born-Oppenheimer approximation, the molecular wave function $\Psi(\mathbf{r}, \mathbf{R}, t)$ is approximated as a single product of a time-independent electronic wave function, $\Phi_j(\mathbf{r}; \mathbf{R})$, that is the $j$th solution of the time-independent electronic Schrödinger equation and a corresponding time-dependent nuclear wave function, $\chi_j(\mathbf{R}, t)$:
$$ \Psi(\mathbf{r}, \mathbf{R}, t) \approx \Psi_{\text{BO}}(\mathbf{r}, \mathbf{R}, t) = \Phi_j(\mathbf{r}; \mathbf{R}) \chi_j(\mathbf{R}, t) \quad (3) $$

This approximation restricts the dynamics of the nuclear wave function $\chi_j(\mathbf{R}, t)$ to the adiabatic electronic state $J$. In other words, the nuclear wave function evolves adiabatically in a given electronic state, and nothing can trigger a change in the electronic state with the Born–Oppenheimer approximation. To move away from an adiabatic representation and account for nonadiabatic effects, one needs to move to a different representation of the molecular wave function called the Born–Huang expansion. In the Born–Huang representation, the molecular wave function is written as

$$ \Psi(\mathbf{r}, \mathbf{R}, t) = \sum_{J}^{\infty} \Phi_j(\mathbf{r}; \mathbf{R}) \chi_j(\mathbf{R}, t) \quad (4) $$

Hence, the Born–Oppenheimer approximation has been corrected via an infinite sum over all electronic states $J$ leading to an in principle exact expansion of the molecular wave function. By inserting this expression into the molecular Schrödinger equation (eq 1), we can (after left multiplication by $\Phi_I(\mathbf{r}; \mathbf{R})$ and integration over all electronic coordinates $\mathbf{r}$) obtain a set of equations of motion for the nuclear amplitudes ($\chi_I(\mathbf{R}, t)$). (We assume a general form for the nuclear amplitudes that could include any phases picked up during the dynamics as a function of $\mathbf{R}$.)

$$ \begin{aligned} i\hbar \frac{\partial}{\partial t} \chi_I(\mathbf{R}, t) = & \left( \sum_{\nu}^{N_n} \frac{-\hbar^2}{2M_\nu} \nabla_{\mathbf{R}_\nu}^2 + \epsilon_{\text{BO}}^{(I)}(\mathbf{R}) \right) \chi_I(\mathbf{R}, t) \\ & + \sum_J^{\infty} \left( \sum_{\nu}^{N_n} \frac{\hbar^2}{2M_\nu} D_{IJ}(\mathbf{R}) - \sum_{\nu}^{N_n} \frac{\hbar^2}{M_\nu} \mathbf{d}_{IJ}(\mathbf{R}) \cdot \nabla_{\mathbf{R}_\nu} \right. \\ & \left. - \mu_{IJ}(\mathbf{R}) \cdot \mathbf{E}(t) \right) \chi_J(\mathbf{R}, t) \end{aligned} \quad (5) $$

The first two terms on the right-hand side denote the nuclear kinetic and potential energy associated with the adiabatic dynamics of the nuclear wave function in electronic state $I$. We can picture the nuclear amplitude associated with state $I$ to be evolving on the (static) potential energy surface $\epsilon_{\text{BO}}^{(I)}(\mathbf{R}) = \langle \Phi_I(\mathbf{R}) | \hat{H}_{\text{BO}}(\mathbf{R}) | \Phi_I(\mathbf{R}) \rangle_r$. We note that $\langle \cdots \rangle_r$ indicates integration over all electronic coordinates $\mathbf{r}$. The three last terms on the right-hand side of eq 5, neglected within the Born–Oppenheimer adiabatic approximation, are responsible for the possible transfer of nuclear amplitude between electronic states: $D_{IJ}(\mathbf{R}) = \langle \Phi_I(\mathbf{R}) | \nabla_{\mathbf{R}}^2 | \Phi_J(\mathbf{R}) \rangle_r$ represents the second-order nonadiabatic couplings, $d_{IJ}(\mathbf{R}) = \langle \Phi_I(\mathbf{R}) | \nabla_{\mathbf{R}} | \Phi_J(\mathbf{R}) \rangle_r$ represents the first-order nonadiabatic coupling vectors, and the final term determines the light–matter interaction with the external time-dependent electric field within the dipole approximation. The (de)excitation process by an electric field $\mathbf{E}(t)$ is mediated by the transition dipole moment $\mu_{IJ}(\mathbf{R}) = \langle \Phi_I(\mathbf{R}) | \hat{\mu}(\mathbf{R}) | \Phi_J(\mathbf{R}) \rangle_r$ of the molecule. We note that $d_{IJ}(\mathbf{R}) = 0$ for real electronic wave functions, and $D_{IJ}(\mathbf{R})$ and $\mu_{IJ}(\mathbf{R})$ are not necessarily zero.

## 2.2. Exact Factorization Picture of Nonadiabatic Dynamics.
Instead of expressing the molecular wave function as an infinite sum over time-independent electronic eigenfunctions and time-dependent nuclear amplitudes, the exact factorization proposes a subtle alternative:

$$ \Psi(\mathbf{r}, \mathbf{R}, t) = \Phi(\mathbf{r}, t; \mathbf{R}) \chi(\mathbf{R}, t) \quad (6) $$

This exact representation of the molecular wave function relies on a single product composed of a nuclear wave function $\chi(\mathbf{R}, t)$ and an electronic wave function $\Phi(\mathbf{r}, t; \mathbf{R})$, both time-dependent. Importantly, the electronic wave function does still depend parametrically on $\mathbf{R}$. Using the partial normalization condition, $\int d\mathbf{r}|\Phi(\mathbf{r}, t; \mathbf{R})|^2 = 1$ $\forall \mathbf{R}, t$, ensures that $|\chi(\mathbf{R}, t)|^2$ will reproduce the nuclear density as obtained from $\Psi(\mathbf{r}, \mathbf{R}, t)$ at all times. Both $\chi(\mathbf{R}, t)$ and $\Phi(\mathbf{r}, t; \mathbf{R})$ are uniquely defined up to a phase factor of $\exp[(i/\hbar)\theta(\mathbf{R}, t)]$ for any choice of the gauge function $\theta(\mathbf{R}, t)$. In general, $\theta(\mathbf{R}, t)$ does not have any physical meaning in the sense that observable properties do not depend on $\theta(\mathbf{R}, t)$. Note that $\Psi(\mathbf{r}, \mathbf{R}, t) = \chi(\mathbf{R}, t) \Phi(\mathbf{r}, t; \mathbf{R}) = \tilde{\chi}(\mathbf{R}, t) \tilde{\Phi}(\mathbf{r}, t; \mathbf{R})$ if one defines $\tilde{\chi}(\mathbf{R}, t) = \exp[-(i/\hbar)\theta(\mathbf{R}, t)]\chi(\mathbf{R}, t)$ and $\tilde{\Phi}(\mathbf{r}, t; \mathbf{R}) = \exp[(i/\hbar)\theta(\mathbf{R}, t)]\Phi(\mathbf{r}, t; \mathbf{R})$. The partial normalization condition also holds for the gauge-transformed electronic wave function. Once this gauge freedom is fixed by imposing a choice of phase, eq 6 is unique.

Because now both electronic and nuclear wave functions depend on time, we get (after insertion into the time-dependent Schrödinger equation) a set of coupled equations of motion:

$$ \begin{aligned} i\hbar \frac{\partial}{\partial t} \chi(\mathbf{R}, t) = & \left[ \sum_{\nu}^{N_n} \frac{[-i\hbar \nabla_{\mathbf{R}_\nu} + \mathbf{A}_\nu(\mathbf{R}, t)]^2}{2M_\nu} + \epsilon(\mathbf{R}, t) \right. \\ & \left. + v_{\text{int}}(\mathbf{R}, t) \right] \chi(\mathbf{R}, t) \\ & + i\hbar \frac{\partial}{\partial t} \Phi(\mathbf{r}, t; \mathbf{R}) = [\hat{H}_{\text{BO}}(\mathbf{r}, \mathbf{R}) + \hat{V}(\mathbf{r}, \mathbf{R}, t) \\ & + \hat{U}_{\text{en}}[\Phi, \chi](\mathbf{R}, t) - \epsilon(\mathbf{R}, t) - v_{\text{int}}(\mathbf{R}, t)] \Phi(\mathbf{r}, t; \mathbf{R}) \end{aligned} \quad (7) $$

These equations of motion introduce two new potentials within the framework of the exact factorization: the time-dependent vector potential, $\mathbf{A}_\nu(\mathbf{R}, t) = \langle \Phi(\mathbf{r}; \mathbf{R}) | -i\hbar \nabla_{\mathbf{R}_\nu} \Phi(\mathbf{r}; \mathbf{R}) \rangle_r$, and the time-dependent potential energy surface consisting of the two scalar potentials, $\epsilon(\mathbf{R}, t) = \langle \Phi(\mathbf{r}; \mathbf{R}) | \hat{H}_{\text{BO}}(\mathbf{R}) + \hat{U}_{\text{en}}[\Phi, \chi](\mathbf{R}, t) - i\hbar \partial_t | \Phi(\mathbf{r}; \mathbf{R}) \rangle_r$ and $v_{\text{int}}(\mathbf{R}, t) = \langle \Phi(\mathbf{r}; \mathbf{R}) | \hat{V}(\mathbf{R}, t) | \Phi(\mathbf{r}; \mathbf{R}) \rangle_r$. The electron–nuclear coupling operator, $\hat{U}_{\text{en}}[\Phi, \chi](\mathbf{R}, t)$, is explicitly dependent on the nuclear wave function $\chi(\mathbf{R}, t)$ and through the time-dependent vector potential, $\mathbf{A}_\nu(\mathbf{R}, t)$, also implicitly contains the electronic wave function $\Phi(\mathbf{r}, t; \mathbf{R})$:

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

It is worth mentioning here that the nuclear momentum field computed from the molecular wave function $\Psi(\mathbf{r}, \mathbf{R}, t)$ can be decomposed as the sum of a curl-free contribution, which is related to the phase of the nuclear wave function $\chi(\mathbf{R}, t)$ and (in general) a non-irrotational contribution, which is given by the time-dependent vector potential $\mathbf{A}_\nu(\mathbf{R}, t)$. The relationship between the nuclear momentum field and the vector potential will be used in Appendix B, where we illustrate the procedure employed to introduce a trajectory-based description of the nuclear dynamics within the framework of the exact factorization. In particular, we clarify the difference between
quantum$^{35–37}$ and classical trajectories$^{38,39}$ used in the numerical studies proposed in section 3.2.

## 2.3. Presentation of the Two-State Two-Dimensional Model. 2.3.1. Computational Details.
We used for this study a 2D two-state molecular model. In the diabatic representation, the general form of the Hamiltonian is given by

$$ \mathbf{H}^{\text{(dia)}}(\mathbf{R}) = \hat{T}_{\text{n}}(\mathbf{R}) + \begin{pmatrix} V_{11}(\mathbf{R}) & V_{12}(\mathbf{R}) \\ V_{12}(\mathbf{R}) & V_{22}(\mathbf{R}) \end{pmatrix} \quad (9) $$

with the following diabatic electronic energies

$$ \begin{aligned} V_{11}(\mathbf{R}) &= \frac{k_X}{2}(X - X_1)^2 + \frac{k_Y}{2}Y^2 \\ V_{22}(\mathbf{R}) &= \frac{k_X}{2}(X - X_2)^2 + \frac{k_Y}{2}Y^2 + \Delta \\ V_{12}(\mathbf{R}) &= \gamma Y \exp(-\alpha(X - X_3)^2) \exp(-\beta Y^2) \end{aligned} \quad (10) $$

using $\mathbf{R} = (X, Y)$. The parameters $k_X = 0.02 E_h a_0^{-2}$, $k_Y = 0.1 E_h a_0^{-2}$, $\Delta = 0.01 E_h$, $\gamma = 0.01 E_h a_0^{-1}$, $\alpha = 3a_0^{-2}$, $\beta = 1.5a_0^{-2}$, $M_X = 20000.0m_e$, and $M_Y = 6667.0m_e$ are based on refs 34 and 40, and we set here $X_1 = 6a_0$, $X_2 = 2a_0$, and $X_3 = 3.875a_0$. The two states are taken to be within the singlet manifold, so henceforth the adiabatic states will be labeled $S_0$ (ground electronic state) and $S_1$ (first excited electronic state).

We propagated nuclear wavepackets with numerically exact quantum dynamics, using for the initial state a Gaussian function with widths $\sigma_X = 0.15a_0$ and $\sigma_Y = 0.197a_0$ initialized in the adiabatic ground electronic state with zero initial nuclear momentum. The nuclear wave function is initially positioned at the Franck–Condon point $\mathbf{R}_{\text{init}} = (2.0, 0.0)a_0$. We note that because of the weak anharmonicity, the Franck–Condon point in the adiabatic representation coincides with the minimum of $V_{22}(\mathbf{R})$.

The external time-dependent electric field of a laser pulse under study here is given by

$$ \mathbf{E}(t) = \boldsymbol{\epsilon}^{\lambda} B_{0} \exp\left(-\frac{(t - t_{0})^{2}}{T^{2}}\right) \left(-2\frac{(t - t_{0})}{T^{2}}\sin(\omega t) + \omega \cos(\omega t)\right) \quad (11) $$

with amplitude $B_0 = 0.065\hbar(ea_0)^{-1}$, centered at $t_0 = 350.0\hbar E_h^{-1}$, with duration $T = 141.421356 \hbar E_h^{-1}$ and frequency $\omega = 0.15E_h \hbar^{-1}$. $\boldsymbol{\epsilon}^{\lambda}$ is the polarization vector, set to $\frac{1}{\sqrt{2}}(1.0, 1.0)$ for all calculations. (The parameters of the external time-dependent electric field were chosen to resemble a typical femtosecond pulse in terms of amplitude and duration, and the frequency was adjusted to be in resonance with the $S_0 \to S_1$ transition at the Franck–Condon point.) According to Maxwell’s equations, the external electric field is the time derivative of a (purely time-dependent, in the long-wavelength approximation) vector potential. Equation 11 is derived from such a vector potential assuming a Gaussian-shaped laser pulse, as discussed in ref 41. The transition dipole moment in the diabatic representation is chosen to be $\boldsymbol{\mu}_{12}(\mathbf{R}) = (f_X(X - X_0), f_Y(Y - Y_0))$, with $f_X = 0.2e$, $f_Y = 2.0e$, $X_0 = -1.0a_0$, and $Y_0 = -0.5a_0$. In the Condon approximation, $\boldsymbol{\mu}_{12}(\mathbf{R} = \mathbf{R}_{\text{init}})$ is used. The diagonal elements of the dipole operator are set to zero.

The full time-dependent Schrödinger equation is solved numerically in the diabatic representation employing a split-operator formalism$^{42,43}$ with a time step of $0.01\hbar E_h^{-1}$. A spatial grid of 800 points per coordinate is used over the ranges of $X \in [0.0, 8.0]a_0$ and $Y \in [-2.0, 2.0]a_0$. Diabatic quantities are consequently transformed to give all of the respective exact-factorization quantities of interest, namely, the time-dependent vector potential and the time-dependent potential energy surface. The expressions used for this transformation are given in Appendix A. To avoid numerical issues in the conversion, the time-dependent vector potential and the time-dependent potential energy surface are calculated only for regions of the nuclear configuration space where the nuclear density is larger than $10^{-8}$.

Phase factor $\theta(\mathbf{R}, t)$ is absorbed here into phase $S(\mathbf{R}, t)$ of the nuclear wave function, $\chi(\mathbf{R}, t) = |\chi(\mathbf{R}, t)| \exp[(i/\hbar)S(\mathbf{R}, t)]$, and the gauge is chosen so that the nuclear wave function is real and non-negative, i.e., $\chi(\mathbf{R}, t) = |\chi(\mathbf{R}, t)| \forall \mathbf{R}, t$. The exact nuclear density is determined as the sum of the squared moduli of the diabatic nuclear wave functions, and its positive square root yields the nuclear wave function $\chi(\mathbf{R}, t)$ in the chosen gauge.

## 2.3.2. Some Considerations for the Proposed Model.
In this work, we propose to simulate explicitly all of the steps in a photochemical process for the model system discussed above, as schematically represented in Figure 1. The model system

<figure><img src="image_8.png" /><figcaption>**Figure 1.** Schematic representation of the in silico photochemical experiment discussed in this work. The adiabatic potential energy surfaces for the ground ($S_0$) and excited ($S_1$) electronic states of the model are represented with a color map that shows positive electronic energies in purple and negative values in yellow/orange tones. The black circles and arrows indicate the different steps of the photochemical experiment: (1) initial state of the system in the Franck–Condon region, (2) photoexcitation from $S_0$ to $S_1$ triggered by a laser pulse, (3) relaxation through a conical intersection, and (4) formation of photoproducts.</figcaption></figure>

consists of two 2D parabolas shifted in the $X$ direction and in energy. At time $t = 0$, our molecular system is in its ground vibrational state, in the Franck–Condon region of the electronic ground state (point (1) in Figure 1). Then, a part of the initial wave function is excited to the $S_1$ excited electronic state via the coupling of the molecule with an explicit ultrashort laser pulse whose frequency is in resonance with the $S_0$-to-$S_1$ transition (point (2) in Figure 1). This photoexcitation generates a nuclear wavepacket in the excited electronic state, which will relax toward the conical intersection—a point of configuration space where the ground and first electronic states are degenerate—and funnel through it (point (3) in Figure 1). At this point, the nuclear wavepacket
experiences branching, preserving some of its amplitude in the excited electronic state while a good part of it has transferred back to the ground state as a result of nonadiabatic effects, evolving now as a photoproduct on the ground electronic state (point (4) in Figure 1).

We stress at this stage that the schematic representation of the photochemical experiment depicted in Figure 1 is deeply rooted in a (post-) Born–Oppenheimer picture. In section 3.1, we will discuss the details of the dynamics in terms of time-dependent nuclear wavepackets and static potential energy surfaces, thus adopting a Born–Oppenheimer vocabulary. In section 3.2, we will revisit this analysis from the perspective of the exact factorization, where the overall dynamics is examined in terms of a time-dependent nuclear wavepacket evolving according to a single time-dependent vector and scalar potential.

Another aspect that requires a comment at this stage is the coupling of an external time-dependent electric field, here a laser pulse, to a molecule. As described above, we use a semiclassical approach to couple an external time-dependent electric field to the molecular dipole operator. When introducing a given basis for our electronic states, once more preserving here a Born–Oppenheimer/Huang picture, we obtain couplings between the time-dependent electric field and the transition dipole moment between the pair of electronic states considered (here S₀ and S₁), μ₁₂(R). The magnitude and direction of the transition dipole moment depend on the nuclear position, as depicted by the color map and white arrows in the upper panel of Figure 2. This R dependence of the transition dipole moment implies that, within the long-wavelength approximation, the time-dependent electric field cannot be considered to be always aligned with the transition dipole moment (gray arrows in Figure 2, symbolizing the electric-field polarization vector), and one has to take the scalar product between the two quantities. The situation in which we account for the explicit R dependence of the transition dipole moment will now be referred to as non-Condon. The Condon approximation proposes to consider the transition dipole moment as a constant, set to its value at the Franck–Condon point, μ₁₂(R₃₄) (lower panel of Figure 2). Hence, only within the Condon approximation could we consider that the time-dependent electric field is always polarized along the transition dipole moment for all nuclear configurations (which is not the case in the present work).

In the non-Condon case, the strength of the coupling between the molecular system and the time-dependent electric field, given by

$$ \mathbf{V}^{(\text{dia})}(\mathbf{R}, t) = \begin{pmatrix} 0 & -\mu_{12}(\mathbf{R}) \cdot \mathbf{E}(t) \\ -\mu_{12}(\mathbf{R}) \cdot \mathbf{E}(t) & 0 \end{pmatrix} \quad (12) $$

depends on R as well. Conversely, in the Condon approximation the relative orientation of the transition dipole moment and the polarization vector of the time-dependent electric field is constant, leading to an R-independent coupling term in the Hamiltonian.

An additional topic that our work addresses is the representation of an in silico photochemical experiment using only the concept of nuclear trajectories, which will allow us to open a broader discussion of the steps toward excited-state molecular dynamics simulations. In section 3.2, we study the use of nuclear trajectories for nonadiabatic dynamics in the context of the exact factorization since well-

<figure><img src="image_9.png" /><figcaption>Figure 2. Representation of the transition dipole moment (white arrows give its direction and color map its intensity in ea₀) in the non-Condon case (upper panel) and the Condon case (lower panel) around the Franck–Condon region. The gray arrows indicate the polarization vector of the time-dependent electric field. The black contour lines show the nuclear density of the initial wave function at time t = 0.</figcaption></figure>

defined nuclear forces can be identified and provide a clear distinction between classical and quantum trajectories. For the quantum and classical trajectory dynamics, the 1000 initial nuclear positions are sampled from the probability density given by the initial nuclear wave function. Special care needs to be taken when selecting the initial nuclear momenta. While for classical trajectories nuclear positions and momenta can be regarded as independent variables, this is no longer the case for quantum trajectories because their nuclear position and momentum at each time step are related. As shown in Appendix B, the nuclear momentum at time t and position R is given by the time-dependent vector potential A(R, t). Thus, the 1000 initial nuclear momenta for the sampled initial nuclear positions R₀ are given by A(R₀, 0). For classical trajectories, the initial nuclear momenta are sampled from the momentum probability distribution computed from the Wigner transform of the initial nuclear wave function (in position representation). For the propagation of the trajectories, a time step of 10ħE⁻¹h is used, and one uses the exact time-dependent vector and scalar potentials as obtained from the quantum dynamics. The time-dependent potentials are not (numerically) available when the nuclear density is small (<10⁻⁸). Furthermore, classical trajectories are propagated according to the force computed as the negative of the gradient of the time-dependent scalar potential. Therefore, numerical errors due to the calculation of the nuclear gradients can cause the trajectories to move in a region of space where the time-dependent potential is not available. Those trajectories need to be removed from the ensemble. We
<figure><img src="image_10.png" /><figcaption>**Figure 3.** Snapshots of the nonadiabatic quantum dynamics at times 0, 10.2, and 31.5 fs (from left to right). The color maps indicate the (adiabatic) electronic energies of the S₁ and S₀ electronic states (see labels) in hartree (Eₕ). The adiabatic contribution to the nuclear density in each electronic state for the three snapshots is indicated by red (lκₕ₀(R, t)²) and purple (lκₕ₁(R, t)²) contour lines. The Franck–Condon point (FC) is indicated by a magenta cross, and the location of the conical intersection (CI) is indicated by a burgundy triangle.</figcaption></figure>

found that an energy-conservation criterion can be used to remove such unstable trajectories. To this end, we imposed that the classical energy should be conserved within a deviation of 0.01Eₕ from the initial value at the end of the laser pulse. A maximum of 51 classical trajectories have been excluded on the basis of this criterion in any of our simulations.

### 3. RESULTS AND DISCUSSION

#### 3.1. Photochemical Experiment in the Born–Oppenheimer Picture. 3.1.1. Nuclear Dynamics.
In a Born–Oppenheimer picture, the overall dynamics of an electronically excited molecule is analyzed by following the evolution of nuclear wavepackets on the potential energy surfaces (as depicted schematically in Figure 1). At time $t = 0$, the nuclear wave function is a stationary state in the adiabatic ground state S₀ (first column in Figure 3). The interaction with the time-dependent external field, consisting of a laser pulse, induces the electronic excitation of part of the original wave function to the first excited state S₁ (second column in Figure 3). The state created in S₁ is a nuclear wavepacket, which evolves and decays

back to the ground state through the conical intersection (third column in Figure 3).

First, let us investigate the quantum dynamics obtained in the non-Condon case (upper two rows of panels in Figure 3). The color map shows the potential energy of the S₁ (top row) and S₀ (second row from top) states, and the locations of two critical points of the potential energy surfaces—the Franck–Condon point (FC) and the conical intersection (CI)—are marked throughout the plots with a cross and a triangle for reference. Because we are now focusing on the Born–Oppenheimer picture of a photochemical process, the potential energy surfaces of the two adiabatic states do not move or change during the dynamics but can rather be seen as the electronic landscape on the support of which the nuclear wavepackets evolve. The nuclear densities associated with the S₀ and S₁ states are superimposed onto the respective potentials in Figure 3 and are indicated as red and purple contour lines, respectively. At time $t = 0$, the complete nuclear density is found in the S₀ state without any contribution from S₁. During the photoexcitation by the laser pulse ($t = 10.2$ fs in Figure 3; the time evolution of the laser field is shown in
Figure 4), nuclear amplitude is transferred to the excited electronic state. In addition, it can be observed that the nuclear

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

**Figure 4.** (Top) Time trace of the ground-state ($S_0$) population for the non-Condon (purple) and Condon (magenta) quantum dynamics simulations and the strength of the electric component of the laser pulse (gray line) in $E_h/ea_0$. (Bottom) Time trace of the transmission probability through an ideal barrier at $X > 3.86a_0$ for the non-Condon (purple lines) and Condon simulations (magenta lines). Solid lines indicate the total ($S_0 + S_1$) probability, and the dashed (dotted) lines give the $S_0$ ($S_1$) population contribution.

contribution remaining in the ground electronic state is also affected by the action of the laser pulse. Once the short laser pulse is over, the excited portion of the nuclear wavepacket relaxes on the $S_1$ potential energy surface, and at $t = 31.5$ fs, it reaches the conical intersection (last column in Figure 3). At this precise moment, part of the nuclear density transfers to the $S_0$ electronic state due to the influence of nonadiabatic effects. Interestingly, the portion of the nuclear wavepacket that remained in the Franck-Condon region in $S_0$ formed a nodal line as a result of the photoexcitation process. This nodal line is somehow tilted, that is, not parallel to the Y axis (red contour lines in the second row, last column panel of Figure 3).

Within the Condon approximation (bottom two rows in Figure 3), the potential energy surfaces remain identical to the non-Condon case but the transition dipole moment and, as explained above, the interaction Hamiltonian are no longer R-dependent. The overall dynamics is very similar to that observed in the non-Condon case, with the main difference being that the nodal line observed after the laser pulse on the ground-state nuclear wavepacket appears to be parallel to the Y axis rather than being tilted as in the non-Condon case. This feature can be clearly observed at times $t = 10.2$ and $31.5$ fs in

the $S_0$ portion of the nuclear density (bottom panels of Figure 3).

It is worth commenting further on the development of the nodal line in the ground-state portion of the nuclear wavefunction after the laser pulse. This nodal line suggests the formation of a higher vibrational eigenstate in the electronic ground state in the Franck-Condon region. The asymmetry of the nuclear density in $S_0$ observed at $t = 31.5$ fs results from the anharmonic nature of the $S_0$ potential energy surface for the low vibrational states. From a more general perspective, this observation highlights something to keep in mind when analyzing a photochemical experiment because the outcomes of a light-induced process—photoproducts or hot ground-state molecules—are usually thought to originate from the dynamics in the excited electronic states following photoexcitation. Our observation underlines the possibility of forming vibrationally excited molecules without nonadiabatic processes simply as a result of the coupling with a laser pulse. However, nonadiabatic dynamics simulations are often initiated directly from the formed nuclear wavepacket in a given excited electronic state, neglecting the remaining ground-state nuclear amplitude. Although that is often an adequate approximation, the simple example presented here could indicate that the influence of a laser pulse on the contribution remaining in the ground electronic state may not always be negligible.

### 3.1.2. Analysis of the Nuclear Dynamics.
A more quantitative analysis of the photochemical dynamics described in section 3.1.1 can be obtained by monitoring the electronic population evolution and the formation of photoproducts (Figure 4).

In the top panel of Figure 4, we present the time evolution of the ground-state population for both the Condon and the non-Condon quantum dynamics simulations with, for reference, the strength of the external electric field over time. Initially, the population is fully in the ground state. Because of the coupling between the molecule and the laser pulse, the population starts to be transferred to the excited state just before 5 fs. Maximum values of around 64% (non-Condon) and 67% (Condon) of the population are excited to $S_1$ at $t = 10.2$ fs, shortly after the laser pulse reaches its maximum intensity. The ground-state population subsequently plateaus at 48% (non-Condon) and 46% (Condon). After 26 fs, the excited-state population starts to decay back to the ground state as the $S_1$ nuclear wavepacket reaches the conical intersection (section 3.1.1). This ground-state population reaches a final plateau at about 98% after about 39 fs. Overall, the evolution of the ground-state population is very similar in the non-Condon and Condon quantum dynamics simulations, with only minor quantitative differences emerging.

It is important at this stage to stress that the electronic population dynamics reported in the top panel of Figure 4 and discussed above is a representation-dependent quantity. In other words, this quantity is not (strictly speaking) an observable, and the assignment of an adiabatic electronic state is intrinsically linked to the Born-Oppenheimer picture. Hence, we complement our earlier observations with an analysis of the evolution of the nuclear density in the configuration space. To this end, we calculate the transmission probability through an ideal barrier parallel to the Y axis and containing the conical intersection along the reaction coordinate (along X). We propose that $X_{CI}$, i.e., the position of the conical intersection in the X coordinate, delimits two regions: the photoreactant region ($X < X_{CI}$) and the
<figure><img src="image_12.png" /><figcaption>**Figure 5.** Schematic representation of the full time-dependent potential energy surface at three different times during the dynamics: 0, 10.2, and 31.5 fs. The positions of the Franck-Condon (magenta cross) point and conical intersection (burgundy triangle) are marked for reference. The color bar is given in hartrees ($E_h$) and we note that its negative range has been extended for the $t = 31.5$ fs snapshot. The TDPES is represented only in the regions of nuclear configuration space where the total nuclear density is $10^{-8}$ or larger.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 6.** Time-dependent potential energy surface at three different times (0, 10.2, and 31.5 fs (from left to right)) for the non-Condon (top row) and Condon quantum dynamics (bottom row). The color bar is given in hartrees ($E_h$). The nuclear density is superimposed (black contour lines). The positions of the Franck-Condon (magenta cross) point and conical intersection (burgundy triangle) are marked for reference.</figcaption></figure>

photoproduct region ($X > X_{CI}$). Therefore, we can estimate the formation of photoproducts by integrating the nuclear density over all values of Y and for $X > 3.86a_0$ (that is, $X > X_{CI}$). In the case of a real photochemical experiment, the formation of photoproducts could be probed,⁶ and the passage through a conical intersection can open up pathways in the ground electronic state that would have been inaccessible by thermal evolution.

The transmission probability for the non-Condon and Condon cases is shown in the bottom panel of Figure 4. Within the Condon approximation, 54.2% of the population is transferred to the photoproduct region. In the non-Condon case, this value is slightly lower, reaching 52.0%. The transmission probability can be decomposed into ground- and excited-state contributions, going back to a Born-Oppenheimer picture. We observe that a majority of the photoproducts are formed in their ground states, 50.2% (non-Condon) and 52.3% (Condon). Hence, the slight differences between non-Condon and Condon excitations observed in the representation-dependent population dynamics are reflected in this “representation-free” observable.

**3.2. Photochemical Experiment in the Exact-Factorization Picture.** Now that we have discussed our in silico

photochemical experiment from a Born-Oppenheimer perspective, we propose in the following section to reinterpret our findings using the tools offered by the exact-factorization picture. We recall that this representation proposes to replace the concepts of static potential energy surfaces associated with adiabatic electronic states and their nonadiabatic couplings with a single time-dependent potential energy surface (TDPES) and a time-dependent vector potential (TDVP). With the exact factorization, we move away from the representation of multiple electronic states visited by time-dependent nuclear wave functions and focus solely on a single nuclear wave function evolving under the influence of the single TDPES and TDVP.

For the non-Condon dynamics, Figure 5 schematically depicts the behavior of the TDPES for the three same times along the dynamics as in the previous section, i.e., $t = 0, 10.2, 31.5$ fs. At time $t = 0$, the TDPES exhibits a single well in the Franck-Condon region, basically reproducing the shape of the $S_0$ potential energy surface in this region (Figure 1). When the laser pulse is interacting with the molecule ($t = 10.2$ fs in Figure 5), the minimum in the TDPES decreases in energy and becomes asymmetrical around the Franck-Condon point. This distortion is due to the fact that the TDPES includes the
<figure><img src="image_14.png" /><figcaption>**Figure 7.** Time-dependent vector potential at three different times (0, 10.2, and 31.5 fs (from left to right)) for the non-Condon (top row) and Condon quantum dynamics (bottom row). The color map indicates the absolute value of the time-dependent vector potential in $ħ/a_0$, and the gray unit vectors show the orientation of the vector potential. The nuclear density is superimposed (black contour lines). The positions of the Franck–Condon (magenta cross) point and conical intersection (burgundy triangle) are marked for reference.</figcaption></figure>

effect of the time-dependent external field via $v(\mathbf{R}, t)$ (section 2.2). At later times ($t = 31.5$ fs in Figure 5), the TDPES develops a step separating the two portions of the nuclear wavepacket: one portion is localized in the Franck–Condon region, and the second one can be found in the vicinity of the conical intersection. In addition, the nodal line observed for the nuclear wave function in the Franck–Condon region is reflected in the TDPES as a potential barrier, reaching high positive and negative values. In the following section, we will provide a deeper analysis of the TDPES and its companion, the TDVP, for the studied photochemical experiment.

### 3.2.1. *Nuclear and Electronic Dynamics.*
As in section 3.1, we propose here to compare the quantum dynamics obtained with and without invoking the Condon approximation, but this time from an exact-factorization perspective.

Let us first focus on the evolution of the TDPES. In Figure 6, the full TDPES is plotted as a color map with superimposed black contour lines indicating the areas where the (full) nuclear density is mainly localized. At $t = 0$, both the non-Condon and the Condon simulations show an identical picture, where the TDPES shows a minimum around the center of the initial nuclear wave function and is curved upward to higher energies toward the borders, mainly reproducing the shape of the $S_0$ potential energy surface (as discussed for Figure 5). When the laser pulse reaches its maximum intensity (at around $t = 10.2$ fs), the TDPES is lower in energy in the area just around the center of the nuclear wave function in comparison to the previously shown time step, where the yellow/orange areas correspond to negative energy values. At this time step, there appears to be no significant differences between the TDPES computed with and without the Condon approximation. Later in time, when the laser pulse is over (Figure 6, $t = 31.5$ fs), a step appears within the TDPES: the portion of the TDPES with $X > 3a_0$ is significantly lower than that at $X < 3a_0$.

Furthermore, a marked difference between the non-Condon and Condon calculations emerges when looking at the TDPES in the region where the nodal line appears on the nuclear wave function around the FC position. In the non-Condon case, a localized peak can be observed between the two portions of the nuclear wave function, whereas a barrier forms within the Condon approximation, almost parallel to the Y axis and stretching all through the TDPES. It is worth stressing that the interesting features of the TDPES at $t = 31.5$ fs are all localized around the FC point and are caused by the formation of an eigenstate in this portion of the configuration. The TDPES in the region where a conical intersection would be observed in the Born–Oppenheimer picture is blatantly featureless, basically leading the nuclear wavepacket toward the photo-product region ($X > X_{CI}$).

The TDPES is one of the time-dependent quantities of interest within the exact factorization. Another key ingredient of this formalism is the TDVP. We recall that the TDVP is equivalent to the nuclear momentum field within the chosen gauge, as detailed in section 2.2. Figure 7 shows the TDVP for the three critical times of the dynamics, with and without the Condon approximation. At time $t = 0$, the magnitude of the TDVP is very small. This can be understood from the definition of the initial condition for our dynamics, which is the ground vibrational state for the ground electronic state, at the FC point. From a quantum trajectory perspective, the momentum field corresponding to a nuclear (real) eigenstate would be zero everywhere.⁴⁴ During the excitation, the TDVP increases in magnitude and triggers the dynamics of the nuclear wavepacket. Conversely with respect to our earlier observations with the TDPES, the TDVP has already developed some differences at $t = 10.12$ fs depending on the use of the Condon approximation. In the non-Condon simulation, the TDVP shows a non-negligible contribution
<figure><img src="image_15.png" /><figcaption>**Figure 8.** Curl of the TDVP at *t* = 31.5 fs for the non-Condon (top row) and Condon quantum dynamics (bottom row). The curl of the TDVP is given for the FC region *X* < 2.8*a*₀ (left panels) as well as the CI and photoproduct regions *X* > 2.8*a*₀ (right panels). The positions of the Franck-Condon (magenta cross) point and conical intersection (burgundy triangle) are marked for reference.</figcaption></figure>

along the Y direction, whereas within the Condon approximation, such a contribution is basically zero and the vector potential has components only along the X axis (see unit arrows in Figure 7). In both cases, however, one can observe an abrupt change in direction of the TDVP along the X-axis at around *X* = 1.9*a*₀ for all values of *Y*. At *t* = 31.5 fs, the TDVP develops two clear portions. In the region of *X* > 3*a*₀, the TDVP has an overall larger magnitude. In this region, the nuclear wavepacket is indeed pushed toward larger *X* values following the slope of the TDPES. Interestingly and consistent with our observation on the TDPES, this portion of the TDVP is smooth and does not reveal any features that would testify to a specific electron-nuclear coupling. This observation is particularly important if one considers that, in the Born-Oppenheimer picture, the nuclear wavepacket is passing through a conical intersection at this specific time (Figure 3). Such a behavior of the TDPES and TDVP simply highlights that the coupled electron-nuclear dynamics presented here should be seen as a non-event, with a molecule simply relaxing in energy, driven by the TDPES and TDVP. However, a Born-Oppenheimer picture enforces the description of this process with adiabatic electronic states, which are not representative of single electronic character. As a result, the *S*₁ nuclear wavepacket in the Born-Oppenheimer picture suffers nearly singular nonadiabatic coupling to transfer to the *S*₀ electronic state, hence preserving its electronic character. Conical intersections and singular nonadiabatic coupling terms therefore emanate from the intrinsic limitations posed by the adiabatic representation of the Born-Oppenheimer picture, and the exact factorization naturally remedies this issue by eliminating the notion of electronic states altogether.

Around the (FC) position, the TDVP has significantly smaller values at *t* = 31.5 fs, highlighting the rather stationary nature of the portion of the nuclear wave function localized in this region (Figure 7). In the region around *X* = 2.1*a*₀ where the nuclear density splits, an intense localized peak appears on the TDVP in the non-Condon case while the Condon simulation leads to the appearance of an intense “delocalized”

barrier parallel to the Y axis. In the non-Condon case, the unit vectors of the TDVP betray the fact that the nuclear wavepacket is somehow rotating around *X* = 2.1, *Y* = 0*a*₀.

The TDVP is in general not irrotational,³⁰ and the calculation of its line integral along a closed loop in nuclear (R) space yields a nonzero quantity. Such a quantity is a geometric phase whose value depends on the integration path. This dependence on the integration path means that this geometric phase is not quantized but is instead the case for the topological Berry phase arising in the post-Born-Oppenheimer picture in the presence of a conical intersection. In Figure 8, we prove numerically that the curl of the TDVP, i.e., $[\text{curl } A(X, Y)]_Z = \partial_X A_Y(X, Y) - \partial_Y A_X(X, Y)$, is nonzero at time $t = 31.5$ fs in the non-Condon (upper panels) and Condon (lower panels) cases. Figure 8 gives the magnitude of the curl of the TDVP, which has only a Z component orthogonal to the X, Y plane in two dimensions. Because of the large difference in the curl magnitude, the plot for the FC region ($X < 2.8a_0$, left panels) is separated from that for the CI region ($X > 2.8a_0$, right panels). In line with our previous observations, the curl of the TDVP highlights the main differences between the non-Condon and Condon dynamics in the FC region. In the region where the ground-state wavepacket is located at the end of the laser pulse, the curl of the TDVP is nonzero only around $X = 2.1$, $Y = 0a_0$ in the non-Condon case. The “rotating” behavior of the TDVP is encoded in the change of sign of its curl. (We checked numerically that this observed structure is robust (and not a numerical artifact) by varying the cutoff parameter for the calculation of the TDVP from $10^{-6}$ to $10^{-10}$.) The observed structure in the non-Condon case is lost in the Condon simulation, where the curl of the TDVP is nonzero and always positive only along a line at $X = 2.1a_0$ parallel to the Y axis. The two simulations yield similar results in the CI region. In particular, the curl of the TDVP is zero at the CI. The change in curl sign between positive and negative values of *Y* attests from the overall spreading of the nuclear amplitude in this region.
To further unravel the differences observed in the TDPES and TDVP for the non-Condon versus Condon treatment of quantum dynamics, we decompose the TDPES into $\epsilon(\mathbf{R}, t)$ and $v_{\text{int}}(\mathbf{R}, t)$ at 2 times close to the maximum in the field envelop. The 8.95 fs time corresponds to a snapshot before the electric field of the laser pulse reaches a maximum while the other snapshot, $t = 9.20$ fs, is just after the maximum is passed. Figure 9 presents a close-up view of the TDPES and TDVP around the FC region for these two specific times. Looking only at the contribution to the TDPES coming from the coupling of the laser pulse to the molecule ($v_{\text{int}}(\mathbf{R}, t)$ in eq 7 ($v_{\text{int}}$ top panels in Figure 9)), we observe in the non-Condon case that an anisotropy along the Y direction is created, which is absent in the Condon approximation. This anisotropy is also visible in the term $\epsilon(\mathbf{R}, t)$ (TDPES without $v_{\text{int}}$ middle panels in Figure 9), where for instance the negative contributions of the TDPES are not symmetric along the Y direction in the non-Condon case. The variation of the TDPES caused by the interferences between the two components of the nuclear wave function is reminiscent of the observations of interference in nonadiabatic processes.³⁹ Zooming in on the TDVP (bottom panels of Figure 9) highlights some additional interesting features. The TDVP already exhibits strong contributions along the Y-direction in the early stage of the dynamics in the non-Condon case, when the laser pulse is present, while in the Condon case no Y-contribution can be observed. The line along which the X-component of the TDVP changes sign is not parallel to the Y-axis in the non-Condon simulation, which is opposite to the Condon simulation.

The observations presented in this section show that the exact-factorization formalism offers an alternative perspective of a photochemical process to the conventional Born–Oppenheimer picture. The contribution to the molecular Hamiltonian coming from the light–matter interaction is fully encoded in the TDPES, and the evolution of the nuclear wavepacket is driven by the (single) TDPES and TDVP, which in turn encode the effect of the electronic dynamics on the nuclei. In the particular gauge chosen in this work, the TVDP is identified as the nuclear momentum field. Furthermore, we highlighted the combined effects of the TDPES and TDVP in inducing the splitting of the nuclear wavepacket into two portions in the Franck–Condon region, resembling the formation of a higher vibrational state of the electronic ground-state potential. Instead, the TDPES and TDVP are smooth where, in the Born–Oppenheimer picture, the nuclear wavepacket passes through a conical intersection and suffers singular nonadiabatic couplings.

The exact factorization and its time-dependent potential also highlight another interesting fact about photochemistry and photophysics. The presence of potential energy surfaces in a Born–Oppenheimer picture invites us to think of photochemical processes and chemical reactions in the ground electronic state in similar ways: by looking at these static potentials and possibly looking at their critical points and how they connect via minimum-energy paths. However, it is crucial to realize that such concepts are often not sufficient or valid in understanding a given photochemical process.²¹ In other words, a molecule does not have a defined photochemical reactivity per se, but such a photochemical reactivity depends on the initialization of the process (type of photoexcitation) and the following out-of-equilibrium dynamics on the coupled potential energy surfaces. The exact factorization highlights this fact by exhibiting different time-dependent potentials even

<figure><img src="image_16.png" /><figcaption>**Figure 9.** Time-dependent potential energy surface and vector potential during the laser pulse, for the non-Condon (top row of each panel) and Condon (bottom row of each panel) dynamics. Top panel: $v_{\text{int}}$, contribution to the time-dependent potential energy surface coming from an external potential (here the interaction</figcaption></figure>
**Figure 9.** continued

between the laser pulse and the molecule). Middle panel: time-
dependent potential energy surface without the external potential.
Bottom panel: time-dependent vector potential, where the color map
indicates the absolute value and the black unit vectors show the
orientation. The Franck–Condon point is indicated by a magenta
cross in all plots. Color bars for the top and middle panels are given in
hartrees (E<sub>h</sub>), and those of the lower panel are given in ℏ/a<sub>0</sub>.

for subtly different photochemical dynamics (here, simply
employing or not employing the Condon approximation).

**3.2.2. Further Analysis of a Photochemical Process with Quantum and Classical Trajectories.** Since the exact factorization encodes all of the coupled electron–nuclear dynamics in the TDPES and TDVP, this framework naturally lends itself to the propagation of nuclear trajectories, which can be used to further analyze and approximate the nuclear quantum dynamics without the need to introduce additional approximations. The challenge of propagating trajectories in a Born–Oppenheimer picture to describe nonadiabatic dynamics is notorious and caused by the presence of the lack of well-defined nuclear forces due to the presence of different (sometimes coupled) potential energy surfaces.<sup>10,38</sup>

As discussed in Appendix B, quantum trajectories can be propagated easily within the exact factorization by using the TDVP as a nuclear momentum field, integrating eq 25. Although the initial nuclear positions can be randomly sampled from the initial nuclear density, the initial nuclear momenta have to be chosen as $A(\mathbf{R}_0, t = 0)$ (where $\mathbf{R}_0$ stands for the whole set of initial positions) because positions and momenta are not independent variables. We propagated 1000 quantum trajectories using the TDVP calculated at all times from the quantum dynamics simulations, without and with the Condon approximation.

As expected from their definition, the quantum trajectories
closely follow the nuclear density of the quantum wavepacket
at all three times as previously discussed (Figure 10). We show
the results only for the non-Condon dynamics because the
Condon dynamics are very similar. Several observations can be
made by looking at the time evolution of the quantum
trajectories in the non-Condon case (movies in the Supporting
Information). First, we can see that the motion of the quantum
trajectories is dominated by an evolution in the X direction.
Second, the trajectories forming the two split portions of the
nuclear wavepacket that remain in the vicinity of the FC region
appear to rotate around a point close to the initial position of

the nuclear wave function (but not exactly at the FC point due
to the anharmonicity of the potential energy surface in this
region). Such behaviors were predicted earlier on the basis of
the analysis of the shape of the TDVP (section 3.2.1).
Conversely, within the Condon approximation the quantum
trajectories representing the nuclear wavepacket around the FC
region do not significantly move once the splitting is complete.
This behavior further supports the stationary nature of the
driving wave function in this region and attests to the fact that
the Y contribution of the momentum field is basically zero at
later times in this region (as observed in Figure 7).

The distribution of the nuclear momenta for each quantum
trajectory of the swarm highlights striking differences between
the non-Condon and Condon dynamics at all times (Figure
11). At *t* = 0, the TDVP is very small in magnitude, i.e., close
to zero, which is why the distribution of initial momenta
appears to be highly localized around (0.0, 0.0). Subsequently,
the quantum trajectories acquire larger momenta and start
spreading to finally be distributed in three distinct areas by *t* =
31.5 fs in the non-Condon case (top right panel in Figure 11).
The first group of trajectories shows large (positive) values for
*P*<sub>X</sub> between 60 and 80 *ħ*a<sub>0</sub><sup>-1</sup> and represents the photoproduct
trajectories. Two groups of trajectories with comparably
smaller *P*<sub>X</sub> values spread along two lines, reaching a value of
-4*ħ*a<sub>0</sub><sup>-1</sup> or 4*ħ*a<sub>0</sub><sup>-1</sup> for *P*<sub>Y</sub>. These last two groups depict the
slow quantum trajectories in the region where a stationary state
is formed. The distribution of nuclear momenta for the
dynamics conducted within the Condon approximation looks
significantly different from those of the non-Condon case
described above (bottom row in Figure 11). At *t* = 10.2 fs, all
of the trajectories have a dominant *X* contribution to their
nuclear momenta. At *t* = 31.5 fs, a partitioning of the
trajectories in momentum space is visible, as observed in the
non-Condon case, forming two main groups: trajectories with
a large *P*<sub>X</sub> value and those with smaller momenta in the *X*
direction. We note that the group with a smaller *P*<sub>X</sub>
contribution appears to be further split at the same position
in *X* as in the non-Condon case. We also notice that a small
portion of the trajectories with a very large *P*<sub>X</sub> contribution at *t*
= 31.5 fs also starts to spread along *P*<sub>Y</sub>.

The fate of the quantum trajectories over time is made
clearer by plotting their traces in time and space (Figure 12).
We note that the quantum trajectories presented here were
initialized from an equally spaced grid to improve the clarity of
the plots in Figure 12. The trace representation of the quantum
trajectories highlights the presence of two groups of quantum
trajectories: those that remain that evolve around the FC

<figure><img src="image_17.png" /><figcaption>**Figure 10.** Positions of the quantum trajectories at three different times (0, 10.2, and 31.5 fs (from left to right)) for non-Condon dynamics. The full nuclear density from the quantum dynamics is superimposed (black contour lines). The positions of the Franck–Condon (magenta cross) point and conical intersection (burgundy triangle) are marked for reference.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 11.** Components of the nuclear momentum in X ($P_X$) and Y ($P_Y$) for each quantum trajectory at three different times (0, 10.2, and 31.5 fs (from left to right)) for non-Condon (top row) and Condon (bottom row) dynamics. The label a.u. stands for $\hbar/a_0$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 12.** Traces of quantum trajectories in space and time for the non-Condon dynamics. These quantum trajectories were specially initialized on an equally spaced grid in configuration space to enhance clarity. The left panel shows in yellow/orange the trajectories remaining in the FC region. The right panel depicts in pink/purple the trajectories moving away from the FC region to reach the photoproduct region. The same trajectories are shown in both panels.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 13.** Positions of the classical trajectories at three different times (0, 10.2, and 31.5 fs (from left to right)) for non-Condon dynamics. The full nuclear density from the quantum dynamics is superimposed (black contour lines). The positions of the Franck-Condon (magenta cross) point and conical intersection (burgundy triangle) are marked for reference.</figcaption></figure>

region (left panel of Figure 12) and those that leave that translate the formation of photoproducts (right panel of Figure 12). Interestingly, it appears from this analysis that quantum trajectories with a larger X value in the FC region are more

likely to follow the photoproduct path than those started at $t = 0$ fs at a smaller value of X. This trace representation for the quantum trajectories also offers an opportunity to show that they do not cross in configuration space.⁴⁵ This noncrossing
<figure><img src="image_21.png" /><figcaption>**Figure 14.** Components of the nuclear momentum in X ($P_X$) and Y ($P_Y$) for each classical trajectory at three different times (0, 10.2, and 31.5 fs (from left to right)) for non-Condon (top row) and Condon (bottom row) dynamics. The label a.u. stands for $\hbar/a_0$.</figcaption></figure>

rule is a strict requirement for quantum trajectories and a direct consequence of the single-valuedness of the nuclear wave function.

An interesting feature of quantum trajectories is that they can easily be transformed into their more classical analogues by neglecting the so-called quantum potential, which acts as a nonlocal glue for the trajectories (Appendix B). We can therefore propagate classical trajectories by using the TDPES and TDVP computed from the quantum wavepacket dynamics, but with the removal of the contribution from the quantum potential $Q_{\text{pot}}(\mathbf{R}, t)$. Hence, the classical trajectories, initialized from 1000 nuclear positions and momenta sampled from a Wigner distribution, come directly from the integration of eqs 27 and 28.

The positions of these classical trajectories at $t = 0, 10.2, 31.5$ fs show rather good agreement with the nuclear wavepacket (Figure 13), even if it is clear by comparison with Figure 10 that the classical trajectories spread more in both the X and Y directions over time. (We note that the same initial positions were used for the classical and quantum trajectories.) Looking at the full evolution of the classical trajectories over time further reinforces this observation (movies in the Supporting Information), with trajectories exhibiting important oscillations along the Y coordinate.

The projection of the classical trajectories in nuclear momentum space for the three selected times (Figure 14) is strikingly different from that of the quantum trajectories discussed earlier (Figure 11). The distribution of nuclear momenta for the classical trajectories is spread at all times and does not exhibit the peculiar structure observed for the quantum trajectories. The magnitude of the nuclear momenta is also significantly larger in the X direction and, more specifically, in the Y direction. This observation explains the oscillating behavior in the Y coordinate described above for the classical trajectories. At $t = 31.5$ fs, the momenta distribution separates into only two regions (Figure 14): one part with

momenta distributed around $P_X = 0$ and the other part with momenta distributed around $P_X = 70\hbar a_0^{-1}$. Interestingly, the distribution of classical nuclear momenta in the non-Condon case closely resembles that in the Condon case, as if the removal of the quantum potential washed out the fine differences observed during the formation of the stationary state.

We conclude this analysis by calculating the transmission toward photoproducts, as performed in section 3.1.2 but here based on the distribution of the quantum and classical trajectories over time. To this end, we simply count the trajectories with $X > 3.86a_0$ (as described in section 3.1.2) for both quantum and classical trajectories and in non-Condon and Condon dynamics. Quantum trajectories appear to slightly underestimate the reference transmission probability (from the quantum dynamics) by around 0.01 for the non-Condon dynamics and 0.02 for the Condon dynamics (top panel of Figure 15). The classical trajectories overestimate the transmission probabilities in both non-Condon and Condon cases by 0.04 and 0.05, respectively (bottom panel of Figure 15). This deviation is marginal, and both quantum and classical trajectory-based dynamics reproduce reasonably well the qualitative evolution of the nuclear density while providing quantitatively good estimates for the transmission probability, an observable that would, for a real molecule, connect to the formation of photoproducts and thus to the quantum yield of a photochemical reaction.

One important aspect that needs to be stressed at this point is the fact that both quantum and classical trajectories benefitted from the definition of the time-dependent potentials used to propagate them thanks to the formalism of the exact factorization (and within our choice of gauge). More specifically, simulating a full photochemical experiment as done here using quantum or classical trajectories within a Born-Oppenheimer picture would have made the simulations dramatically more complex. Transfers, hops, averaging, or
<figure><img src="image_22.png" /><figcaption>Figure 15. Time trace of the transmission probabilities through the line defined by X > 3.86a₀ obtained for the quantum (top panel, red lines) and classical trajectories (bottom panel, orange lines) in the non-Condon and Condon cases. The results from quantum dynamics (QD) are given for reference (non-Condon with purple lines and Condon with magenta lines).</figcaption></figure>

spawns would have to be invoked to allow the trajectories to visit different electronic states due to the influence of a laser pulse or nonadiabatic transitions. The exact factorization simplifies that by providing the single TDPES and TDVP.

## 4. CONCLUSIONS

We introduced in this work a comparison between the Born–Oppenheimer and exact-factorization picture of an in silico photochemical experiment, from the initial photoexcitation with an ultrashort laser pulse to the formation of photo-products. Our extensive analysis of the subsequent photo-dynamics revealed how to picture a photochemical reaction in the exact factorization formalism and how its quantities, namely, the time-dependent vector and scalar potentials, behave when an external laser pulse shakes the electronic wave function and triggers the formation of a nuclear wavepacket. As a side product of the excitation with a laser pulse, we could also observe the formation of a stationary state in the Franck–Condon region and the corresponding behavior of the time-dependent potentials. The Condon approximation can alter the dynamics of the nuclear wavepacket, even if the formation of photoproducts, which is a representation-free quantity, does not suffer from this approximation in the presented model. Interestingly, the passage through a conical intersection in the Born–Oppenheimer picture becomes a nonevent in the exact factorization: the singularities and degeneracies at conical intersections are in stark contrast with the featureless time-dependent potentials of the exact factorization. Finally, we highlighted one of the exciting

features of the exact factorization in the context of a full photochemical process: the possibility to naturally introduce trajectories as an approximation for the nuclear dynamics. The propagation of these trajectories is trivial because the exact factorization has only a single time-dependent vector and scalar potential, meaning that no hops or spawns are required to describe regions of strong nonadiabaticity. Distinguishing between classical and quantum trajectories, we discussed the suitability of using trajectories to simulate photochemical processes with the exact factorization as well as the importance of properly selecting initial conditions.

## ■ APPENDIX A: EXPRESSING THE TDVP AND TDPES WITHIN A DIABATIC BASIS

Quantum mechanics allows one to switch from the Dirac representation $|\psi\rangle$ of a state to the wave function representation $\psi(\mathbf{r})$ with the operation $\langle\mathbf{r}|\psi\rangle = \psi(\mathbf{r})$. The same holds true for the exact-factorization quantities: the electronic-state vector at time $t$, $|\Phi(t; \mathbf{R})\rangle$, which parametrically depends on $\mathbf{R}$, is related to the conditional wave function of eq 6 as $\langle\mathbf{r}|\Phi(t; \mathbf{R})\rangle = \Phi(\mathbf{r}, t; \mathbf{R})$. Since the electronic Hamiltonian (eq 9) is given in the diabatic basis, we represent in our calculations the state $|\Phi(t; \mathbf{R})\rangle$ in such a basis and identify the expansion coefficients (or diabatic nuclear amplitudes) as $\chi_l^{(d)}(\mathbf{R}, t)$. All quantities analyzed in this work, namely, the nuclear density, the TDPES, and the TDVP, are obtained by integrating over the electronic degrees of freedom, thus by summing over the diabatic states. We give their expressions below.

The nuclear density $|\chi(\mathbf{R}, t)|^2$ is simply

$$|\chi(\mathbf{R}, t)|^2 = \sum_l |\chi_l^{(d)}(\mathbf{R}, t)|^2 \qquad (13)$$

The explicit expression of the TDPES is

$$\epsilon(\mathbf{R}, t) = \langle\Phi(t; \mathbf{R})|\hat{H}_{\text{BO}}(\mathbf{R})|\Phi(t; \mathbf{R})\rangle_{\mathbf{r}} \\ + \sum_{\nu=1}^{N_{\nu}} \left[ \frac{\hbar^2}{2M_{\nu}} \langle\nabla_{\mathbf{R}_{\nu}} \Phi(t; \mathbf{R}) | \nabla_{\mathbf{R}_{\nu}} \Phi(t; \mathbf{R})_{\mathbf{r}} \rangle - \frac{\mathbf{A}_{\nu}^2(\mathbf{R}, t)}{2M_{\nu}} \right] \\ + \langle\Phi(t; \mathbf{R})| - i\hbar\partial_{\nu}|\Phi(t; \mathbf{R})\rangle_{\mathbf{r}} \qquad (14)$$

and can thus be decomposed as the sum of four terms $\epsilon(\mathbf{R}, t) = \epsilon_{\text{GI1}}(\mathbf{R}, t) + \epsilon_{\text{GI2}}(\mathbf{R}, t) + \epsilon_{\text{GI3}}(\mathbf{R}, t) + \epsilon_{\text{GD}}(\mathbf{R}, t)$, with the first three being gauge-invariant (GI) and the last one being gauge-dependent (GD) under the (gauge) transformations discussed in section 2.2. Transforming the integration over $\mathbf{r}$ to the sum over the diabatic states, considering that the nonadiabatic coupling vectors are identically zero in this basis, and using the notation for the matrix elements of the electronic Hamiltonian introduced in eq 9, the terms in the expression of the TDPES are

$$\epsilon_{\text{GI1}}(\mathbf{R}, t) = \sum_{l,k} \frac{\chi_k^{(d)*}(\mathbf{R}, t) \chi_l^{(d)}(\mathbf{R}, t)}{|\chi(\mathbf{R}, t)|^2} V_{kl}(\mathbf{R}) \qquad (15)$$

$$\epsilon_{\text{GI2}}(\mathbf{R}, t) = \sum_{\nu} \frac{\hbar^2}{2M_{\nu}} \sum_{l} \left| \nabla_{\mathbf{R}_{\nu}} \frac{\chi_l^{(d)}(\mathbf{R}, t)}{\chi(\mathbf{R}, t)} \right|^2 \qquad (16)$$
$$ \epsilon_{\text{GI3}}(\mathbf{R}, t) = -\sum_{\nu} \frac{\mathbf{A}_{\nu}^2(\mathbf{R}, t)}{2M_{\nu}} \qquad (17) $$

$$ \epsilon_{\text{GD}}(\mathbf{R}, t) = -i\hbar \sum_l \frac{\chi_l^{(d)*}(\mathbf{R}, t)}{\chi^*(\mathbf{R}, t)} \partial_l \frac{\chi_l^{(d)}(\mathbf{R}, t)}{\chi(\mathbf{R}, t)} \qquad (18) $$

Similarly, the TDVP is given by

$$ \mathbf{A}_{\nu}(\mathbf{R}, t) = -i\hbar \sum_l \frac{\chi_l^{(d)*}(\mathbf{R}, t)}{\chi^*(\mathbf{R}, t)} \nabla_{\mathbf{R}_{\nu}} \frac{\chi_l^{(d)}(\mathbf{R}, t)}{\chi(\mathbf{R}, t)} \qquad (19) $$

While the TDPES and TDVP can be expressed in the diabatic basis, as done above, it is critical to note that they do not depend on any particular choice of electronic representation.

## ■ APPENDIX B: TRAJECTORY-BASED SOLUTION OF THE NUCLEAR EXACT FACTORIZATION EQUATION

Because the coupled nuclear and electronic equations of motion give rise to a single time-dependent vector potential and a single time-dependent potential energy surface in the exact factorization, one intriguing question is whether one could simplify the nuclear dynamics by representing the nuclear probability density with trajectories evolving on the basis of the two potentials. In previous work,³⁵,⁴³ it was shown that by inserting the polar form of the nuclear wave function, $\chi(\mathbf{R}, t) = |\chi(\mathbf{R}, t)|\exp(\frac{i}{\hbar}S(\mathbf{R}, t))$, into the time-dependent Schrödinger equation, one can derive an evolution equation for the phase that can furthermore be identified as a (nuclear) Hamilton–Jacobi equation (with $\nabla_{\mathbf{R}_{\nu}} S(\mathbf{R}, t) = \tilde{\mathbf{P}}_{\nu}(\mathbf{R}, t)):

$$ \begin{aligned} -\frac{\partial}{\partial t} S(\mathbf{R}, t) = & H_n(\tilde{\mathbf{P}}, \mathbf{R}, t) \\ = & \sum_{\nu} \frac{|\tilde{\mathbf{P}}_{\nu}(\mathbf{R}, t) + \mathbf{A}_{\nu}(\mathbf{R}, t)|^2}{2M_{\nu}} + \epsilon(\mathbf{R}, t) + v_{\text{int}}(\mathbf{R}, t) \\ & + \sum_{\nu} \frac{-\hbar^2 \nabla_{\mathbf{R}_{\nu}}^2 |\chi(\mathbf{R}, t)|}{2M_{\nu} |\chi(\mathbf{R}, t)|} \end{aligned} \qquad (20) $$

The last term on the right-hand side is the so-called quantum potential, $Q_{\text{pot}}(\mathbf{R}, t)$, which is responsible for an important portion of the nuclear quantum effects. The term containing the nuclear momentum field, $\mathbf{P}_{\nu}(\mathbf{R}, t) = \tilde{\mathbf{P}}_{\nu}(\mathbf{R}, t) + \mathbf{A}_{\nu}(\mathbf{R}, t)$, is the nuclear kinetic energy, and the sum $\epsilon(\mathbf{R}, t) + v_{\text{int}}(\mathbf{R}, t) + Q_{\text{pot}}(\mathbf{R}, t)$ represents the potential energy. As discussed in an earlier work,⁴³ this Hamilton–Jacobi partial differential equation can be solved by using the method of characteristics. This results in Hamilton-like evolution equations (which are ordinary differential equations to be solved for an a priori infinite number of initial conditions)

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

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

or equivalently

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

$$ \dot{\mathbf{P}}_{\nu}(t) = -\nabla_{\mathbf{R}_{\nu}} H_n(\mathbf{P}(t), \mathbf{R}(t), t) + \dot{\mathbf{A}}_{\nu}(\mathbf{R}(t), t) \qquad (24) $$

for the trajectories $\mathbf{R}(t) \equiv \{\mathbf{R}_{\nu}(t)\}_{\nu=1,...,N_{\nu}}$ with corresponding momenta $\mathbf{P}(t) \equiv \{\mathbf{P}_{\nu}(t)\}_{\nu=1,...,N_{\nu}}$.

Enforcing our choice of gauge for this work ($S(\mathbf{R}, t) = 0$, see section 2.2) in the characteristic representation of the nuclear dynamics, we find that eqs 21 and 22 simplify to

$$ \dot{\mathbf{R}}_{\nu}(t) = \frac{\mathbf{A}_{\nu}(\mathbf{R}(t), t)}{M_{\nu}} \qquad (25) $$

$$ \tilde{\mathbf{P}}_{\nu}(t) = 0 \qquad (26) $$

proving that the time-dependent vector potential fully accounts for the nuclear momenta of the trajectories. It is worth noting that, within the exact factorization, the nuclear momentum field computed with $\Psi(\mathbf{r}, \mathbf{R}, t)$ is the sum of a curl-free term—the gradient of $S(\mathbf{R}, t)$—and a (in general not irrotational) contribution due to the time-dependent vector potential, as shown in ref 38. With our choice of gauge, the curl-free contribution from the gradient of $S(\mathbf{R}, t)$ becomes zero. In Appendix C, we prove analytically the identity given in eq 26 and show that, in the chosen gauge, $\nabla_{\mathbf{R}_{\nu}} H_n(\tilde{\mathbf{P}}(t), \mathbf{R}(t), t) = 0$.

Propagating trajectories solely using either eq 25 or coupled eqs 23 and 24 yields what we refer to as quantum trajectories. In particular, eq 25 corresponds to the Bohmian definition of the velocity field for quantum-mechanical particles described by a real wave function in the presence of an external vector potential.⁴⁶ By neglecting the quantum potential in $H_n(\mathbf{P}(t), \mathbf{R}(t), t)$, the nuclear Hamiltonian becomes classical, i.e., $H_n^{\text{cl}}(\mathbf{P}(t), \mathbf{R}(t), t)$. This classical approximation is used when integrating eqs 23 and 24 to generate classical trajectories. The corresponding evolution equations for the classical trajectories then read

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

$$ \dot{\mathbf{P}}_{\nu}(t) = -\nabla_{\mathbf{R}_{\nu}} H_n^{\text{cl}}(\mathbf{P}(t), \mathbf{R}(t), t) + \dot{\mathbf{A}}_{\nu}(\mathbf{R}(t), t) \qquad (28) $$

## ■ APPENDIX C: NUCLEAR FORCES FOR THE QUANTUM TRAJECTORIES

In this appendix, we prove eq 26 of Appendix B by demonstrating that, within the chosen gauge, $S(\mathbf{R}, t) = 0$, $\nabla_{\mathbf{R}_{\nu}} H_n(\tilde{\mathbf{P}}, \mathbf{R}, t) = 0$.

The nuclear Hamiltonian identified in the Hamilton–Jacobi equation in eq 20 reads

$$ H_n(\mathbf{R}, t) = \sum_{\nu} \frac{\mathbf{A}_{\nu}(\mathbf{R}, t)^2}{2M_{\nu}} + \epsilon(\mathbf{R}, t) + v_{\text{int}}(\mathbf{R}, t) \\ + \sum_{\nu} \frac{-\hbar^2 \nabla_{\mathbf{R}_{\nu}}^2 |\chi(\mathbf{R}, t)|}{2M_{\nu} |\chi(\mathbf{R}, t)|} \qquad (29) $$

when imposing the gauge as $\mathbf{P} \approx 0$. From the definition of the time-dependent scalar potential $\epsilon(\mathbf{R}, t)$ given in section 2.2 and in eq 14, we identify $\epsilon_{\text{BO}}(\mathbf{R}, t) = \epsilon_{\text{GI1}}(\mathbf{R}, t)$ and $\epsilon_{\text{bBO}}(\mathbf{R}, t) = \epsilon_{\text{GI2}}(\mathbf{R}, t)$ such that the full nuclear Hamiltonian reads

$$ H_n(\mathbf{R}, t) = \epsilon_{\text{BO}}(\mathbf{R}, t) + \epsilon_{\text{bBO}}(\mathbf{R}, t) + \epsilon_{\text{GD}}(\mathbf{R}, t) \\ + v_{\text{int}}(\mathbf{R}, t) + Q_{\text{pot}}(\mathbf{R}, t) \qquad (30) $$
It is easy to prove that while $\epsilon_{BO}(\mathbf{R}, t)$ (Born–Oppenheimer term) and $\epsilon_{bBO}(\mathbf{R}, t)$ (beyond-Born–Oppenheimer term) are gauge-invariant contributions because they remain unchanged if the electronic wave function is modified by a phase factor that depends only on $\mathbf{R}$ and $t$, $\epsilon_{GD}(\mathbf{R}, t)$ is gauge-dependent (GD).³⁸ Furthermore, it might seem that in the chosen gauge the nuclear Hamiltonian is purely potential energy. This is not actually the case because it has been shown in refs 29 and 47 that $\epsilon_{bBO}(\mathbf{R}, t)$ contributes to the nuclear kinetic energy despite appearing as a scalar potential included in the expression of $\epsilon(\mathbf{R}, t)$.

With the aim of manipulating previous equations, we use the definition of the electronic wave function as $\Phi(\mathbf{r}, t; \mathbf{R}) = \frac{\Psi(\mathbf{r}, \mathbf{R}, t)}{|\chi(\mathbf{R}, t)|}$. This expression is inserted into the definitions of the $\epsilon_{BO}(\mathbf{R}, t)$ and $\epsilon_{GD}(\mathbf{R}, t)$ potentials. In addition, we write explicitly the expectation values as integrals over the electronic configuration space such that

$$ \epsilon_{BO}(\mathbf{R}, t) = \frac{1}{|\chi(\mathbf{R}, t)|^2} \int d\mathbf{r} \, \Psi^*(\mathbf{r}, \mathbf{R}, t) \hat{H}_{BO}(\mathbf{r}, \mathbf{R}) \Psi(\mathbf{r}, \mathbf{R}, t) \quad (31) $$

$$ \begin{aligned} \epsilon_{GD}(\mathbf{R}, t) &= \int d\mathbf{r} \frac{\Psi^*(\mathbf{r}, \mathbf{R}, t)}{|\chi(\mathbf{R}, t)|} \left(-i\hbar \frac{\partial}{\partial t}\right) \frac{\Psi(\mathbf{r}, \mathbf{R}, t)}{|\chi(\mathbf{R}, t)|} \\ &= \int d\mathbf{r} \frac{\Psi^*(\mathbf{r}, \mathbf{R}, t)}{|\chi(\mathbf{R}, t)|} (-i\hbar) \frac{\left(\frac{\partial}{\partial t}\Psi(\mathbf{r}, \mathbf{R}, t)\right)|\chi(\mathbf{R}, t)| - \Psi(\mathbf{r}, \mathbf{R}, t) \frac{\partial}{\partial t}|\chi(\mathbf{R}, t)|}{|\chi(\mathbf{R}, t)|^2} \\ &= \frac{1}{|\chi(\mathbf{R}, t)|^2} \int d\mathbf{r} \Psi^*(\mathbf{r}, \mathbf{R}, t) \left(-i\hbar \frac{\partial}{\partial t}\right) \Psi(\mathbf{r}, \mathbf{R}, t) \end{aligned} \quad (32) $$

The second term in $\epsilon_{GD}(\mathbf{R}, t)$, being purely imaginary, is clearly identically zero because the whole expression is purely real.

The contribution to the time-dependent scalar potential arising from the external field, i.e., $v_{int}(\mathbf{R}, t)$, can be rewritten as well as

$$ v_{int}(\mathbf{R}, t) = \frac{1}{|\chi(\mathbf{R}, t)|^2} \int d\mathbf{r} \, \Psi^*(\mathbf{r}, \mathbf{R}, t) \hat{V}(\mathbf{r}, \mathbf{R}, t) \Psi(\mathbf{r}, \mathbf{R}, t) \quad (33) $$

The “kinetic” part of the time-dependent scalar potential, $\epsilon_{bBO}(\mathbf{R}, t)$, becomes

$$ \begin{aligned} \epsilon_{bBO}(\mathbf{R}, t) &= \sum_{\nu} \frac{\hbar^2}{2M_\nu} \int d\mathbf{r} \, (\nabla_{\mathbf{R}_\nu} \Phi(\mathbf{r}, t; \mathbf{R}))^* (\nabla_{\mathbf{R}_\nu} \Phi(\mathbf{r}, t; \mathbf{R})) \\ &= \sum_{\nu} \frac{\hbar^2}{2M_\nu} \left[ \nabla_{\mathbf{R}_\nu} \int d\mathbf{r} \, \Phi^*(\mathbf{r}, t; \mathbf{R}) \nabla_{\mathbf{R}_\nu} \Phi(\mathbf{r}, t; \mathbf{R}) \right. \\ & \qquad \left. - \int d\mathbf{r} \, \Phi^*(\mathbf{r}, t; \mathbf{R}) \nabla_{\mathbf{R}_\nu}^2 \Phi(\mathbf{r}, t; \mathbf{R}) \right] \end{aligned} \quad (34) $$

Here, we used the chain rule for the derivative with respect to $\mathbf{R}_\nu$, namely, $\nabla_{\mathbf{R}_\nu}$. Applying once more the chain rule to the first term on the right-hand side of eq 34 yields

$$ \int d\mathbf{r} \, \Phi^*(\mathbf{r}, t; \mathbf{R}) \nabla_{\mathbf{R}_\nu} \Phi(\mathbf{r}, t; \mathbf{R}) = \nabla_{\mathbf{R}_\nu} \int d\mathbf{r} \, \Phi^*(\mathbf{r}, t; \mathbf{R}) \Phi(\mathbf{r}, t; \mathbf{R}) \\ - \int d\mathbf{r} \, (\nabla_{\mathbf{R}_\nu} \Phi(\mathbf{r}, t; \mathbf{R}))^* \Phi(\mathbf{r}, t; \mathbf{R}) \quad (35) $$

Using the partial normalization condition, $\int d\mathbf{r}|\Phi(\mathbf{r}, t; \mathbf{R})|^2 = 1$ $\forall \mathbf{R}$, $t$, to show that the first term on the right-hand side of eq 35 is identically zero implies that

$$ \begin{aligned} & \int d\mathbf{r} \, \Phi^*(\mathbf{r}, t; \mathbf{R}) \nabla_{\mathbf{R}_\nu} \Phi(\mathbf{r}, t; \mathbf{R}) \\ &= - \int d\mathbf{r} \, (\nabla_{\mathbf{R}_\nu} \Phi(\mathbf{r}, t; \mathbf{R}))^* \Phi(\mathbf{r}, t; \mathbf{R}) \end{aligned} \quad (36) $$

proving that eq 36 is purely imaginary. However, $\epsilon_{bBO}(\mathbf{R}, t)$ is real overall, which implies that such a purely imaginary contribution, from the first integral in its definition (eq 34), should be canceled out exactly by the imaginary part of its second term. As a result, we can rewrite eq 34 as

$$ \epsilon_{bBO}(\mathbf{R}, t) = -\Re \left[ \sum_{\nu} \frac{\hbar^2}{2M_\nu} \int d\mathbf{r} \, \Phi^*(\mathbf{r}, t; \mathbf{R}) \nabla_{\mathbf{R}_\nu}^2 \Phi(\mathbf{r}, t; \mathbf{R}) \right] \quad (37) $$

Summarizing the results obtained so far, we can rewrite the nuclear Hamiltonian as

$$ \begin{aligned} H_n(\mathbf{R}, t) &= \frac{1}{|\chi(\mathbf{R}, t)|^2} \int d\mathbf{r} \, \Psi^*(\mathbf{r}, \mathbf{R}, t) \left(-i\hbar \frac{\partial}{\partial t} + \hat{H}_{BO}(\mathbf{r}, \mathbf{R})\right) \\ & \quad + \hat{V}(\mathbf{r}, \mathbf{R}, t) \left(\Psi(\mathbf{r}, \mathbf{R}, t)\right) \\ & \quad - \Re \left[ \sum_{\nu} \frac{\hbar^2}{2M_\nu} \int d\mathbf{r} \, \Phi^*(\mathbf{r}, t; \mathbf{R}) \nabla_{\mathbf{R}_\nu}^2 \Phi(\mathbf{r}, t; \mathbf{R}) \right] \\ & \quad + Q_{\text{pot}}(\mathbf{R}, t) \end{aligned} \quad (38) $$

Note that the first integral of eq 38 can be simplified using the time-dependent Schrödinger equation. In fact, the operator in parentheses equals the negative of the nuclear kinetic energy, i.e., $(-i\hbar\frac{\partial}{\partial t} + \hat{H}_{BO}(\mathbf{r}, \mathbf{R}) + \hat{V}(\mathbf{r}, \mathbf{R}, t)) = -\hat{T}_n(\mathbf{R})$, and the whole integral is the negative of the expectation value of $\hat{T}_n(\mathbf{R})$ over the molecular wave function. Hence, we can now define the condition for (the nuclear gradient of) $H_n(\mathbf{R}, t)$ to be equal to zero:

$$ \begin{aligned} -\Re \left[ \sum_{\nu} \frac{\hbar^2}{2M_\nu} \int d\mathbf{r} \, \Phi^*(\mathbf{r}, t; \mathbf{R}) \nabla_{\mathbf{R}_\nu}^2 \Phi(\mathbf{r}, t; \mathbf{R}) \right] \\ \doteq \sum_{\nu} \frac{-\hbar^2}{2M_\nu} \int d\mathbf{r} \, \frac{\Psi^*(\mathbf{r}, \mathbf{R}, t) \nabla_{\mathbf{R}_\nu}^2 \Psi(\mathbf{r}, \mathbf{R}, t)}{|\chi(\mathbf{R}, t)|^2} \\ - Q_{\text{pot}}(\mathbf{R}, t) \end{aligned} \quad (39) $$

We use once again the definition of $\Phi(\mathbf{r}, t; \mathbf{R})$ in terms of the molecular and nuclear wave functions, and we apply $\nabla_{\mathbf{R}_\nu}^2$ on the left-hand side of eq 39 such that
$$
\begin{align}
& -\sum_{\nu} \frac{\hbar^2}{2M_{\nu}} \int dr \, \Phi^*(\mathbf{r}, t; \mathbf{R}) \nabla_{\mathbf{R}_{\nu}}^2 \Phi(\mathbf{r}, t; \mathbf{R}) \nonumber \\
& = -\sum_{\nu} \frac{\hbar^2}{2M_{\nu}} \int dr \frac{\Psi^*(\mathbf{r}, \mathbf{R}, t)}{|\chi(\mathbf{R}, t)|} \left( \frac{\nabla_{\mathbf{R}_{\nu}}^2 \Psi(\mathbf{r}, \mathbf{R}, t)}{|\chi(\mathbf{R}, t)|} \right. \nonumber \\
& \qquad \left. - \frac{2 \nabla_{\mathbf{R}_{\nu}} \Psi(\mathbf{r}, \mathbf{R}, t) \nabla_{\mathbf{R}} |\chi(\mathbf{R}, t)|}{|\chi(\mathbf{R}, t)|^2} \right. \nonumber \\
& \qquad \left. - \frac{\Psi(\mathbf{r}, \mathbf{R}, t) \nabla_{\mathbf{R}_{\nu}}^2 |\chi(\mathbf{R}, t)|}{|\chi(\mathbf{R}, t)|^2} + \frac{2 \Psi(\mathbf{r}, \mathbf{R}, t) (\nabla_{\mathbf{R}_{\nu}} |\chi(\mathbf{R}, t)|)^2}{|\chi(\mathbf{R}, t)|^3} \right) \nonumber \\
& = \sum_{\nu} \frac{\hbar^2}{2M_{\nu}} \left( - \int dr \frac{\Psi^*(\mathbf{r}, \mathbf{R}, t) \nabla_{\mathbf{R}_{\nu}}^2 \Psi(\mathbf{r}, \mathbf{R}, t)}{|\chi(\mathbf{R}, t)|^2} \right. \nonumber \\
& \qquad \left. + \frac{\nabla_{\mathbf{R}_{\nu}}^2 |\chi(\mathbf{R}, t)|}{|\chi(\mathbf{R}, t)|} \right. \nonumber \\
& \qquad \left. + 2 \int dr \frac{\Psi^*(\mathbf{r}, \mathbf{R}, t) \nabla_{\mathbf{R}_{\nu}} \Psi(\mathbf{r}, \mathbf{R}, t)}{|\chi(\mathbf{R}, t)|^2} \frac{\nabla_{\mathbf{R}_{\nu}} |\chi(\mathbf{R}, t)|}{|\chi(\mathbf{R}, t)|} \right. \nonumber \\
& \qquad \left. - 2 \frac{(\nabla_{\mathbf{R}_{\nu}} |\chi(\mathbf{R}, t)|)^2}{|\chi(\mathbf{R}, t)|^2} \right) \tag{40}
\end{align}
$$

Therefore, we see that the condition given by eq 39 would be
fulfilled if, by taking the real part of eq 40, the last two terms
would cancel out.

To show that, let us investigate the real part of the
contribution depending on the full molecular wave function,
where we use $\Psi(\mathbf{r}, \mathbf{R}, t) = \Phi(\mathbf{r}, t; \mathbf{R})|\chi(\mathbf{R}, t)|$, namely,

$$
\begin{align}
& \mathcal{R} \left[ \int dr \frac{\Psi^*(\mathbf{r}, \mathbf{R}, t) \nabla_{\mathbf{R}} \Psi(\mathbf{r}, \mathbf{R}, t)}{|\chi(\mathbf{R}, t)|} \right. \nonumber \\
& \qquad \left. + \mathcal{R} \left[ \int dr \Phi^*(\mathbf{r}, t; \mathbf{R}) \nabla_{\mathbf{R}} \Phi(\mathbf{r}, t; \mathbf{R}) \right] \right. \nonumber \\
& \qquad \left. + \mathcal{R} \left[ \int dr \Phi^*(\mathbf{r}, t; \mathbf{R}) \Phi(\mathbf{r}, t; \mathbf{R}) \frac{\nabla_{\mathbf{R}} |\chi(\mathbf{R}, t)|}{|\chi(\mathbf{R}, t)|} \right] \right] \tag{41}
\end{align}
$$

The first term on the right-hand side of eq 41 is identically zero because it is purely imaginary (eq 36); the integral in the second term yields 1 thanks to the partial normalization condition. As a result, we have

$$
\mathcal{R} \left[ \int dr \frac{\Psi^*(\mathbf{r}, \mathbf{R}, t) \nabla_{\mathbf{R}} \Psi(\mathbf{r}, \mathbf{R}, t)}{|\chi(\mathbf{R}, t)|^2} \frac{\nabla_{\mathbf{R}} |\chi(\mathbf{R}, t)|}{|\chi(\mathbf{R}, t)|} \right] = \frac{(\nabla_{\mathbf{R}} |\chi(\mathbf{R}, t)|)^2}{|\chi(\mathbf{R}, t)|^2} \quad (42)
$$

proving that when taking the real part of the last two terms of eq 40, they exactly cancel out. Therefore, the condition posed by eq 39 is fulfilled, as stated above.

Finally we can show that the nuclear Hamiltonian of eq 38 is
zero in the chosen gauge

$$
H_n(\mathbf{R}, t) = \frac{\langle \Psi(\mathbf{R}, t) | -i\hbar \frac{d}{\partial \mathbf{r}} + \hat{H}_{\text{BO}}(\mathbf{R}) + \hat{V}(\mathbf{R}, t) + \hat{T}_n(\mathbf{R}) | \Psi(\mathbf{R}, t) \rangle_{\mathbf{r}}}{|\chi(\mathbf{R}, t)|^2} \\
- Q_{\text{pot}}(\mathbf{R}, t) + Q_{\text{pot}}(\mathbf{R}, t) = 0 \qquad (43)
$$

meaning that its gradient is zero as well.

■ ASSOCIATED CONTENT

Supporting Information

The Supporting Information is available free of charge at
https://pubs.acs.org/doi/10.1021/acs.jpca.1c09604.

Movie of the classical and quantum trajectories without
invoking the Condon approximation (MP4)

Movie of the classical and quantum trajectories with
invoking the Condon approximation (MP4)

Movie of the classical and quantum trajectories without invoking the Condon approximation (MP4)

Movie of the classical and quantum trajectories with
invoking the Condon approximation (MP4)

Description of the movies (PDF)

■ AUTHOR INFORMATION

Corresponding Authors

Basile F. E. Curchod - Department of Chemistry, Durham University, Durham DH1 3LE, United Kingdom;

orcid.org/0000-0002-1705-473X;

Email: basile.f.curchod@durham.ac.uk

Federica Agostini - Université Paris-Saclay, CNRS, Institut de Chimie Physique UMR8000, 91405 Orsay, France;

orcid.org/0000-0003-2951-4964;

Email: federica.agostini@universite-paris-saclay.fr

Author

Lea M. Ibele - Department of Chemistry, Durham University,
Durham DH1 3LE, United Kingdom

Complete contact information is available at:
https://pubs.acs.org/10.1021/acs.jpca.1c09604

Notes

The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement no. 803718, project SINDAM). L.M.I. acknowledges the EPSRC for an EPSRC Doctoral Studentship (EP/R513039/1). The authors acknowledge financial support from CNRS via the International Emerging Action.

■ REFERENCES

(1) Born, M.; Oppenheimer, R. Zur Quantentheorie der Molekeln. *Ann. Phys.* **1927**, 389, 457–484.

(2) Tully, J. C. Perspective on Zur Quantentheorie der Molekeln. *Theoretical Chemistry Accounts*; Springer: Berlin, 2000; Vol. 103, pp 173–176.

(3) Zewail, A. H. Laser Femtochemistry. *Science* **1988**, 242, 1645–1653.

(4) Kobayashi, Y.; Chang, K. F.; Zeng, T.; Neumark, D. M.; Leone, S. R. Direct mapping of curve-crossing dynamics in IBr by attosecond transient absorption spectroscopy. *Science* **2019**, 365, 79–83.

(5) Ibele, L. M.; Curchod, B. F. E. A molecular perspective on Tully models for nonadiabatic dynamics. *Phys. Chem. Chem. Phys.* **2020**, 22, 15183–15196.

(6) Rose, T. S.; Rosker, M. J.; Zewail, A. H. Femtosecond real-time probing of reactions. IV. The reactions of alkali halides. *J. Chem. Phys.* **1989**, 91, 7415–7436.

(7) Hahn, S.; Stock, G. Quantum-mechanical modeling of the femtosecond isomerization in rhodopsin. *J. Phys. Chem. B* **2000**, 104, 1146–1149.
(8) Böckmann, M.; Doltsinis, N. L.; Marx, D. Unraveling a Chemically Enhanced Photoswitch: Bridged Azobenzene. *Angew. Chem., Int. Ed.* **2010**, *49*, 3382–3384.

(9) Domcke, W.; Yarkony, D. R. Role of conical intersections in molecular spectroscopy and photoinduced chemical dynamics. *Annu. Rev. Phys. Chem.* **2012**, *63*, 325–352.

(10) Agostini, F.; Curchod, B. F. E. Different flavors of nonadiabatic molecular dynamics. *WIREs Comput. Mol. Sci.* **2019**, *9*, No. e1417.

(11) Born, M.; Huang, K. *Dynamical Theory of Crystal Lattices*; Clarendon: Oxford, 1954.

(12) Baer, M. *Beyond Born-Oppenheimer: Electronic Nonadiabatic Coupling Terms and Conical Intersections*; John Wiley & Sons, Inc.: Hoboken, NJ, 2006.

(13) Yarkony, D. R. Diabolical conical intersections. *Rev. Mod. Phys.* **1996**, *68*, 985.

(14) Domcke, W.; Yarkony, D.; Köppel, H., Eds.; *Conical Intersections: Electronic Structure, Dynamics & Spectroscopy*; World Scientific Publishing Company Inc.: Singapore, 2004; Vol. 15.

(15) Berry, M. V.; Wilkinson, M. Diabolical points in the spectra of triangles. *Proc. R. Soc. London Ser. A* **1984**, 15–43.

(16) Mead, C. A.; Truhlar, D. G. On the determination of Born-Oppenheimer nuclear motion wave functions including complications due to conical intersections and identical nuclei. *J. Chem. Phys.* **1979**, *70*, 2284–2296.

(17) Gueye, M.; Paolino, M.; Gindensperger, E.; Haacke, S.; Olivucci, M.; Léonard, J. Vibrational coherence and quantum yield of retinal-chromophore-inspired molecular switches. *Faraday Discuss.* **2020**, *221*, 299–321.

(18) Cattaneo, P.; Persico, M. Diabatic and adiabatic potential-energy surfaces for azomethane photochemistry. *Theor. Chem. Acc.* **2000**, *103*, 390–398.

(19) Cusati, T.; Granucci, G.; Persico, M. Photodynamics and Time-Resolved Fluorescence of Azobenzene in Solution: A Mixed Quantum-Classical Simulation. *J. Am. Chem. Soc.* **2011**, *133*, 5109–5123.

(20) Fregoni, J.; Granucci, G.; Persico, M.; Corni, S. Strong Coupling with Light Enhances the Photoisomerization Quantum Yield of Azobenzene. *Chem.* **2020**, *6*, 250–265.

(21) Levine, B. G.; Martínez, T. J. Isomerization through conical intersections. *Annu. Rev. Phys. Chem.* **2007**, *58*, 613–634.

(22) Mignolet, B.; Curchod, B. F. E.; Martínez, T. J. Rich Athermal Ground-State Chemistry Triggered by Dynamics through a Conical Intersection. *Angew. Chem.* **2016**, *128*, 15217–15220.

(23) Pathak, S.; Ibele, L. M.; Boll, R.; Callegari, C.; Demidovich, A.; Erk, B.; Feifel, R.; Forbes, R.; Di Fraia, M.; Giannessi, L.; et al. Tracking the ultraviolet-induced photochemistry of thiophenone during and after ultrafast ring opening. *Nat. Chem.* **2020**, *12*, 795–800.

(24) Crespo-Otero, R.; Barbatti, M. Recent Advances and Perspectives on Nonadiabatic Mixed Quantum-Classical Dynamics. *Chem. Rev.* **2018**, *118*, 7026–7068.

(25) Curchod, B. F. E.; Martínez, T. J. Ab initio nonadiabatic quantum molecular dynamics. *Chem. Rev.* **2018**, *118*, 3305–3336.

(26) Agostini, F.; Gross, E. K. U. Ultrafast dynamics with the exact factorization. *Eur. Phys. J. B* **2021**, *94*, 179.

(27) González, L.; Lindh, R. *Quantum Chemistry and Dynamics of Excited States: Methods and Applications*; John Wiley & Sons: Chichester, 2020.

(28) Abedi, A.; Maitra, N.; Gross, E. Exact Factorization of the Time-Dependent Electron-Nuclear Wave Function. *Phys. Rev. Lett.* **2010**, *105*, 123002.

(29) Abedi, A.; Maitra, N. T.; Gross, E. Correlated electron-nuclear dynamics: Exact factorization of the molecular wavefunction. *J. Chem. Phys.* **2012**, *137*, 22A530.

(30) Curchod, B. F. E.; Agostini, F. On the dynamics through a conical intersection. *J. Phys. Chem. Lett.* **2017**, *8*, 831–837.

(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.

(32) Min, S. K.; Abedi, A.; Kim, K. S.; Gross, E. K. U. Is the Molecular Berry Phase an Artifact of the Born-Oppenheimer Approximation? *Phys. Rev. Lett.* **2014**, *113*, 263004.

(33) Requist, R.; Tandetzky, F.; Gross, E. Molecular geometric phase from the exact electron-nuclear factorization. *Phys. Rev. A* **2016**, *93*, 042108.

(34) Agostini, F.; Curchod, B. F. E. When the exact factorization meets conical intersections. *Eur. Phys. J. B* **2018**, *91*, 141.

(35) 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.

(36) Schild, A. Electronic quantum trajectories with quantum nuclei. *arXiv:2109.13632 [physics.chem-ph]*, 2021, Vol. 54.

(37) Suzuki, Y.; Watanabe, K. Bohmian mechanics in the exact factorization of electron-nuclear wave functions. *Phys. Rev. A* **2016**, *94*, 032517.

(38) Agostini, F.; Abedi, A.; Suzuki, Y.; Min, S. K.; Maitra, N. T.; Gross, E. The exact forces on classical nuclei in non-adiabatic charge transfer. *J. Chem. Phys.* **2015**, *142*, 084303.

(39) 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.

(40) Ferretti, A.; Granucci, G.; Lami, A.; Persico, M.; Villani, G. Quantum mechanical and semiclassical dynamics at a conical intersection. *J. Chem. Phys.* **1996**, *104*, 5517.

(41) Mignolet, B.; Curchod, B. F. E. Excited-State Molecular Dynamics Triggered by Light Pulses—Ab Initio Multiple Spawning vs Trajectory Surface Hopping. *J. Phys. Chem. A* **2019**, *123*, 3582–3591.

(42) Feit, M. D.; Fleck, J. A.; Steiger, A. Solution of the Schrödinger equation by a spectral method. *J. Comput. Phys.* **1982**, *47*, 412–433.

(43) Agostini, F.; Tavernelli, I.; Ciccotti, G. Nuclear quantum effects in electronic (non) adiabatic dynamics. *Eur. Phys. J. B* **2018**, *91*, 139.

(44) Wyatt, R. E. *Quantum Dynamics with Trajectories: Introduction to Quantum Hydrodynamics*; Interdisciplinary Applied Mathematics; Springer: New York, 2005.

(45) Oriols, X.; Martín, F.; Suñé, J. Implications of noncrossing properties of Bohm trajectories in one-dimensional tunneling configurations. *Phys. Rev. A* **1996**, *54*, 2594–2604.

(46) Holland, P. R. *The Quantum Theory of Motion: An Account of the de Broglie-Bohm Causal Interpretation of Quantum Mechanics*; Cambridge University Press, 1993.

(47) Agostini, F.; Abedi, A.; Gross, E. K. U. Classical nuclear motion coupled to electronic non-adiabatic transitions. *J. Chem. Phys.* **2014**, *141*, 214101.