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

HAL
open science

# Exploring Exact-Factorization-Based Trajectories for Low-Energy Dynamics near a Conical Intersection

Lea M. Ibele, Federica Agostini

► To cite this version:

Lea M. Ibele, Federica Agostini. Exploring Exact-Factorization-Based Trajectories for Low-Energy Dynamics near a Conical Intersection. Journal of Physical Chemistry A, 2024, 128 (18), pp.3672-3684. 10.1021/acs.jpca.4c00555. hal-05287018

HAL Id: hal-05287018

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

Submitted on 27 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Exploring exact-factorization-based
trajectories for low-energy dynamics near a
conical intersection

Lea M. Ibele* and Federica Agostini

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

E-mail: lea-maria.ibele@universite-paris-saclay.fr

## Abstract

We study low-energy dynamics generated by a two-dimensional two-state Jahn-Teller Hamiltonian in the vicinity of a conical intersection using quantum wavepacket and trajectories dynamics. Recently, these dynamics were studied by comparing the adiabatic representation and the exact factorization, with the purpose to highlight the different nature of topological- and geometric-phase effects arising in the two theoretical representation of the same problem. Here, we employ the exact factorization to understand how to model accurately low-energy dynamics in the vicinity of a conical intersection using an approximate description of the nuclear motion that uses trajectories. We find that, since nonadiabatic effects are weak but non-negligible, the trajectory-based description that invokes the classical approximation struggles to capture the correct behavior.
# 1 Introduction

Conical intersections, points of degeneracy between two adiabatic potential energy surfaces are ubiquitous in molecular physics and photochemistry. Due to the strong nonadiabatic electron-nuclear coupling in their vicinity, they act as funnels for electronically excited population to decay to its ground state and therefore strongly influence the photoinduced processes in molecules.¹⁻⁶ If those degeneracies occur between the excited state and the ground state, they can also influence low-energy processes and induce a nonadiabatic behavior of the dynamics of a system evolving mainly in the ground state.⁷⁻¹¹ In the past, such influences have been postulated and discussed for transition states where higher lying electronically excited states are involved.¹²,¹³ As conical intersections represent points of degeneracies between adiabatic states and, thus, of a singular nonadiabatic coupling, it has been postulated that they can induce a *molecular geometric phase* on the dynamics.¹⁴⁻²⁴ This means that a nuclear wavepacket picks up a phase when encircling a conical intersection, and if parts of the nuclear wavepacket pass the conical intersection on different sides, they acquire opposite phases leading to destructive interference upon recombination.

This discussion intrinsically assumes an adiabatic representation of the photochemical or low-energy process. In this representation, the molecular wavefunction is expressed in a basis of adiabatic electronic states, obtained from the electronic Hamiltonian, with the corresponding nuclear coefficients. Since the electronic Hamiltonian depends on the nuclear positions, the adiabatic states and the corresponding eigenvalues, i.e. the adiabatic potential energy surfaces, acquire a characteristic dependence on the nuclear positions. This dependence produces degeneracies, singularities and, thus, conical intersections, which become signatures of the adiabatic representation.

Recently, an experiment mimicking a nuclear wavepacket encircling a conical intersection has been performed for a two-dimensional Jahn-Teller Hamiltonian using a mixed qudit-boson simulator.²⁵ Destructive interference between the wavepacket parts has been reported leading to a nodal-like feature in the nuclear density which has been proposed as a direct
experimental observation of the geometric phase. We studied<sup>26</sup> the same system through numerical quantum dynamics simulations and analyzed the dynamics comparing the adiabatic representation with the exact factorization. Our study showed that in the adiabatic representation the geometric phase is actually topological, as it is well-known in the presence of conical intersections.<sup>14,27</sup> However, such a topological feature is not robust under changes of representation as it disappears when the exact factorization is employed for the analysis. A physical observable, coined *dynamics-induced geometric phase*, can be introduced instead, independent of the representation which is, therefore, unrelated to conical intersections.<sup>26,28,29</sup>

The exact factorization is a representation of the molecular wavefunction that writes it as a single product of time-dependent electronic and nuclear wavefunction. This formalism<sup>30-33</sup> gives rise to single, time-dependent scalar and vector potentials governing the nuclear wavefunction evolution – “single” and time-dependent, as opposed to multiple static potential energy surfaces and nonadiabatic couplings arising in the more commonly-used adiabatic representation. Thus, it has been shown in the past that these potentials allow one to introduce in a natural way the concept of classical-like force to drive the nuclear dynamics and, thus, are well-suited for the development of trajectory-based algorithms.<sup>34-37</sup> In this work, our aim is to investigate the time-dependent scalar and vector potentials of the low-energy dynamics around a conical intersection of Refs. [25,26], in relation to various trajectory schemes based on the exact factorization. Such an analysis will allow us to understand if and how one can circumvent including the geometric phase of the adiabatic representation when an approximate description of the dynamics using nuclear trajectories is employed.

In this paper, Section 2 is devoted to the presentation of the exact factorization and of
the trajectory schemes that will be used for the analysis of the dynamics. Section 3 intro-
duces the model Hamiltonian and provides information on the simulations. In Section 4,
we report our results and our interpretation of the trajectories dynamics. Using the time-
dependent potentials obtained from the numerically exact quantum dynamics, we first look
at quantum trajectories, then we study the influence of the classical approximation on the trajectories' dynamics, in order to understand the effect of the different terms of the time-dependent potentials. Based on this analysis, we investigate why the coupled-trajectory mixed-quantum-classical (CT-MQC) algorithm$^{36,38,39}$ fails to reproduce the correct dynamics and we compare it to adiabatic-like dynamics using, as done previously, quantum and classical trajectories. We conclude in Section 5 by summarizing the shortcomings of the approximate CT-MQC scheme for such a dynamical problem.

## 2 Methods

The time evolution of the electron-nuclear wavefunction describing the state of a molecule, $\Psi(r, R, t)$ is governed by the time-dependent Schrödinger equation

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

with the full, molecular Hamiltonian $\hat{H}(r, R)$ consisting of the sum of the nuclear kinetic energy operator $\hat{T}_n(R)$ and the electronic, Born-Oppenheimer Hamiltonian $\hat{H}_{BO}(r, R)$, consisting of the electronic kinetic energy and of the interaction potential. $r$ and $R$ are collective variables for the electronic and nuclear coordinates of the molecule, respectively.

### 2.1 Exact Factorization of the Molecular Wavefunction

The exact factorization$^{30-33}$ proposes to express the molecular wavefunction as a single product of time-dependent nuclear $\chi(R, t)$ and electronic $\Phi_R(r, t)$ wavefunctions

$$\Psi(r, R, t) = \Phi_R(r, t)\chi(R, t) \quad (2)$$

It is important to note that the electronic wavefunction depends parametrically on the nuclear coordinates, denoted by the subscript $R$. Once the partial normalization condition is
imposed, i.e. $\langle \Phi_R(t) | \Phi_R(t) \rangle_r = 1 \forall R, t$ (the notation $\langle \cdot \rangle_r$ signifies integration over the whole space of electronic coordinates), it can be proven that the single product form (2) is exact and unique up to a gauge-like transformation.³⁰,³³ In the current work we fix the gauge so that the nuclear wavefunction is real and non-negative, i.e. $\chi(R, t) = |\chi(R, t)| \forall R, t$. This choice of gauge allows for nuclear wavefunctions that are reminiscent of vibrational excited states, i.e., with nodes, because they yield (gauge-independent) nuclear probability densities that can be zero.

The time-evolution of the molecular wavefunction can be obtained within the exact factorization upon inserting the factored form (2) into the time-dependent Schrödinger equation (1). Using the partial normalization condition, we can derive coupled equations of motion for the nuclear and electronic wavefunctions:

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

$$i\hbar \frac{\partial}{\partial t} \Phi_R(\mathbf{r}, t) = \left[ \hat{H}_{\text{BO}} + \hat{U}_{\text{en}}[\Phi, \chi](\mathbf{R}, t) - \epsilon(\mathbf{R}, t) \right] \Phi_R(\mathbf{r}, t) \quad (4)$$

where $\mu$ denotes the $N_n$ nuclei. The evolution equation for the nuclear wavefunction couples to the electronic wavefunction implicitly through the two new potentials, the time-dependent vector potential (TDVP) $A_{\mu}(\mathbf{R}, t)$,

$$A_{\mu}(\mathbf{R}, t) = \langle \Phi_{\mathbf{R}}(t) | -i\hbar\nabla_{\mu}\Phi_{\mathbf{R}}(t) \rangle_{\mathbf{r}} \quad (5)$$

and the time-dependent potential energy surface (TDPES) $\epsilon(\mathbf{R}, t)$

