Simulations of photoinduced processes with the exact factorization: state of the art and perspectives

Abstract

This perspective offers an overview of the applications of the exact factorization of the electron–nuclear wavefunction to the domain of theoretical photochemistry, where the aim is to gain insights into the ultrafast dynamics of molecular systems via simulations of their excited-state dynamics beyond the Born–Oppenheimer approximation. The exact factorization offers an alternative viewpoint to the Born–Huang representation for the interpretation of dynamical processes involving the electronic ground and excited states as well as their coupling through the nuclear motion. Therefore, the formalism has been used to derive algorithms for quantum molecular-dynamics simulations where the nuclear motion is treated using trajectories and the electrons are treated quantum mechanically. These algorithms have the characteristic features of being based on coupled and on auxiliary trajectories, and have shown excellent performance in describing a variety of excited-state processes, as this perspective illustrates. We conclude with a discussion on the authors' point of view on the future of the exact factorization.

---

1 Introduction

Photochemical and photophysical phenomena in molecular systems occur around us at every moment. Light absorption initiates the chemical reactions at the basis of photosynthesis in plants,¹ activates the isomerization of the retinal chromophore in rhodopsin culminating in animal vision,² and converts provitamin D into vitamin D in our body.³ For technological and medical applications, the goal of understanding these and similar phenomena is the design of systems that function as well as nature itself, to produce photovoltaic current upon absorption of light,⁴ to label biological cells carrying diseases to the human body,⁵ or to produce light in a clean and efficient way via electroluminescence in organic light emitting diodes.⁶ Understanding and controlling the behavior of matter in all these situations require a deep knowledge of the fundamental quantum-mechanical laws governing the complex interplay of the molecules with light and with their environment at the microscopic scale over time scales ranging from femtoseconds to nanoseconds.

In this rich landscape of molecular processes and systems, the field of quantum molecular dynamics has sparkled in recent years.⁷ Great progress has been achieved at the more fundamental, theoretical level as well as from the point of view of the algorithms/software developments and applications. Many novel ideas have been proposed, motivated by the challenges encountered when simulating complex organic materials^4,6,8–10 and biological systems,^3,9,11,12 or when investigating phenomena over long time scales to look at slow processes such as fluorescence^13,14 or vibrational relaxation, or when creating hybrid light–matter states in the strong coupling regime in optical and plasmonic microcavities.^15–19

Simulating photochemical and photophysical processes in molecular systems requires to describe the interplay of electronic and nuclear dynamics, possibly with the explicit inclusion of light, and, in particular, accounting for excited-state effects, i.e., the so-called nonadiabatic or beyond-Born–Oppenheimer effects.²⁰ For excited-state processes, on-the-fly molecular-dynamics simulations are perhaps the most widely used approaches, allowing one to access complex systems over reasonably long time scales up to few tens of picoseconds. Aiming to solve the time-dependent Schrödinger equation and to approximate the fundamental quantity of interest, i.e., the time-dependent molecular wavefunction,²¹ these methods rely on the support of various kinds of trajectories to mimic the nuclear dynamics under the effect of the electronic ground state and the (necessary) excited states. Ranging from quantum to classical, the nuclear evolution can be treated using:

• trajectory-basis functions, generally used to evolve Gaussian wavepackets by integrating ordinary differential equations that in some cases closely resemble Hamilton's equations,^22–33

• quantum,^34–42 coupled,^43–47 auxiliary^48,49 trajectories, retaining fully or partially the quantum character of the overall dynamics,

• semiclassical methods combined with the path-integral formulation of the quantum-mechanical propagator,^11,50–55

• independent classical trajectories,^9,47,56–84 allowing one to access efficiently large systems of “experimental complexity”, even if, sometimes, at the cost a losing accuracy.

In on-the-fly simulations, the other aspect of the problem, i.e., the electrons, enters the nuclear evolution in a somehow static way, as the necessary properties, such as energies, gradients, nonadiabatic couplings, spin–orbit couplings, or transition dipole moments, are determined (on-the-fly) along the trajectories at the visited nuclear geometries. In addition, using an ab initio or a semi-empirical or an analytical description of the electronic Hamiltonian, the most suitable electronic representation for the problem at hand and for the methodology chosen for the simulations needs to be chosen, e.g., the adiabatic vs. the diabatic representation.

It is worth mentioning here that other classes of methodologies exist to simulate nonadiabatic dynamics, such as approaches based on the density-matrix formalism^10,85–93 or on quantum wavepackets propagation,^94–99 but will not be covered in this perspective.

As it is clear from the above discussion, the field of quantum molecular dynamics is extremely vast. Consequently, this perspective will only discuss the impact that theoretical developments based on the exact factorization of the electron–nuclear wavefunction^100,101 had in this domain so far, and will present our viewpoint on some promising avenues for the future of the exact factorization.

In molecular systems, the standard theoretical construction to go beyond the Born–Oppenheimer approximation¹⁰² and to account for non-adiabatic effects arising from the coupled electron–nuclear nature of these systems, is the Born–Huang representation¹⁰³ of the molecular wavefunction.²¹ This is an exact representation of the electron–nuclear wavefunction as a linear combination of the so-called adiabatic states, i.e., the eigenstates of the electronic Hamiltonian at each nuclear geometry. An alternative to the Born–Huang representation is the exact factorization, which expresses the electron–nuclear wavefunction as the product of a marginal nuclear amplitude and a conditional electronic amplitude, parametrically dependent on the nuclear configuration.¹⁰⁰ The exact factorization can be extended to any multicomponent many-body wavefunction, since no assumption is made on the physical properties of the subsystems to identify the marginal and the conditional amplitudes, e.g., the small electron–nuclear mass ratio usually invoked in the Born–Oppenheimer approximation.

The exact factorization depicts a photochemical reaction in terms of the dynamics of nuclei and electrons, offering a reformulation of the quantum dynamics of a molecule that is free from concepts such as static potential energy surfaces, conical intersections,^104–107 and electronic transitions.⁵⁶ Specifically, it was formulated in 2010 by Gross and co-workers to tackle the problem of describing the nonadiabatic dynamics of electrons and nuclei in a molecule (H₂⁺ in that case) under the effect of an external laser pulse.^100,108 This work proposed to analyze and to simulate the coupled dynamics of electrons and nuclei in nonadiabatic conditions employing the time-dependent potential energy surface and the time-dependent vector potential, which are concepts arising purely from the dynamics of the electrons even in the absence of external time-dependent fields. These time-dependent potentials provided an original tool to rethink our way of looking at nonadiabatic processes: in the Born–Huang framework, nuclear wavepackets (blue and pink shaded areas in Fig. 1 represent the modulus squared of the wavepackets) evolve under the effect of static electronic potential energy surfaces (blue and pink curves in Fig. 1 represent the energies of the ground state and the first excited state as functions of a one-dimensional nuclear coordinate) and exchange amplitude in the regions where the energies are close, thus accounting for electronic nonadiabatic transitions; in the eye of the exact factorization, the evolution of a single nuclear wavefunction (purple-colored areas in Fig. 1 represent the modulus squared of the wavefunction) is driven by a single time-dependent (scalar and vector) potential incorporating the dynamical effects of the electrons (the purple curve in Fig. 1, i.e., the time-dependent potential energy surface, changes in time but its shape is reminiscent of the shapes of the blue and pink curves of the upper panels).

Fig. 1 Schematic comparison of the “standard” way of looking at the steps of a photochemical “experiment” and the exact-factorization way. Top panels: Pink and blue solid lines indicate the potential energy curves of the ground and excited state as a function of a one-dimensional nuclear coordinate. The Gaussian-like shapes represent the nuclear densities associated to the ground state (pink) and excited state (blue) evolving on the support of the respective potential energy curve. Bottom panels: The solid purple line shows the time-dependent potential energy curve as a function of a one-dimensional nuclear coordinate. The purple coloured area indicates the nuclear density that evolves on the support of the time-dependent potential energy curve. From left to right: Initially the system is in its ground state (Initialisation), absorption of a photon excites part of the system to the excited state, creating a wavepacket (Excitation), this wavepacket evolves first adiabatically on the excited state (Evolution) until it reaches a conical intersection where the population is transfered back to the ground state (Deactivation) where it can finally form photoproducts by equilibrating in another region of configuration space (Photoproducts).

The introduction of the time-dependent potentials in the framework of the exact factorization, along with a natural analogy with classical electromagnetism, lends itself to the introduction of the concept of a classical force,^109,110 that is uniquely defined in the framework of the exact factorization. The concept of a classical force goes hands in hands with that of trajectories:¹¹¹ in this way, the formalism has been employed successfully and largely in the domain of quantum molecular dynamics and of on-the-fly simulations.^48,112–114 Perhaps, this is the field where the exact factorization made the most progress and proposed the most interesting developments, some of which will be discussed in this perspective.

The time-dependent version of the exact factorization followed previous work based on the stationary Schrödinger equation.^115–119 We refer the reader to ref. 119–128 for an overview of the topics tackled by various authors employing the time-independent formulation of the exact factorization. In addition, alternative viewpoints have been proposed to study electron–nuclear,^{126,129–134} electron–electron^135–139 and photon–electron–nuclear^140–145 systems.

In this perspective, we will focus on the theoretical and on the algorithmic aspects of the exact factorization of the electron–nuclear wavefunction. The formulation of the theory will be briefly recalled in Section 2 and we will describe the coupled-trajectory and auxiliary-trajectory algorithms derived from it to simulate nonadiabatic processes in Section 3. Some applications of the exact factorization and of these algorithms are presented in Section 4. Section 5 presents a brief overview of recent work on the dynamics of photon–electron–nuclear systems employing the exact-factorization framework. Our conclusions are discussed in Section 6.

2 The exact factorization framework

A non-relativistic quantum-mechanical system comprised of two sets of interacting particles, with coordinates r and R, can be described by the Hamiltonian Ĥ(r, R), which can be written as

Ĥ(r, R) = T̂(R) + Ĥ_BO(r, R) (1)

We indicate with T̂(R) the kinetic energy of the particles identified by the coordinates R, in our case the nuclei, and with Ĥ_BO(r, R) the sum of the remaining terms, i.e., the kinetic energy of the particles identified by the coordinates r, in our case the electrons, and all the interactions. In the electron–nuclear problem, Ĥ_BO(r, R) is usually referred to as Born–Oppenheimer (BO) Hamiltonian, which explains our choice of labelling this term with the subscript “BO”. When describing the interaction of a molecule with (classical) light, as in the case of a laser pulse^146–149 or a continuous-wave laser,^150–156 an additional time-dependent term may appear in the Hamiltonian, which can be expressed as V̂(r, R, t) = −[small mu, Greek, circumflex](r, R)·E(t) in the dipole approximation. In this case, the interaction potential V̂(r, R, t) depends on [small mu, Greek, circumflex](r, R), the molecule electric dipole moment, and on E(t), the external time-dependent electric field. In some situations,^157–167 the relativistic correction due to the geometry-dependent spin–orbit coupling Ĥ_SOC(x, R) can be included in the system Hamiltonian, where the variable x = [r,σ] accounts for the electronic spin σ as well as for the electronic positions r.

The time-dependent Schrödinger equation (TDSE) with Ĥ(r, R) dictates the evolution of the wavefunction Ψ(r, R, t) as

The exact factorization of the time-dependent wavefunction has been proposed by Gross and coworkers^100,108 in the form

Ψ(r, R, t) = χ(R, t)Φ(r, t; R) (3)

by expressing the full wavefunction as the product of a marginal amplitude, i.e., χ(R, t), and a conditional amplitude, i.e., Φ(r, t; R), that parametrically depends on the coordinates R. The interpretation of |Ψ(r, R, t)|² as a joint probability density allows one to easily identify |χ(R, t)|² as the marginal probability density of R, with |Φ(r, t; R)|² the conditional probability density, whose normalization over r has to be imposed ∀R, t. The TDSE combined with such a partial normalization condition yieldswhich are the coupled evolution equations for χ(R, t) and Φ(r, t; R), respectively. The index ν is used to label the nuclei.

Eqn (4) is itself a TDSE where the coupling to the dynamics of Φ(r, t; R) is expressed in terms of a time-dependent vector potential (TDVP)

A_ν(R, t) = 〈Φ(t; R)| − iħ∇_νΦ(t; R)〉_r (6)

and a time-dependent potential energy surface (TDPES)In the definitions of the TDVP and TDPES, we indicate with the symbol 〈·〉_r an integration over r, and we use the convention that such variable does not appear explicitly in the term in the bracket as it is integrated over.

The TDVP encodes information about the momentum field of the nuclei, since

where χ(R, t) = exp[(i/ħ)S(R, t)]|χ(R, t)|.

Note that the product form of the wavefunction Ψ(r, R, t) is invariant under the phase transformations [capital Phi, Greek, tilde](r, t; R) = exp[(i/ħ)θ(R, t)]Φ(r, t; R) and [small chi, Greek, tilde](R, t) = exp[(−i/ħ)θ(R, t)]χ(R, t). Thus, under these transformations, the TDVP and TDPES transform as well, as standard gauge potentials, namely Ã_ν(R, t) = A_ν(R, t) + ∇_νθ(R, t) and [small epsilon, Greek, tilde](R, t) = ε(R, t) + ∂_tθ(R, t), and eqn (4) and (5) are form-invariant.

Eqn (5) yields the evolution of Φ(r, t; R), where the coupling to the dynamics of eqn (4) is provided by

The exact factorization of the time-dependent wavefunction has been employed extensively in the field of quantum molecular dynamics as a tool to interpret and to simulate the coupled electron–nuclear dynamics in nonadiabatic conditions, i.e., beyond the BO approximation. In the Born–Huang representation of the coupled electron–nuclear problem, the molecular wavefunction is written as a linear combination of adiabatic electronic states φ_l(r; R), that are the eigenstates of the BO Hamiltonian, with coefficients χ_l(R, t) usually referred to as nuclear amplitudes, namelyIf in the exact-factorization form of the molecular wavefunction, the electronic conditional amplitude is expressed as a Born–Huang-like expansion,the relation χ_l(R, t) = χ(R, t)C_l(R, t) follows. This relation is quite interesting as it relates clearly the exact factorization to the Born–Huang representation. First, it shows that information is not lost when expressing the molecular wavefunction as a single product rather than as a sum of products, since the nuclear amplitudes can be reconstructed using the nuclear wavefunction and the electronic coefficients. Second, it introduces the concept of spatially-resolved electronic coefficients, i.e., C_l(R, t), yielding the contributions of the electronic states at each point in configuration space R (note that since the adiabatic states are orthonormal, the normalization condition on Φ(r, t; R) implies that ). Third, it allows to somehow partition the nuclear wavefunction, which does not contain any information about the adiabatic states because it is a purely nuclear quantity, into adiabatic contributions.

While attempts have been made to solve exactly eqn (4) and (5),^168,169 the most interesting developments of the exact factorization of the electron–nuclear wavefunction focused on introducing the quantum-classical perspective within this framework, such that the nuclear dynamics is ultimately approximated using classical-like trajectories while the electronic dynamics is treated quantum mechanically.^{39,41,43,112,170–176} The quantum-classical perspective naturally emerges in the exact factorization, since the classical limit of eqn (4) simply means to interpret the TDVP and TDPES as standard (classical) electromagnetic potentials producing a classical force on the nuclei.¹¹¹ The main efforts in this context have been devoted, then, to the calculation and approximation of the TDVP and of the TDPES from the solution of the quantum-mechanical electronic eqn (5). To this end, extensive work has been devoted to analyze the TDVP and TDPES in model situations, such as in photo-activated processes,^109,110,177 in relaxation dynamics through conical intersections,^{147,178–180} in the presence of quantum interferences,¹⁸¹ under the effect of an external time-dependent field,^{100,108,147,150,182–185} or including spin–orbit coupling.^113,186 Those numerical analyses, together with theoretical developments focusing on understanding the behavior of the equations in some limiting cases, e.g., the adiabatic limit,¹⁸⁷ or on resolution strategies of the nuclear equation, e.g., using semiclassical techniques¹⁸⁸ or the method of characteristics,^111,169 yielded various trajectory-based schemes readily applicable to on-the-fly molecular-dynamics simulations, as we will present now.

