Full multiple spawning (FMS) offers an exciting framework for the development of strategies to simulate the excited-state dynamics of molecular systems. FMS proposes to depict the dynamics of nuclear wavepackets by using a growing set of traveling multidimensional Gaussian functions called *trajectory basis functions* (TBFs). Perhaps the most recognized method emanating from FMS is the so-called *ab initio* multiple spawning (AIMS). In AIMS, the couplings between TBFs—in principle exact in FMS—are approximated to allow for the on-the-fly evaluation of required electronic-structure quantities. In addition, AIMS proposes to neglect the so-called second-order nonadiabatic couplings and the diagonal Born–Oppenheimer corrections. While AIMS has been applied successfully to simulate the nonadiabatic dynamics of numerous complex molecules, the direct influence of these missing or approximated terms on the nonadiabatic dynamics when approaching and crossing a conical intersection remains unknown to date. It is also unclear how AIMS could incorporate geometric-phase effects in the vicinity of a conical intersection. In this work, we assess the performance of AIMS in describing the nonadiabatic dynamics through a conical intersection for three two-dimensional, two-state systems that mimic the excited-state dynamics of bis(methylene)adamantyl, butatriene cation, and pyrazine. The population traces and nuclear density dynamics are compared with numerically exact quantum dynamics and trajectory surface hopping results. We find that AIMS offers a qualitatively correct description of the dynamics through a conical intersection for the three model systems. However, any attempt at improving the AIMS results by accounting for the originally neglected second-order nonadiabatic contributions appears to be stymied by the hermiticity requirement of the AIMS Hamiltonian and the independent first-generation approximation.

## I. INTRODUCTION

Simulating the dynamics of a molecule following photoexcitation in an excited electronic state is a dantesque task due to the breakdown of the Born–Oppenheimer approximation, as it implies that one should explicitly account for nonadiabatic effects resulting from the coupling between electronic motion and nuclear motion.^{1,2} An accurate description of such processes would require an exact solution of the molecular time-dependent Schrödinger equation (TDSE), a task only achievable for the smallest molecular systems, unfortunately. As a result, simulating the excited-state dynamics of molecular systems requires an approximate solution of the TDSE.

One family of approximate methods, coined mixed quantum/classical, proposes to propagate the nuclear degrees of freedom classically while preserving a quantum treatment of the electrons.^{3} One of the most famous members of the mixed-quantum/classical family is trajectory surface hopping (TSH),^{4,5} where the dynamics of nuclear wavepackets are approximated by swarms of independent classical trajectories that can jump between adiabatic potential energy surfaces (PESs) when they encounter regions of strong nonadiabaticity. Mixed-quantum classical methods are appealing from a computational perspective, but they can sometimes suffer from their inherent independent trajectory approximation.

Another family of methods for nonadiabatic dynamics emerged from the idea of expanding the nuclear wavefunctions into a basis of coupled moving multidimensional Gaussian functions, the so-called trajectory basis functions (TBFs).^{6–9} In stark contrast with mixed-quantum/classical methods, these Gaussian-based strategies can be derived from first-principles and have an exact limit. By approximating the couplings between TBFs, Gaussian-based approaches can be made compatible with on-the-fly dynamics and become an alternative to TSH for simulating the photodynamics of molecular systems. Different families of Gaussian-based techniques emerge depending on the approximations made to the couplings between TBFs and the TBF dynamics itself. Variational multiconfigurational Gaussian (vMCG) uses a Dirac–Frenkel variational principle for the time evolution of the TBFs and relies on a diabatic representation for the electronic states.^{10,11} Multiconfigurational Ehrenfest (MCE) employs classically evolving TBFs that follow a time-dependent Ehrenfest (mean-field) potential.^{12,13} Full multiple spawning (FMS) and *ab initio* multiple spawning (AIMS) propagate the TBFs classically along adiabatic PESs and allow the number of TBFs to grow during the dynamics for an accurate description of nonadiabatic events.^{14,15} In the following, we will focus our attention on the FMS and AIMS strategies.

FMS offers an in principle exact strategy to simulate the excited-state dynamics of molecular systems. However, the calculation of coupling terms between TBFs hampers the use of this method for the study of molecular systems (as explained in more detail in Sec. II). Two key approximations can be applied to the couplings between TBFs in FMS, leading to the AIMS method that is compatible with on-the-fly dynamics. The first approximation—the saddle-point approximation (SPA)—proposes to approximate the electronic-structure quantity appearing in the coupling elements between TBFs. The second approximation—the independent first-generation approximation (IFGA)—suggests neglecting the couplings between all the parent TBFs, that is, the set of TBFs used to describe the initial state of the nonadiabatic dynamics. The SPA and IFGA will be discussed in detail in Sec. II, together with the spawning algorithm in AIMS.

Earlier works have focused on understanding the influence of the IFGA and the SPA in the excited-state dynamics of simple model systems.^{16–19} Nevertheless, a gray area still remains around additional contributions to the coupled electron-nuclear dynamics that are typically neglected in any practical implementations of AIMS: the second-order nonadiabatic couplings (NACs) and the diagonal Born–Oppenheimer corrections (DBOCs). As AIMS proposes to perform nonadiabatic dynamics in the adiabatic representation of the electronic states, one may also inquire about the inclusion of geometric-phase (GP) effects in its formalism. Geometric (or Berry) phase effects can emerge when a nuclear wavepacket evolving on adiabatic PESs encircles a conical intersection (see Sec. II).^{20} The inclusion of DBOCs in the dynamics of the TBFs in AIMS has been discussed in an earlier work,^{21} and the question of how the lack of GP might affect the AIMS dynamics has already been raised in the past.^{22–24} No work has, however, tested the quality of the AIMS approximations to describe the nonadiabatic dynamics through conical intersections and the influence of the missing NACs, DBOCs, and GP terms mentioned above.

In the following, we propose to closely investigate the quality of the coupling elements between TBFs in AIMS in the context of a dynamics through a conical intersection. In particular, we focus our attention on the neglected contributions to these coupling terms in the original formulation of AIMS. We propose a thorough discussion of the original coupling elements in AIMS and how one could include the missing NACs, DBOCs, and geometric phase while preserving the hermiticity of the AIMS Hamiltonian (Sec. II). To highlight the quality of the AIMS dynamics, we propose a numerical comparison between AIMS, numerically exact quantum dynamics (QD), and TSH for a series of challenging two-state two-dimensional model systems parameterized to reproduce key features of the excited-state dynamics of bis(methylene)adamantyl (BMA), butatriene cation, and pyrazine.^{20,25} We compare the excited-state population traces and the (reconstructed) ground-state nuclear densities for the different nonadiabatic methods (Secs. IV A and IV B). While preserving the philosophy of the AIMS approximations, we test the inclusion of NACs, DBOCs, and GP effects in the coupling terms (Sec. IV C), before questioning the use of Born–Huang (BH) vs Born–Oppenheimer PESs for the dynamics of the TBFs (Sec. IV D). In short, we show that the original AIMS formalism can qualitatively capture well the excited-state dynamics on these model systems and that including additional terms challenges the existing practical approximations of AIMS.