$$\epsilon(\mathbf{R}, t) = \langle \Phi_R(t) | \hat{H}_{\text{BO}} + \hat{U}_{\text{en}}[\Phi, \chi](\mathbf{R}, t) - i\hbar \frac{\partial}{\partial t} | \Phi_R(t) \rangle_r \quad (6)$$

The electronic evolution equation is coupled to the nuclear wavefunction through the electron-nuclear coupling operator, $\hat{U}_{\text{en}}[\Phi, \chi](\mathbf{R}, t)$, which explicitly depends on the nuclear wavefunc-
tion and also implicitly contains the electronic wavefunction through the TDVP:

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

Note that for conciseness the $(\mathbf{R}, t)$ dependencies of $\chi$ and $\mathbf{A}_{\mu}$ have been dropped.

As introduced above, the exact factorization of the molecular wavefunction is unique up to a gauge transformation,

$$ \Psi(\mathbf{r}, \mathbf{R}, t) = \Phi_{\mathbf{R}}(\mathbf{r}, t)\chi(\mathbf{R}, t) = \Phi_{\mathbf{R}}(\mathbf{r}, t)e^{-i/\hbar\theta(\mathbf{R}, t)}\chi(\mathbf{R}, t)e^{i/\hbar\theta(\mathbf{R}, t)} = \tilde{\Phi}_{\mathbf{R}}(\mathbf{r}, t)\tilde{\chi}(\mathbf{R}, t) \quad (8) $$

for a real gauge function $\theta(\mathbf{R}, t)$. While the nuclear and electronic equations of motion, Eqs. (3) and (4), respectively, are form-invariant under a gauge transformation, the TDVP and TDPES transform according to

$$ \tilde{\epsilon}(\mathbf{R}, t) = \epsilon(\mathbf{R}, t) + \frac{\partial}{\partial t}\theta(\mathbf{R}, t) \quad (9) $$

$$ \tilde{\mathbf{A}}_{\mu}(\mathbf{R}, t) = \mathbf{A}_{\mu}(\mathbf{R}, t) + \nabla_{\mu}\theta(\mathbf{R}, t) \quad (10) $$

The TDVP is inherently gauge dependent and the TDPES can be decomposed into two gauge-invariant $(\epsilon_{\text{BO}}, \epsilon_{\text{geo}})$ and one gauge-dependent contributions $(\epsilon_{\text{GD}})$:

$$ \epsilon_{\text{BO}} = \langle \Phi_{\mathbf{R}}(t) | \hat{H}_{\text{BO}}(\mathbf{R}) | \Phi_{\mathbf{R}}(t) \rangle_{\mathbf{r}} \quad (11) $$

$$ \epsilon_{\text{geo}} = \sum_{\mu} \frac{\hbar^2 \langle \nabla_{\mu} \Phi_{\mathbf{R}}(t) | \nabla_{\mu} \Phi_{\mathbf{R}}(t) \rangle_{\mathbf{r}}}{2M_{\mu}} - \sum_{\mu} \frac{|\mathbf{A}_{\mu}|^2}{2M_{\mu}} \quad (12) $$

$$ \epsilon_{\text{GD}} = \langle \Phi_{\mathbf{R}}(t) | -i\hbar \frac{\partial}{\partial t} | \Phi_{\mathbf{R}}(t) \rangle_{\mathbf{r}} \quad (13) $$

The first of the two gauge-invariant contributions, $\epsilon_{\text{BO}}$, arises from the Born-Oppenheimer Hamiltonian and yields a mean-field-like contribution to the TDPES. The second term, $\epsilon_{\text{geo}}$, accounts for the spatial variation of the electronic wavefunction in the parametric space $\mathbf{R}$,
as it is evident from its definition (12). Therefore, it accounts for nonadiabatic effects in the TDPES, which become important when $\Phi_R(r, t)$ changes abruptly as function of $R$, as for instance – but not only – in the vicinity of the avoided crossings or the conical intersections of the adiabatic representation. we interpret it as somehow responsible for nonadiabatic effects in the TDPES. In addition, if an analogy is drawn between TDPES/TDVP and classical electromagnetism, it has been argued that the gradient of $\epsilon_{geo}$ gives rise to a geometric contribution to the electric-field-like effect of the electrons acting on the nuclei.²⁹ Note that, in the literature,⁴⁰ an extended divergence relation between the (first-order) nonadiabatic coupling vectors and the (second-order) nonadiabatic couplings has been identified in a adiabatic representation, namely $-\langle\phi_R^J|\nabla_\mu^2\phi_R^J\rangle_r = \sum_{K \neq J} |\langle\phi_R^K|\nabla_\mu\phi_R^J\rangle_r|^2$, leading to $-\nabla_\mu \cdot \langle\phi_R^J|\nabla_\mu\phi_R^J\rangle_r + \langle\nabla_\mu\phi_R^J|\nabla_\mu\phi_R^J\rangle_r = \sum_{K \neq J} |\langle\phi_R^K|\nabla_\mu\phi_R^J\rangle_r|^2$ and to $\langle\nabla_\mu\phi_R^J|\nabla_\mu\phi_R^J\rangle_r = \sum_{K \neq J} |\langle\phi_R^K|\nabla_\mu\phi_R^J\rangle_r|^2$ in the case of real adiabatic eigenstates $|\phi_R^J\rangle$ of the electronic Hamiltonian $\hat{H}_{BO}$. This last relation looks surprisingly similar to Eq. (12), even though it is not directly applicable to the time-dependent electronic wavefunction of the exact factorization, since Eq. (12) is non-zero in general.

## 2.2 Trajectory formalisms from the exact factorization

As ensured by the partial normalization condition, at all times $|\chi(R, t)|^2$ reproduces the nuclear density obtained from $\Psi(r, R, t)$. Thus, we can equally express the nuclear wavefunction as

$$ \chi(R, t) = \sqrt{\langle\Psi(R, t)|\Psi(R, t)\rangle_r} \exp\left(\frac{i}{\hbar}S(R, t)\right) = |\chi(R, t)| \exp\left(\frac{i}{\hbar}S(R, t)\right) \quad (14) $$

Inserting this polar form of the nuclear wavefunction into the time-dependent Schrödinger equation allows us to derive an evolution equation for the phase, as shown in Refs. [34,35]. This can be identified as a nuclear Hamilton-Jacobi equation using as canonical momentum
$\nabla_{\mu}S(\mathbf{R}, t) = \tilde{\mathbf{P}}_{\mu}(\mathbf{R}, t):$

$$ -\frac{\partial}{\partial t} S(\mathbf{R}, t) = H_n(\tilde{\mathbf{P}}, \mathbf{R}, t) = \sum_{\mu} \frac{[\tilde{\mathbf{P}}_{\mu}(\mathbf{R}, t) + \mathbf{A}_{\mu}(\mathbf{R}, t)]^2}{2M_{\mu}} + \epsilon(\mathbf{R}, t) + Q_{\text{pot}}(\mathbf{R}, t) \quad (15) $$

The quantum potential, $Q_{\text{pot}}(\mathbf{R}, t) = \sum_{\mu} \frac{-\hbar^2}{2M_{\mu}} \frac{\nabla_{\mu}^2 |\chi(\mathbf{R}, t)|}{|\chi(\mathbf{R}, t)|}$, introduces most of the nuclear quantum effects, while the other contributing terms are the nuclear kinetic energy with the nuclear mechanical momentum field $\mathbf{P}_{\mu}(\mathbf{R}, t) = \tilde{\mathbf{P}}_{\mu}(\mathbf{R}, t) + \mathbf{A}_{\mu}(\mathbf{R}, t)$ and the “classical” potential energy $\epsilon(\mathbf{R}, t)$.

Using the method of characteristics to solve this Hamilton-Jacobi partial differential equations yields Hamilton-like evolution equations,

$$ \dot{\mathbf{R}}_{\mu}(t) = \frac{\tilde{\mathbf{P}}_{\mu}(t) + \mathbf{A}_{\mu}(\mathbf{R}(t), t)}{M_{\mu}} \quad (16) $$

$$ \dot{\tilde{\mathbf{P}}}_{\mu}(t) = -\nabla_{\mu} H_n(\tilde{\mathbf{P}}(t), \mathbf{R}(t), t) \quad (17) $$

where the canonical momentum appears, which are equivalent to

$$ \dot{\mathbf{R}}_{\mu}(t) = \frac{\mathbf{P}_{\mu}(t)}{M_{\mu}} \quad (18) $$

$$ \dot{\mathbf{P}}_{\mu}(t) = -\nabla_{\mu} H_n(\mathbf{P}(t), \mathbf{R}(t), t) + \dot{\mathbf{A}}_{\mu}(\mathbf{R}(t), t) \quad (19) $$

using the mechanical momentum. Introducing the characteristics, i.e. a set of ordinary differential equations that always satisfy the original Hamilton-Jacobi equation, allows us to define the trajectories $\mathbf{R}(t) \equiv \{\mathbf{R}_{\mu}(t)\}_{\mu=1,...,N_n}$ and their corresponding momenta $\mathbf{P}(t) \equiv \{\mathbf{P}_{\mu}(t)\}_{\mu=1,...,N_n}$.