3 Trajectory-based methods for nonadiabatic dynamics in molecules

The exact factorization of the electron–nuclear wavefunction Ψ(r, R, t), especially thanks to the TDVP and TDPES, has been remarkably useful for understanding how the quantum-classical perspective is introduced starting from the quantum formulation of the dynamical problem, once eqn (4) and (5) are given. However, accurately and efficiently calculating the electronic TDVP and TDPES in the course of a quantum molecular-dynamics simulation requires to develop algorithms^{43,44,46,48,49,150,171,186,189,190} and software^{112,191–194} exploiting the strengths of current quantum-chemistry techniques for on-the-fly dynamics. With this idea in mind, let us present and compare the various quantum-classical schemes that have been introduced in recent years to simulate nonadiabatic (photochemical) processes.

In general, the quantum-classical perspective on the exact factorization requires to replace the concept of a quantum nuclear wavefunction with an ensemble of trajectories. Formally, this is done by replacing in the equations R, a 3N_n-dimensional vector, with N_n the number of nuclei of the system, with the symbol R^α(t). Here, α = 1, …, N_tr labels the trajectories, that have to be many in order to reproduce the delocalization of the nuclei in configuration space, and for each α and at each time t, R^α(t) is a 3N_n-dimensional vector. The trajectory-based nuclear dynamics can be simply summarized using Hamilton's equations

Ṗ^α_ν(t) = F^α_ν(t) (13)

where we introduced here, for every trajectory α and for every nucleus ν, the concepts of velocity, i.e., Ṙ^α_ν(t), of momentum, i.e., P^α_ν(t), and of force, i.e., F^α_ν(t). In particular, P^α_ν(t) and F^α_ν(t) are expressed in terms of the TDVP and of the TDPES.

The trajectories can be assimilated to a moving grid: while in the quantum treatment the value of a function f(R, t) can be determined at time t at any point R, in the quantum-classical treatment only the values f(R^α(t), t) are accessible. Clearly, for a very large number of trajectories N_tr, information is not lost when going from the quantum to the quantum-classical treatment. It is important to note that, in order to evaluate how functions of the type f(R^α(t), t) evolve, only total time-derivatives can be computed instead of partial time-derivatives.

Aiming to develop on-the-fly procedures to solve the exact-factorization equations by exploiting quantum chemistry, the electronic wavefunction is expanded in the adiabatic basis, as shown in eqn (11). Therefore, the expression is inserted in eqn (5) to derive a set of evolution equations for the expansion coefficients C_l(R, t). Following from the above observations, when the idea of trajectories is introduced, the electronic eqn (5) is affected as well, and the evolution of the coefficients is ultimately expressed as a total time-derivative Ċ_l(R^α(t),t) = Ċ^α_l(t).

The key quantity arising in such a quantum-classical formulation of the electron–nuclear dynamics is the so-called quantum momentum. In eqn (5), the operator Û[Φ, χ] depends on the nuclear wavefunction, and, when it is expressed in polar form in terms of its modulus |χ| (or ) and phase S, one gets

where the term in square brackets is the classical momentum P_ν = M_νṘ_ν, that depends on the TDVP as anticipated above, and it is summed to the quantum momentum , which appears as a purely imaginary correction. The quantum momentum encloses information about the spatial delocalization of the nuclear density, and this information has to be recovered even when trajectories are used to mimic the nuclear dynamics. Therefore, the non-local character of the quantum nuclear density is encoded in the quantum momentum, that is why the algorithms derived from the exact factorization require either a coupled-trajectory scheme^43,44 or an auxiliary-trajectory scheme.^42,48

In general, the exact-factorization-based algorithms derived so far can be viewed as variations of the Ehrenfest scheme or of the surface-hopping method, where the effect of the quantum momentum is included. The advantage of these particular formulations is twofold: first, since Ehrenfest and surface hopping are the most widely-used approaches for simulations of photochemical processes, they are well-understood by the community and are implemented in many codes; second, it is easy to show how the corrections related to the inclusion of the quantum momentum cure the major drawbacks of both approaches, i.e., the mean-field character of Ehrenfest¹⁷⁴ and the overcoherence problem of surface hopping.^59,66,195 To demonstrate this second point, we find instructive to show here, for a simple model case, the performance of Ehrenfest dynamics and of surface hopping. Similarly to the scheme of Fig. 1, we use a one-dimensional two-electronic-state model defined by the green and orange potential energy curves reported in the upper panels of Fig. 2. The dynamics is initiated in the excited state with a Gaussian wavepacket centered at R = 2 bohr, whose density is represented as light-blue curves in both upper panels of the figure. The dynamics proceeds with such a photo-excited wavepacket moving towards the right and crossing the region where the energies of the ground and of the excited states are close. There, a nonadiabatic event takes place and the wavepacket branches, following partially the ground-state potential and partially the excited-state potential. This behavior can be verified at t = 12 fs, when the magenta nuclear density presents two peaks, attesting to the fact that the potential energy curves with different slopes drive the dynamics. Later on, at t = 24 fs, the nuclear density, represented in purple, is completely delocalized, with the portion on the right associated to the ground state while the portion on the left crosses once again the nonadiabatic region after being reflected by the barrier on the far right of the excited state potential curve.

Fig. 2 Left panels: Ehrenfest dynamics (using 1000 trajectories). Right panels: Surface hopping dynamics (using 2000 trajectories). In the upper panels, the ground-state (dark-green line) and the excited-state (orange line) potential energy curves are shown, along with the quantum nuclear density at three times along the simulated dynamics compared to the Ehrenfest or surface hopping trajectories at the same times. In the lower panels, the population of the excited state is shown as gray curves, superimposed to the population estimated during Ehrenfest dynamics (left) or during surface hopping dynamics (right).

The mean-field character of the Ehrenfest approach is not able, by construction, to capture with trajectories the dynamics just described, characterized by a final density that proceeds along diverging paths. This is very clearly represented by the circles in the upper left panel of Fig. 2, that show the positions of the Ehrenfest trajectories at the same time snapshots of the density, along with the potential energy felt by those trajectories, whose gradient is used to compute the force. The absence of the splitting in the distribution of trajectories yields wrong nuclear dynamics, and, consequently, it is not surprising that the population of the excited state (lower left panel of Fig. 2) at the end of the simulated dynamics is not correctly reproduced by Ehrenfest. The surface-hopping approach cures the mean-field problem of Ehrenfest. Specifically, the trajectories follow at all times either one or the other potential energy, that is why the circles on the right panel of Fig. 2 are distributed always either on the ground-state or on the excited-state energy curve. Surface-hopping trajectories correctly reproduce the splitting of the nuclear density since a hopping algorithm allows them to change instantaneously the potential that drives their dynamics (so-called, active state), with high probability in the region of strong nonadiabaticity. Nonetheless, the disconnect between the potential (or the force) that drives the dynamics of the trajectories and the evolution of the electronic coefficients (similarly to the Ċ^α_l(t) above) can be source of disagreement between the two ways of estimating the electronic populations: one way is to count the fraction of trajectories in each state, while the other is via the coefficients. A surface-hopping procedure is internally consistent if these two ways of calculating the electronic populations yield the same or a similar result, as it should. When this is not the case, the algorithm suffers from overcoherence. This is clearly shown at the final times of the dynamics simulated with the fewest-switches surface-hopping algorithm⁵⁶ (without decoherence corrections) and reported in the lower right panel of Fig. 2, where the purple-dashed and blue-continuous curves diverge. The small discrepancy between the exact curve (gray) and the curve calculated from the fraction of trajectories (blue) is probably due to the limited statistics, as only 2000 trajectories are used in this simulation. Note that many authors have devoted their efforts towards understanding and curing the problem of overcoherence in surface hopping^{45,59,66,72,196–202} (as well as in Ehrenfest^75,203), since usually a surface-hopping simulation needs to account for decoherence corrections. Here, however, we show in Fig. 2 results obtained without including decoherence simply to explain the issue and to justify the need to go beyond the original algorithm.

Various algorithms based on the exact factorization are capable to circumvent naturally at the same time the mean-field and overcoherence issues of Ehrenfest and surface hopping thanks to the presence of the quantum momentum. In the past few years, different groups introduced different flavours of quantum-classical algorithms: the coupled-trajectory mixed quantum-classical (CTMQC) scheme was developed by Gross and coworkers in 2015,⁴³ and as the name says, it employs an ensemble of coupled trajectories to mimic the nuclear dynamics; later on, the method surface hopping based on the exact factorization (SHXF) was introduced by Min and coworkers in 2018,⁴⁸ aiming to use a surface-hopping scheme to mitigate the problems related to the large computational cost of CTMQC by evolving independent trajectories and adopting auxiliary trajectories when necessary; with a similar purpose, in 2021, Pieroni and Agostini developed the coupled-trajectory Tully surface hopping (CTTSH),⁴⁴ which uses the idea of surface hopping within a coupled-trajectory procedure; Ha and Min proposed in 2022 an independent-trajectory version of CTMQC, named Ehrenfest dynamics based on the exact factorization (EhXF),⁴⁹ which requires auxiliary trajectories similarly to SHXF; aiming to alleviate the issue of non-conservation of the total energy in CTMQC, Villaseco Arribas and Maitra introduced in 2023 a variation of the algorithm, i.e., CTMQC-E,^190,204 by imposing energy conservation over the ensemble of trajectories, following a similar idea that in EhXF is imposed at the single-trajectory level; in a similar spirit, Blumberger and coworkers developed in 2023 CTMQC-(E)DI,⁴⁶ where a double-intercept (DI) idea is introduced in CTMQC and in CTMQC-E to cure numerical instabilities encountered in CTMQC(-E) when calculating the quantum momentum and, thus, greatly improving energy and norm conservation; in 2024 Maitra and coworkers, inspired by the work of Martens,²⁰¹ developed the quantum-trajectory surface-hopping based on the exact factorization (QTSH-XF),⁴² where a phase-space approach combined with the exact factorization cures the frustrated hops and velocity rescaling issues that SHXF inherits from surface hopping; an assessment of various exact-factorization-based algorithms, i.e., those based on auxiliary trajectories, along with their implementation in the Libra package²⁰⁵ has been recently presented by Han and Akimov.¹⁹⁴ Currently, to the best of our knowledge, the available open-source codes able to perform exact-factorization-based calculations are PyUNIxMD,¹⁹³ developed in the group of Min and interfaced with various quantum-chemistry packages, G-CTMQC,²⁰⁶ developed by our group and interfaced with the QuantumModelLib library of analytical potentials,²⁰⁷ and Libra,¹⁹⁴ which for the moment only allows calculations with model potentials; CTTSH has been recently implemented in MOPAC-PI.^208,209

As an example of the structure of the equations defining the CTMQC algorithm^43,173 as a variation of the Ehrenfest scheme, let us show the following electronic and nuclear (force) equations

Ċ^α_l(t)|_CTMQC = Ċ^α_l(t)|_Ehr + Ċ^α_l(t)|_qm (15)

F^α_ν(t)|_CTMQC = F^α_ν(t)|_Ehr + F^α_ν(t)|_qm (16)

Here, the first terms on the right-hand side are Ehrenfest-like terms while the additional terms depend on the quantum momentum. Similarly, CTTSH⁴⁴ can be summarized as

Ċ^α_l(t)|_CTTSH = Ċ^α_l(t)|_TSH + Ċ^α_l(t)|_qm (17)

F^α_ν(t)|_CTTSH = F^α_ν(t)|_TSH (18)

where Ċ^α_l(t)|_TSH = Ċ^α_l(t)|_Ehr, as shown previously,¹⁷⁴ and F^α_ν(t)|_TSH is simply an adiabatic force calculated as the gradient of the electronic active state that is chosen stochastically at each time according to the fewest-switches algorithm.⁵⁶

To give an idea of the performance of the exact factorization, we show in Fig. 3 how CTMQC works in the situation described in Fig. 2. The delocalization of the distribution of trajectories is correctly reproduced by CTMQC, especially as it shown by the circles at the final simulated time t = 24 fs (green circles). Such a delocalization seems to appear more slowly than in the quantum dynamics, as the distribution of trajectories at t = 12 fs (blue circles) suggests. However, it is interesting to note how a single potential is able to yield different forces in different portions of space so as to induce the splitting of the distribution of trajectories (that is not possible to achieve in Ehrenfest while in surface hopping it is only possible by introducing the hopping idea). In this simple case, the TDPES can be calculated exactly, and it is shown in Fig. 3 as crosses at times t = 0 (light-blue crosses), t = 12 fs (magenta crosses) and t = 24 fs (purple crosses). Only at the intermediate time shown in the figure, the distribution of CTMQC trajectories along the potential that drives their dynamics, i.e., the blue circles, does not follow the TDPES, i.e., the magenta crosses, confirming, as observed previously, that CTMQC trajectories reproduce the splitting of the nuclear density with some delay. Nonetheless, such a delay does not affect significantly the population of the excited state (lower panel in Fig. 3) which perfectly agrees with the reference (gray line) all along the dynamics.

Fig. 3 Same as in Fig. 2 but for the CTMQC algorithm (using 1000 trajectories).

3.1 Exact-factorization-based surface-hopping approaches

Among the various trajectory-based algorithms derived from the exact factorization and listed in the previous section, those that adopt a surface-hopping flavor have shown to be the most promising for molecular applications in photochemistry. Surface hopping is a widely-used, numerically-efficient, (almost) black-box and largely-implemented^{64,193,205,206,208,210–212} algorithm, thus it is easy to use as a starting point for new developments, either theoretical or computational. In addition, some analyses^109,110 performed on the features of the TDPES and of the TDVP as well as their role in the dynamics draw interesting connections with surface hopping, motivating our efforts for combining it with the exact-factorization idea. Therefore, CTTSH and SHXF stand out as alternative choices to (standard) surface hopping, which ultimately has to be always cured with at least one of the many decoherence-correction schemes.^59,66

The main reason for moving away from CTMQC, which can be considered as the algorithm rigorously derived from the exact factorization, and for turning towards surface-hopping-based approaches is the computational advantage of the latter. Specifically,

• CTMQC requires the calculation of the nonadiabatic coupling vectors at all times. They are 3N_n-dimensional vectors for each pair of adiabatic states, and account for the spatial variation of the electronic states as functions of the nuclear coordinates. Their calculation is computationally expensive, but in surface hopping (i) only their projections along the nuclear velocity, which is a scalar quantity, is needed at all times, and (ii) depending on the implementations, they are calculated in full dimensionality only after a trajectory hop occurs.

• CTMQC potentially yields unstable long-time dynamics if the trajectories cross multiple times (the same) regions of configuration space with strong nonadiabaticity^44,213 or in the case of quantum interferences.^181,214 This happens because the algorithm, ultimately, produces a nuclear dynamics driven by an approximate TDPES. Such TDPES becomes noisy in the situations just mentioned, but in surface hopping this problem is naturally circumvented because the nuclear dynamics takes place on the adiabatic surfaces.

• CTMQC does not guarantee energy conservation, either along a single trajectory or over the ensemble of coupled trajectories, which is mainly due to the approximate nature of the quantum momentum used in the equations. This problem is partially alleviated in CTMQC-E by modifying the expression of the nuclear forces by imposing energy conservation a posteriori over the ensemble of trajectories. In surface hopping, instead, energy conservation is natural when the trajectories evolve adiabatically and is enforced on each trajectory when a hop occurs.

Both surface-hopping-based algorithms, i.e., CTTSH and SHXF, have been employed for calculations of molecular properties, as will be discussed in Section 4.2. The basic equations are the same, as both methods follow the same electronic evolution equation (eqn (17)), nuclear propagation scheme (eqn (18)) and fewest-switches hopping probability,^44,56 but differ in the way they reconstruct the nuclear density for computing the quantum momentum (eqn (14)).