## II. THEORY

This section introduces the standard working equations for the full and *ab initio* multiple spawning starting from the time-dependent molecular Schrödinger equation. We then present the different approximations introduced by AIMS and the challenge in improving such approximations while preserving the hermiticity of the Hamiltonian matrix. We then show, for a two-state system, how to include second-order nonadiabatic couplings and the DBOCs. Finally, we present a set of transformed equations for the AIMS matrix elements that incorporate GP effects.

### A. Full and *ab initio* multiple spawning

The time evolution of a molecular wavefunction Ψ(**r**, **R**, *t*) is obtained by solving the time-dependent Schrödinger equation

with the electronic Hamiltonian $H\u0302el(r,R)$ and the nuclear kinetic energy operator $T\u0302n(R)=\u2212\u2211\rho 3N12M\rho \u2202R\rho 2$ (with $\u2202R\rho 2\u2261\u22022\u2202R\rho 2$). **R** and **r** are collective variables for the coordinates of the *N* nuclei and *N*_{el} electrons forming the molecule, respectively, and *ρ* is used as an index for the 3*N* nuclear coordinates. Note that we use atomic units throughout this work.

The molecular wavefunction Ψ(**r**, **R**, *t*) can be expressed within the Born–Huang ansatz,^{26} which uses the eigenfunctions of $H\u0302el(r,R)$, the so-called electronic wavefunctions ${\Phi J(r;R)}J=1\u221e$ with *J* being a label for electronic states, as a basis,

The time-dependent expansion coefficients, ${\chi J(R,t)}J=1\u221e$, can be seen as nuclear wavefunctions, each one associated with a given electronic state *J*. In this work, we focus on real electronic wavefunctions.

Inserting the Born–Huang ansatz in the time-dependent Schrödinger equation, left multiplying with $\Phi I*(r;R)$, and integration over the electronic coordinates **r** lead to a set of coupled equations of motion for the nuclear wavefunctions,

The first two terms on the right-hand side are related to the adiabatic propagation of the nuclear wavefunction in electronic state *I*, with $EIel(R)$ being the adiabatic electronic energy for state *I* (an eigenvalue of the electronic Hamiltonian) often called potential energy surface. The two last terms are related to the nonadiabatic coupling terms, responsible for the coupled electron-nuclear dynamics, mediated by the first-order nonadiabatic coupling vectors (NACVs) $dIJ(R)=\u3008\Phi I|\u2202R|\Phi J\u3009r$ and the second-order nonadiabatic couplings (NACs), $DIJ(R)=\u3008\Phi I|\u2202R2|\Phi J\u3009r$. The notation ⟨⋯⟩_{r} is used to symbolize an integration over all electronic coordinates. For real electronic eigenfunctions, the diagonal NACVs, **d**_{II}(**R**), are equal to zero. The diagonal NACs, *D*_{II}(**R**), are non-zero and often called diagonal Born–Oppenheimer corrections (DBOCs).

In the full multiple spawning (FMS) framework, the Born–Huang ansatz is employed and each nuclear wavefunction is expressed as a linear combination of multidimensional Gaussian functions, the so-called trajectory basis functions (TBFs),^{14,15,27}

The TBF *m* evolving in electronic state *J*, $\chi \u0303m(J)R;R\u0304m(J)(t),P\u0304m(J)(t),\gamma \u0304m(J)(t),\alpha $, is characterized by the central positions $R\u0304m(J)(t)$ and momenta $P\u0304m(J)(t)$ of the nuclei, a phase $\gamma \u0304m(J)(t)$, and a width ** α**. $CmJ(t)$ is the complex (nuclear) coefficient associated with the TBF

*m*evolving in electronic state

*J*.

As indicated by its name, a TBF is not static but evolves on a given potential energy surface, offering an adequate support to describe the dynamics of the nuclear wavefunctions—in other words, the TBFs can be seen as a moving grid to represent the nuclear wavefunctions. In FMS, the TBFs evolve according to classical equations of motion, and their widths are kept frozen. More importantly, it can be seen in Eq. (4) that the number of TBFs, $NTBFsJ(t)$, is time-dependent. This time dependence comes from the fact that any TBF in the dynamics has the possibility to create new TBFs in the other electronic states when it enters regions of strong nonadiabaticity. These *spawning* events allow for a smooth transfer of nuclear amplitude between the TBFs as a result of nonadiabatic couplings. More details on the spawning algorithm and the dynamics of the TBFs can be found in Refs. 27–29. To initiate an FMS dynamics, we need to project the initial nuclear wavefunction/wavepacket onto a certain number of TBFs. These TBFs, called *parent* TBFs, are coupled from time *t* = 0 and will evolve and spawn *child* TBFs. One common approximation to this picture is to consider that the parent TBFs are uncoupled, which is justified by the fact that the initial nuclear wavepacket will spread rapidly at the beginning of the dynamics. This so-called independent first-generation approximation (IFGA) is employed in practical applications of FMS and AIMS. A detailed discussion of these approximations can be found in Refs. 19, 27, and 29.

To better understand the couplings between TBFs and how TBFs can exchange amplitude, we can insert the FMS ansatz [Eq. (4)] into the time-dependent Schrödinger equation [Eq. (1)]. We obtain, after left multiplication by $\chi \u0303k(I)R;R\u0304k(I)(t),P\u0304k(I)(t),\gamma \u0304k(I)(t),\alpha \Phi I(r;R)*$ and integration over all electronic and nuclear coordinates, a set of coupled equations of motion for the complex coefficients,

Here, **S**_{II} is an overlap matrix between TBFs with elements $SkmII=\u27e8\chi \u0303k(I)|\chi \u0303m(I)\u27e9R$. **H**_{II} and **H**_{IJ} are intra- and interstate Hamiltonian matrices that will be discussed in detail in the following Sections. The Hamiltonian matrix elements are crucial components of the FMS and AIMS methods, as they couple the TBFs and offer exciting opportunities to simplify the FMS equations of motion. As such, they will be at the center of our attention for the remaining parts of this section.

### B. Original version of the *ab initio* multiple spawning matrix elements

We derived above the equations of motion for the (nuclear) complex coefficients in FMS [Eq. (5)]. In this part, we will focus on the Hamiltonian matrix elements in FMS and how they can be approximated—taking us to the AIMS method.

Let us start by describing the two types of Hamiltonian matrix elements, the intrastate (coupled TBFs are in the same electronic state) and interstate (coupled TBFs are in different electronic states) couplings. Following our notation above, the intrastate couplings are characterized by *I* = *J* ∀*k*, *m*, while the interstate couplings are given by *I* ≠ *J* ∀*k*, *m*.

An example of an *intrastate* coupling would be the Hamiltonian matrix element between two TBFs, *k* and *m*, that are evolving in the same electronic state *J*,