In this work, we chose a gauge where the nuclear wavefunction is real and non-negative,
thus, $S(\mathbf{R}, t) = 0$ which yields the evolution equations

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

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

Equation (21) clearly shows that the nuclear momentum field is the TDVP. These evolution equations have been derived by simply introducing the characteristics to solve the nuclear Hamilton-Jacobi equation, Eq. (15), without invoking the classical approximation yet. Consequently, the “trajectories” emerging as solutions of the evolution equations (21) and (20) are *quantum* trajectories. For these quantum trajectories the TDVP is the only quantity driving the dynamics.³⁷ Note that, in general, when working using the exact factorization, our choices of gauges are usually made with the scope to either convey a direct interpretation of the analyzed quantities, e.g., the TDVP in our gauge fully encodes the nuclear momentum field, or simplify the numerical calculations, e.g., the TDVP only will be used to propagate the quantum trajectories in the numerical studied reported in Section 4.

Once we apply the classical approximation, i.e. neglecting the quantum potential in Eq. (15), the evolution equations become

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

$$ \dot{\mathbf{P}}_{\mu}(t) = -\nabla_{\mu} H_{\mathrm{n}}^{\mathrm{cl}}(\mathbf{P}(t), \mathbf{R}(t), t) + \dot{\mathbf{A}}_{\mu}(\mathbf{R}(t), t) \qquad (23) $$

with classical Hamiltonian

$$ H_{\mathrm{n}}^{\mathrm{cl}}(\mathbf{P}(t), \mathbf{R}(t), t) = \sum_{\mu} \frac{\mathbf{P}_{\mu}(\mathbf{R}, t)^2}{2M_{\mu}} + \epsilon(\mathbf{R}, t) \qquad (24) $$

Having found evolution equations for both quantum and classical trajectories, we only need to discuss how the initial conditions are selected. For this, it is important to note that for
the quantum trajectories, position and momenta are not independent variables. Therefore, the initial positions, for both quantum and classical trajectories, can be sampled stochastically from the initial nuclear probability distribution. While for the classical trajectories the momenta can equally be stochastically and independently sampled from a momentum probability distribution, e.g. calculated from the Wigner transform of the initial nuclear wavefunction, for the quantum trajectories the initial momenta need to be selected as $A_μ(R_0, t = 0)$ for an ensemble of initial positions $R_0$.

Analyzing the nuclear dynamics using these quantum and classical trajectories allows us to focus purely on the accuracy of trajectories in capturing the simulated dynamics. The potential driving the dynamics are numerically exact, as they are determined from the solution of the full time-dependent Schrödinger equation. The exact factorization is a suitable formalism for such an analysis because it is not based on the adiabatic representation and inherently couples the nuclear dynamics on different electronic states.

## 2.3 Coupled-trajectory mixed quantum-classical approach

Having introduced trajectory formalisms to propagate nuclear trajectories using the exact TDVP and TDPES, we briefly introduce in this section a formalism that approximates the TDVP and TDPES from the electronic wavefunction, which allows one to simulate nonadiabatic dynamics on-the-fly: the so-called Coupled-Trajectory Mixed Quantum-Classical (CT-MQC) algorithm. We outline here just the main ideas of this formalism and refer the interested reader to more in depth descriptions in literature.³⁶,⁴¹