In CTTSH the quantum momentum is computed by approximating the nuclear density as a sum of Gaussians centered at the position of the trajectories and by adjusting this analytical expression by finding a modified expression guaranteeing that when averaged over the ensemble of coupled trajectories, the populations of the adiabatic states do not change if the nonadiabatic couplings are zero.^{112,173,177,209,215} The modified expression is used by default, unless it diverges numerically, in which case it is replaced by the analytical expression.

SHXF, instead, is an independent-trajectory approach, thus the quantum momentum has to be constructed with the support of auxiliary trajectories.⁴⁸ In this case, frozen Gaussians are centered at the positions of these auxiliary – non-physical – trajectories such that the quantum momentum can be approximated locally using an analogous expression as the analytical expression used in CTTSH. SHXF is less computationally demanding than CTTSH since each physical trajectory can be propagated independently form the others, and the auxiliary trajectories are associated only to the electronic states that are coupled nonadiabatically to the active state. The auxiliary trajectories are launched when non-active electronic states become populated, with initial positions chosen to be the same as the position of the physical trajectory and with initial velocities determined by imposing that each auxiliary trajectory has the same kinetic or the same total energy of the physical trajectory. After they are created and associated to a certain electronic state, the auxiliary trajectories change positions with a velocity that at each time step is calculated by rescaling in all directions the velocity of the physical trajectory. The rescaling factor is determined by imposing energy conservation using the potential of the corresponding electronic state at the position of the physical trajectory. This procedure avoids additional electronic-structure calculations, namely the energy and the gradient of the electronic state(s) at the position of the corresponding auxiliary trajectories.

In Fig. 4 and 5, we summarize the steps of CTTSH and of SHXF. Fig. 4 shows that the procedure starts at time t = 0 with the sampling of the initial conditions (ICs) for the trajectories and the initialization of the electronic coefficients. Then, the electronic structure is computed at the positions of all trajectories, a task that can be performed in parallel, thus simultaneously for each trajectory. Note that in the scheme in Fig. 4 a thick gray arrow indicates a set of tasks, one for each trajectory, that can be performed in parallel. In CTTSH, the quantum momentum is calculated only once using information collected for all trajectories, therefore this task cannot be parallelized (it is indicated as a thin gray arrow). After the calculation of the quantum momentum, the nuclear trajectories and the electronic coefficients are propagated in time and the probability for a surface hop is determined, eventually imposing energy conservation if the hop is successful; these tasks can be parallelized over the trajectories. The procedure is repeated from the electronic-structure calculations at the new trajectories positions and iterated until the final time t_max of the simulation. Fig. 5 shows that all tasks can be performed in parallel, thus all steps in the scheme representing the algorithm are indicated as thin gray arrows. In addition, the calculation of the quantum momentum in SHXF is performed only after the auxiliary trajectories are created on the non-active states that are nonadiabatically coupled to the active state.

Fig. 4 Schematic representation of the CTTSH algorithm, with the thick gray arrows indicating the tasks that can be performed in parallel and with the thin gray arrow indicating the calculation of the quantum momentum that cannot be parallelized over the trajectories.

Fig. 5 Schematic representation of the SHXF algorithm, with all the gray arrows (5 in this case, so as to indicate the propagation of 5 physical trajectories) being represented as thin lines, since all the associated tasks can be performed trivially in parallel.

While the computational cost of SHXF is the same as “standard” surface hopping, CTTSH is more expensive. From the electronic-structure perspective, in addition to what is done in SHXF, CTTSH requires calculating the gradients of all electronic states involved in the dynamics (which is similar to ab initio multiple spawning²⁷ after the trajectories start spawning, if the stochastic-selection idea^216–218 is not applied, or to the augmented-fewest-switches surface hopping¹⁹⁷). Note that, even though the calculation of the quantum momentum scales as N_tr², it does not require additional electronic-structure information.

In ref. 209 we compared CTTSH to surface hopping and to a “synchronized” version of surface hopping using MOPAC-PI.²⁰⁸ This code uses a semi-empirical Floating Occupation Molecular Orbitals-Configuration Interaction (FOMO-CI) electronic structure approach, which is computationally very efficient, thus the time-limiting step in the overall dynamics is not the electronic structure, as it would be when using ab initio methods. The comparison was performed on the trans-azobenzene photodynamics using two electronic states, with ensembles of 10, 20 and 30 trajectories propagated for 50 fs (500 time steps) with computations carried out on a machine equipped with an Intel Xeon CPU E5-2450 0 with 32 CPUs, 2 sockets, 8 cores per package, and 2 threads per core. In surface hopping, the total time of the simulation was calculated as the sum of the time used in each trajectory. For CTTSH and “synchronized” surface hopping, the trajectories run synchronized in parallel using the same Message Passing Interface (MPI) strategy, but in the “synchronized” version of surface hopping, the quantum momentum is not calculated. The most significant time-consuming factor in CTTSH is the communication among processes necessary for the computation of quantum momentum. In Table 1 it is showed that the time spent on communication and synchronizing trajectories is significant, however, we expect that if a more expensive electronic-structure method was employed, the higher computational time required to perform CTTSH simulations would be less pronounced.

Table 1 Total time required (in seconds) to propagate 10, 20 and 30 trajectories with CTTSH and Tully SUrface Hopping (TSH) for 50 fs (500 time steps)

N tr	CTTSH	TSH	S-TSH^a
a S-TSH refers for the “synchronized” version of SH (see text).
10	3870.4	1705.4	2224.6
20	13 640.2	4239.2	7842.8
30	19 815.9	8899.7	11 452.5

---

Apart from an efficient parallelization for the simultaneous propagation of CTTSH trajectories, there are indeed some ways to make the algorithm affordable for molecular studies: using a local diabatization scheme for the evolution (see Section 3.2), such that the nonadiabatic couplings are not needed for the propagation;^194,208 using an efficient approach for the electronic structure information alternative to ab initio methods, as done with semi-empirical models or by building on-the-fly a database of nuclear configurations where the electronic structure is explicitly computed or by employing linear vibronic coupling models or with machine learning.^{208,210–212}

3.2 General observations on nonadiabatic dynamics with trajectories

It is important to mention here that, in general, the outcome of a nonadiabatic quantum molecular-dynamics simulation strongly depends not only on the algorithm used to propagate nuclear and electronic dynamics, but also (i) on the choice of the initial conditions, (ii) on the representation used for the electronic problem, and, perhaps mainly, (iii) on the electronic-structure approach employed to determine electronic properties on-the-fly along the trajectories.²¹⁹ Since the application of the exact factorization in the domain of photochemistry is independent of point (iii), below, we will briefly discuss only points (i) and (ii), as they have been analyzed precisely in the framework of the exact factorization.^{113,145,150,180} In relation to point (iii) above, some ideas have been developed in the context of the time-independent formulation of the exact factorization by Requist and Gross in ref. 120. There, it was proposed that the electronic equation, equivalent to eqn (5) but for the stationary solutions of the electron–nuclear Hamiltonian, is expressed and solved using density functional theory^220,221 and some test studies were reported.¹²¹ In addition, Min and coworkers recently recast the exact factorization in such a way that a trajectory-based scheme can be made operational in combination with real-time time-dependent density functional theory.¹⁸⁹

An adequate preparation of the initial state of the system, point (i) above, depends on the process under investigation or on the experiment that one wants to simulate, and, thus, a general “recipe” for the correct procedure does not exist. Our group addressed this issue in earlier work focusing independently either on the nuclear¹⁸⁰ or on the electronic¹⁴⁵ aspects within the exact-factorization framework, but the problem has been discussed in the literature by other authors.^222–228 The choice of initial conditions, namely the nuclear positions and momenta together with the electronic coefficients/populations at the initial time, can affect the calculated dynamics and observables, like photoproducts, lifetimes, and quantum yields, in many different ways. The dependence on the initial nuclear distribution might result more severe for the exact-factorization-based coupled-trajectory methods than for other approaches, since properly reconstructing the nuclear density along the dynamics is crucial to compute the quantum momentum that enters the expression of the classical force, in CTMQC (16), and of the evolution equation of the electronic coefficients, in CTMQC (15) and in CTTSH (17).

On the one hand, in recent work on the photo-isomerization dynamics of trans-azobenzene using CTTSH (see Section 4.2), it was necessary to freeze some high-energy internal vibrations while performing the initial harmonic Wigner sampling to avoid creating highly energetic trajectories: due to the limited number of affordable CTTSH trajectories (100 and 150 coupled trajectories were used) in such a large nuclear configuration space (of 72 dimensions), the energetic trajectories were rapidly separating in space and needed to be decoupled from the rest of the ensemble for numerical convenience, essentially reducing CTTSH to standard surface hopping. With a similar idea in mind, in order to avoid ensembles, or bundles, of trajectories that dramatically diverge in configuration space, the independent bundle approximation (IBA) was proposed using CTMQC.⁴⁴ The IBA consists in organizing N_tr trajectories in n bundles each containing n_tr trajectories, such that N_tr = n × n_tr, where the n_tr trajectories are grouped in the bundles according to their initial kinetic energy (or potential energy, or total energy, etc.). On the other hand, the initialization of the electronic dynamics in a pure adiabatic state is the commonly-accepted practice for photochemical studies, since with a femtosecond laser pulse it is possible to target the excitation of a low-lying electronic state. Nonetheless, photo-ionization with an attosecond laser pulse, broadly distributed in energy, is capable of creating a coherent superposition of electronic states (see Section 4.3). For applications in attochemistry, then, the way of creating such a superposition should be adapted to the used simulation method as well as to the observables under investigation, as discussed in ref. 145.

Point (ii) above concerns the theoretical representation adopted to describe the electronic dynamics. When aiming to simulate the ultrafast relaxation process of a photo-excited system via conical intersections or avoided crossings between electronic states of the same spin multiplicity, i.e., internal conversions, the Hamiltonian given in eqn (1) is able to fully capture this dynamics. Therefore, usually, the adiabatic representation is employed, formed by the eigenstates of the electronic Hamiltonian, i.e., Ĥ_BO(r, R) in eqn (1). In this representation, a peculiar parametric dependence on the nuclear configuration R appears, as consequence of the dependence of Ĥ_BO on R, and stimulates steadily interests in various communities due to its potential relation with geometric phases^229–239 (see Section 4.1). The form of the electronic Hamiltonian Ĥ_BO is diagonal in this representation and the coupling between pairs of electronic states is encoded in the nonadiabatic coupling vectors – connecting two electronic states via the nuclear displacement operator. Due to the very localized nature of these nonadiabatic coupling vectors, local diabatization schemes are often employed to solve the electronic evolution equation locally in the diabatic basis, designed to capture even with a finite integration time step the localized degeneracy region between two electronic states along the nuclear evolution. When describing processes as intersystem crossings involving electronic states with different spin multiplicity, the spin–orbit coupling Ĥ_SOC(x, R), with x = r, σ, has to be included in the molecular Hamiltonian to open spin-forbidden relaxation pathways upon photo-excitation. In such situations, the question presents itself as to whether the spin-diabatic or the spin-adiabatic representation is most suited for the simulations.^{113,158,186,240} The former is composed by the eigenstates of Ĥ_BO, and the pairs of electronic states are coupled either via the nonadiabatic coupling vectors or via the spin–orbit couplings, depending on the spin multiplicity; the latter is formed by the eigenstates of Ĥ_BO + Ĥ_SOC, and the couplings between pairs of states becomes fully of nonadiabatic character. Also in this case, the literature presents various options depending on the algorithm and on the code employed for the simulations. Dynamics in the presence of spin–orbit coupling was analyzed based on CTMQC^113,186 and pointed out that, while at the exact-factorization level there is no difference between the two representations, as expected, the approximations underlying the CTMQC algorithm perform better in the spin-adiabatic basis than in the spin-diabatic basis. A very similar issue is encountered when an external laser field is added to the molecular Hamiltonian, i.e., Ĥ_BO(r, R) + Ĥ_ext(r, R, t), since now the ideas of field-adiabatic or field-diabatic bases can be introduced.²⁴¹ Also in this case, authors discussed this point in relation to various trajectory-based algorithms for nonadiabatic dynamics, and in the context of the exact factorization only the particular case of a periodic drive was considered, and CTMQC was combined with the Floquet formalism,¹⁵⁰ designed to solve a TDSE with a Hamiltonian that is periodic in time.

4 Studies of electron–nuclear dynamics based on the exact factorization

This section of the perspective provides an overview of the variety of dynamical problems that have been tackled using the exact factorization. In the following sections, we will demonstrate the flexibility of the formalism in describing different effects and different properties, often providing new insights into the studied processes or new understanding of the theory.

First, in Section 4.1, we will report on a dynamics in the vicinity of a conical intersection in a Jahn–Teller model, recently engineered in a quantum simulator.²³⁹ In this study, it was discussed the fundamental difference between the concept of topological phases due to the presence of conical intersections that arise in the Born–Huang representation of the dynamics and the concept of dynamics-induced geometric phase that is independent of the underlying theoretical representation used to analyze the dynamics. Then, in Section 4.2, we report on various applications of exact-factorization-based algorithms to study the dynamics of molecules induced by the creation of a vibrational wavepacket produced by photoexcitation. In Section 4.2.1, we describe in some detail the work of our group that employs a coupled-trajectory scheme, in particular, the photo-isomerization process of trans-azobenzene upon nπ* and ππ* excitation with CTTSH in full dimensionality using the semi-empirical FOMO-CI electronic structure method with the reparametrized semi-empirical electronic-structure model AM1. The results of CTTSH dynamics are compared to surface hopping w/o accounting for decoherence corrections. In Section 4.2.2, an overview of the auxiliary-trajectory scheme SHXF is provided, highlighting the works of Min and of Maitra on a large variety of molecular systems and properties. Finally, in Section 4.3, we discuss a recent analysis reported on the different strategies to initialize the dynamics simulated by various trajectory-based procedures for excited-state dynamics in a model potential upon creation of a coherent superposition of electronic states,¹⁴⁵ usually referred to as electronic wavepacket. Such an electronic wavepacket can be created by an attosecond laser pulse, and the ensuing dynamics is usually simulated using the (multi-trajectory) Ehrenfest method, whose performance in comparison to CTMQC have been discussed.

4.1 Low-energy dynamics at a conical intersection

As discussed in the previous sections, the exact factorization provides a fundamentally different picture of photochemistry if compared to the Born–Huang representation, as it introduces a single time-dependent potential governing the nuclear evolution, in the form of a TDPES and a TDVP, bypassing concepts such as multiple static adiabatic potential energy surfaces, conical intersections^{104–106,234,236,239,242–244} and nonadiabatic couplings.²⁴⁵

When a dynamical process is studied in the Born–Huang representation, conical intersections appear ubiquitous in molecules: since the electronic Hamiltonian depends on the nuclear positions, its eigenstates, i.e., the adiabatic states, and its eigenvalues, the adiabatic potential energy surfaces, are functions of the nuclear positions as well. Conical intersections are regions of degeneracy of the adiabatic energies, where the correlation between electrons and nuclei manifests a singular behavior resulting in the divergence of the nonadiabatic couplings. In this picture, conical intersections give rise to topological phenomena of great interest for photochemistry and for low-energy molecular collisions or dissociations,^{231–235,238,246,247} which are related to geometric phases and often referred to as molecular Aharonov-Bohm effects.^237,248,249

In the light of the broad interest of the physical-chemistry and chemical-physics communities in conical intersections and geometric phases, several studies employing the exact factorization^{147,178,179,250,251} have been conducted on prototypical situations of coupled electron–nuclear dynamics in the vicinity of conical intersections. Overall, these works concluded that in the studied cases, the time-dependent potentials do not reflect the singular behavior of the adiabatic potential energy surfaces and nonadiabatic couplings at the positions of the conical intersections. Hence, the signatures of conical intersections (seem to) disappear in the eye of the exact factorization. These studies raised the question as to whether it is possible to identify observable effects related to conical intersections without relying on the Born–Huang, or any other theoretical, representation, since physical observables are independent of the representation. This question is strictly related to a longstanding debate about the possibility of providing experimental evidence of effects directly related to geometric phases in molecules.^{232,233,235,238,239,246,247}