The terms on the right-hand side encode the contributions of the nuclear kinetic energy operator, the electronic energy, and the DBOC introduced earlier. An *interstate* coupling is characterized by a matrix element where the two TBFs, *k* and *m*, are evolving in different electronic states, *I* and *J*,

The right-hand side of Eq. (7) shows that the non-Born–Oppenheimer coupling between two TBFs is mediated by the first-order NACVs and the second-order NACs. It is important to realize at this stage that the Hamiltonian matrix elements given above are formally exact, and the integrations over the nuclear coordinates imply that the required electronic-structure quantities—electronic energies, NACVs, DBOCs, and NACs—should be known over the full nuclear configuration space. This requirement dramatically limits the applicability of FMS to molecular systems, and a key strategy to alleviate this issue is to apply an approximation to the matrix elements containing these electronic-structure quantities. This approximation will lead to the AIMS method, which is an approximate strategy within the FMS framework.

Within the AIMS strategy, the **R**-dependence of the electronic energies and the NACVs is approximated by Taylor-expanding these terms around the centroid position of the two TBFs under consideration. In AIMS, only the contribution to zeroth order is preserved, leading to the so-called saddle-point approximation to order zero (SPA0). In addition, AIMS neglects the NACs and the DBOCs, that is, all *D*_{IJ}(**R**) and *D*_{II}(**R**) terms, due to their small size—we will discuss this issue in a following part of this section. We note that neglecting the DBOCs and NACs is an approximation also employed in the practical implementation of *ab initio* MCE and multiple cloning^{13} within the adiabatic representation.^{30} As a result, the intrastate Hamiltonian matrix elements within the AIMS strategy reduce to

while the interstate ones read

where, $EJel(R\u0304km)$ is the electronic energy for state *J* evaluated at $R\u0304km$, the centroid position of the TBFs *k* and *m*. The NACV contribution, $dIJ\rho (R\u0304km)$, follows the same logic. One can appreciate at this stage that the SPA0 addresses elegantly the complexity linked to the **R**-dependence of the electronic energies and the NACVs: the Hamiltonian matrix elements can now be simply calculated by performing at each time step a single point calculation at the centroid position between each pair of TBFs, allowing us to form the different matrices in Eq. (5) and propagate on-the-fly the complex coefficients on the support of the TBFs. In other words, AIMS is a strategy fully compliant with (on-the-fly) *ab initio* molecular dynamics and has been used to simulate the nonadiabatic quantum dynamics of molecular systems in their full dimensionality (see Ref. 27 for a list of examples).

An important consideration has been swept under the carpet in our presentation of the AIMS approximations above: neglecting the NACs in the interstate matrix elements should break the hermiticity of the Hamiltonian, that is, $HIJ\u2260HJI*$.^{31} This is a result of the NACVs being anti-Hermitian, **d**_{IJ}(**R**) = −**d**_{JI}(**R**), and would have the dramatic consequence that AIMS dynamics should not conserve norm. Interestingly, the use of the SPA0 compensates for the neglect of the NACs.^{27} Looking closely at the definition of the AIMS interstate Hamiltonian matrix elements in the SPA0 [Eq. (9)], one can see that the AIMS Hamiltonian *is* Hermitian as $dIJ(R\u0304km)=\u2212dJI(R\u0304km)$ and $\u3008\chi \u0303k(I)|\u2207R|\chi \u0303m(J)\u3009R=\u2212\u3008\chi \u0303m(J)|\u2207R|\chi \u0303k(I)\u3009R$. This compensation between the SPA0 and the neglect of the NACs explains the Hermitian nature of the AIMS Hamiltonian and its norm-conserving dynamics. However, this observation also indicates that one needs to be extremely careful when trying to modify or improve the interstate Hamiltonian matrix elements in AIMS if one wants to preserve its Hermitian nature, as exemplified in Sec. II C.

### C. Issues when moving to the SPA1 in AIMS

Improving the quality of the AIMS matrix elements beyond the SPA0 appears to be a trivial task: one can simply include higher-order terms in the Taylor expansion discussed above. We will show in the following that care has to be taken with the interstate couplings.

Applying the SPA1 to the intrastate Hamiltonian matrix elements is straightforward as it only requires the addition of the first-order term of the Taylor expansion for the electronic energy—the nuclear gradient—evaluated at the centroid position of the two TBFs considered. The general expression for the intrastate couplings is

where the evaluation of the last term on the right-hand side only requires the calculation of a nuclear gradient of the electronic energy for state *J* at the centroid position.

Deriving the SPA1 for the interstate Hamiltonian matrix elements is substantially more challenging for two reasons. The first issue comes from the fact that adding the first-order term of the Taylor expansion for these matrix elements implies the evaluation of the Jacobian of the NACVs—a quantity not commonly available in electronic-structure packages. The Jacobian of the NACVs, $JdIJ(R)$, is defined as

(The explicit form of the Jacobian for the NACVs for a two-state two-dimensional case is provided in the supplementary material.) Considering that we can obtain $JdIJ(R)$, the interstate couplings within the SPA1 would be given by

with $JdIJ(R\u0304km)\rho \rho \u2032$ being the element *ρρ*′ of the Jacobian matrix [see Eq. (11)]. Without considering the difficulties in obtaining the Jacobian, Eq. (12) reveals the second issue mentioned above: it is non-Hermitian. Upon expanding the integrand of the Gaussian integral in the last term of the right-hand side of Eq. (12) (see the supplementary material for the explicit expressions), it becomes clear that most resulting contributions are Hermitian, with the exception of the term $\u2211\rho 3N12M\rho JdIJ(R\u0304km)\rho \rho \u27e8\chi \u0303k(I)|\chi \u0303m(J)\u27e9R$. This term is anti-Hermitian as $JdIJ(R\u0304km)\rho \rho =\u2212JdJI(R\u0304km)\rho \rho $ and $\u27e8\chi \u0303k(I)|\chi \u0303m(J)\u27e9R=\u27e8\chi \u0303m(J)|\chi \u0303k(I)\u27e9R$. The reader is referred to the supplementary material for an explicit discussion of all the contributing terms for a two-level two-dimensional system.

From the above analysis of the SPA1 interstate couplings, it becomes clear that the original approximations in AIMS place this method in a sweet spot where the SPA0 compensates for the loss of hermiticity caused by neglecting the NACs in the interstate couplings. We may now wonder whether including back the NACs could resolve the hermiticity issue observed within the SPA1.

### D. Including the NACs and DBOCs in AIMS—Two-state systems

Based on the findings above that the SPA1 breaks the hermiticity of the AIMS interstate couplings, we wish to investigate whether including back the (long-neglected) NACs could lead to Hamiltonian matrix elements that are Hermitian. We take this opportunity to also discuss the inclusion of the DBOCs in the interstate matrix elements, as they originate from the second-order couplings. We note that discussions of the role of NACs and DBOCs in AIMS have been proposed in Refs. 21 and 32, for example.

To simplify our analysis, we propose here to restrict ourselves to a two-state model, which will allow us to simplify the expressions for the NACs and DBOCs. We start by expanding the NACs into a symmetrized form,