CT-MQC has been derived³⁸ starting from the evolution equation of the phase of the nuclear wavefunction (Eq. (15)) following the procedure detailed in Section 2.2 to derive Eqs. (20) and (21). However, in CT-MQC the gauge is fixed so that $\epsilon(R(t), t) + \sum_{\mu'} \dot{R}_{\mu'}(t) \cdot A_{\mu'}(R(t), t) = 0$. To evaluate the TDPES $\epsilon(R(t), t)$ and TDVP $A_μ(R(t), t)$ along each trajectory $R(t)$, one has to solve the electronic evolution equation, Eq. (4), along the same trajectory.
In order to make the CT-MQC formalism applicable for molecules, the electronic wave-
function $\Phi_{\mathbf{R}(t)}(\mathbf{r}, t)$ is expressed in a basis of adiabatic electronic eigenfunctions $\phi_{\mathbf{R}(t)}^J(\mathbf{r})$, i.e.
a Born-Huang-like expansion

$$
\Phi_{\mathbf{R}(t)}(\mathbf{r}, t) = \sum_{J} C_J(\mathbf{R}(t), t) \phi_{\mathbf{R}(t)}^{J}(\mathbf{r}) \quad (25)
$$

The time-dependence is collected in the coefficients so that the electronic eigenfunctions, $\phi_{\mathbf{R}(t)}^J(\mathbf{r})$, only contain an implicit time-dependence through the parametric dependence on the nuclear positions. The adiabatic potential energy surfaces $\epsilon_{\text{BO}}^J(\mathbf{R}(t))$, i.e. the energy eigenvalues, as well as the nonadiabatic coupling vectors $\mathbf{d}_{\mu,JI}(\mathbf{R}(t)) = \langle \phi_{\mathbf{R}(t)}^J | \nabla_\mu \phi_{\mathbf{R}(t)}^I \rangle_\mathbf{r}$ are introduced. In general, these quantities shall be obtained through electronic-structure calculations performed on-the-fly: they are determined on a specific nuclear configuration encountered along the trajectory and, therefore, they are inherently dependent on the nuclear motion.

Gathering the time-dependence of the electronic wavefunction into the expansion coef-
ficients allows us to replace the evolution equation for the electronic wavefunction with a
set of ordinary differential equations for the evolution of the coefficients $\dot{C}_J(\mathbf{R}(t), t)$. Some
approximations are then introduced, along with the choice of the gauge. The main approxi-
mations are: (i) in the definition of the electron-nuclear coupling operator, Eq. (7), the first
term is neglected within CT-MQC,⁴² as well as its contribution in the expression of the TD-
PES, i.e. $\epsilon_{\text{geo}}$, to maintain the gauge invariance; (ii) the spatial variation of the coefficients
is expressed in the form $\nabla_{\mu} C_J(\mathbf{R}(t), t) = \frac{i}{\hbar} \mathbf{f}_{\mu}^{J}(\mathbf{R}(t), t) C_J(\mathbf{R}(t), t)$, with the force for state J
accumulated along the trajectory

$$
\mathbf{f}_{\mu}^{J}(\mathbf{R}(t), t) = - \int^{t} \nabla_{\mu} \epsilon_{\text{BO}}^{J}(\mathbf{R}(\tau)) d\tau; \qquad (26)
$$

(iii) the TDVP along a trajectory, using the Born-Huang-like expansion of the electronic
wavefunction in its definition (5), is approximated as

$$A_{\mu}(\mathbf{R}(t), t) = \sum_{J} |C_{J}(\mathbf{R}(t), t)|^{2} \mathbf{f}_{\mu}^{J}(\mathbf{R}(t), t), \qquad (27)$$

arguing that the accumulated force, which is integrated over the trajectory and thus accumulates over time, is expected to dominate over the neglected contribution that depends on spatially localized nonadiabatic coupling vectors $d_{\mu,JI}(\mathbf{R}(t))$; (iv) the quantum momentum, i.e. the signature quantity of CT-MQC that naturally emerges within the exact factorization

$$Q_{\mu}(\mathbf{R}(t), t) = \frac{-\hbar \nabla_{\mu} |\chi(\mathbf{R}(t), t)|}{|\chi(\mathbf{R}(t), t)|} \qquad (28)$$

is approximated by expressing the nuclear density as a sum of frozen Gaussians centered at the positions of the trajectories, as discussed in detail elsewhere.³⁶,⁴³ The quantum momentum has a form that is reminiscent of the quantum potential, which, we remind, contains the second spatial derivative of the modulus of the nuclear wavefunction. While in CT-MQC the quantum potential is neglected, the quantum momentum is not, and this does not alter the classical nature of the trajectories. The quantum momentum couples the trajectories and it induces decoherence effects.³⁶

With all these approximations, we arrive at the final evolution equation for the electronic coefficients along a trajectory $\mathbf{R}(t)$

$$\dot{C}_J(\mathbf{R}(t), t) = -\frac{i}{\hbar} \epsilon_{\text{BO}}^J(\mathbf{R}(t), t) C_J(\mathbf{R}(t), t) - \sum_I \sum_{\mu} \dot{\mathbf{R}}(t) d_{\mu,JI}(\mathbf{R}(t), t) C_I(\mathbf{R}(t), t) \\ + \sum_{\mu} \frac{Q_{\mu}(\mathbf{R}(t), t)}{\hbar M_{\mu}} (\mathbf{f}_{\mu}^J(\mathbf{R}(t), t) - A_{\mu}(\mathbf{R}(t), t)) C_J(\mathbf{R}(t), t) \qquad (29)$$

Recalling the characteristic equation (21), we can define the classical forces for the trajectories, $\dot{\mathbf{P}}_{\mu}(t) = \mathbf{F}(\mathbf{R}(t), t) = \frac{d\mathbf{A}_{\mu}(\mathbf{R}(t), t)}{dt}$; taking the total time derivative of the approximated
form of $A_μ(R(t), t)$ allows us to arrive at an equivalent expression for the forces

$$ 
\begin{aligned}
F_{\mu}(\mathbf{R}(t), t) = & -\sum_{J} |C_J \mathbf{R}(t), t|^2 \nabla_{\mu} \epsilon_{\text{BO}}^{J}(\mathbf{R}(t), t) \\
& -\sum_{J,I} \overline{C_J(\mathbf{R}(t), t)} C_I(\mathbf{R}(t), t) (\epsilon_{\text{BO}}^{J}(\mathbf{R}(t), t) - \epsilon_{\text{BO}}^{I}(\mathbf{R}(t), t)) \mathbf{d}_{\mu,JI}(\mathbf{R}(t), t) \\
& +\sum_{J} |C_J(\mathbf{R}(t), t)|^2 \left( \sum_{\mu'} \frac{2\mathbf{Q}_{\mu'}(\mathbf{R}(t), t)}{\hbar M_{\mu'}} \cdot \mathbf{f}_{\mu'}^{J}(\mathbf{R}(t), t) \right) (\mathbf{f}_{\mu}^{J}(\mathbf{R}(t), t) - \mathbf{f}_{\mu'}^{J}(\mathbf{R}(t), t))
\end{aligned}
\quad (30) $$

with $\overline{C_J}$ the complex conjugated of $C_J$. In both expressions for the evolution of the coefficients and for the force, the first two terms are identical to standard Ehrenfest dynamics; the new terms containing the quantum momentum couple the trajectories through its dependence on the nuclear density. Note that the form of the CT-QMC equations has been derived in a particular gauge, thus, it would be difficult to disentangle the effect of the approximation from that of the chosen gauge. Indeed, it is possible to derive an alternative algorithm in a different gauge so as to improve, potentially, the performance of the algorithm. However, we are not aware that such route has been investigated.

## 3 Computational details

The low-energy dynamics studied in this work are initiated in the ground-state of a $E \otimes e$ Jahn-Teller model which was recently explored by Valahu et al.²⁵ The electronic Hamiltonian in the diabatic basis $\{|\phi^1\rangle, |\phi^2\rangle\}$ reads

$$ H_{\text{el}}(X, Y) = \begin{pmatrix} \frac{\omega}{2}(X^2 + Y^2) + \kappa X & \kappa Y \\ \kappa Y & \frac{\omega}{2}(X^2 + Y^2) - \kappa X \end{pmatrix} \quad (31) $$

where $\kappa = 2\pi\text{kHz}$, $\omega = 4\pi/3\text{kHz}$. $R = (X, Y)$ are the nuclear coordinates. The nuclear kinetic energy operator is $\hat{T}_n = \omega(\hat{P}_X^2 + \hat{P}_Y^2)/2$. At time $t = 0$ the function $\chi(R, 0)$ is a real
Gaussian centered at $\mathbf{R}_0 = (-2.5, 0)$ with variance 1 to ensure $|\chi(\mathbf{R}, 0)|^2 \approx |\chi_{S_0}(\mathbf{R}, 0)|^2$. Exact propagation is performed in the diabatic basis using the split-operator technique$^{44}$ for 1.5 ms with a time step of 10 ns. A grid consisting of 501 grid points in X and Y is used in the range $X, Y \in [-8, 8]$. The TDVP and TDPES are reconstructed from the diabatic quantities as detailed in Ref. [37]. We note that the original formulation of this model Hamiltonian was used in a quantum simulator, and therefore the units and scales do not reflect the model of a molecule.

Based on the TDPES and TDVP, 500 quantum and classical trajectories are propagated with initial positions sampled from the initial probability density of the nuclear wavefunction. For the classical trajectories, positions and momenta are treated as independent variables and the momenta are sampled from the momentum probability distribution from the Wigner transform of the initial nuclear wavefunction. The momenta of the quantum trajectories at a certain time $t$ and position $\mathbf{R}$ are given by the TDVP $A(\mathbf{R}, t)$. Consequently, the initial momenta corresponding to initial positions $\mathbf{R}_0$ are obtained as $A(\mathbf{R}_0, 0)$. A time step of 1 $\mu$s is used and the stored TDPES and TDVP are used at each time step.

For CT-MQC, the 500 Wigner sampled trajectories (positions and momenta) are propagated using the G-CTMQC software.$^{45}$ A nuclear and electronic integration time step of 10 ns is used.

It is worth mentioning here that, since the model Hamiltonian of Eq. (31) is given in the diabatic basis, one can work with the exact factorization bypassing the use of a conditional electronic wavefunction that depends on the electronic positions $\mathbf{r}$. Specifically, a conditional electronic state vector can be introduced in $|\Psi(\mathbf{R}, t)\rangle = \chi(\mathbf{R}, t)|\Phi_{\mathbf{R}}(t)\rangle$ where $\chi(\mathbf{R}, t)$ is the same as in Eq. (2), while $\langle\mathbf{r}|\Psi(\mathbf{R}, t)\rangle = \Psi(\mathbf{r}, \mathbf{R}, t)$ and $\langle\mathbf{r}|\Phi_{\mathbf{R}}(t)\rangle = \Phi_{\mathbf{R}}(\mathbf{r}, t)$. Furthermore, in the calculations, the conditional electronic state vector is fully represented using its diabatic
$C_J^{(d)}(\mathbf{R}, t)$ or adiabatic $C_J(\mathbf{R}, t)$ components, namely

$$|\Phi_{\mathbf{R}}(t)\rangle = \begin{pmatrix} C_1^{(d)}(\mathbf{R}, t) \\ C_2^{(d)}(\mathbf{R}, t) \end{pmatrix} = \begin{pmatrix} C_1(\mathbf{R}, t) \\ C_2(\mathbf{R}, t) \end{pmatrix} \qquad (32)$$

where the former are used in the quantum-dynamical calculations and the latter in CT-MQC. In Ref. [37] the expressions of the TDPES and TDVP are given using the diabatic representation and the coefficients $C_J^{(d)}(\mathbf{R}, t)$.

# 4 Results

We discuss in this section the results of the quantum wavepacket dynamics from the model Hamiltonian of Eq. (31) adopting the adiabatic representation. In Section 4.1, instead, we will discuss the dynamics in the eye of the exact factorization, but we will relate our analysis to the dynamics described here.

<figure><img src="image_3.png" /><figcaption>Figure 1: Dynamics in the adiabatic representation. Left panel: populations of the adiabatic states $S_0$ (brown) and $S_1$ (green) as functions of time. Right panels, top: adiabatic excited-state potential energy surface represented as colormaps with the black contour lines representing the excited-state contribution to the total density at the four representative snapshots indicated in the panels. Right panels, bottom: adiabatic ground-state potential energy surface represented as colormaps with the black contour lines representing the ground-state contribution to the total density at the same four representative snapshots.</figcaption></figure>

As discussed in Section 3, the simulated low-energy dynamics starts with the nuclear
density centered around $\mathbf{R}_0 = (-2.5, 0)$, thus $X < 0, Y = 0$, and evolves towards the conical intersection at $\mathbf{R}_{CI} = (0, 0)$, reaching the region $X > 0, Y = 0$. Such low-energy dynamics remain mostly in the adiabatic ground state $S_0$, as confirmed by the populations of the two electronic states plotted as functions of time in the left panel of Fig. 1. However, during the whole simulation, we observe that up to 5% of population is transferred to the excited state $S_1$, confirming that nonadiabatic effects are weak, but not negligible. This is due to the fact that the conical intersection coupling between the electronic states is singular, and manifests itself as a “reverted funnelling” process from $S_0$ to $S_1$.

The adiabatic potential energy surfaces are shown in Fig. 1 in the right panels, with $S_1$ in reported at the top and $S_0$ at the bottom as colormaps. The ground state presents a Mexican-hat-like shape with a cusp at the center at the position of the conical intersection, while the excited state is basically a well that terminates at the bottom with a cone shape whose vertex is degenerate with the cusp of the ground state: there, the typical double-cone shape at the conical intersection appears clearly. The four panels on the right of Fig. 1 show the same (static) potentials but at the different time snapshots we observe the ground-state and the excited-state contributions to the nuclear density that evolve from $X < 0$ towards $X > 0$ remaining symmetric around $Y = 0$.

The adiabatic nuclear densities, after having been initialized around $\mathbf{R}_0$ with the majority of the population in the ground state, move towards the conical intersection. At $t = 0.5$ ms, the ground-state density transfers some population to the excited state, and density appears in the direct vicinity of the conical intersection. Note that ground-state density is completely absent at $\mathbf{R}_{CI}$ because it has been funnelled to the excited state. In moving towards $X > 0, Y = 0$, at $t = 1.0$ ms, the ground-state density branches in the positive and negative directions along $Y$, remaining (numerically) zero for $X > 0, Y = 0$. This behavior is the manifestation of destructive interference during the evolution of the $S_0$ wavepacket encircling the conical intersection, but, as observed in our previous study,²⁶ it does not affect the full nuclear density, only its ground-state component. At the conical intersection, non-zero
density appears only in the excited state. Finally, at $t = 1.45$ ms, while the dynamics progresses towards $X > 0$, we note that the region $X > 0, Y = 0$ is only populated by nuclear density in the excited state.

## 4.1 Nuclear dynamics through quantum and classical trajectories

<figure><img src="image_4.png" /><figcaption>Figure 2: The nuclear density at four snapshots along the dynamics (colour map and contour lines). The green dots indicate the positions of the quantum trajectories at the same time step. $R_{CI}$ is indicated by the red cross.</figcaption></figure>

In order to investigate the low-energy dynamics around the region where the conical intersection is located in the adiabatic representation, we first look at the evolution of the full nuclear density. In Fig. 2, we show the total nuclear density (in contrast to the adiabatic ground state density discussed before) at four snapshots throughout the dynamics. At $t = 0$ ms, the dynamics is initialized centered around $R_0 = (-2.5, 0)$. At $t = 0.5$ ms, the nuclear density reaches the point of the conical intersection at $R_{CI} = (0, 0)$, indicated by a red cross in all plots. At $t = 1.0$ ms, the wavepacket evolves towards positive $X$ values, splitting in $Y$ direction symmetrically around $Y = 0$. Finally, at $t = 1.45$ ms, two parts of the wavepacket have formed that evolve at positive $X$ values, while a diffuse part of the wavepacket remains spread over $X < 0$.

Subsequently, we compare the wavepacket evolution with the dynamics of the quantum trajectories. Quantum trajectories are by construction expected to reproduce closely the quantum dynamics, and in the limit of an infinite amount of trajectories, they would converge to the same dynamics. The quantum trajectories at all snapshots in Fig. 2 follow closely
<figure><img src="image_5.png" /><figcaption>Figure 3: Snapshots of the positions of the quantum trajectories indicated by green dots propagated along the full nonadiabatic dynamics. The overall time traces of the quantum trajectories are shown in each plot with the time indicated by the colour bar. $R_{CI}$ is indicated by the red cross.</figcaption></figure>

the evolution of the nuclear density. The most interesting observation is that along and in close vicinity of $Y = 0$, there are indeed several quantum trajectories evolving past $X = 0$. This clearly indicates the presence of nuclear density in this region of configuration space, which supports our findings of Ref. [26] against the formation of a zero-density line, which would appear as a consequence of the effect of topological-phase interferences in relation to the conical intersection.

The behaviour of the quantum trajectories can be clearly understood by plotting the traces of the trajectory and specific positions at the same four snapshots discussed above (Fig. 3). The colour of the trace indicates the time. The trajectories, randomly sampled from the initial density, are initially spread within [-4 : 0] in the X direction and [-2 : 2] in the Y direction. Once they evolve towards $R_{CI}$ ($t = 0.5$ ms), there is a clear trend of the majority of trajectories avoiding $X > 0, Y = 0$ visible through the deviation of the traces towards positive and negative Y values. However, a small number evolves unfaced through the region of $R_{CI}$, very close to and at $Y = 0$, without any deviation in the Y direction; these trajectories keep evolving solely in the X-direction throughout the dynamics, and the time traces validates this observation. Thus, in general, there is a tendency of avoiding the configuration space around $Y = 0$ for $X > 0$ resulting in a low density of trajectories in the vicinity of this region. However, this region is not completely unaccessible, as shown
by the non-negligible number of the trajectories evolving there. These quantum trajectories strongly underline the fact that in the region of configuration space around $Y = 0$ for $X > 0$ the nuclear density is non-zero. This is of particular interest, since, under the effect of a quantized topological phase, one would expect destructive interference in the quantum wavepacket, leading to a zero-density line in the nuclear density and, subsequently, in the trajectory distribution. However, this is clearly not the case.

While the quantum trajectories reproduce all the significant features of the nuclear-density dynamics, we can evaluate the effect that the classical approximation, i.e. neglecting the quantum potential, has on these complex dynamics. This topic has been previously investigated for adiabatic dynamics as well as excited state dynamics.³⁷,⁴⁶ However, to the best of our knowledge, few studies have looked into this topic for low-energy dynamics where nonadiabatic effects still come into play in the vicinity of the region where the conical intersection appears in the adiabatic representation. Note that a major part of simulations of molecular adiabatic and nonadiabatic dynamics rely on a classical trajectory framework. Therefore, it is highly relevant to understand the direct limitations of this approximation and estimate if it is appropriate to describe low-energy dynamics at conical intersections.

Figure 4 depicts the time traces and position of classical trajectories propagated based on the exact TDVP and TDPES but without the quantum potential. The positions of the trajectories give an overall similar picture as the quantum trajectories and thus, reasonably well reproduce the quantum dynamics. Nonetheless, when looking at the time traces, qualitative differences arise: The classical trajectories completely avoid the region of $Y \in [-0.1, 0.1]$ for $X > 0$. The majority of trajectories approaches this region at different times later during the dynamics, but these trajectories are seemingly repelled. Especially at $t = 1.45$ ms, a large number of trajectories have approached $Y = 0$ for $X > 0$, but their momentum is being reverted in $Y$ direction. We note that the traces of the classical trajectories show a rather “chaotic” behavior. This can be attributed to numerical effects arising from the finite spatial and temporal resolution of the TDVP and TDPES that are employed in their propagation.
<figure><img src="image_6.png" /><figcaption>Figure 4: Snapshots of the positions of classical trajectories indicated by green dots propa-
gated along the full nonadiabatic dynamics. The overall time traces of the classical trajec-
tories are shown in each plot with the time indicated by the colour bar. **R<sub>CI</sub>** is indicated by
the red cross.</figcaption></figure>

In order to explain this behaviour and its origin, we turn to the time-dependent potentials
of the exact factorization, the TDPES and TDVP. It is important to recall the difference
between the quantum and classical trajectory frameworks. As explained in Section 2.2, the
TDVP is the nuclear momentum field and is used in our quantum trajectory calculations
to drive the evolution of the trajectories. For the classical trajectories, the TDPES and
its gradients are used to reconstruct the nuclear momenta and, thus, drive the trajectories.
Figure 5 shows the TDPES and TDVP at a representative time step *t* = 1.0 ms. The TDPES
(left plot in Fig. 5) forms a strong barrier along *X* > 0, *Y* = 0 that has some finite width
in the *Y* direction for all *X* > 0. This barrier persists in time and it is responsible for
the classical trajectories not being able to access the region *X* > 0, *Y* = 0. In comparison,
the TDVP does not show a particular feature along *Y* = 0. Due to the symmetry of the
potential and the initial wavefunction, the *Y* component of the TDVP is zero along *Y* = 0
which explains the motion of the quantum trajectories in the vicinity of *Y* = 0 that move
strictly along the *X* direction all along the dynamics.

To understand in more detail the barrier formed in the TDPES, we plot in Fig. 6 the
three contributions to the TDPES at the representative time step *t* = 1.0 ms as given in
Eqs. (11), (12), and (13): The Born-Oppenheimer term, *ϵ*<sub>BO</sub>, in the left panel, the geometric
term, *ϵ*<sub>geo</sub>, in the central panel, and the gauge-dependent term *ϵ*<sub>GD</sub>, in the right panel. We
<figure><img src="image_7.png" /><figcaption>Figure 5: Left panel: The time-dependent potential energy surface (TDPES) at $t = 1.0$ ms is indicated by the colour map. Right panel: The absolute value of the time-dependent vector potential (TDVP) at $t = 1.0$ ms is indicated by the colour map and its directionality is given by a unit vector. In both panels, the nuclear density at the same time step is indicated by the contour lines and $R_{CI}$ is indicated by the red cross.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: The three terms contributing to the time-dependent potential energy surface (TD-PES) at time $t = 1.0$ ms. In all panels, the nuclear density at the same time step is indicated by the contour lines and $R_{CI}$ is indicated by the red cross. Left panel: Born-Oppenheimer term $\epsilon_{BO}$, central panel: geometric term $\epsilon_{geo}$, right panel: gauge-dependent term $\epsilon_{GD}$.</figcaption></figure>
observe a different manifestation of this barrier in all three terms. The term $\epsilon_{BO}$ encodes, as expected from its definition, the overall shape of the TDPES that closely resembles the potential energy surface of the adiabatic ground state (see Fig. 1 for comparison), where the main part of the dynamics takes place. If one were to decompose into contributions from the adiabatic ground state and from the adiabatic excited state the shape of $\epsilon_{BO}$ at $t = 1.0$ ms, one would associate only the region around $R_{CI}$ and along $X > 0, Y = 0$ to the excited state, because in this region the ground-state wavepacket disappears, as it is transferred to the excited state in a “reverted funnelling” process. This yields a small barrier that is closely localized around $Y = 0$ for $X > 0$. The geometric contribution, $\epsilon_{geo}$, has a very peculiar behavior, as it is almost zero in the majority of the configuration space. However, in the region close to $R_{CI}$, where one expects mixing of the two adiabatic states involved, $\epsilon_{geo}$ contributes significantly to the TDPES with a barrier. Most importantly, this is the main contribution to the observed wide barrier along $Y = 0$ for $X > 0$ in the full TDPES of Fig. 5 (left panel). Thus, the depletion in nuclear density at long times at $X > 0, Y = 0$ can be traced back to be largely caused by this contribution. Finally, the gauge-dependent contribution, $\epsilon_{GD}$, oscillates over the configuration space and seems to form a small barrier which confines the density in X and Y direction and prevents it from delocalizing further in configuration space. Additionally, it also forms a barrier around $R_{CI}$.

Following from these observations of the TDPES, the question arises as to how these contributions affect the classical trajectories over the full course of the dynamics. In the approximate trajectory schemes derived from the exact factorization,³⁶,⁴³,⁴⁷ some of these or similar contributions to the dynamics are neglected or approximated. Therefore, by providing a deeper understanding of the effect that each term has on a classical trajectory-based dynamics, we can offer insight into which approximations deserve particular care.

All these terms are time dependent, so looking at just one specific snapshot, although informative, is insufficient to estimate entirely its contribution and importance for the evolution of classical trajectories. Therefore, we repeated the same dynamics for classical trajectories
<figure><img src="image_9.png" /><figcaption>Figure 7: Snapshots of the positions of classical trajectories at time $t = 1.0$ ms indicated by green dots using different terms of the TDPES. Left panel: Full TDPES, central panel: partial TDPES constructed with only Born-Oppenheimer and gauge-dependent terms, right panel: only Born-Oppenheimer term. The overall time traces of the classical trajectories are shown in each plot with the time indicated by the colour bar. $R_{CI}$ is indicated by the red cross.</figcaption></figure>

progressively removing terms of the TDPES. In Fig. 7, we show the snapshot of the trajectories at time $t = 1.0$ ms as well as the time traces of each trajectory. The left panel reports again, for comparison, the dynamics on the full TDPES that clearly shows the trajectories avoiding a large region around $X > 0, Y = 0$. In the central panel of Fig. 7, we show the dynamics without the geometric contribution $\epsilon_{geo}$. From the positions of the trajectories, it can already be seen that the dynamics are sped up with respect to the full TDPES. The most significant feature becomes clear when looking at the time traces: the trajectories still mostly avoid the region $X > 0, Y = 0$ but they approach significantly closer the line $Y = 0$ and several individual trajectories even cross through $Y = 0$. This brings the dynamics qualitatively closer to what is observed in the quantum trajectories, although they now overestimate how close the trajectories can approach $Y = 0$. We can therefore argue that the absence of classical trajectories at $X > 0, Y = 0$ is majorly caused by the $\epsilon_{geo}$ contribution. It appears that the quantum potential counteracts the effect of $\epsilon_{geo}$, thus, without the quantum potential, one gets spurious zero-density region at $X > 0, Y = 0$. However, $\epsilon_{geo}$ cannot simply be neglected as it is crucial to describe the nuclear density/trajectories surrounding $X > 0, Y = 0$. When also the gauge-dependent term, $\epsilon_{GD}$, is removed, i.e. using only $\epsilon_{BO}$,
(Fig. 7, right panel) the trajectories approach closely $R_{CI}$ and there is almost no repulsion visible in the region surrounding $R_{CI}$. The trajectories still avoid a very small portion of space around $X > 0, Y = 0$ and the dynamics are even further sped up with the trajectories already moving back towards $X < 0$ at 1.0 ms.

Summarizing, we see substantial differences between the quantum and classical trajectories in particular when focusing on the region of low nuclear density along $X > 0, Y = 0$. Both sets of trajectories are able to capture the small nuclear density in this region, but the classical trajectories overestimate this effect.

The TDPES shows a fairly wide, finite barrier along $Y = 0$ for $X > 0$, which is predominantly caused by the $\epsilon_{geo}$ contribution. As this barrier is finite in height, it can be expected that in the limit of an infinite number of trajectories, this region of configuration space will be populated.

Recalling the above presentation of the wavepacket dynamics decomposed into ground-state and excited-state contributions, we observe that the density transferred through the conical intersection to the first excited state is the density that evolves from the negative to the positive $X$ “undisturbed”, i.e. with only a negligible momentum contribution in $Y$ direction, in close vicinity of $Y = 0$. As the quantum trajectories are propagated in the exact-factorization framework, i.e. only based on the TDPES and the TDVP, such a direct relation to the adiabatic contributions is not possible. However, we can assume that the trajectories accessing the configuration space around $X > 0, Y = 0$ can be “attributed” to the adiabatic excited state: as the nuclear density evolves towards positive $X$ with a contribution of the momentum field only along $X$, as seen from the wavepacket propagation, so the quantum trajectories evolving in this part of configuration space have zero momentum along $Y$. During the wavepacket dynamics, up to 5% of population is accessing the excited state, which is qualitatively well reproduced by the fraction of quantum trajectories evolving in $X > 0, Y = 0$, while largely underestimated by the classical trajectories. Thus, it can be argued that the quantum potential is essential to ensure a proper description of these
dynamics in the trajectory framework.

## 4.2 CT-MQC for a low-energy process around $R_{CI}$

As we have seen the classical trajectories being strongly influenced by the terms of the TD-PES that are included in the dynamics, we explore in this section the fully approximate CT-MQC approach, where the TDPES is constructed on-the-fly for each trajectory from the adiabatic quantities. This approach is mainly targeted for the use on molecular systems. Its performance has been estimated for molecular systems³⁶ and for a large variety of low dimensional model systems recreating the deactivation through a conical intersection.⁴⁸⁻⁵⁰ However, so far, its capabilities of describing low-energy dynamics which exhibit some nonadiabatic behaviour has not been tested. This is not a pure theoretical curiosity but also a relevant problem for real molecular systems.⁸,¹²,¹³

In Fig. 8 we depict, as done previously, the positions of the trajectories at four snapshots of the dynamics overlaid with the time traces of all the trajectories. The trajectories are not showing any repulsion along $X > 0, Y = 0$ and they evolve smoothly surrounding $R_{CI}$ without avoiding the region $X > 0, Y = 0$. This behavior can be attributed to the fact that the simulation is performed in the adiabatic basis (the algorithm, and thus the underlying approximations, is derived in the adiabatic representation). First, we observe from the calculation of the populations of the adiabatic states (not shown) that the trajectories do not feel the nonadiabatic effect in the vicinity of the CI and all the coefficients of the adiabatic excited state remain zero throughout the dynamics. Thus, since in CT-MQC one reconstructs the TDPES and the TDVP through the adiabatic quantities, the equations of motion derived in Section 2.3 show that in the case of zero population transfer, the force acting on the trajectories reverts back to the negative gradient of the adiabatic energy. This implies that the trajectories are classically evolving only according to the adiabatic ground state. Due to the shape of the adiabatic ground state PES, the trajectories evolve surrounding $R_{CI}$ but there is no hint towards an energy barrier for $X > 0, Y = 0$ as in the TDPES. This
analysis explains why CT-MQC classical trajectories behave qualitatively differently from
the classical trajectories propagated using the full TDPES. On the other hand, also those
classical trajectories did not feel the nonadiabatic effect encoded in the TDPES through ϵ<sub>geo</sub>:
none of the 500 trajectories of Fig. 4 felt such a nonadiabatic effect, if we assume that it is
only localized in the region X > 0, Y = 0 where ϵ<sub>geo</sub> is strongly non-zero.

<figure><img src="image_10.png" /><figcaption>Figure 8: Snapshots of the positions of classical CT-MQC trajectories indicated by green dots. The overall time traces of the CT-MQC trajectories are shown in each plot with the time indicated by the colour bar. $R_{CI}$ is indicated by the red cross.</figcaption></figure>

In addition, the trajectories are evolving significantly faster than the quantum and clas-
sical trajectories propagated on the TDPES. Specifically, after 1.45 ms the majority of the
trajectories is back in the area of initialization. In the classical trajectories (Section 4.1), ne-
glecting the gauge-dependent and geometric terms of the TDPES caused a significantly faster
evolution of the trajectories. The same is true for the CT-MQC trajectories, where the these
terms are in part neglected (geometric term) and in part approximated (gauge-dependent).

In order to differentiate the effects of the classical approximation and of the reconstruction
of the TDPES and TDVP using adiabatic quantities, we perform numerically exact quantum
dynamics simulations on the adiabatic ground state potential uncoupled from the excited
state. These dynamics can be subsequently analyzed further using quantum and classical
trajectories, in analogy to Section 4.1.

In Fig. 9, we show this ground-state evolution of the nuclear density. The initial evo-
lution up to 0.5 ms looks very similar to the fully nonadiabatic case, with the wavepacket
approaching $R_{CI}$ and starting to get diverted towards positive and negative Y values. After
1 ms, the wavepacket has split into two equal parts in positive and negative Y, but these parts recombine in the positive X region showing constructive interference along Y = 0. The area around $R_{CI}$ is less densely populated than the $X > 0, Y = 0$ region, however, there is still a significant population there. Finally, after 1.45 ms, the wavepacket starts moving back towards $X < 0$, with the larger amount of density centered along $Y = 0$ for $X > 0$ and with a smaller contribution, that split towards positive and negative Y values, that moves back towards $X < 0, Y = 0$.

<figure><img src="image_11.png" /><figcaption>Figure 9: The nuclear density at four snapshots along the adiabatic ground-state dynamics (colour map and contour lines). The green dots indicate the positions of the quantum trajectories at the same time step. $R_{CI}$ is indicated by the red cross.</figcaption></figure>

The quantum trajectories of Fig. 9 follow closely this ground-state adiabatic evolution of the nuclear density. Most importantly, the trajectories evolve in the vicinity of $R_{CI}$. At $t = 1.0$ ms, for instance, when the nuclear density is reduced in that region of configuration space, there is still a non-negligible number of trajectories present. The classical trajectories, as depicted in Fig. 10, are distributed similarly to the quantum nuclear density. Following the time traces, it can be seen that the area around $R_{CI}$ is not entirely avoided, however, there are not many trajectories passing through, which is in line with the reduced nuclear density observed in the quantum propagation. In general, we note that in the adiabatic propagation, the classical approximation does not significantly affect the dynamics. Both, qualitatively and quantitatively, the propagation of the classical trajectories agrees well with the quantum trajectories as well as the wavepacket propagation.

Comparing the classical, adiabatic trajectories with the CT-MQC trajectories, we observe
<figure><img src="image_12.png" /><figcaption>Figure 10: Snapshots of the positions of classical trajectories indicated by green dots prop-
agated along adiabatic ground-state dynamics. The overall time traces of the classical tra-
jectories are shown in each plot with the time indicated by the colour bar. $R_{CI}$ is indicated
by the red cross.</figcaption></figure>

that the overall dynamics shows several similarities. In particular, the trajectories evolve in
a circular-like shape around $R_{CI}$, while slightly avoiding the region close to $R_{CI}$ for $X > 0$.
There are two main differences that the trajectories exhibit: 1) as demonstrated by the
different time traces reported in Figs. 8 and 10, the classical, adiabatic trajectories do not
form an elliptic smooth shape around $R_{CI}$, as done, instead, by the CT-MQC trajectories;
2) the dynamics are significantly faster with CT-MQC, as CT-MQC trajectories return to
the region of configuration space where they were initialized already after 1.45 ms, while the
classical, adiabatic trajectories are still predominantly located at $X > 0$. The faster time
scale, however, can be explained by the approximations underlying the CT-MQC algorithm,
that neglect the geometric term of the TDPES: the results presented in Fig. 7 show that
removing this contribution from the TDPES leads to a significant speed up in the trajectories
dynamics.

Overall, we conclude that the classical approximation has to be used delicately to describe
this low-energy dynamics, and, due to the missing description of nonadiabatic effects, the
approximated CT-MQC dynamics is reverted to a purely adiabatic process unaffected by the
presence of the conical intersection.
# 5 Conclusion

In this work, we investigated different trajectory frameworks to describe low-energy dynamics in the vicinity of a conical intersection (that appears in the adiabatic representation of the dynamics). The particularly challenging nature of this process lies within the fact that nonadiabatic effects are weak and very localized in configuration space, but significantly alter the observed dynamics. The comparison between quantum wavepacket dynamics and quantum trajectories propagated on the TDPES (and with the quantum potential) confirms previous observations on the absence of topological phase effects, that would manifest themselves with a zero-density line in $X > 0, Y = 0$ in the present study, due to destructive interferences when two branches of the wavepackets/trajectories move around the position where the conical intersection is located in the adiabatic representation. Nonetheless, quantum trajectories are able to accurately capture the depletion of density in the same region where the zero-density line would be expected. At the quantum-dynamical level, the symmetry of the model Hamiltonian and of the initial nuclear density studied in this work might naturally indicate that a zero-density line for $X > 0, Y = 0$ can only arise in the adiabatic ground-state component and not in the full density. While in this perspective our observations might appear trivial, an important point of our analysis is the fact that the only (physical) observable quantity is full nuclear density, whose dynamical behavior is solely determined by the TDPES and the TDVP.

We tested the performance of the classical approximation by employing two frameworks of classical trajectories, and both have difficulties to accurately reproduce nonadiabatic effects. The classical trajectories propagated with the exact TDPES (without the quantum potential) actually completely avoid the region of configuration space $X > 0, Y = 0$, thus manifesting a behavior reminiscent of the destructive interference related to the topological phase. This occurs because the nonadiabatic effects that counteract the formation of the zero-density line are underestimated. Thus, in this dynamics, it is the quantum potential that ensures a proper description of the nonadiabatic nature of the dynamics. Dissecting the contributions
of the different terms of the TDPES has shown that the $\epsilon_{geo}$ term as well contributes to a proper account of nonadiabatic effects. The classical trajectories propagated according to the CT-MQC algorithm and, thus, according to an approximated TDPES, do not exhibit any nonadiabatic effects. Therefore, the dynamics reduces to a fully classical propagation along the adiabatic ground state. These dynamics entirely fail to reproduce the influence of the conical intersection: there is no depletion of density in the region $X > 0, Y = 0$ once the trajectories pass the conical intersection.

Finally, we propagated quantum dynamics adiabatically on the ground state and obtained the numerically exact TDPES. The trajectories evolved according to such a TDPES show, indeed, a similar behavior as the CT-MQC trajectories. However, the lack of the geometric contribution to the TDPES and the approximate nature of the gauge-dependent term in CT-MQC, together with the complete absence of nonadiabaticity, leads to a significantly faster dynamics.

Overall, this analysis through trajectories highlights an interesting point: nonadiabatic effects in the vicinity of a conical intersection induce what we can refer to as “partial destructive interference”, manifesting itself with a depletion of nuclear density/trajectories, and, at the same time, this effect counteracts the formation of a zero-density line that would be the signature of a topological phase.

## Acknowledgement

This work was supported by the ANR Q-DeLight project, Grant No. ANR-20-CE29-0014 of the French Agence Nationale de la Recherche.

## References

(1) Baer, M. *Beyond Born-Oppenheimer: Electronic nonadiabatic coupling terms and conical intersections*; John Wiley & Sons, Inc., 2006.
(2) Yarkony, D. R. Diabolical conical intersections. *Rev. Mod. Phys.* **1996**, *68*, 985.

(3) Domcke, W.; Yarkony, D.; Köppel, H. *Conical intersections: electronic structure, dynamics & spectroscopy*; World Scientific, 2004; Vol. 15.

(4) Boeije, Y.; Olivucci, M. From a one-mode to a multi-mode understanding of conical intersection mediated ultrafast organic photochemical reactions. *Chem. Soc. Rev.* **2023**,

(5) Xie, C.; Malbon, C. L.; Guo, H.; Yarkony, D. R. Up to a Sign. The Insidious Effects of Energetically Inaccessible Conical Intersections on Unimolecular Reactions. *Acc. Chem. Res.* **52**, 501–509.

(6) Domcke, W.; Yarkony, D. Role of Conical Intersections in Molecular Spectroscopy and Photoinduced Chemical Dynamics. *Annu. Rev. Phys. Chem.* **2012**, *63*, 325.

(7) Ryabinkin, I. G.; Izmaylov, A. F. Geometric phase effects in dynamics near conical intersections: Symmetry breaking and spatial localization. *Phys. Rev. Lett.* **2013**, *111*, 220406.

(8) Kendrick, B.; Pack, R. T. Geometric phase effects in H+ O2 scattering. I. Surface function solutions in the presence of a conical intersection. *J. Chem. Phys.* **1996**, *104*, 7475–7501.

(9) Bernardi, F.; De, S.; Olivucci, M.; Robb, M. A. The mechanism of ground-state-forbidden photochemical pericyclic reactions: evidence for real conical intersections. *J. Am. Chem. Soc.* **1990**, *112*, 1737–1744.

(10) Celani, P.; Bernardi, F.; Olivucci, M.; Robb, M. A. Conical intersection mechanism for photochemical ring opening in benzospiropyran compounds. *J. Am. Chem. Soc.* **1997**, *119*, 10815–10820.

(11) Xie, C.; Malbon, C. L.; Yarkony, D. R.; Xie, D.; Guo, H. Signatures of a Conical
Intersection in Adiabatic Dissociation on the Ground Electronic State. J. Am. Chem. Soc. **2018**, *140*, 1986–1989.

(12) Desouter-Lecomte, M.; Sannen, C.; Lorquet, J.-C. A transition state theory of nonadiabatic unimolecular reactions controlled by a conical intersection. Application to the C2H+ 4 ion. J. Chem. Phys. **1983**, *79*, 894–904.

(13) Carpenter, B. K.; Ellison, G. B.; Nimlos, M. R.; Scheer, A. M. A conical intersection influences the ground state rearrangement of fulvene to benzene. J. Phys. Chem. A **2022**, *126*, 1429–1447.

(14) Mead, C. A. The molecular Aharonov-Bohm effect in bound states. Chem. Phys. **1980**, *49*, 23–32.

(15) Hazra, J.; Kendrick, B. K.; Balakrishnan, N. Geometric phase effects in ultracold hydrogen exchange reaction. J. Phys. B: Ato. Mol. Opt. Phys. **2016**, *49*, 194004.

(16) Kendrick, B. K.; Hazra, J.; Balakrishnan, N. Geometric Phase Appears in the Ultracold Hydrogen Exchange Reaction. Phys. Rev. Lett. **2015**, *115*, 153201.

(17) Ribeiro, R. F.; Yuen-Zhou, J. Vibronic Ground-State Degeneracies and the Berry Phase: A Continuous Symmetry Perspective. J. Phys. Chem. Lett. **2018**, *9*, 242–247.

(18) Juanes-Marcos, J. C.; Althorpe, S. C.; Wrede, E. Theoretical Study of Geometric Phase Effects in the Hydrogen-Exchange Reaction. Science **2005**, *309*, 1227–1230.

(19) Yuan, D.; Guan, Y.; Chen, W.; Zhao, H.; Yu, S.; Luo, C.; Tan, Y.; Xie, T.; Wang, X.; Sun, Z.; Zhang, D. H.; Yang, X. Observation of the geometric phase effect in the H + HD → H₂+D reaction. Science **2018**, *362*, 1289–1293.

(20) Yuan, D.; Huang, Y.; Chen, W.; Zhao, H.; Yu, S.; Luo, C.; Tan, Y.; Wang, S.; Wang, X.; Sun, Z.; Yang, X. Observation of the geometric phase effect in the H+HD→H2+D reaction below the conical intersection. Nat. Commun. **2020**, *11*, 3640.
(21) Yuan, D.; Yu, S.; Chen, W.; Sang, J.; Luo, C.; Wang, T.; Xu, X.; Casavecchia, P.; Wang, X.; Sun, Z.; Zhang, D. H.; Yang, X. Direct observation of forward-scattering oscillations in the H+HD→H2+D reaction. *Nature Chem.* **2018**, *10*, 653–658.

(22) Xie, Y.; Zhao, H.; Wang, Y.; Huang, Y.; Wang, T.; Xu, X.; Xiao, C.; Sun, Z.; Zhang, D. H.; Yang, X. Quantum interference in H + HD → H₂+ D between direct abstraction and roaming insertion pathways. *Science* **2020**, *368*, 767–771.

(23) Daoud, H.; Joubert-Doriol, L.; Izmaylov, A. F.; Miller, R. D. Exploring vibrational ladder climbing in vibronic coupling models: Toward experimental observation of a geometric phase signature of a conical intersection. *Chem. Phys.* **2018**, *515*, 28–35.

(24) Xie, C.; Kendrick, B. K.; Yarkony, D. R.; Guo, H. Constructive and Destructive Interference in Nonadiabatic Tunneling via Conical Intersections. *J. Chem. Theory Comput.* **2017**, *13*, 1902–1910.

(25) Valahu, C. H.; Olaya-Agudelo, V. C.; MacDonell, R. J.; Navickas, T.; Rao, A. D.; Millican, M. J.; Pérez-Sánchez, J. B.; Yuen-Zhou, J.; Biercuk, M. J.; Hempel, C.; Tan, T. R.; Kassal, I. Direct observation of geometric phase in dynamics around a conical intersection. *Nature Chem.* **2023**, *15*, 1503–1508.

(26) Ibele, L. M.; Sangiogo Gil, E.; Curchod, B. F.; Agostini, F. On the Nature of Geometric and Topological Phases in the Presence of Conical Intersections. *J. Phys. Chem. Lett.* **2023**, *14*, 11625–11631.

(27) Mead, D. G., C. Alden; Truhlar 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.

(28) Martinazzo, R.; Burghardt, I. Dynamics of the molecular geometric phase. *arXiv preprint arXiv:2312.02823* **2023**,
(29) Requist, R.; Tandetzky, F.; Gross, E. Molecular geometric phase from the exact electron-nuclear factorization. Phys. Rev. A **2016**, *93*, 042108.

(30) Abedi, A.; Maitra, N. T.; Gross, E. K. Exact factorization of the time-dependent electron-nuclear wave function. Phys. Rev. Lett. **2010**, *105*, 123002.

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

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

(33) Agostini, F.; Curchod, B. F. E. Different flavors of nonadiabatic molecular dynamics. Wiley Interdiscip. Rev. Comput. Mol. Sci. **2019**, *0*, e1417.

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

(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) Min, S. K.; Agostini, F.; Tavernelli, I.; Gross, E. K. Ab initio nonadiabatic dynamics with coupled trajectories: A rigorous approach to quantum (de) coherence. J. Phys. Chem. Lett. **2017**, *8*, 3048-3055.

(37) Ibele, L. M.; Curchod, B. F.; Agostini, F. A photochemical reaction in different theoretical representations. J. Phys. Chem. A **2022**, *126*, 1263-1281.

(38) Min, S. K.; Agostini, F.; Gross, E. K. Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes. Phys. Rev. Lett. **2015**, *115*, 073001.
(39) Curchod, B. F.; Agostini, F.; Tavernelli, I. CT-MQC-a coupled-trajectory mixed quantum/classical method including nonadiabatic quantum coherence effects. Eur. Phys. J. B **2018**, 91, 1-12.

(40) Meek, G. A.; Levine, B. G. Wave function continuity and the diagonal Born-Oppenheimer correction at conical intersections. J. Chem. Phys. **2016**, 144, 184109.

(41) Talotta, F.; Morisset, S.; Rougeau, N.; Lauvergnat, D.; Agostini, F. Internal conversion and intersystem crossing with the exact factorization. J. Chem. Theory Comput. **2020**, 16, 4833-4848.

(42) Eich, F. G.; Agostini, F. The adiabatic limit of the exact factorization of the electron-nuclear wave function. J. Chem. Phys. **2016**, 145.

(43) Pieroni, C.; Sangiogo Gil, E.; Ibele, L. M.; Persico, M.; Granucci, G.; Agostini, F. Investigating the photodynamics of trans-azobenzene with coupled trajectories. J. Chem Theory Comput.

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

(45) Agostini, F.; Marsili, E.; Talotta, F.; Pieroni, C. G-CTQMC. https://gitlab.com/agostini.work/g-ctmqc, 2021.

(46) Dupuy, L.; Talotta, F.; Agostini, F.; Lauvergnat, D.; Poirier, B.; Scribano, Y. Adiabatic and Nonadiabatic Dynamics with Interacting Quantum Trajectories. J. Chem. Theory Comput. **2022**, 18, 6447-6462.

(47) Ha, J.-K.; Lee, I. S.; Min, S. K. Surface hopping dynamics beyond nonadiabatic couplings for quantum coherence. J. Phys. Chem. Lett. **2018**, 9, 1097-1104.

(48) Talotta, F.; Lauvergnat, D.; Agostini, F. Describing the photo-isomerization of a retinal
chromophore model with coupled and quantum trajectories. J. Chem. Phys. **2022**, *156*,
184104.

(49) Pieroni, C.; Marsili, E.; Lauvergnat, D.; Agostini, F. Relaxation dynamics through a conical intersection: Quantum and quantum-classical studies. J. Chem. Phys. **2021**, *154*, 034104.

(50) Marsili, E.; Olivucci, M.; Lauvergnat, D.; Agostini, F. Quantum and quantum-classical studies of the photoisomerization of a retinal chromophore model. J. Chem Theory Comput. **2020**, *16*, 6032–6048.