Particular cases of geometric phases, which are those usually encountered in molecular processes within the Born–Huang representation, are the topological phases. When an electronic adiabatic eigenstate is slowly varied along a closed path in nuclear configuration space encircling the conical intersection, the corresponding wavefunction picks up a phase, in addition to the dynamical phase arising from the Schrödinger equation: if the phase depends on the geometry of the path, it is called geometric, whereas if it depends on the winding number of the path around the conical intersection, it is called topological. For a real-valued adiabatic wavefunction, the topological phase is γ = 0 or π and the phase factor is e^iγ = ±1, yielding a change of sign of the electronic wavefunction when transported along the path. This idea of adiabatic transport of the electronic wavefunction in nuclear space implies the BO approximation. In an explicitly dynamical formulation, the BO approximation states that the molecular wavefunction is

Ψ_BO(r,R, t) = ϕ^BO(r; R)χ^BO(R, t) (19)

thus only one electronic state ϕ^BO(r; R), e.g., the ground state, is considered in the dynamics. On the other hand, the exact factorization yields the molecular wavefunction as in eqn (11), i.e., Ψ(r,R, t) = Φ(r, t; R)χ(R, t), which is formally similar to eqn (19), as it is a single product of a nuclear wavefunction and of an electronic wavefunction with a parametric dependence on R. In the BO approximation, one usually introduces a static potential energy surface and static vector potential, or derivative coupling, namelywhose expressions are reminiscent of the TDPES (7) and of the TDVP (6). The static vector potential of eqn (20) has been discussed extensively in the literature,^248,249 since it appears in connection to conical intersections and geometric phases, as it will be discussed below, and it arises in the so-called extended BO equations.^252–256 The latter formulation was proposed in order to circumvent the fact that nonadiabatic problems in the Born–Huang framework require solving a set of coupled equations. Instead, the extended BO equations produce single-potential Schroödinger-like equations with the effect of the nonadiabatic coupling vectors reformulated in terms of a vector potential. This procedure is applicable only under certain conditions, i.e., the vector matrix of nonadiabatic coupling vectors needs to be expressed as a product of a vector function and an antisymmetric scalar matrix. Instead, the TDVP has more general validity, even in the formulation of the exact factorization for the stationary states, which is where the extended BO equations were applied. Nonetheless, we might expect strong similarities between the extended BO equations, when valid, and the nuclear equation(s) of the exact factorization.

The topological phase within the BO approximation γ^BO(Γ) is defined as the circulation of the corresponding vector potential along the path Γ that encloses the conical intersection,

An analogous quantity can be introduced using the TDVP, yieldingWhile γ^BO(Γ) = 0 or π is a topological phase, and is thus an intrinsic property of the Hamiltonian of the system, γ(Γ, t) depends on the dynamics that the system undergoes and is a geometric quantity. In addition, since the TDVP is related to the nuclear momentum field of the molecular system as shown in eqn (8), its circulation is gauge-independent and representation-independent. Therefore, eqn (22) is a robust physical observable. Note that, in addition, the curl of the TDVP is a physical observable as well, since it is the curl of the nuclear momentum field.†

Preliminary studies on geometric phases in the exact factorization were conducted by Gross and coworkers employing its static formulation and drew different conclusions depending on the systems that were considered.^124,125,257 In ref. 124, for instance, it was shown that, in the studied model for proton coupled-electron transfer, the concept of topological phase was merely a consequence of the BO approximation. Thus, it completely disappeared in the framework of the exact factorization, only to reappear as a topological quantity in the limit of infinite nuclear mass, i.e., the BO limit. Instead, in ref. 125, a model of a pseudorotating triatomic molecule²⁵⁸ was used to show that the topological phase reverts to a geometric quantity if the system manifests a degenerate ground state, even within the exact factorization. These studies suggest, then, that the topological phase of the Born–Huang representation becomes geometric in particular cases related to degeneracies of the full Hamiltonian, while neither topological nor geometric phases appear otherwise. Later on, to shed some light on the time-dependent picture of geometric phases in the exact factorization and to relate theoretical observations to experiments, our group investigated the dynamics generated by the E⊗e Jahn–Teller Hamiltonian recently engineered in a trapped-ion quantum simulator.²³⁹ The experiment allowed to directly image the evolution of the nuclear density at the conical intersection in order to detect directly the signature of the conical intersection. The working hypothesis was that, since the real-valued electronic wavefunction changes sign when transported around the conical intersection, in order to yield a single-valued molecular wavefunction, the nuclear wavefunction in eqn (19) needs to change sign as well. Therefore, the nuclear wavefunction has to develop a node or a zero-density line, which are indications of destructive interferences caused by the geometric phase. However, this hypothesis relies on the BO approximation and, as shown in ref. 250, it cannot be verified when employing the exact factorization.

In ref. 250 and 251, the exact quantum dynamics was simulated using the same Hamiltonian as in ref. 239 and compared to the BO approximation, by calculating and analyzing the nuclear density, the TDPES and the TDVP of the exact factorization. Overall, the observations allowed to discard all traces of singularities or destructive interferences related to topological-phase effects. As an example, let us discuss here the TDPES in Fig. 6, which is shown as colour map at three times along the dynamics. First, in the left panel, the nuclear density (depicted as black contour lines) approaches the location of the conical intersection (indicated by a red cross at the origin X, Y = 0, 0) at 0.24 ms; then, in the central panel, the density has reached the location of the conical intersection at 0.69 ms; finally, at 1.35 ms, the density has passed the origin X, Y = 0, 0 and exhibits a depletion along the line X > 0, Y = 0. In particular at this last time, it is interesting to observe that the TDPES develops a barrier-like feature along Y = 0 for X > 0. However, while the height and width of this barrier change over time, it remains finite and does not diverge. Thus, this barrier induces a depletion of the nuclear density along Y = 0 but cannot cause a node in the corresponding nuclear wavefunction. Indeed, a finite (strictly larger than 10⁻⁵), yet small, amount of nuclear density is always present in the region X > 0, Y = 0, attesting to the fact that there is no nodal line formed and thus, no topological-phase effect induced in the nuclear dynamics.

Fig. 6 The color maps show the TDPES at different times along the simulated dynamics of the E⊗e Jahn–Teller model of ref. 239. The contour lines show the nuclear density at the same times. The cross at X, Y = 0, 0 indicates the position of the conical intersection. Reprinted with permission from J. Phys. Chem. Lett., 2023, 14 (51), 11625–11631. Copyright 2023 American Chemical Society.

As indicated in eqn (7), the geometric phase arising in the exact factorization is calculated using the TDVP. While the TDVP is a gauge-dependent quantity, its curl and its circulation are not, thus they are physical observables. In Fig. 7, we show the curl of the TDVP at three different times along the simulated dynamics, and we compare it to the curl of the nonadiabatic coupling vectors of the Born–Huang representation. Note that in the particular case of a two-state system with real adiabatic eigenstates, γ^BO(Γ) of eqn (21) can be calculated by replacing A^BO_ν(R) with the nonadiabatic coupling vectors. Therefore, we find instructive to compare the two quantities, i.e., the curl of the TDVP and the curl of the coupling, that give rise to the geometric phases. In Fig. 7 it is clear that the two quantities are qualitatively different, one being smooth everywhere in space and the other presenting a particular behavior only at the location of the conical intersection. Furthermore, we observe that the TDVP is not irrotational and, thus, cannot be eliminated by a suitable choice of gauge. It is clearly dependent on the dynamics of the electron–nuclear system while the nonadiabatic coupling vector is an intrinsic property of the system Hamiltonian itself.

Fig. 7 The color maps show the curl of the TDVP at different times along the simulated dynamics of the E⊗e Jahn–Teller model of ref. 239 and the curl of the nonadiabatic coupling vectors (NACV) for the same model. The circles show the paths along which the line integrals of eqn (21) and (22) are performed. Reprinted with permission from J. Phys. Chem. Lett., 2023, 14 (51), 11625–11631. Copyright 2023 American Chemical Society.

The geometric phases, γ^BO(Γ) and γ(Γ,t), have been calculated from eqn (21) and (22), respectively, and are reported in Table 2. We used different paths Γ, namely Γ_CI1 and Γ_CI2 encircling the position of the conical intersection, and Γ_X and Γ_Y that do not encircle the position of the conical intersection. The paths are shown in Fig. 7. Table 2 shows that in the adiabatic representation, namely for γ^BO(Γ), a value close to π is obtained only for Γ_CI1 and Γ_CI2, and that we find a value much smaller, i.e., close to zero, along the other paths. Using the TDVP, all the values change with time and with the path. Therefore, our numerical calculations show that the result of the integral in eqn (22) depends on the behavior and on the dynamics of the TDVP, i.e., the smooth function to be integrated, as well as on the geometry of the chosen path: γ(Γ, t) is a geometric phase. Instead, the integral of the nonadiabatic coupling vectors used to determine eqn (21) over the same paths chosen for eqn (22) only depends on the winding number of the path(s) around the point of singularity of the function to be integrated: γ^BO(Γ) is a topological phase. Thus, in ref. 250 we introduced the concept of dynamics-induced geometric phase to indicate γ(Γ, t), a concept that arises in the exact factorization. The dynamics-induced geometric phase is not a topological feature but rather caused and governed by the dynamics of the system, in stark contrast with the topological, path-independent phase that arises only in the adiabatic representation.

Table 2 Geometric phase γ(Γ, t) at three times throughout the dynamics (first three columns) along four different paths, encircling the point representing the location of the conical intersection (Γ_CI1, Γ_CI2) and without encircling this point (Γ_X, Γ_Y). The last column reports the values of γ^BO(Γ) using the same paths

	0.24 ms	0.69 ms	1.35 ms	Adiabatic
Γ CI1	−0.0209	0.0265	0.0506	3.12
Γ CI2	−0.794	2.05	−0.00973	3.12
Γ X	0.00269	0.0624	0.000178	−0.00281
Γ Y	−0.758	0.157	0.188	0.000581

---

4.2 Dynamics induced upon creation of a vibrational wavepacket

In Section 4.2.1, we report on our recent study of the photodynamics of trans-azobenzene using CTTSH, whereas in Section 4.2.2 we present some examples of the work by Min and by Maitra with SHXF.

4.2.1 Coupled-trajectory methodologies. The photoisomerization mechanism of trans-azobenzene was investigated using the CTTSH approach²⁰⁹ in comparison to the Tully surface-hopping (TSH) scheme, in order to assess the performance of the implementation of CTTSH in MOPAC-PI.²⁰⁸ This study represented the first application of CTTSH to a molecular system in full dimensionality, as previous applications of CTTSH have been limited to model systems of one or two dimensions, typically comprising only two electronic states. While the theoretical foundations of the CTTSH approach have been previously discussed for low-dimensional model systems, the algorithm required customization to address numerical and computational challenges encountered in simulations involving systems with a large number of nuclear degrees of freedom.

In ref. 209, the photoisomerization dynamics was investigated upon nπ* and ππ* excitation, describing and comparing a two-state dynamics and a three-state dynamics, respectively. The semi-empirical FOMO-CI method was used to evaluate the electronic properties required for the on-the-fly dynamics, entering the evolution equation of the electronic coefficients (17) and the expression of the force needed for the nuclei classical propagation (18). The semi-empirical FOMO-CI method demonstrates high computational efficiency and reasonable accuracy when properly parameterized, making it suitable for computing multiple trajectories for medium to large molecular systems over time scales of up to a few picoseconds.

In CTTSH the force for the propagation of the nuclear trajectories is the same as in TSH, whereas CTTSH includes an additional term in the propagation of electronic coefficients that depends on quantum momentum, as shown in eqn (17). TSH is known to suffer from a systematic issue termed “overcoherence”, arising from the disconnect between the evolution of electrons and nuclei. Throughout time, the nuclei evolve on a single adiabatic potential energy surface and are allowed to hop to another surface at any time based on a stochastic algorithm. Therefore, on the one hand, the nuclear dynamics is governed by a single electronic states, while, on the other hand, the electronic evolution remains always coherent along a trajectory and yields a superposition of adiabatic states along that trajectory. CTTSH, thanks to the coupled-trajectory term of eqn (17), improves the description of quantum decoherence compared to TSH as it accounts in the electronic evolution equation for information about the spatial delocalization of the trajectories. This information is encoded in quantum momentum (14), which is related to the spatial variation of the nuclear density. Then, the nuclear density needs to be constructed from the distribution of trajectories, and this is done in CTTSH, as in CTMQC, with a sum of normalized Gaussians centered at the positions of trajectories. Therefore, at a given time t, the positions of all the trajectories need to be “shared” among the trajectories to construct the nuclear density and, thus, the trajectories cannot be evolved independently. From a practical standpoint, this is the main drawback of CTTSH when compared to standard TSH. For an efficient performance of CTTSH, a parallelized implementation using MPI was employed in ref. 209 (see Table 1). Despite the additional computation required for the quantum momentum, the runtime of our parallelized implementation is similar to that of a standard TSH approach, as all other calculations are performed in parallel. Furthermore, a local diabatization algorithm was utilized within the framework of the Runge–Kutta integrator to bypass the computation of the expensive and sometimes unavailable nonadiabatic coupling vectors.

This work on the photodynamics of trans-azobenzene demonstrated that, in CTTSH, the accuracy of the approximation of the nuclear density is crucial for a proper description of the quantum momentum and, consequently, for a correct account of decoherence effects in the electronic dynamics. Indeed, reconstructing the nuclear density for a polyatomic molecule is much more challenging than for a few-dimensional system. Therefore, we proposed four variants of the CTTSH method depending on the set of nuclear coordinates considered to construct the quantum momentum and depending on the approaches used to determine the widths of the Gaussian functions. The implementations of CTTSH can be summarized as follows:

1. CTTSH-FGAC (CTTSH with Frozen Gaussians and All Coordinates): in this approach, we utilize frozen Gaussians and all coordinates contribute to the calculation of the quantum momentum.

2. CTTSH-TGAC (CTTSH with Thawed Gaussians and All Coordinates): in this approach, we employ Gaussians whose widths change over time, and all coordinates contribute to the calculation of the quantum momentum.

3. CTTSH-FGLA (CTTSH with Frozen Gaussians and List of Active Atoms): in this approach, we utilize frozen Gaussians and only “active atoms”, whose list is provided as an input, contribute to the calculation of the quantum momentum.

4. CTTSH-TGLA (CTTSH with Thawed Gaussians and List of Active Atoms): in this approach, we utilize thawed Gaussians and only “active atoms” contribute to the calculation of the quantum momentum.

To assess the performance of these four flavours of CTTSH implementations, we compared our numerical results with TSH and TSH corrected with the inclusion of overlap-based decoherence corrections²⁰⁰ (TSH-ODC).