Inserting the resolution of the identity, ∑_{K}|Φ_{K}⟩⟨Φ_{K}| = **1**, into the last term of the right-hand side of Eq. (13), we obtain

For a two-state system, we can derive from Eq. (14) a simplified expression for the DBOCs,

and the NACs,

Hence, both the DBOCs and the NACs can be readily obtained from the NACVs for a two-state model.

We can now use the definitions of DBOCs and NACs provided by Eqs. (15) and (16) in the AIMS matrix elements within the SPA0. The additional term in any matrix elements incorporating the DBOC reads, within the SPA0, as

while the one containing the NACs is

In light of our discussion of Eq. (12), we know that the interstate coupling term, Eq. (18), is anti-Hermitian, as $JdIJ(R\u0304km)\rho \rho =\u2212JdJI(R\u0304km)\rho \rho $, while $\u27e8\chi \u0303k(I)|\chi \u0303m(J)\u27e9R=\u27e8\chi \u0303m(J)|\chi \u0303k(I)\u27e9R$. However, this anti-Hermitian contribution coming from Eq. (18) is exactly the same as the anti-Hermitian contribution obtained from the expansion of the last term of the right-hand side of Eq. (12) (interstate coupling within the SPA1). Hence, if one combines the SPA1 for the NACVs with the SPA0 for the NACs, the two anti-Hermitian parts cancel out, providing an overall Hermitian Hamiltonian. We note that expanding the interstate matrix elements to different orders of the SPA allows us to include all terms that contain the NACVs up to their first-order derivatives.

In summary, we can propose the AIMS Hamiltonian matrix elements for a two-state system using a combination of the SPA0 and SPA1 for both the intrastate couplings,

and the interstate couplings,

We denote this combination of approximations for the Hamiltonian matrix elements *SPA*1′. The SPA1′ is the simplest approximation to the FMS Hamiltonian matrix elements that allows for the inclusion of the NACs and DBOCs while preserving hermiticity. We can now deal with the last missing contribution emanating from the inherent use of the adiabatic representation for the electronic eigenstates in AIMS—the geometric phase.

### E. On the inclusion of geometric-phase effects in *ab initio* multiple spawning

A GP arises in the vicinity of conical intersections when using the adiabatic representation of the electronic states, that is, the eigenfunctions of $H\u0302el(r,R)$, ${\Phi J(r;R)}J=1\u221e$, which parametrically depend on the nuclear coordinates. (See the Appendix for a more extended discussion on this topic.) When following a closed path in (parametric) nuclear coordinates, an adiabatic electronic wavefunction will acquire a phase. This phase will be exactly *π* if the path encircles a conical intersection, leading to a change in sign of the (real) electronic wavefunction that breaks its single-valuedness. To counteract this effect and ensure that the molecular wavefunction is singled valued, one can introduce a position-dependent phase factor, $ei\Omega J(R)$, for the corresponding nuclear wavefunction that will also change sign when encircling a conical intersection. Based on the seminal work by Mead and Truhlar^{33} on the topic, this phase factor will be the same for both states considered in a two-state system and can be chosen to be the mixing angle between the diabatic states used in the unitary transformation between diabatic and adiabatic states. Hence, we can chose Ω_{1}(**R**) = Ω_{2}(**R**) = Ω(**R**) for a two-state system. To perform quantum dynamics in the adiabatic representation, it can be more convenient to perform a transformation of the adiabatic molecular Hamiltonian, $H\u0302mol(r,R)$, to include the effect of a GP while keeping the standard nuclear wavefunctions.^{25} As a result, the molecular Hamiltonian for a two-state system becomes $H\u0302mol(r,R)\u2192H\u0302molGP(r,R)=ei\Omega (R)H\u0302mol(r,R)e\u2212i\Omega (R)$. As stated above, GP Ω(**R**) can be related to the mixing angle, *θ*(**R**), that mediates the diabatic-to-adiabatic transformation for a two-level system. As a result, knowing that ∇_{R}*θ*(**R**) = **d**_{12}(**R**) allows us to choose the GP such that ∇_{R}Ω(**R**) = **d**_{12}(**R**).^{20,33–37} The interested reader is referred to the Appendix for a detailed discussion on this choice of GP.

Coming back to FMS and AIMS, one can evaluate the Hamiltonian matrix elements for $H\u0302molGP(r,R)$, where the action of the Laplacian in the diagonal elements and ∇_{R} in the off-diagonal elements gives rise to additional terms. These terms contain ∇_{R}Ω(**R**) and thus are connected to **d**_{12}(**R**) (see the Appendix for a complete derivation of the GP effects in FMS and AIMS). Expanding all AIMS matrix elements as done in the previous parts of this work and including all terms up to the first-order derivative of **d**_{IJ}(**R**), we obtain the following expression for the intrastate GP Hamiltonian matrix elements:

where the last three terms of the right-hand side arise from the inclusion of the GP. For the interstate GP Hamiltonian matrix elements, we find

Evaluating the AIMS matrix elements within SPA1′ for the GP-transformed Hamiltonian $H\u0302molGP(r,R)$ results in intrastate and interstate couplings containing both contributions from the DBOCs and the NACs.

Armed with equations for the AIMS Hamiltonian matrix elements accounting for NACs, DBOCs, and GP effects while preserving the hermiticity of the Hamiltonian matrix, we propose to unravel the effects of these different contributions by studying the nonadiabatic dynamics of a nuclear wavepacket evolving through a conical intersection for two-dimensional, two-state model systems.

## III. COMPUTATIONAL DETAILS

We used for this study a series of two-dimensional, two-state linear vibronic coupling (LVC) models. In the diabatic representation, the general form of the Hamiltonian is given by

with the following diabatic electronic energies:

with **R** = (*X*, *Y*). The parameters *a*, *ω*_{1}, *ω*_{2}, Δ, and *c* were obtained from Ref. 25 and chosen to create three models representative of the conical intersection for the molecules bis(methylene)adamantyl (BMA), butatriene cation, and pyrazine. The model potentials describe, in a diabatic picture, two parabolas displaced in the *X* direction and shifted in energy by the parameter Δ. Figure 1 depicts the general shape of the adiabatic potential energy surfaces for the different molecules. The ground electronic state is indicated by green contour lines, and the excited state is indicated with purple ones. We note that in the case of BMA, Δ = 0.0 a.u., which means that the model is symmetric with a CI at **R**_{CI} = (*X*_{CI}, *Y*_{CI}) = (0.0, 0.0) (see Fig. 1). In the following, units for the nuclear coordinates are given in bohr.

We simulated the nonadiabatic passage of a nuclear wavepacket through a conical intersection using three different methods: AIMS, numerically exact QD, and TSH.