The minimum energy path on the S₁ potential energy surface of azobenzene connects the trans (and cis) geometries to the S₁/S₀ conical intersection, located at approximately 95° of torsion of the CNNC dihedral. Consequently, nonadiabatic transitions to S₀ in isolated trans-azobenzene occur predominantly in this region or slightly before, within a CNNC range of 90 to 120°. During the decay to S₀, this molecule can isomerize to cis-azobenzene. The partial‡ photoisomerization quantum yields Φ_nπ* calculated using all four CTTSH methods, TSH, and TSH-ODC for the trajectories starting in the nπ* state are shown in the second column of Table 3. The calculated Φ_nπ* values are consistent across methods, are in agreement with each other and in good agreement with experimental results, which fall within the range of 20–32%.^259–266 It is worth noting that trajectories with faster decay predominantly remain in the trans configuration after relaxing to the ground state, whereas slower trajectories have a higher probability of undergoing photoisomerization. In the third column of Table 3, we report the values of the partial photoisomerization quantum yields Φ_ππ* upon ππ* excitation. The lower photoisomerization quantum yield after ππ* excitation has been previously discussed.^267–269 It is related to an exception to Kasha's rule, resulting from a competition between “reactive” and “unreactive” internal conversion at the S₁/S₀ conical intersection, i.e., between internal conversion with photoisomerization or without photoisomerization, respectively. Both processes require a certain degree of progress along the CNNC torsional coordinate. However, the unreactive S₁ → S₀ internal conversion process can occur earlier (i.e., farther from the 90° midpoint of the torsional pathway) if more vibrational energy is available, as is typically the case with ππ* excitation. TSH-ODC, CTTSH-FGAC, and CTTSH-TGAC approaches presented photoisomerization quantum yields for ππ* in a very good agreement with the ones reported experimentally (0.09–0.16), whereas TSH, CTTSH-FGLA, and CTTSH-TGLA were unable to do so. These last three approaches showed higher or equivalent photoisomerization quantum yields for ππ* and nπ* excitation. However, it is important to emphasize that only approximately 40% of the trajectories for the TSH and CTTSH-TGLA methods successfully reached the ground state by the end of our simulations. Therefore, it is challenging to attribute a “precise” quantum yield in such cases. In addition, our simulations show that when using only the CNNC atoms for the calculations of the quantum momentum, the number of hops and back-hops between S₁ and S₂ increases if compared to the other CTTSH implementations. Thus, it is our interpretation that the non-satisfactory account of decoherence in the thawed-Gaussian method combined with the increased number of hops and back-hops in the list-of-active-atoms option bring the trajectories to explore regions of configuration space where the molecule isomerizes with higher probability.

Table 3 Partial photoisomerization quantum yields (Φ) calculated with our four implementations of CTTSH, with TSH and TSH-ODC. The second column reports the results for the dynamics initiated in the nπ* state and the third column reports the results for the ππ* excitation

Method	Φ nπ*	Φ ππ*
CTTSH-FGAC	0.22 ± 0.04	0.11 ± 0.04
CTTSH-FGLA	0.19 ± 0.03	0.19 ± 0.05
CTTSH-TGAC	0.23 ± 0.04	0.13 ± 0.04
CTTSH-TGLA	0.20 ± 0.04	0.38 ± 0.08
TSH	0.18 ± 0.04	0.17 ± 0.05
TSH-0DC	0.22 ± 0.03	0.13 ± 0.04

---

Fig. 8 shows the populations of the ground state S₀ (in red) and of the excited state S₁ (in blue) as functions of time for the simulation starting after the nπ* excitation. The four panels refer to the four flavors of CTTSH implementations presented above. Fig. 9 reports the analogous quantities calculated using TSH and TSH-ODC. In Fig. 8 and 9, the thick lines represent the populations calculated as the fraction of trajectories in a given electronic state m, namely as the ratio F_m(t) = N_m(t)/N_tr, where N_m(t) is the number of trajectories for which the active state at time t is m, and N_tr is the total number of trajectories; the thin lines represent the electronic population estimated as the average P_m(t) over the trajectories of |C^α_m(t)|², where C^α_m(t) are the electronic coefficients along the trajectory α. In general F_m(t) ≠ P_m(t), and the difference is particularly severe if decoherence effects are important and not accounted for, like observed for TSH in Fig. 9. The decoherence corrections used in TSH-ODC are designed to restore the internal consistency of the TSH procedure in an ad hoc manner, resulting in F_m(t) ≃ P_m(t), as seen in Fig. 9. Instead, we note a remarkably good agreement across implementations of CTTSH with respect to the electronic populations evaluated as F_m(t) or as P_m(t). Internal consistency, thus F_m(t) ≃ P_m(t), is correctly achieved with CTTSH when utilizing frozen Gaussians (fixed widths) all along the simulated dynamics, instead, the agreement between F_m(t) and P_m(t) decreased over long times when thawed Gaussians (variable widths) are employed to calculate quantum momentum. This happens because, as the trajectories become delocalized in space, the widths of the individual Gaussians used for reconstructing the nuclear density increase. Consequently, the quantum momentum, which is evaluated as the spatial derivative of the density, decreases, leading to a reduced effect of decoherence. Nevertheless, we believe that by significantly increasing the number of coupled trajectories, we can attain an accurate representation of the nuclear distribution in configuration space without requiring a substantial increase in Gaussian width. In an ideal scenario, each trajectory would be strongly coupled to several others, making the nuclear density almost insensitive to variations in Gaussian widths within a reasonable range. However, achieving this would necessitate a much larger number of trajectories than the 100–150 (total number of trajectories used in this study) that can currently be accommodated on a single computing node, with the required number increasing with the number of coordinates.

Fig. 8 Populations of the S₀ state (in red) and of the S₁ state (in blue) as functions of time obtained from the calculations of the photodynamics of trans-azobenzene upon nπ* excitation using the four implementations of CTTSH. Thick lines refer to the fraction of trajectories F_m(t) and thin lines refer to the average over the trajectories of the electronic populations P_m(t). Reprinted with permission from J. Chem. Theor. Comput., 2024, 20 (2), 580–596. Copyright 2024 American Chemical Society.

Fig. 9 Same as in Fig. 8 but for TSH and TSH-ODC. Reprinted with permission from J. Chem. Theor. Comput., 2024, 20 (2), 580–596. Copyright 2024 American Chemical Society.

4.2.2 Auxiliary-trajectory methodologies. As presented in Section 3.1, an efficient technique to combine the surface-hopping algorithm with the exact factorization has been proposed by Min and coworkers and dubbed SHXF, i.e., surface hopping based on the exact factorization.⁴⁸ The main idea at the basis of SHXF is to use the surface-hopping procedure to propagate independent nuclear trajectories adiabatically, with the fewest-switches probability for the hops, and, at the same, to “spawn” auxiliary trajectories on the electronic states that are nonadiabatically coupled to the active/force state. The auxiliary trajectories are only used for the calculation of the quantum momentum (14), similarly to what is done in the original CTMQC algorithm: the quantum momentum induces quantum decoherence effects in the electronic evolution equation, thus bypassing the overcoherence problem of standard fewest-switches surface hopping.^59,66 SHXF is an interesting alternative to CTMQC and CTTSH, as it is computationally less demanding, due to the fact that for each physical trajectory, only a few auxiliary trajectories need to be generated in the vicinity of a region of strong coupling as long as such a region remains well localized in space. SHXF has been applied to a large variety of molecular systems to predict their nonadiabatic behavior upon photoexcitation. Min and co-workers studied the photoinduced isomerization of the trans-penta-2,4-dieniminium cation²⁷⁰ and the ring-opening in cyclohexa-1,3-diene²⁷¹ to calculate their excited-state lifetime, isomerization quantum yield and branching ratio, they worked on a class of molecular rotary motors that utilize the fulgide motif^272,273 to determine their isomerization quantum yields, they studied the effect of the solvent on the isomerization quantum yield of light-driven rotary motors,²⁷⁴ and they simulated the primary steps of the photodynamics of cyclopropanone²⁷⁵ to compute the dissociation quantum yields of the various photoproducts. These simulations have been performed using the PyUNIxMD code,^193,276 mainly in combination with spin-restricted ensemble-referenced Kohn–Sham method,^270,271 that allows for a proper characterization of the conical intersections with the ground state. It is worth mentioning, however, that Min proposed recently an interesting generalization of the exact factorization and of SHXF such that the electronic evolution equation accommodates for a hermitian formulation in terms of the electronic density. With this proposition, the photodynamics of ethylene was studied employing real-time time-dependent density functional theory.¹⁸⁹

In order to assess the performance of various surface-hopping procedures, including SHXF, Maitra and coworkers presented an in-depth analysis of the photodynamics of ethylene, fulvene and the methaniminium cation employing ab initio multiple spawning as benchmark.²⁷⁷ Specifically, the work analyzed various decoherence-correction schemes, which are found to operate in very different ways on the individual trajectories while yielding comparable results when averaged over the trajectories. In addition, the choice of the velocity-rescaling algorithm as well as the influence of nuclear time step for the propagation were studied as they were found to have an equally, if not more, important role on electronic and nuclear properties compared to the decoherence correction. Finally, Maitra and coworkers showed that the new mechanism for electronic transitions given by the quantum momentum in the exact-factorization based methods, and in particular SHXF, goes beyond a mere decoherence correction. The quantum-momentum induced effects were shown to be essential to capture the dynamics of the uracil cation as it is driven by a three-state conical intersection. This cannot be correctly captured using corrections to the electronic coefficients that involve only a pair of states (active and non active state).^114,278

4.3 Dynamics induced upon creation of an electronic wavepacket

The creation of an electronic wavepacket can be achieved following the interaction of a molecule with an attosecond laser pulse^279–282 that ionizes the molecule. Such a photoionization is able to trigger an ultrafast process that implies electronic charge redistribution coupled to the nuclear motion. The ultimate goal of the emerging field of attochemistry⁷¹ is to be able to control the evolution of such a wavepacket that implies the control of chemical reactivity, which is potentially achievable by tuning the electronic coherences among the states which are coupled via a second probe pulse or via the nuclear dynamics.

An accurate description of electron–nuclear correlation effects, including coherence and decoherence, is essential for understanding and controlling the products of a photochemical reaction. However, the theoretical description of these processes presents several well-known challenges, and in ref. 145 we focused on the problem of the initialization of the electronic dynamics in trajectory-based simulations. The initial conditions for the classical-like nuclear dynamics {R^α_ν, P^α_ν} are often sampled from the (harmonic) Wigner distribution of the ground state nuclear wavepacket assuming vertical excitation or from a classical Boltzmann distribution at some given temperature.¹⁸⁰ However, several options can be proposed for the initialization of the electronic dynamics via the electronic coefficients {C^α_k(0)} when more than one electronic state needs to be populated to initiate the nonadiabatic simulations. In ref. 145, we defined and distinguished two situations, a pure state, where all trajectories carry the same set of electronic coefficients at the initial time, and a mixed state, where each trajectory is allowed to carry a different initial superposition of electronic states in such a way that the ensemble-average electronic populations at the initial time reproduces the “expected” population distribution of the molecular wavefunction after the pulse, i.e., . Our work aimed to assess the performance of Ehrenfest and CTMQC in describing the nonadiabatic dynamics initiated by the creation of an electronic wavepacket distinguishing the two cases of a pure state and of a mixed state. We used a model system for our calculations, described by a two-electronic-state one-nuclear-mode Hamiltonian that consists, in the diabatic representation, of two coupled one-dimensional harmonic oscillators displaced in position and energy. The initial condition for the quantum dynamics is chosen as a coherent superposition of 80% and 20% of the ground state and of the excited state, respectively.

For such a model system, reference results can be provided by solving exactly the TDSE. Turning to the trajectory-based calculations, in the pure state, the initial electronic state for all trajectories is described by and whereas in the mixed state 80% of the trajectories are initialized by {C^(α)₀(0)} = 1 and the remaining 20% by {C^(β)₀(0)} = 1. Here, the index α runs over 80% of the trajectories, and β runs over the remaining 20%.

Fig. 10 shows the electronic excited state population and the electronic coherence as functions of time. Surprisingly, we observe that starting with a pure state in Ehrenfest dynamics yields negligible net population transfer with the ground and excited state remaining coherent throughout the dynamics. This behavior is in stark contrast with the expected result predicted by the quantum dynamics (QD) simulation. On the other hand, CTMQC captures both the population and coherence behavior in very good agreement with the reference. Note that all methods miss the recoherence just after 90 fs that, unlike the other recoherences, occurs far from the nonadiabatic coupling region. Details on this mechanism are discussed in ref. 283.

Fig. 10 Upper panel: Indicator of coherence as function of time. Lower panel: Population of the excited state as function of time. In black, the reference results of the quantum dynamics (QD) are reported. Trajectory-based results are indicated in green (for Ehrenfest with pure-state initialization), in blue (for Ehrenfest with mixed-state initialization) in red (for CTMQC with pure-state initialization), and in purple (for CTMQC with mixed-state initialization). Reprinted from J. Chem. Phys., 2024, 160, 054102, with the permission of AIP Publishing.

Let us focus on the first decoherence event occurring within the first 20 fs, during which the excited-state nuclear wavepacket moves away from the trapped ground state wavepacket before entering a region of strong nonadiabatic coupling. This situation represents the first step after a molecule gets excited via the laser pulse. Fig. 11 shows time snapshots of the exact nuclear density and the spatially-resolved energies (upper panel) and excited state populations (lower panel) during this event. We observe that Ehrenfest trajectories evolve in a mean-field surface which has 80% ground state character, and their electronic populations (and coherences) do not change away from regions of strong nonadiabaticity. Conversely, the quantum momentum term accounted for in the CTMQC equation of motion (15) induces both a branching of the electronic coefficients and a splitting of the trajectory distribution achieving decoherence. We now direct our attention to the mixed-state results, where both Ehrenfest and CTMQC capture the periodic population transfers, with Ehrenfest deviating from the reference at ∼120 fs. In terms of coherences, both methods cannot account for the initial decoherence, because in a mixed state the coherences are zero initially. However, after this first decoherence event, CTMQC captures the behavior of the indicator of coherence quite well, even though it misses the large recoherence at 90 fs as pointed out earlier; instead, Ehrenfest remains overcoherent after the first increase of coherence at ∼25 fs.

Fig. 11 Upper panel: Time snapshots of the exact density (orange), adiabatic potential energy curves (red lines), nonadiabatic coupling (magenta regions), TDPES (black), distribution of CTMQC trajectories starting in a pure ensemble (blue dots) and of Ehrenfest trajectories starting in a pure ensemble (green dots). Lower panel: Excited state populations as function of R or R^(α)(t) for QD (black line), CTMQC starting in a pure ensemble (blue dots) and Ehrenfest starting in a pure ensemble (green dots). Reprinted from J. Chem. Phys., 2024, 160, 054102, with the permission of AIP Publishing.

In summary, this work pointed out that the problem of the initialization of the electronic dynamics is critical when trajectory-based methods are used to simulate the dynamics triggered by an ultrashort laser pulse that creates an electronic wavepacket. We presented two strategies: a pure state and a mixed state. We tested these concepts in the simulation of the nonadiabatic dynamics of two-state electronic wavepacket in a one-dimensional model system. The pure state seems the natural choice for both Ehrenfest and CTMQC since all trajectories are associated to the same coherent superposition of electronic states. However, Ehrenfest performs very poorly, and if the mixed state is used instead, the population behavior greatly improves but the coherences, and consequently coherence-dependent observables such as the electronic current-density,¹⁴⁵ will be wrong from the very beginning. On the other hand, CTMQC starting from a pure state accurately captures both populations and coherences behavior. The quantum momentum is essential to capture those events at the spatially-resolved level, which is key to predict the right dynamics when regions of nonadiabatic electron–nuclear coupling are encountered.

5 Photon–electron–nuclear systems: the future of the exact factorization?

Little over a decade ago, the field of polaritonic chemistry started to emerge as an innovative way to steer photochemical reactions by creating hybrid light–matter states, i.e., the polaritons.^284–286 The strong-coupling regime between the molecular excitations and confined light can be achieved in optical and plasmonic cavities as it scales inversely with the volume of the cavity and is proportional to the dipole moment of the molecule.^17,287–304 The modifications of the photochemical landscape underlying a cavity-free molecular process are achieved when the bare electronic energies are altered by the coupling to the photons and produce the polaritonic states, whose energies, in turn, are affected by the number of molecules in the cavity. Theoretical modeling of polaritonic chemistry clearly requires a quantum-mechanical treatment of the coupled photon–electron–nuclear (PEN) problem and the extension of state-of-the-art simulation techniques for molecular dynamics to treat light as an intrinsic variable of the process rather than as an external field.

The fields of physical chemistry and chemical physics have sparkled with ideas for novel developments to study processes in the strong light–matter coupling regime, reformulating in cavities concepts such as the BO approximation^305–307 and nonadiabatic dynamics.^{15,18,19,144,302,303,308–310} In this respect, the exact factorization appears to be an extremely well adapted formalism to gain new insights into the correlated PEN dynamics and to develop simulation algorithms, in the spirit of similar efforts using surface hopping and the Ehrenfest scheme.

Maitra and coworkers,¹⁴¹ as well as Tokatly and coworkers,¹⁴⁰ laid the basis for the conception of the exact-factorization formalism to treat the PEN wavefunction. Specifically, the factored form of the wavefunction Ψ(r, R, t) = χ(R, t)Φ(r, t; R) has been generalized and used to analyze various kinds of problems, like photon–electron dynamics, with Ψ(r,q,t) = χ(q,t)Φ(r, t; R) interpreting q as the photon displacement and r as the electronic coordinate, or the full PEN dynamics, by interpreting R as the nuclear coordinates and introducing the polaritonic coordinates r, q. Maitra's preliminary work demonstrates that different factorizations yield complementary information on the fundamental interpretation of the PEN problem,^{142–144,311} as well as, perhaps most importantly, on the directions to follow to introduce approximations for an efficient computational treatment. For instance, employing the factorization Ψ(r, q, t) = χ(q, t)Φ(r, t; q), it was observed that even though the free-photon Hamiltonian is harmonic in the photon displacement coordinate and the light–matter coupling is bilinear in the nonrelativistic dipole approximation, the q-TDPES is far from harmonic. Therefore, a quasi-classical trajectory-based representation of the photon displacements and conjugated momenta should be carefully assessed and does not guarantee reproducing correctly the quantum dynamics. Specifically, the limitations of the multi-trajectory Ehrenfest approach have been discussed in the light of the different forces guiding the photon-displacement dynamics within such a mean-field theory and the exact factorization.¹⁴⁴ Furthermore, the factorization Ψ(r, R, q, t) = χ(R, t)Φ(r, q, t; R) was used to introduced the TDPES as an alternative tool to the polaritonic surfaces to interpret how photochemical reactivity is affected by the cavity in a model proton coupled-electron transfer.^142,311 Such an exact factorization naturally inspires simulations based on the coupled-trajectory schemes presented in previous sections, similarly to what the community already proposed using, for instance, surface hopping³⁰⁸ or ab initio multiple spawning³¹² and in the spirit of the Floquet-based CTMQC¹⁵⁰ method discussed briefly in Section 3.2. Nonetheless, this possibility has not been yet investigated.

Despite Maitra's observations on the difficulties in capturing correctly the photon dynamics using classical Ehrenfest trajectories, the exact factorization has the potential to shed light into alternative approximations of the PEN quantum dynamics where the photon displacements and momenta are treated with trajectories, perhaps within a CTMQC-like scheme. In the presence of many cavity modes¹⁴³ or even when envisaging to employ the formalism of the exact factorization to study the vibrational strong-coupling regime,^313–315 such a treatment might represent the only valuable route for efficient and accurate simulations, at the cost of rethinking the underlying approximations of CTMQC and similar algorithms. Developments in this direction are currently ongoing in our group and preliminary results have been reported elsewhere.³¹⁶

6 Conclusions

This Perspective provided an overview of the application of the exact factorization of the electron–nuclear wavefunction in the domain of photochemistry. After a brief presentation of the fundamental theory at the quantum-mechanical level, we introduced the concept of classical-like trajectories to be used to approximate the nuclear dynamics coupled to the quantum electronic dynamics and to derive algorithms for quantum molecular dynamics. Key concepts arising in the exact factorization are the time-dependent vector potential and time-dependent potential energy surface that completely account for the effect of the electronic ground and excited states on the nuclear dynamics. These time-dependent potentials can be interpreted classically such that the concept of classical nuclear force naturally arises. Another key ingredient of the exact-factorization equations is the quantum momentum, that tracks the spatial delocalization of the nuclear density, or equivalently of the distribution of classical trajectories, such that quantum-decoherence effects can be reproduced in the coupled-trajectory or auxiliary-trajectory schemes derived from the exact factorization. In this respect, we discussed previous work where the time-dependent potentials and the quantum momentum have shown their power in interpreting and simulating dynamical processes in the presence of strong nonadiabatic effects.

We concluded this perspective with some observations on the future of the exact factorization. In particular, we discussed some recent work in the field of polaritonic chemistry based on the introduction of the photon–electron–nuclear wavefunction together with some avenues for the extension to photochemistry in cavities of the algorithms discussed in the course of this perspective.

Data availability

No primary research results, software or code have been included and no new data were generated or analysed as part of this review.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the ANR Q-DeLight project, grant no. ANR-20-CE29-0014 of the French Agence Nationale de la Recherche, by a public grant from the Laboratoire d'Excellence Physics Atoms Light Matter (LabEx PALM) overseen by the French National Research Agency (ANR) as part of the “Investissements d'Avenir” program (reference: ANR-10-LABX-0039-PALM), and by the Chateaubriand Fellowship from the Mission Scientifique et Technologique of the Embassy of France in the United States.

Notes and references