In the QD simulations, the initial state was given by a Gaussian nuclear wavepacket with widths $\sigma X=2\omega 1$ and $\sigma Y=2\omega 2$ initialized in the adiabatic excited electronic state with zero initial nuclear momenta. The widths (*σ*_{X}, *σ*_{Y}) were set to (16.07, 17.30) for BMA, (14.47, 24.43) for the butatriene cation, and (23.41, 21.86) for pyrazine. Following earlier work by Izmaylov and co-workers,^{25,37} the nuclear wavepacket was initially positioned on the Franck–Condon point of the corresponding full-dimensional model, which gives (15.53, 0.0) for BMA, (−2.08, 0.0) for the butatriene cation, and (5.10, 0.0) for pyrazine. The full time-dependent Schrödinger equation was solved numerically in the diabatic representation employing a split-operator formalism.^{38,39}

The TSH dynamics were carried out with the LVC interface of the SHARC 2.1 program.^{40,41} For consistency, we compared the dynamics obtained from the LVC interface and the analytical potential routine of SHARC and obtained identical results. A nuclear time step of 0.05 fs, that is, ∼2 atomic time units (atu), was used for the classical propagation of the nuclei. The energy-based decoherence correction was used with the standard parameter *a* = 0.1 hartree. Only the component of the nuclear velocities parallel to the NACVs was rescaled after a hop or reflected after a frustrated hop.

The AIMS dynamics were performed with a modified version of the FMS90 code implemented in MOLPRO.^{42} The additional correction terms presented in Sec. II were included into the overlap module of FMS90. A time step of 1 atu was used for the propagation of the TBFs, which was reduced to 0.25 atu in regions of strong nonadiabaticity. The threshold to enter the spawning mode (absolute value of the nonadiabatic coupling) was set to 0.0001 a.u.^{−1}, with a minimum population of 0.001 required for a TBF to spawn. The widths for each TBF were chosen to be the same as those of the initial nuclear wavepacket in the quantum dynamics, and all AIMS simulations presented employed the IFGA.

The AIMS and TSH dynamics used the same 2000 initial conditions, obtained by randomly sampling the Wigner distribution corresponding to the initial Gaussian nuclear wavepacket used in each QD simulations. We note that the additional terms included in AIMS and discussed in Sec. IV C caused some numerical instabilities for the AIMS dynamics of a few initial conditions. As a result, the number of initial conditions used for the different dynamics discussed in Sec. IV C varies between 1980 and 2000.

## IV. RESULTS AND DISCUSSION

### A. Excited-state population decay

Let us first start by investigating the excited-state population decay for the three different models, comparing the QD results to the decays predicted by TSH and AIMS within the SPA0 [Eqs. (8) and (9)].

The BMA model is set up to recreate the diabatic trapping occurring in the excited-state dynamics of the full molecule. As a result, the numerically exact QD dynamics exhibits an almost complete depletion of the initial excited-state population within the first 250 atu [black line in Fig. 2(a)], rapidly followed by a complete revival of the excited-state population after 625 atu. The excited-state population then decays again back to the ground electronic state after around 1000 atu. TSH reproduces quite closely the QD results for the excited-state population [green line in Fig. 2(a)]. Deviations are observed after the first decay (at around 250 atu), where TSH does not fully decay the excited-state population. Another marked difference can be observed during the revival of the excited-state population (at around 750 atu), where TSH predicts that nearly 20% of the population remains in the ground state. These results are entirely consistent with the analysis of Izmaylov and co-workers, who rationalized the overall good performance of TSH for dynamics around conical intersections by a compensation of errors between neglecting DBOCs and GP effects.^{43,44} The population trace obtained with AIMS-SPA0 for BMA [purple line in Fig. 2(a)] is almost identical to that obtained with QD, with only a small deviation observed from the reference during the revival of the excited-state population. This deviation is smaller than what is observed for TSH, with only 10% of population left on the ground state for AIMS-SPA0.

The second model is based on the excited-state dynamics of the butatriene cation. The excited-state population described by QD suffers an immediate decay, leaving ∼15% of population in the excited state after 350 atu of dynamics [black line in Fig. 2(b)]. A brief increase in this population is observed at around 550 atu, before it carries on its decay. A sudden repopulation of the excited state is, however, observed between 1500 and 2000 atu from around 5% to over 30%. The time-trace of the excited-state population is reproduced qualitatively by both TSH and AIMS-SPA0 [green and purple lines in Fig. 2(b)]. The overall shape of the population trace obtained with AIMS is closer to the QD reference than that of TSH, but AIMS-SPA0 underestimates the overall decay of population. The TSH population trace matches that of QD between 850 and 1600 atu. Perhaps more noticeably, both approximate methods miss the strong repopulation of the excited state after 1500 atu.

The third and last model studied is based on the photodynamics of pyrazine. The QD population trace shows an initial decay starting after ∼400 atu of dynamics until the excited-state population reaches ∼15% of population after 1000 atu [black line in Fig. 2(c)]. While the population decay dramatically slows down after this first transfer, the excited-state population decreases to less than 5% after 2250 atu of dynamics. A significant repopulation of the excited state occurs at this time, leading to a maximum of excited-state population at around 3200 atu, before a decrease is again observed. This rather complex population trace is quite well mimicked by TSH and AIMS-SPA0 [green and purple lines in Fig. 2(c)]. Both approximate methods describe the initial fast decay of excited-state population, with AIMS transferring slightly less population to the ground state than TSH. The plateau is then reached with, for both TSH and AIMS, a slightly higher excited-state population than what QD predicts (20% with TSH and 23% with AIMS). In contrast with the dynamics of the butatriene cation described above, both AIMS-SPA0 and TSH manage to describe the repopulation of the excited state after 2250 atu. The overall shape of the repopulation (and then decay) observed during the QD is better reproduced by AIMS, while the variation of population depicted by TSH is smaller (overestimation before the repopulation and underestimation after it).

Overall, the TSH and AIMS-SPA0 agree surprisingly well with the QD reference for these three non-trivial dynamics through a conical intersection. The excellent agreement between the AIMS and QD results obtained for BMA comes somehow as a surprise as it is important to remember that the AIMS-SPA0 dynamics presented here uses the SPA0 and the IFGA. The latter was originally justified for the study of molecules in their full dimensionality, but in resonance with previous work,^{19,45} the IFGA appears to perform well for models with a small number of nuclear degrees of freedom. While these approximations can explain, in part, the deviations in population between AIMS and QD observed for the two other models, it is crucial to realize that the AIMS dynamics presented here do not account for GP effects and neglect the DBOCs and NACs. As observed by Izmaylov for mixed-quantum classical methods,^{43,44} the somewhat encouraging results obtained for the three models with AIMS-SPA0 seem to imply that AIMS benefits from a similar cancelation of errors as TSH—the lack of GP effects in AIMS is compensated by the neglect of DBOCs. We note that the BMA model is characterized by a rather asymmetric DBOC in the branching space, while the models for the butatriene cation and pyrazine show a more symmetric DBOC.^{25,43} These observations and the results above suggest that the approximations of AIMS-SPA0 may perform best, in general, for molecules with less symmetric DBOCs.