1. A. Chenu and G. D. Scholes, Ann. Rev. Phys. Chem., 2015, 66, 69–96.
2. S. Gozem, H. L. Luk, I. Schapiro and M. Olivucci, Chem. Rev., 2017, 117, 13502–13565.
3. J. W. Snyder Jr., B. F. E. Curchod and T. J. Martínez, J. Phys. Chem. Lett., 2016, 7, 2444–2449.
4. H. Tamura and I. Burghardt, J. Am. Chem. Soc., 2013, 135, 16364–16367.
5. M. Gao, F. Yu, C. Lv, J. Choo and L. Chen, Chem. Soc. Rev., 2017, 46, 2237–2271.
6. Y. Olivier, J.-C. Sancho-Garcia, L. Muccioli, G. D'Avino and D. Beljonne, J. Phys. Chem. Lett., 2018, 9, 6149–6163.
7. L. González and R. Lindh, Quantum Chemistry and Dynamics of Excited States: Methods and Applications, John Wiley & Sons, Chichester, 2020.
8. H. Oberhofer, K. Reuter and J. Blumberger, Chem. Rev., 2017, 117, 10319–10357.
9. K. Takatsuka, Bull. Chem. Soc. Jpn., 2021, 94, 421–1477.
10. A. Kelly, Faraday Discuss., 2020, 221, 547–563.
11. G. Tao and W. H. Miller, J. Phys. Chem. Lett., 2010, 1, 891–894.
12. L.-M. Ibele, P. A. Sanchez-Murcia, S. Mai, J. J. Nogueira and L. González, J. Phys. Chem. Lett., 2020, 11, 7483–7488.
13. S. Mukherjee, M. Kar, M. Bhati, X. Gao and M. Barbatti, Theor. Chem. Acc., 2023, 142, 85.
14. M. Pérez-Escribano, J. J. G. Granucci and D. Escudero, J. Chem. Phys., 2023, 158, 124104.
15. J. Fregoni, G. Granucci, E. Coccia, M. Persico and S. Corni, Nat. Commun., 2018, 9, 4688.
16. M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel and A. Rubio, Nat. Rev. Chem., 2018, 2, 0118.
17. J. Feist, J. Galego and F. J. Garcia-Vidal, ACS Photonics, 2018, 5, 205–216.
18. H. L. Luk, J. Feist, J. J. Toppari and G. Groenhof, J. Chem. Theory Comput., 2017, 13, 4324–4335.
19. B. Gu, Y. Gu, V. Y. Chernyak and S. Mukamel, Acc. Chem. Res., 2023, 56, 2753–2762.
20. F. Agostini and B. F. E. Curchod, Philos. Trans. R. Soc., A, 2022, 380, 20200375.
21. F. Agostini and B. F. E. Curchod, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2019, 9, e1417.
22. G. A. Worth and I. Burghardt, Chem. Phys. Lett., 2003, 368, 502–508.
23. T. Ma, M. Bonfanti, P. Eisenbrandt, R. Martinazzo and I. Burghardt, J. Chem. Phys., 2018, 149, 244107.
24. G. Richings, I. Polyak, K. E. Spinlove, G. A. Worth, I. Burghardt and B. Lasorne, Int. Rev. Phys. Chem., 2015, 34, 269–308.
25. B. Lasorne, M. A. Robb and G. A. Worth, Phys. Chem. Chem. Phys., 2007, 9, 3210–3227.
26. G. A. Worth, M. A. Robb and I. Burghardt, Faraday Discuss., 2004, 127, 307–323.
27. B. F. E. Curchod and T. J. Martínez, Chem. Rev., 2018, 118, 3305–3336.
28. B. Mignolet and B. F. E. Curchod, J. Chem. Phys., 2018, 148, 134110.
29. D. V. Makhov, W. J. Glover, T. J. Martínez and D. V. Shalashilin, J. Chem. Phys., 2014, 141, 054110.
30. D. V. Shalashilin, Faraday Discuss., 2011, 153, 105–116.
31. A. J. Jenkins, K. E. S. M. Vacher and G. A. W. M. A. Robb, J. Chem. Phys., 2018, 149, 094108.
32. M. Vacher, M. J. Bearpark and M. A. Robb, Theor. Chem. Acc., 2016, 135, 187.
33. V. M. Freixas, A. J. White, T. Nelson, H. Song, D. V. Makhov, D. Shalashilin, S. Fernandez-Alberti and S. Tretiak, J. Phys. Chem. Lett., 2021, 12, 2970–2982.
34. R. E. Wyatt, Quantum dynamics with trajectories: introduction to quantum hydrodynamics, Springer Science & Business Media, 2005, vol. 28.
35. R. E. Wyatt and E. R. Bittner, J. Chem. Phys., 2000, 113, 8898–8907.
36. B. Poirier, Chem. Phys., 2010, 370, 4–14.
37. B. F. E. Curchod, I. Tavernelli and U. Rothlisberger, Chimia, 2012, 4, 174.
38. B. Gu and S. Garashchuk, J. Phys. Chem. A, 2016, 120, 3023–3031.
39. B. Gu and I. Franco, J. Chem. Phys., 2017, 146, 194104.
40. G. Albereda and I. Tavernelli, Quantum Chemistry and Dynamics of Excited States: Methods and Applications, John Wiley & Sons, 2021.
41. L. Dupuy, F. Talotta, F. Agostini, D. Lauvergnat, B. Poirier and Y. Scribano, J. Chem. Theory Comput., 2022, 18, 6447–6462.
42. L. Dupuy, A. Rikus and N. T. Maitra, J. Phys. Chem. Lett., 2024, 15, 2643–2649.
43. S. K. Min, F. Agostini and E. K. U. Gross, Phys. Rev. Lett., 2015, 115, 073001.
44. C. Pieroni and F. Agostini, J. Chem. Theory Comput., 2021, 17, 5969.
45. C. C. Martens, J. Phys. Chem. Lett., 2016, 7, 2610–2615.
46. A. Dines, M. Ellis and J. Blumberger, J. Chem. Phys., 2023, 159, 234118.
47. S. Gómez, E. Spinlove and G. Worth, Phys. Chem. Chem. Phys., 2023, 26, 1829–1844.
48. J.-K. Ha, I. S. Lee and S. K. Min, J. Phys. Chem. Lett., 2018, 9, 1097–1104.
49. J.-K. Ha and S. K. Min, J. Chem. Phys., 2022, 156, 174109.
50. W. H. Miller, J. Chem. Phys., 2012, 136, 210901.
51. G. Trenins and J. O. Richardson, J. Chem. Phys., 2022, 156, 174115.
52. R. Conte and M. Ceotto, Quantum Chemistry and Dynamics of Excited States: Methods and Applications, John Wiley & Sons, 2021.
53. A. Scheidegger, N. V. Golubev and J. J. L. Vaní'ček, Chimia, 2023, 77, 201.
54. S. Malpathak, M. S. Church and N. Ananth, J. Phys. Chem. A, 2022, 126, 6359–6375.
55. S. Kundu and N. Makri, Annu. Rev. Phys. Chem., 2022, 73, 349–375.
56. J. C. Tully, J. Chem. Phys., 1990, 93, 1061.
57. J. C. Tully, Faraday Discuss., 1998, 110, 407.
58. J. R. Mannouch and J. O. Richardson, J. Chem. Phys., 2023, 158, 104111.
59. M. Persico, G. Granucci and D. Accomasso, Comprehensive Computational Chemistry, Elsevier, Oxford, 1st edn, 2024, pp. 273–292.
60. A. Akimov, Comprehensive Computational Chemistry, Elsevier, 2022.
61. M. Barbatti, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 620–633.
62. T. R. Nelson, S. Fernandez-Alberti, A. E. Roitberg and S. Tretiak, J. Chem. Phys., 2013, 138, 224111.
63. M. Richter, P. Marquetand, J. González-Vázquez, I. Sola and L. González, J. Chem. Theory Comput., 2011, 7, 1253–1258.
64. S. Mai, P. Marquetand and L. González, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1370.
65. R. Crespo-Otero and M. Barbatti, Chem. Rev., 2018, 118, 7026–7068.
66. J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang and N. Bellonzi, Annu. Rev. Phys. Chem., 2016, 67, 387–417.
67. J. Suchan, J. Janos and P. Slavícěk, J. Chem. Theory Comput., 2020, 16, 5809–5820.
68. A. K. Belyaev, C. Lasser and G. Trigila, J. Chem. Phys., 2014, 140, 224108.
69. D. A. Fedorov and B. G. Levine, J. Phys. Chem. Lett., 2019, 10, 4542–4548.
70. L. Wang, A. Akimov and O. V. Prezhdo, J. Phys. Chem. Lett., 2016, 7, 2100–2112.
71. I. C. D. Merritt, D. Jacquemin and M. Vacher, J. Phys. Chem. Lett., 2021, 12, 8404–8415.
72. Y. Shu and D. G. Truhlar, J. Chem. Theory Comput., 2023, 19, 380–395.
73. T. R. Nelson, S. Fernandez-Alberti and S. Tretiak, Nat. Comput. Sci., 2022, 2, 689–692.
74. J. Qiu, X. Bai and L. Wang, J. Phys. Chem. Lett., 2018, 9, 4319–4325.
75. J. Xu and L. Wang, J. Phys. Chem. Lett., 2020, 11, 8283–8291.
76. C. Stopera, T. V. Grimes, P. M. McLaurin, A. Privett and J. A. Morales, Adv. Quantum Chem., 2013, 66, 113–194.
77. S. M. Parker and C. J. Schiltz, J. Chem. Phys., 2020, 153, 174109.
78. S. Ghosh, S. Giannini, K. Lively and J. Blumberger, Faraday Discuss., 2020, 221, 501–525.
79. X. Yang, M. Manathunga, S. Gozem, J. Léonard, T. Andruniów and M. Olivucci, Nat. Chem., 2022, 14, 441–449.
80. A. Kelly, J. Chem. Phys., 2019, 150, 204107.
81. J. Spencer, L. Scalfi, A. Carof and J. Blumberger, Faraday Discuss., 2016, 195, 215–236.
82. S. J. Cotton and W. H. Miller, J. Chem. Phys., 2019, 150, 104101.
83. J. F. K. Halliday and E. Artacho, SciPost Phys., 2022, 12, 020.
84. P. Chakraborty and S. Matsika, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2024, 14, e1715.
85. S. Bonella, G. Ciccotti and R. Kapral, Chem. Phys. Lett., 2010, 484, 399–404.
86. W. C. Pfalzgraff, A. Montoya-Castillo, A. Kelly and T. E. Markland, J. Chem. Phys., 2019, 150, 244109.
87. E. R. Bittner and P. J. Rossky, J. Chem. Phys., 1995, 103, 8130–8143.
88. H. Kim, A. Nassimi and R. Kapral, J. Chem. Phys., 2008, 129, 084102.
89. D. F. Coker and L. Xiao, J. Chem. Phys., 1995, 102, 496–510.
90. S. Bonella and D. F. Coker, J. Chem. Phys., 2005, 122, 194102.
91. F. Gay-Balmaz and C. Tronci, J. Phys. A: Math. Theor., 2023, 56, 144002.
92. A. J. Dunnett and A. W. Chin, Entropy, 2021, 23, 77.
93. A. Chin, E. Mangaud, V. Chevet, O. Atabek and M. Desouter-Lecomte, Chem. Phys., 2019, 525, 110392.
94. H.-D. Meyer, U. Manthe and L. S. Cederbaum, Chem. Phys. Lett., 1990, 165, 73–78.
95. A. J. M. H. Beck, G. A. Worth and H.-D. Meyer, Phys. Rep., 2000, 324, 1–105.
96. G. A. Worth, H.-D. Meyer, H. Koeppel, L. S. Cederbaum and I. Burghardt, Int. Rev. Phys. Chem., 2008, 27, 569.
97. H.-D. Meyer, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 351–374.
98. D. Lauvergnat and A. Nauts, Spectrochim. Acta, Part A, 2014, 119, 18–25.
99. S. Sasmal and O. Vendrell, J. Chem. Phys., 2022, 157, 134102.
100. A. Abedi, N. T. Maitra and E. K. U. Gross, Phys. Rev. Lett., 2010, 105, 123002.
101. F. Agostini and E. K. U. Gross, Eur. Phys. J. B, 2021, 94, 179.
102. M. Born and R. Oppenheimer, Ann. Phys., 1927, 84, 457–484.
103. M. Born and K. Huang, The Dynamical Theory of Crystal Lattices, Clarendon Press, Oxford, 1954.
104. G. A. Worth and L. S. Cederbaum, Annu. Rev. Phys. Chem., 2004, 55, 127–158.
105. W. Domcke and D. R. Yarkony, Annu. Rev. Phys. Chem., 2012, 63, 325–352.
106. S. Matsika and P. Krause, Annu. Rev. Phys. Chem., 2011, 62, 621.
107. S. Matsika and D. R. Yarkony, J. Am. Chem. Soc., 2003, 125, 10672–10676.
108. A. Abedi, N. T. Maitra and E. K. U. Gross, J. Chem. Phys., 2012, 137, 22A530.
109. A. Abedi, F. Agostini, Y. Suzuki and E. K. U. Gross, Phys. Rev. Lett., 2013, 110, 263001.
110. F. Agostini, A. Abedi, Y. Suzuki, S. K. Min, N. T. Maitra and E. K. U. Gross, J. Chem. Phys., 2015, 142, 084303.
111. F. Agostini, I. Tavernelli and G. Ciccotti, Eur. Phys. J. B, 2018, 91, 139.
112. S. K. Min, F. Agostini, I. Tavernelli and E. K. U. Gross, J. Phys. Chem. Lett., 2017, 8, 3048–3055.
113. F. Talotta, S. Morisset, N. Rougeau, D. Lauvergnat and F. Agostini, Phys. Rev. Lett., 2020, 124, 033001.
114. P. Vindel-Zandbergen, S. Matsika and N. T. Maitra, J. Phys. Chem. Lett., 2022, 13, 1785–1790.
115. G. Hunter, Int. J. Quantum Chem., 1975, 9, 237–242.
116. D. M. Bishop and L. M. Cheung, Chem. Phys. Lett., 1977, 50, 172–174.
117. G. Hunter, Int. J. Quantum Chem., 1981, 19, 755–761.
118. R. Lefebvre, J. Chem. Phys., 2015, 142, 214105.
119. N. I. Gidopoulos and E. K. U. Gross, Philos. Trans. R. Soc., A, 2014, 372, 20130059.
120. R. Requist and E. K. U. Gross, Phys. Rev. Lett., 2019, 117, 193001.
121. C. Li, R. Requist and E. K. U. Gross, J. Chem. Phys., 2018, 148, 084110.
122. L. Lacombe and N. T. Maitra, Phys. Rev. Lett., 2020, 124, 206401.
123. E. K. U. G. Ryan Requist, Phys. Rev. Lett., 2019, 127, 116401.
124. S. K. Min, A. Abedi, K. S. Kim and E. K. U. Gross, Phys. Rev. Lett., 2014, 113, 263004.
125. R. Requist, F. Tandetzky and E. K. U. Gross, Phys. Rev. A, 2016, 93, 042108.
126. R. Requist, C. R. Proetto and E. K. U. Gross, Phys. Rev. B, 2019, 99, 165136.
127. L. S. Cederbaum, J. Chem. Phys., 2013, 138, 224110.
128. L. S. Cederbaum, Chem. Phys., 2015, 457, 129.
129. Y. Suzuki, A. Abedi, N. T. Maitra, K. Yamashita and E. K. U. Gross, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 89, 040501(R).
130. A. Schild, F. Agostini and E. K. U. Gross, J. Phys. Chem. A, 2016, 120, 3316.
131. A. Scherrer, F. Agostini, D. Sebastiani, E. K. U. Gross and R. Vuilleumier, Phys. Rev. X, 2017, 7, 031035.
132. C. Li, R. Requist and E. K. U. Gross, Phys. Rev. Lett., 2022, 128, 113001.
133. R. Martinazzo and I. Burghardt, Phys. Rev. Lett., 2022, 128, 206002.
134. R. Martinazzo and I. Burghardt, arXiv, 2023, preprint, arXiv:2310.08766v1 [physics.chem-ph] DOI:10.48550/arXiv.2310.08766.
135. A. Schild and E. K. U. Gross, Phys. Rev. Lett., 2017, 118, 163202.
136. J. Kocák and A. Schild, Phys. Rev. Res., 2020, 2, 043365.
137. S. Giarrusso, P. Gori-Giorgi and F. Agostini, ChemPhysChem, 2024, 25, e202400127.
138. S. Giarrusso and F. Agostini, arXiv, 2024, preprint, arXiv:2408.06038v3 [physics.chem-ph] DOI:10.48550/arXiv.2408.06038.
139. J. Kock, E. Kraisler and A. Schild, J. Phys. Chem. Lett., 2021, 12, 3204–3209.
140. A. Abedi, E. Khosravi and I. Tokatly, Eur. Phys. J. B, 2018, 91, 194.
141. N. M. Hoffmann, H. Appel, A. Rubio and N. Maitra, Eur. Phys. J. B, 2018, 91, 180.
142. L. Lacombe, N. M. Hoffmann and N. T. Maitra, Phys. Rev. Lett., 2019, 123, 083201.
143. N. M. Hoffmann, L. Lacombe, A. Rubio and N. T. Maitra, J. Chem. Phys., 2020, 153, 104103.
144. B. Rosenzweig, N. M. Hoffmann, L. Lacombe and N. T. Maitra, J. Chem. Phys., 2022, 156, 054101.
145. E. Villaseco Arribas, N. T. Maitra and F. Agostini, J. Chem. Phys., 2024, 160, 054102.
146. B. Mignolet and B. F. E. Curchod, J. Phys. Chem. A, 2019, 123, 3582.
147. L. M. Ibele, B. F. E. Curchod and F. Agostini, J. Phys. Chem. A, 2022, 126, 1263–1281.
148. T. Fiedlschuster, J. Handt, E. K. U. Gross and R. Schmidt, Phys. Rev. A, 2017, 95, 063424.
149. T. Fiedlschuster, J. Handt and R. Schmidt, Phys. Rev. A, 2016, 93, 053409.
150. M. Schirò, F. G. Eich and F. Agostini, J. Chem. Phys., 2021, 154, 114101.
151. G. Jolicard and O. Atabek, Phys. Rev. A: At., Mol., Opt. Phys., 1992, 46, 5845–5855.
152. X. He, O. Atabek and A. Giusti-Suzor, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 5586–5594.
153. A. Giusti-Suzor and F. H. Mies, Phys. Rev. Lett., 1992, 68, 38638729.
154. L. J. Frasinski, J. Plumridge, J. H. Posthumus, K. Codling, P. F. Taday, E. J. Divall and A. J. Langley, Phys. Rev. Lett., 2001, 86, 2541–2544.
155. P. H. Bucksbaum, A. Zavriyev, H. G. Muller and D. W. Schumacher, Phys. Rev. Lett., 1990, 64, 1883–1886.
156. A. Giusti-Suzo, X. He, O. Atabek and F. H. Mies, Phys. Rev. Lett., 1990, 64, 515–518.
157. F. Talotta, S. Morisset, N. Rougeau, D. Lauvergnat and F. Agostini, J. Phys. Chem. A, 2021, 125, 6075.
158. B. F. E. Curchod, C. Rauer, P. Marquetand, L. González and T. J. Martínez, J. Chem. Phys., 2016, 144, 101102.
159. L. Favero, G. Granucci and M. Persico, Phys. Chem. Chem. Phys., 2013, 15, 20651.
160. J. S. Garcáa, F. Talotta, F. Alary, I. M. Dixon, J.-L. Heully and M. Boggio-Pasqua, Molecules, 2017, 22, 1667.
161. A. J. Atkins, F. Talotta, L. Freitag, M. Boggio-Pasqua and L. González, J. Chem. Theory Comput., 2017, 13, 4123–4145.
162. F. Talotta, M. Boggio-Pasqua and L. González, Chem. – Eur. J., 2020, 26, 11522–11528.
163. F. Talotta, J.-L. Heully, F. Alary, I. M. Dixon, L. González and M. Boggio-Pasqua, J. Chem. Theory Comput., 2017, 13, 6120–6130.
164. I. Tavernelli, B. F. E. Curchod and U. Rothlisberger, Chem. Phys., 2011, 391, 101–109.
165. M. Richter, P. Marquetand, J. González-Vázquez, I. Sola and L. González, J. Phys. Chem. Lett., 2012, 3, 3090–3095.
166. G. Capano, M. Chergui, U. Rothlisberger, I. Tavernelli and T. J. Penfold, J. Phys. Chem. A, 2014, 118, 9861–9869.
167. G. Cui and W. Thiel, J. Chem. Phys., 2014, 141, 124101.
168. G. H. Gossel, L. Lacombe and N. T. Maitra, J. Chem. Phys., 2019, 150, 154112.
169. F. Talotta, F. Agostini and G. Ciccotti, J. Phys. Chem. A, 2020, 124, 6764–6777.
170. F. Agostini, A. Abedi, Y. Suzuki and E. K. U. Gross, Mol. Phys., 2013, 111, 3625–3640.
171. F. Agostini, A. Abedi and E. K. U. Gross, J. Chem. Phys., 2014, 141, 214101.
172. A. Abedi, F. Agostini and E. K. U. Gross, Europhys. Lett., 2014, 106, 33001.
173. F. Agostini, S. K. Min, A. Abedi and E. K. U. Gross, J. Chem. Theory Comput., 2016, 12, 2127–2143.
174. F. Agostini, Eur. Phys. J. B, 2018, 91, 143.
175. L. M. Ibele, C. Pieroni, F. Talotta, B. F. Curchod, D. Lauvergnat and F. Agostini, Comprehensive Computational Chemistry, Elsevier, Oxford, 1st edn, 2024, pp. 188–211.
176. Y. Suzuki and K. Watanabe, Phys. Rev. A, 2016, 94, 032517.
177. E. V. Arribas, F. Agostini and N. T. Maitra, Molecules, 2022, 27, 4002.
178. B. F. E. Curchod and F. Agostini, J. Phys. Chem. Lett., 2017, 8, 831–837.
179. F. Agostini and B. F. E. Curchod, Eur. Phys. J. B, 2018, 91, 141.
180. C. Pieroni, E. Marsili, D. Lauvergnat and F. Agostini, J. Chem. Phys., 2021, 154, 034104.
181. B. F. E. Curchod, F. Agostini and E. K. U. Gross, J. Chem. Phys., 2016, 145, 034103.
182. E. Khosravi, A. Abedi and N. T. Maitra, Phys. Rev. Lett., 2015, 115, 263002.
183. Y. Suzuki, A. Abedi, N. T. Maitra and E. K. U. Gross, Phys. Chem. Chem. Phys., 2015, 17, 29271–29280.
184. E. Khosravi, A. Abedi, A. Rubio and N. T. Maitra, Phys. Chem. Chem. Phys., 2017, 19, 8269–8281.
185. T. Fiedlschuster, J. Handt, E. K. U. Gross and R. Schmidt, Phys. Rev. A, 2017, 95, 063424.
186. F. Talotta, S. Morisset, N. Rougeau, D. Lauvergnat and F. Agostini, J. Chem. Theory Comput., 2020, 16, 4833–4848.
187. F. G. Eich and F. Agostini, J. Chem. Phys., 2016, 145, 054110.
188. F. Agostini, S. K. Min and E. K. U. Gross, Ann. Phys., 2015, 527, 546–555.
189. D. Han, J.-K. Ha and S. K. Min, J. Chem. Theory Comput., 2023, 19, 2186–2197.
190. E. Villaseco Arribas and N. T. Maitra, J. Chem. Phys., 2023, 158, 161105.
191. F. Agostini, B. F. E. Curchod, R. Vuilleumier, I. Tavernelli and E. K. U. Gross, Handbook of Materials Modeling, Springer, Netherlands, 2018, pp. 1–47.
192. B. F. E. Curchod, F. Agostini and I. Tavernelli, Eur. Phys. J. B, 2018, 91, 168.
193. I. S. Lee, J. Ha, D. Han, T. I. Kim, S. W. Moon and S. K. Min, J. Comput. Chem., 2021, 42, 1755–1766.
194. D. Han and A. V. Akimov, J. Chem. Theory Comput., 2024, 20, 5022–5042.
195. G. Granucci and M. Persico, J. Chem. Phys., 2007, 126, 134114.
196. N. Shenvi, J. E. Subotnik and W. Yang, J. Chem. Phys., 2011, 134, 144102.
197. J. E. Subotnik and N. Shenvi, J. Chem. Phys., 2011, 134, 024105.
198. H. M. Jaeger, S. Fischer and O. V. Prezhdo, J. Chem. Phys., 2012, 137, 22A545.
199. G. Granucci and M. Persico, J. Chem. Phys., 2007, 126, 134114.
200. G. Granucci, M. Persico and A. Zoccante, J. Chem. Phys., 2010, 133, 134111.
201. C. C. Martens, J. Phys. Chem. A, 2019, 123, 1110–1128.
202. C. Zhu, S. Nangia, A. W. Jasper and D. G. Truhlar, J. Chem. Phys., 2004, 121, 7658.
203. M. P. Esch and B. G. Levine, J. Chem. Phys., 2020, 153, 114104.
204. E. Villaseco Arribas, L. M. Ibele, D. Lauvergnat, N. T. Maitra and F. Agostini, J. Chem. Theory Comput., 2023, 19, 7787–7800.
205. A. V. Akimov, J. Comput. Chem., 2016, 37, 1626–1649.
206. F. Agostini, E. Marsili, F. Talotta, C. Pieroni, E. Villaseco Arribas, L. M. Ibele and E. Sangiogo Gil, G-CTMQC, 2024, accessed in June, https://gitlab.com/agostini.work/g-ctmqc.
207. D. Lauvergnat, QuantumModelLib, 2024, accessed in June, https://github.com/lauvergn/QuantumModelLib.
208. G. Granucci, M. Persico, D. Accomasso, E. Sangiogo Gil, S. Corni, J. Fregoni, T. Laino, M. Tesi and A. Toniolo, MOPAC-PI: a program for excited state dynamics simulations based on nonadiabatic trajectories and semiempirical electronic structure calculations, 2024, https://gitlab.com/granucci/mopacpi.git.
209. C. Pieroni, E. Sangiogo Gil, L. M. Ibele, M. Persico, G. Granucci and F. Agostini, J. Chem. Theory Comput., 2024, 20, 580–596.
210. G. A. Worth, K. Giri, G. W. Richings, I. Burghardt, M. H. Beck, A. Jäckle and H.-D. Meyer, The QUANTICS Package, Version 2.0, University College London, London, UK, 2020.
211. S. Mai, D. Avagliano, M. Heindl, P. Marquetand, M. F. S. J. Menger, M. Oppel, F. Plasser, S. Polonius, M. Ruckenbauer, Y. Shu, D. G. Truhlar, L. Zhang, P. Zobel and L. González, SHARC3.0: Surface Hopping Including Arbitrary Couplings – Program Package for Non-Adiabatic Dynamics, 2023https://sharc-md.org/.
212. M. Barbatti, M. Bondanza, R. Crespo-Otero, B. Demoulin, P. O. Dral, G. Granucci, F. Kossoski, H. Lischka, B. Mennucci, S. Mukherjee, M. Pederzoli, M. Persico, M. Pinheiro Jr, J. Pittner, F. Plasser, E. S. Gil and L. Stojanovic, J. Chem. Theory Comput., 2022, 18, 6851.
213. L. M. Ibele and B. F. E. Curchod, Phys. Chem. Chem. Phys., 2020, 22, 15062.
214. G. Granucci and M. Persico, Chem. Phys. Lett., 1995, 246, 228–234.
215. D. L. Audrey Pollien, E. Villaseco Arribas and F. Agostini, Mol. Phys., 2024, e2378960.
216. B. F. E. Curchod, W. J. Glover and T. J. Martínez, J. Phys. Chem. A, 2020, 124, 6133–6143.
217. L. M. Ibele, Y. Lassmann, T. J. Martínez and B. F. E. Curchod, J. Chem. Phys., 2021, 154, 104110.
218. Y. Lassmann and B. F. E. Curchod, J. Chem. Phys., 2021, 154, 211106.
219. J. Janoš and P. Slavíček, J. Chem. Theory Comput., 2023, 19, 8273–8284.
220. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B 864.
221. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A 1133.
222. S. Kube, C. Lasser and M. Weber, J. Comput. Phys., 2009, 228, 1947–1962.
223. M. Barbatti and K. Sen, Int. J. Quantum Chem., 2015, 116, 762–771.
224. T. Plé, S. Huppert, F. Finocchi, P. Depondt and S. Bonella, J. Chem. Phys., 2019, 151, 114114.
225. F. Kröninger, C. Lasser and J. J. L. Vaní’ček, Front. Phys., 2023, 11, 1106324.
226. J. Suchan, D. Hollas, B. F. E. Curchod and P. Slavíček, Faraday Discuss., 2018, 212, 307.
227. C. Pieroni, F. Becuzzi, L. Creatini, G. Granucci and M. Persico, J. Chem. Theory Comput., 2023, 19, 2430–2445.
228. D. Avagliano, E. Lorini and L. González, Philos. Trans. R. Soc., A, 2022, 380, 20200381.
229. H. C. Longuet-Higgins, U. Öpik, M. H. L. Pryce and R. Sack, Proc. R. Soc. A, 1958, 244, 1–16.
230. M. V. Berry, Proc. R. Soc. London, Ser. A, 1984, 392, 45–57.
231. S. C. Althorpe, F. Fernández-Alonso, B. D. Bean, J. D. Ayers, A. E. Pomerantz, R. N. Zare and E. Wrede, Nature, 2002, 416, 67–70.
232. J. C. Juanes-Marcos, S. C. Althorpe and E. Wrede, Science, 2005, 309, 1227–1230.
233. B. K. Kendrick, J. Hazra and N. Balakrishnan, Phys. Rev. Lett., 2015, 115, 153201.
234. C. Xie, C. L. Malbon, D. R. Yarkony, D. Xie and H. Guo, J. Am. Chem. Soc., 2018, 140, 1986–1989.
235. D. Yuan, Y. Guan, W. Chen, H. Zhao, S. Yu, C. Luo, Y. Tan, T. Xie, X. Wang, Z. Sun, D. H. Zhang and X. Yang, Science, 2018, 362, 1289–1293.
236. C. Xie, C. L. Malbon, H. Guo and D. R. Yarkony, Acc. Chem. Res., 2019, 52, 501–509.
237. E. Cohen, H. Larocque, F. Bouchard, F. Nejadsattari, Y. Gefen and E. Karimi, Nat. Rev. Phys., 2019, 1, 437–449.
238. D. Yuan, Y. Huang, W. Chen, H. Zhao, S. Yu, C. Luo, Y. Tan, S. Wang, X. Wang, Z. Sun and X. Yang, Nat. Commun., 2020, 11, 3640.
239. C. H. Valahu, V. C. Olaya-Agudelo, R. J. MacDonell, T. Navickas, A. D. Rao, M. J. Millican, J. B. Pérez-Sánchez, J. Yuen-Zhou, M. J. Biercuk, C. Hempel, T. R. Tan and I. Kassal, Nat. Chem., 2023, 15, 1503–1508.
240. S. Mai, P. Marquetand and L. González, Int. J. Quantum Chem., 2015, 115, 1215–1231.
241. B. Mignolet, B. F. E. Curchod and T. J. Martínez, J. Chem. Phys., 2016, 145, 191104.
242. T. J. Martinez, Nature, 2010, 467, 412–413.
243. W. Domcke and D. Yarkony, Annu. Rev. Phys. Chem., 2012, 63, 325.
244. D. R. Yarkony, Rev. Mod. Phys., 1996, 68, 985.
245. M. Baer, Beyond Born-Oppenheimer: Electronic nonadiabatic coupling terms and conical intersections, John Wiley & Sons, Inc., 2006.
246. J. Hazra, B. K. Kendrick and N. Balakrishnan, J. Phys. B: At. Mol. Phys., 2016, 49, 194004.
247. D. Yuan, S. Yu, W. Chen, J. Sang, C. Luo, T. Wang, X. Xu, P. Casavecchia, X. Wang, Z. Sun, D. H. Zhang and X. Yang, Nat. Chem., 2018, 10, 653–658.
248. C. A. Mead, Chem. Phys., 1980, 49, 23–32.
249. C. A. Mead and D. G. Truhlar, J. Chem. Phys., 1979, 70, 2284.
250. L. M. Ibele, E. Sangiogo Gil, B. F. E. Curchod and F. Agostini, J. Phys. Chem. Lett., 2023, 14, 11625–11631.
251. L. M. Ibele and F. Agostini, J. Phys. Chem. A, 2024, 128, 3672–3684.
252. M. Baer, S. H. Lin, A. Alijah, S. Adhikari and G. D. Billing, Phys. Rev. A: At., Mol., Opt. Phys., 2000, 62, 032506.
253. S. Adhikari, G. D. Billing, A. Alijah, S. H. Lin and M. Baer, Phys. Rev. A: At., Mol., Opt. Phys., 2000, 62, 032507.
254. B. Sarkar and S. Adhikari, J. Chem. Phys., 2006, 124, 074101.
255. B. Sarkar and S. Adhikari, Int. J. Quantum Chem., 2009, 109, 650–667.
256. A. K. Paul, S. Sardar, B. Sarkar and S. Adhikari, J. Chem. Phys., 2009, 131, 124312.
257. R. Requist, C. R. Proetto and E. K. U. Gross, Phys. Rev. A, 2017, 96, 062503.
258. P. B. Allen, A. G. Abanov and R. Requist, Phys. Rev. A: At., Mol., Opt. Phys., 2005, 71, 043203.
259. J. Ronayette, R. Arnaud, P. Lebourgeois and J. Lemaire, Can. J. Chem., 1974, 52, 1848–1857.
260. G. Gauglitz and S. Hubig, J. Photochem., 1985, 30, 121–125.
261. P. Bortolus and S. Monti, J. Phys. Chem., 1979, 83, 648–652.
262. N. Siampiringue, G. Guyot, S. Monti and P. Bortolus, J. Photochem., 1987, 37, 185–188.
263. A. W. Adamson, A. Vogler, H. Kunkely and R. Wachter, J. Am. Chem. Soc., 1978, 100, 1298–1300.
264. H. Rau and E. Lueddecke, J. Am. Chem. Soc., 1982, 104, 1616–1620.
265. H. Rau, J. Photochem., 1984, 26, 221–225.
266. V. Ladányi, P. Dvořák, J. Al Anshori, L. Vetráková, J. Wirz and D. Heger, Photochem. Photobiol. Sci., 2017, 16, 1757–1761.
267. C. Ciminelli, G. Granucci and M. Persico, Chem. – Eur. J., 2004, 10, 2327–2341.
268. A. Nenov, R. Borrego-Varillas, A. Oriana, L. Ganzer, F. Segatta, I. Conti, J. Segarra-Marti, J. Omachi, M. Dapor, S. Taioli, C. Manzoni, S. Mukamel, G. Cerullo and M. Garavelli, J. Chem. Phys. Lett., 2018, 9, 1534–1541.
269. J. K. Yu, C. Bannwarth, R. Liang, E. G. Hohenstein and T. J. Martínez, J. Am. Chem. Soc., 2020, 142, 20680–20690.
270. M. Filatov, S. K. Min and K. S. Kim, J. Chem. Theory Comput., 2018, 14, 4499–4512.
271. M. Filatov, S. K. Min and K. S. Kim, Mol. Phys., 2019, 117, 1128–1141.
272. M. Filatov, M. Paolino, S. K. Min and K. S. Kim, J. Phys. Chem. Lett., 2018, 9, 4995–5001.
273. M. Filatov, M. Paolino, S. K. Min and C. H. Choi, Chem. Commun., 2019, 55, 5247.
274. M. Filatov(Gulak), M. Paolino, D. Kaliakin, M. Olivucci, E. Kraka and S. K. Min, Commun. Phys., 2024, 7, 219.
275. M. Filatov, S. K. Min and C. H. Choi, Phys. Chem. Chem. Phys., 2019, 21, 2489–2498.
276. T. I. Kim, J.-K. Ha and S. K. Min, Top. Curr. Chem., 2022, 380, 1–27.
277. P. Vindel-Zandbergen, L. M. Ibele, J.-K. Ha, S. K. Min, B. F. E. Curchod and N. T. Maitra, J. Chem. Theory Comput., 2021, 17, 3852–3862.
278. E. V. Arribas, P. Vindel-Zandbergen, S. Roya and N. T. Maitra, Phys. Chem. Chem. Phys., 2023, 25, 26380.
279. F. Krausz and M. Ivanov, Rev. Mod. Phys., 2009, 81, 163–234.
280. F. Lépine, M. Y. Ivanov and M. J. J. Vrakking, Nat. Photonics, 2014, 8, 195–204.
281. M. Nisoli, P. Decleva, F. Calegari, A. Palacios and F. Martn, Chem. Rev., 2017, 117, 10760–10825.
282. A. Palacios and F. Martn, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2020, 10, e1430.
283. N. T. Maitra and E. Villaseco Arribas, arXiv, 2024, preprint, arXiv:2405.00649 [physics.chem-ph] DOI:10.48550/arXiv.2405.00649.
284. J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux and T. W. Ebbesen, Angew. Chem., Int. Ed., 2012, 51, 1592.
285. S. Wang, T. Chervy, J. George, J. A. Hutchison, C. Genet and T. W. Ebbesen, J. Phys. Chem. Lett., 2014, 5, 1433.
286. T. W. Ebbesen, Acc. Chem. Res., 2016, 49, 2403.
287. L. A. Martínez-Martínez, M. Du, R. F. Ribeiro, S. Kéa-Cohen and J. Yuen-Zhou, J. Phys. Chem. Lett., 2018, 9, 1951.
288. M. Du, L. A. Martínez-Martínez, R. F. Ribeiro, Z. Hu, V. M. Menon and J. Yuen-Zhou, Chem. Sci., 2018, 9, 6659.
289. L. A. Martínez-Martínez, R. F. Ribeiro, J. Campos-González-Angulo and J. Yuen-Zhou, ACS Photonics, 2018, 5, 167.
290. R. T. Grant, P. Michetti, A. J. Musser, P. Gregoire, T. Virgili, E. Vella, M. Cavazzini, K. Georgiou, F. Galeotti, C. Clark, J. Clark, C. Silva and D. G. Lidzey, Adv. Opt. Mater., 2016, 4, 1615.
291. V. C. Nikolis, A. Mischok, B. Siegmund, J. Kublitski, X. Jia, J. Benduhn, U. Hörmann, D. Neher, M. C. Gather, D. Spoltore and K. Vandewal, Nat. Commun., 2019, 10, 3706.
292. P. Vasa and C. Lienau, ACS Photonics, 2018, 5, 2.
293. D. G. Baranov, M. Wersäll, J. Cuadra, T. J. Antosiewicz and T. Shegai, ACS Photonics, 2018, 5, 24.
294. J. T. Hugall, A. Singh and N. F. van Hulst, ACS Photonics, 2018, 5, 43.
295. E. Orgiu, J. George, J. A. Hutchison, E. Devaux, J. F. Dayen, B. Doudin, F. Stellacci, C. Genet, J. Schachenmayer, C. Genes, G. Pupillo, P. Samori and T. W. Ebbesen, Nat. Mater., 2015, 14, 1123.
296. J. Fregoni, F. J. Garcia-Vidal and J. Feist, ACS Photonics, 2022, 9, 1096–1107.
297. J. Galego, F. J. Garcia-Vidal and J. Feist, Phys. Rev. X, 2015, 5, 041022.
298. J. Galego and F. J. G.-V. J. Feist, Phys. Rev. Lett., 2017, 119, 136001.
299. T. E. Li, A. Nitzan and J. E. Subotnik, J. Chem. Phys., 2021, 154, 094124.
300. T. E. Li, A. Nitzan and J. E. Subotnik, J. Chem. Phys., 2022, 156, 134106.
301. M. R. Fiechter, J. E. Runeson, J. E. Lawrence and J. O. Richardson, J. Phys. Chem. Lett., 2023, 14, 8261–8267.
302. E. Hulkko, S. Pikker, V. Tiainen, R. H. Tichauer, G. Groenhof and J. J. Toppari, J. Chem. Phys., 2021, 154, 154303.
303. R. H. Tichauer, J. Feist and G. Groenhof, J. Chem. Phys., 2021, 154, 104112.
304. R. H. Tichauer, D. Morozov, I. Sokolovskii, J. J. Toppari and G. Groenhof, J. Phys. Chem. Lett., 2022, 13, 6259–6267.
305. J. Flick, H. Appel, M. Ruggenthaler and A. Rubio, J. Chem. Theory Comput., 2017, 13, 1616–1625.
306. M. R. Fiechter and J. O. Richardson, J. Chem. Phys., 2024, 160, 184107.
307. E. W. Fischer and P. Saalfrank, J. Chem. Theory Comput., 2023, 19, 72157229.
308. J. Fregoni, G. Granucci, M. Persico and S. Corni, Chem, 2020, 6, 250–265.
309. O. Vendrell, Phys. Rev. Lett., 2018, 121, 253001.
310. M. Kowalewski, K. Bennett and S. Mukamel, J. Phys. Chem. Lett., 2016, 7, 2050.
311. P. Martinez, B. Rosenzweig, N. M. Hoffmann, L. Lacombe and N. T. Maitra, J. Chem. Phys., 2021, 154, 014102.
312. B. Rana, E. G. Hohenstein and T. J. Martínez, J. Phys. Chem. A, 2023, 128, 139–151.
313. A. Thomas, J. George, A. Shalabney, M. Dryzhakov, S. J. Varma, J. Moran, T. Chervy, X. Zhong, E. Devaux, C. Genet, J. A. Hutchison and T. W. Ebbesen, Angew. Chem., Int. Ed., 2016, 55, 11462–11466.
314. T. E. Lia, J. E. Subotnik and A. Nitzan, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 18324–18331.
315. A. Z. Lieberherr, S. T. E. Furniss, J. E. Lawrence and D. E. Manolopoulos, J. Chem. Phys., 2023, 158, 234106.
316. E. Sangiogo Gil, D. Lauvergnat and F. Agostini, J. Chem. Phys., 2024, 161, 084112.

---

Footnotes

† According to Stokes' theorem, the circulation of the TDVP along Γ is equal to the flux of the curl of the TDVP across a surface whose boundary is Γ. Therefore, this proves again that γ(Γ, t) is a physical gauge-independent observable.

‡ Note that the photosisomerization is termed “partial” because the photoisomerization quantum yield is calculated only for the molecules that reached the ground state by the end of the trajectory.

---