The results reported here are somehow reassuring for the use of AIMS-SPA0 to describe the dynamics through conical intersections for molecular systems. Given that AIMS-SPA0 encodes more nuclear quantum effects than a mixed-quantum/classical method such as TSH, it was not given that AIMS-SPA0 would perform well in nonadiabatic dynamics around conical intersections when NACs and geometric effects are significant. Differences are, however, non-negligible, and we propose in the following to investigate more closely the influence of the different missing terms in the AIMS Hamiltonian matrix elements. We focus our attention on the butatriene cation model [Fig. 2(b)] as it appears to be the most challenging model for AIMS-SPA0.

### B. Ground-state nuclear dynamics for the butatriene cation model

To complement our analysis of the excited-state population dynamics, we propose here to visualize the ground-state nuclear density obtained with the different methods. This quantity is straightforward to access in the QD simulation (we note here that we plot the adiabatic ground-state nuclear density). In AIMS-SPA0, we can easily reconstruct the ground-state nuclear density from the TBFs and their complex coefficients owing to the FMS ansatz discussed above. We note that the use of the IFGA implies that the nuclear density is reconstructed for each AIMS-SPA0 run (each parent TBF with their respective child) individually. The overall AIMS-SPA0 nuclear density is obtained by performing an incoherent average over the nuclear density of each AIMS-SPA0 run. TSH does not offer an ansatz for the nuclear wavefunctions due to its independent trajectory approximation. Hence, we approximate the nuclear density in TSH by broadening each TSH trajectory with a Gaussian function, whose width is chosen to be identical to that of the TBFs in AIMS-SPA0.

We compare the ground-state nuclear density at four different times for the butatriene dynamics—200, 552, 1000, and 2000 atu, indicated by colored stars in Fig. 2(b). These times were selected as they correspond to a specific behavior in the excited-state population trace: the initial decay of the population (200 atu), the maximum of the weak repopulation (552 atu), the beginning of the plateau (1000 atu), and the end of the dynamics after the repopulation (2000 atu). Figure 3 shows the different snapshots of the ground-state nuclear density for the QD, AIMS-SPA0, and TSH dynamics (from top to bottom), reconstructed at these four times. The pink circle in the plots indicates the location of the conical intersection.

At 200 atu, nuclear density starts to appear in the ground electronic state. While the location of the nuclear density given by AIMS-SPA0 and TSH matches that of the QD result, the QD nuclear density shows a nodal region near the position of the conical intersection, which is not reproduced by either AIMS-SPA0 or TSH. The approximations of AIMS-SPA0 can be held responsible for this deviation from the QD result, considering that the AIMS-SPA0 simulation presented invokes the SPA0 and the IFGA (and no geometric effects are included). At later times, the ground-state nuclear density spreads on the ground-state potential energy surface (see Fig. 1), ending up being mostly delocalized by *t* = 2000 atu. While AIMS-SPA0 and TSH reproduce well the overall shape and localization of the nuclear density at most times, it is surprising to see how qualitatively well AIMS-SPA0 captures the main features of the structure of the QD nuclear density, in particular at *t* = 552 atu and *t* = 1000 atu, considering the crude approximation of its Hamiltonian matrix elements and the use of the IFGA. This good match between QD and AIMS-SPA0 is further observed in the actual movies of the evolution of the ground-state nuclear density (see the supplementary material). These movies further highlight that the good agreement between the QD and AIMS-SPA0 nuclear densities is not coincidental at the selected times. Comparing the movies for the AIMS-SPA0 and TSH dynamics offers a clear illustration of the fundamental difference between the two methods. The TSH movie makes it obvious that the motion of the nuclear density is simply correlated to the dynamics of the independent classical trajectories on the ground-state potential energy surface. Conversely, the AIMS-SPA0 movie shows that the complex amplitude is transferred between the TBFs, which serve as moving grid points on which the nuclear density is distributed.

This observation is further supported by comparing the positions of the TSH trajectories with those of the TBF centers in AIMS-SPA0 at the various snapshot times (Fig. 4). In AIMS-SPA0, the TBFs evolve from being mainly around *y* = 0 at 200 atu to mostly everywhere in space at later times—an essential element for the success of the method if one considers that the actual positions of the TBFs are not critical as long as they offer a proper support to describe the nuclear wavepacket dynamics, described in AIMS-SPA0 via the complex coefficients. This coverage of the configuration space by the TBFs is made possible by the spawning algorithm, which increases the number of TBFs continuously during the dynamics. This growing distribution and spreading of the TBFs are in stark contrast with the distribution of the TSH trajectories.

### C. Adding new contributions to the Hamiltonian matrix elements in AIMS

Now that we characterized the nonadiabatic dynamics obtained with AIMS-SPA0 for the different models, we wish to investigate the effect of including additional contributions to the AIMS Hamiltonian matrix elements, focusing, in particular, on some terms that are usually neglected in AIMS—GP effect, DBOCs, and NACs—and discussed in Sec. II.

Let us first start by considering the influence of improving the quality of the intrastate matrix elements in AIMS-SPA0 by moving to the SPA1 [*intrastate*-*SPA*1 in Fig. 5(a), dark purple dashed line], that is, by using Eq. (10) for the intrastate couplings and maintaining the SPA0 [Eq. (9)] for the interstate couplings. No major variations of the population trace can be observed by improving the intrastate couplings. This observation could be rationalized by the fact that, within a given AIMS run, TBFs might separate in phase space sufficiently such that they are not strongly affected by their respective intrastate couplings. An interesting test to validate the importance of the coupling between TBFs on the same state consists in removing all intrastate couplings [pseudo-independent trajectory approximation in AIMS, *pITA*-*AIMS*, in Fig. 5(a), gray dashed line]. Once more, the population trace obtained within this approximation does not significantly deviate from the original AIMS one, validating the weak influence of the intrastate coupling within the current approximations. It is important to note that removing the IFGA may significantly alter this result, as all parent TBFs and their child TBFs would become coupled. Before moving to the next step and including the NACs and GP effects, we note that we tested the bra-ket averaged Taylor (BAT) expansion to order zero proposed in Ref. 13 for the NACVs (instead of the SPA0) and obtained a population trace in very close agreement with AIMS-SPA0.

We are now in a good position to include terms in the Hamiltonian matrix elements that are usually neglected within the conventional AIMS framework based on our analysis in Sec. II. The first step consists in including the NACs in the interstate couplings. As discussed in our analysis in Sec. II, preserving the hermiticity of the Hamiltonian matrix requires that we use the SPA1 for the interstate couplings. The resulting strategy [*AIMS* + *NACs*, gray line in Fig. 5(b)] uses the SPA1 for the intrastate couplings [Eq. (10)] and the interstate couplings [Eq. (20)], including the NACs. AIMS + NACs shows a population transfer initially accelerated in comparison to AIMS-SPA0, but this acceleration rapidly slows down, and a significantly smaller population transfer toward the ground state is observed between 400 and 750 atu. The excited-state population then slowly converges toward that of AIMS-SPA0. Based on our earlier consideration with intrastate-SPA1 and pITA-AIMS, it does not come as a surprise that adding the DBOCs to this dynamics [*AIMS*-*SPA*1′, dark blue dashed line in Fig. 5(b), using Eq. (19) for intrastate couplings and Eq. (20) for interstate couplings] does not alter the population transfer in comparison to AIMS + NACs. Overall, the inclusion of NACs and DBOCs appears to slow down the population transfer for the butatriene cation model—an observation for AIMS in resonance with earlier findings on TSH, where the inclusion of DBOCs was found to slow down significantly the population transfer.^{43} (See Sec. IV D for additional information on the inclusion of DBOCs in AIMS.) The final contribution that we can include is the GP correction within SPA1′ (*GP*-*SPA*1′, palatinate line in Fig. 5) using Eq. (21) for intrastate couplings and Eq. (22) for interstate couplings. The population trace obtained with AIMS under GP-SPA1′ deviates even more from the AIMS-SPA0 (and QD result), exhibiting a dramatic slowdown of the excited-state population transfer. This result is surprising as one would expect, based on earlier works,^{25,43} that including a GP correction in the interstate coupling terms should increase the nonadiabatic transfer of population.

To shed further light on this observation, we propose to focus our attention on a representative parent–child TBF pair and monitor their interstate coupling when approaching the conical intersection.^{46} Following the strength of the interstate coupling for each TBF using AIMS-SPA0 [Fig. 6(b)], one can see how the TBFs meet in configuration space when their coupling strength is maximal—a behavior encoded in the spawning algorithm. When now monitoring the strength of the interstate couplings for the same pair of TBFs but including GP effects [GP + SPA1′, Fig. 6(c)], one can immediately observe that the nonadiabatic strength is almost ten times stronger than without the GP effects. Hence, it appears that, as expected, the nonadiabatic strength is dramatically enhanced when GP effects are included in AIMS. What happens, then, to the population transfer between this pair of TBFs?

The norm of the Hamiltonian matrix element for the interstate coupling for AIMS-SPA0, $Hm,kmIJAIMSSPA0$, or for GP-SPA1′, $Hm,kmGP,IJAIMSGP-SPA1\u2032$, are taken as proxies for the nonadiabatic coupling strength in a situation without or with GP effects. We start by monitoring the excited-state population trace for this pair of TBFs using AIMS-SPA0 and GP + SPA1′ for the matrix elements [Fig. 6(a)]. The population trace using AIMS-SPA0 shows a smooth and monotonic decay. However, upon inclusion of the GP effects, one sees that, while the population starts to decay earlier than in AIMS-SPA0, it then exhibits some rapid oscillations between 1.0 and 0.5, hindering the population transfer. While the TBFs feel a much stronger mutual interstate coupling, the nuclear amplitude oscillates in a Rabi-way between the two TBFs without a stable transfer to the ground-state one. One reason for that could be that there are, at this time, only two TBFs in the dynamics (one on each state) and, therefore, that once transferred to the TBF in the ground state, the amplitude cannot spread toward other ground-state TBFs and remains trapped on the child TBF, still strongly coupled with the parent TBF—leading to an oscillation of the amplitude between the two TBFs, a sort of overcoherent effect due to the locality of the parent/child TBF pair. Including the GP effects only in the interstate couplings [palatinate dashed line in Fig. 6(a)] further enhances these oscillations. Hence, we appear to be facing, in this particular case, one of the limitations of AIMS: the use of the IFGA, which limits the number of coupled TBFs at the early times. To circumvent this issue, it would be necessary to start the dynamics with a large number of coupled TBFs. While this solution would potentially be tractable for a model system, its cost would simply be prohibitive for molecular systems.

We come to the interesting conclusion that including GP effects in AIMS would likely require getting rid of some of the practical approximations (SPA0 and IFGA) that made AIMS a successful strategy to investigate the nonadiabatic dynamics of molecular systems. These tests also reveal that these nonadiabatic processes can be easily described within these strong approximations by relying on cancelation of errors—a reason why TSH performs so well for molecular systems. However, trying to account properly for additional quantum effects such as GP in AIMS would require a dramatic improvement of its underlying approximations.

### D. Born–Oppenheimer or Born–Huang potential energy surfaces in AIMS

Up to this point, we have focused our attention only on the Hamiltonian matrix elements connecting the TBFs. However, one should keep in mind that, while the dynamics of the complex coefficients is obtained via the time-dependent Schrödinger equation [Eq. (5)], the TBFs are propagated classically and meant to follow the dynamics of the nuclear wavepackets. In the standard version of AIMS (AIMS-SPA0), the TBFs evolve on the Born–Oppenheimer potential energy surfaces (BO PESs) given by $EIel(R)$. Considering the inclusion of DBOCs in the Hamiltonian matrix elements of AIMS-SPA0 raises the question of also including these terms in the dynamics of the TBFs. The resulting Born–Huang (BH) PESs are then defined for the electronic state *I* as $EIel(R)+DII(R)$. The BH PESs have been investigated extensively for electron transfer processes in the context of non-Born–Oppenheimer effects.^{47–49} The BO PESs are very similar to the BH PESs for a given molecule as long as electronic states are far apart in energy. However, one can show for a two-level system that the DBOCs are proportional to the squared norm of the NACVs [as discussed earlier; see Eq. (15)]. Hence, the shape of the BH PESs will dramatically differ from that of the BO PESs in the vicinity of a conical intersection, as the NACVs are inversely proportional to the electronic energy gap between the two coupled states. Looking at our model for the butatriene cation, one can see that the typical conical shape of the BO PES [Fig. 7(a)] is completely altered by the inclusion of the DBOCs, leading to BH PESs both displaying a singularity at the exact location of the conical intersection [Fig. 7(b)]. The addition of this singular barrier on both BH PESs is made more visible in Fig. 7(c) by performing a cut on the BH and BO PESs along the *Y* = 0 axis.

This comparison between the shape of BO and BH PESs makes it clear that nonadiabatic dynamics methods based on classical trajectories will be strongly affected by this additional barrier. Earlier work on TSH showed that including the DBOCs in the propagation of the classical trajectories dramatically hampers the nonadiabatic population transfer.^{43,44} Hence, focusing now on AIMS-SPA0, it does not come as a surprise that the TBFs are avoiding the regions of the conical intersection when evolving on the BH PESs [Fig. 8(b)]. Interestingly, though, the width of the TBFs makes that the reconstructed ground-state nuclear density [Fig. 8(a)], combined with the approximate nature of the interstate and intrastate couplings, still displays some nuclear density in the direct vicinity of the conical intersection. As expected from the distribution of the TBFs, the overall excited-state population decay predicted by AIMS-SPA0 when employing BH PESs to propagate the TBFs is dramatically decelerated (light blue line in Fig. 9) in comparison to the AIMS-SPA0 on BO PES (purple line in Fig. 9). This observation is in line with earlier findings on propagating TBFs on BH PESs within the AIMS-SPA0 formalism^{21} and earlier work on the inclusion of DBOCs only in QD and TSH.^{20,25,43,44}

Accounting for GP effects into the AIMS dynamics on BH PESs (at the GP-SPA1′ level; see Sec. IV C) results in a slight speed-up of the excited-state population transfer. This behavior is expected due to the enhanced interstate couplings resulting from the inclusion of the GP (see Fig. 6) and shows that TBFs can still interact next to the conical intersection. However, the classical motion of the TBFs on the BH PESs, combined with the limitations in the description of the interstate and intrastate couplings and the IFGA (described above), does not allow for an overall improved description of the excited-state population transfer. This finding reinforces earlier work on the incorporation of DBOCs into AIMS.^{21,32}

## V. CONCLUSION

In this work, we analyzed the influence of the missing contributions to the Hamiltonian matrix elements in the AIMS method for nonadiabatic dynamics around conical intersections. Based on three two-dimensional model systems, we could show that AIMS—in its original formulation, that is, using the SPA0 and IFGA and neglecting the NACs, DBOCs, and GP corrections—offers a least qualitatively correct description of the dynamics through a conical intersection. The results obtained with AIMS align with what is observed for the mixed-quantum/classical method TSH. As AIMS is not *per se* a mixed-quantum/classical method but is derived from the in-principle exact FMS framework, possible improvements of the method can be envisaged. We derived different sets of equations aiming at improving the Hamiltonian matrix elements coupling the TBFs. We discussed the potential loss of hermiticity for the AIMS Hamiltonian when including either higher-order terms of the SPA or the NACs and GP corrections. Interestingly, adding terms in the interstate or intrastate couplings outside of the realm of the original AIMS approximations does not appear to improve the description of population transfer at a conical intersection. The limited accuracy in the Hamiltonian matrix elements combined with the IFGA seems to be responsible for such shortcomings. Considering that the complexity related to dynamics through a conical intersection is particularly enhanced in the three low-dimensional models presented in this work, the adequate behavior of the original AIMS in such conditions offers an empirical validation for the simulation of nonadiabatic processes in higher dimensions using this method. The shortcomings of AIMS when improving its Hamiltonian matrix elements in the adiabatic representation are also highly stimulating for developing new TBF-based strategies for nonadiabatic dynamics employing dedicated (time-dependent) diabatic electronic quantities.^{21,50,51} The exact factorization of the molecular wavefunction,^{52} which produces smooth time-dependent electronic quantities in the vicinity of a conical intersection,^{53} also offers an exciting venue for the development of coupled-trajectory methods for nonadiabatic dynamics.^{54,55}

## SUPPLEMENTARY MATERIAL

See the supplementary material for a derivation of the explicit form of the SPA1 for the first-order nonadiabatic couplings, an overview of the terms included in the matrix elements for various AIMS flavors, and a movie of the ground-state nuclear density for the butatriene cation model as depicted by AIMS, QD, and TSH.

## ACKNOWLEDGMENTS

This project received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 803718, project SINDAM). L.M.I. acknowledges the EPSRC for an EPSRC Doctoral Studentship (Grant No. EP/R513039/1). The authors thank Federica Agostini and Loïc Joubert-Doriol for fruitful discussions.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: ADIABATIZATION, GEOMETRIC PHASE, AND CONNECTION WITH HAMILTONIAN MATRIX ELEMENTS IN AIMS

In this appendix, we provide additional information on the adiabatization of the model potentials used in the main text and how the geometric phase appears in this context. We also derive the working equations for the AIMS Hamiltonian matrix elements that include geometric phase effects.

The diabatic Hamiltonian of our two-state model system is defined by

Diagonalizing the potential energy matrix to obtain adiabatic electronic states can be done through the unitary transformation matrix^{20,34}

where *θ*(**R**) is the mixing angle, also known as the adiabatic-to-diabatic transformation angle, between the diabatic electronic states *V*_{11}(**R**) and *V*_{22}(**R**). *θ*(**R**) is defined as

The diabatic-to-adiabatic transformation produces adiabatic electronic wavefunctions, |Φ_{J}(**R**)⟩, from the diabatic electronic wavefunctions, $|\Phi Jdia(R)\u3009$, according to

The ground and excited adiabatic electronic energies, $E\u2212el(R)$ and $E+el(R)$, respectively, are obtained as

The nonadiabatic coupling terms can be expressed in terms of derivatives of the mixing angle *θ*(**R**),^{34,36,37}

Closer investigation of *θ*(**R**) shows that upon encircling a conical intersection (a point of degeneracy of the adiabatic electronic eigenstates), the phase of the adiabatic electronic wavefunction will change by a factor of *π*.^{25,33,37,56–58} This factor of *π* means that the (real) adiabatic wavefunctions, given by Eq. (A4), will change their sign by circling around a conical intersection. This phase shift and subsequent flip of sign are a manifestation of the geometric phase.^{56,57} However, the geometric phase leads to a double-valuedness of the electronic wavefunctions. For a proper description of the molecular wavefunction (which is now doubled-valued), one would need to use doubled-valued boundary conditions for the nuclear wavefunctions as well. However, as suggested by Mead and Truhlar,^{33} this can be avoided by, instead, transforming the Hamiltonian as

The choice of *e*^{iθ(R)} for the geometric phase is valid in the two-state models used in this work, but it is not universal. For a detailed discussion of this choice, its validity, and further implications, we refer the interested reader to Refs. 59 and 60.

We will make use of this approach for the inclusion of the geometric phase in AIMS. The molecular Hamiltonian used in FMS, denoted as $H\u0302mol(r,R)$ in the main text, will be transformed into $H\u0302mol(r,R)\u2192H\u0302molGP(r,R)=ei\theta (R)H\u0302mol(r,R)e\u2212i\theta (R)$. Using this transformed Hamiltonian in FMS leads to a new set of equations for the Hamiltonian matrix elements, where for the intrastate couplings, the action of the Laplacian (coming from the nuclear kinetic energy operator) on the phase gives rise to new terms,

The last three terms on the right-hand side are the additional terms arising from the geometric phase.

For the interstate couplings, the presence of the ∇_{R} operator (with the NACV term) gives rise to an additional term (last term on the right-hand side),

As seen above in Eq. (A6), it is possible to connect the derivative of the mixing angle ∇_{R}*θ*(**R**) and the NACVs **d**_{12}(**R**). Thus, we can rewrite Eqs. (A8) and (A9) with the additional contributions now given in terms of the NACVs $d12\rho (R)$,

From this set of equations, we can now derive a practical expression for AIMS by applying a series of approximations. One can expand the electronic energy term and all terms including the NACVs within the SPA to first order (SPA1), while all the DBOC and NAC terms [including $|d12\rho (R)|2$ and $\u2202R\rho dIJ\rho (R)$] are expanded to zeroth order (SPA0). Within this framework, we arrive at the following equation for the intrastate couplings:

and we arrive at the following one for the interstate couplings:

These two equations are those employed in the main text under the name GP-SPA1′. This framework preserves the hermiticity of the AIMS Hamiltonian.

## REFERENCES

We note that this parent–child TBF pair has been randomly chosen from the set of initial conditions and that the vast majority of parent–child TBF pairs exhibits similar behavior.