Electronic decoherence along a single nuclear trajectory

PHYSICAL REVIEW RESEARCH 7, 043075 (2025) Electronic decoherence along a single nuclear trajectory Matisse Wei-Yuan Tu 1,2,*and E. K. U. Gross 1 1Fritz Haber Center for Molecular Dynamics, Institute of Chemistry, The Hebrew University of Jerusalem, Jerusalem 91904, Israel 2Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany (Received 7 February 2025; accepted 17 September 2025; published 17 October 2025) We describe an approach to subsystem decoherence without the usual tracing out of the environment. The subsystem of focus is described entirely by a pure state evolving nonunitarily along a single classical trajectory of its environment. The approach is deduced from the exact factorization framework for arbitrary systems of electrons and nuclei. The nonunitarity of the electronic dynamics arises exclusively from nonadiabatic correlations between electrons and nuclei. We demonstrate that the approach correctly describes the coherence gain and the subsequent decoherence for the example of a nuclear trajectory passing through an avoided crossing, the prototypical case where single-trajectory Ehrenfest dynamics fails to produce decoherence. DOI: 10.1103/nnkr-phm5 I. INTRODUCTION Decoherence is a multifaceted phenomenon that appears in many branches of physics, chemistry, and biology whenever a system consists of interacting distinct subsystems. It represents a major obstacle in the technical realization of a scalable quantum computer. The theoretical description starts from designating one subsystem as the primary focus, called the principal system, while treating the rest as its “environment.” The wave function of the whole composite system satisfies one single time-dependent Schrödinger equation (TDSE). Decoherence results from the entanglement between the principal system and its environment. In the standard treatment [1,2], decoherence effects are obtained by tracing out the environment, which causes the principal system to be in a mixed state, i.e., an ensemble average of pure states, described by a reduced density operator (RDO). In this work, starting from the full TDSE of the complete system, we deduce an alternative way of tackling subsystem decoherence by propagating just one pure state without ensemble averages. This pure state evolves nonunitarily in a norm-conserving fashion along one classical trajectory of its environment. The approach is derived here for arbitrary systems of interacting electrons and nuclei, where the electronic degrees of freedom (DOF) constitute the principal system interacting with its nuclear environment. While this is a specific choice of twocomponent system, we emphasize that the approach presented here is general and can be applied to other composite systems in the same way [3]. A popular family of approaches dealing with the coupled system of electrons and nuclei are the so-called mixed quantum-classical (MQC) approaches where the nuclear mo- *Contact author: [email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. tion is described in terms of classical trajectories while the electronic DOF are treated fully quantum mechanically. It is usually assumed that electronic decoherence cannot be described purely within a single trajectory, without appealing to the finite spatial extension of nuclear wave packets. This holds true for Ehrenfest dynamics, where the electrons evolve unitarily along a single classical nuclear trajectory and decoherence arises only after averaging over multiple trajectories [4,5]. More broadly, many MQC methods aim to incorporate decoherence readily on the single-trajectory level by introducing modifications, often termed decoherence corrections, to the electronic evolution at each time step [6–17]. These corrections produce deviations from the purely unitary Schrödinger dynamics by adjusting the electronic-state coefficients according to a decoherence time parameter, which is typically determined heuristically. A common rationale is the spatial separation of nuclear wave packets on different electronic surfaces: As these wave packets lose overlap, coherence between the corresponding electronic states should diminish. Even energy-based decoherence-time models, which do not explicitly track the wave packet separation, are based on this intuition [9,17]. This prevailing strategy naturally leads to a curiosity: Can electronic decoherence emerge natively from a single classical trajectory, without appealing to the finite spatial extension of nuclear wave packets? The exact factorization (EF) formalism offers a compelling framework in this regard [18–26]. Briefly, the EF framework, without approximation, represents the fully correlated wave function, (r,R,t), of interacting electrons and nuclei as a single product (r,R,t)=χ(R,t)φ(r|t,R),(1) where χ(R,t) is called the nuclear wave function and φ(r|t,R) an electronic wave function, which parametrically depends on (t,R). Here, r= {rj}and R= {Rν}denote the electronic and nuclear positions, respectively. The physical meaning of the two factors becomes immediately clear by remembering that the absolute square |(r,R,t)|2is the joint probability density of finding, at time t, the electrons at rand 2643-1564/2025/7(4)/043075(8) 043075-1 Published by the American Physical Society MATISSE WEI-YUAN TU AND E. K. U. GROSS PHYSICAL REVIEW RESEARCH 7, 043075 (2025) the nuclei at R. Being a joint probability, |(r,R,t)|2can be written as the product of a marginal and a conditional probability, implying that |χ(R,t)|2is the marginal probability density of finding the nuclei at positions R, while |φ(r|t,R)|2 is the conditional probability density of finding the electrons at r, given the nuclei are at Rat time t. The nuclear factor χ(R,t), on one hand, satisfies a standard TDSE with a time-dependent scalar potential ε(R,t) and a time-dependent vector potential A(R,t). The latter has the algebraic structure of a Berry connection but, being formally exact, it is not tied to an adiabatic approximation [24–26]. On the other hand, the electronic factor φ(r|t,R) satisfies a TDSE-like equation of motion (EOM) [18,23], with a non-Hermitian Hamiltonian, i.e., its time evolution is nonunitary but still conserves the norm of φ(r|t,R), namely, dr|φ(r|t,R)|2is normalized to 1 for each fixed (t,R). In fact, in view of the mathematical method of characteristics [27,28], a set of (arbitrarily dense) classical trajectories still gives a formally exact representation of the nuclear wave function χ(R,t). These trajectories satisfy standard classical EOMs [29], and the corresponding exact forces have been identified, both in the quantum case [21] and in the classical limit [30]. Since the exact factorization of Eq. (1) is a formally exact representation of the full electron-nuclear wave function, it should certainly describe electronic decoherence. Indeed, by treating the nuclear dynamics with a swarm of trajectories, and combining it with the electronic EOM of the EF (with further-going approximations for the parametric dependence of the electronic factor on nuclear coordinates), the resulting coupled-trajectory mixed quantumclassical (CTMQC) algorithm was demonstrated to capture electronic decoherence, essentially through the nuclear quantum momentum [20,22,31]. On one hand, the nuclear quantum momentum explicitly derives from the nuclear wave function, it is not intrinsic to a single classical trajectory, which underlies the rise of coupling of trajectories. On the other hand, a numerical difficulty associated with this algorithm is that the coupled nature of the trajectory propagation makes it harder to converge with the number of trajectories. While several simplifications have led to independent-trajectory variants [32–34], these developments are primarily motivated by computational efficiency rather than by the conceptual issue we address here, namely, the nature of decoherence when nuclear dynamics is represented by a single trajectory. This raises two key open questions. First, although the nonunitary yet norm-preserving evolution of the electronic conditional state is an explicit feature of the EF formalism in its fully quantum formulation [18–26], it remains unclear whether this nonunitarity persists after taking the classical limit of the nuclei (cf. [30]). To our knowledge, this has not yet been explicitly demonstrated. Second, and more fundamentally, while nonunitarity is commonly associated with decoherence in the context of statistical mixtures and quantum entanglement with an environment [1,2,35], its implications for pure state, trajectory-based descriptions remain largely unexplored. The aim of this work is to clarify how electronic decoherence can arise within the EF formalism when the nuclear motion is approximated by a single classical trajectory. By isolating and analysing the internal structure of this nonunitary evolution, we seek to identify trajectory-native factors of decoherence. In Sec. II, we review the core EF formalism, with particular focus on the dynamics of the electronic factor under the classical limit of the nuclei. In Sec. III A, we explicitly demonstrate that the nonadiabatic correlations between electrons and nuclei continue to dominate the nonunitarity of the electronic subsystem’s pure-state dynamics, even in the classical limit for the nuclei. In Sec. III B, we analyze a prototypical scenario in which a nuclear trajectory passes through an avoided crossing, comparing three approaches: (1) the full exact quantum solution, (2) mean-field Ehrenfest dynamics, and (3) the single-trajectory representation of EF. These comparisons reveal that the deorthogonalization processes underlying the nonunitarity of the pure electronic state dynamics can indeed account for decoherence at the level of a single classical trajectory. Broader implications and concluding remarks are provided in Sec. IV. II. THE SINGLE-TRAJECTORY LIMIT OF THE EF We start from the Hamiltonian of a general system of interacting electrons and nuclei, given by H=Tn+HBO(R),(2) where Tnis the nuclear kinetic energy operator and HBO(R) is the so-called Born-Oppenheimer (BO) Hamiltonian, HBO(R)=Te+Vee +Ven (R)+Vnn(R),(3) which consists of the electronic kinetic energy Teand the bare Coulomb interactions between electrons and nuclei. HBO(R) acts on the electronic Hilbert space alone, and it depends parametrically on the nuclear positions R. The eigenvalues and eigenfunctions of the BO Hamiltonian are obtained by solving HBO(R)|ϕk(R)=εk(R)|ϕk(R),(4) at each fixed nuclear configuration R. The full quantum state |(t)of the interacting system of electrons and nuclei is subject to the TDSE i|˙ (t)=H|(t). Atomic (Hartree) units are used throughout. In terms of the solutions of Eq. (4), the dynamical processes governed by Hcan be divided into adiabatic processes, where the expansion of |(t)in the complete system of BO states is dominated by a single BO state and nonadiabatic processes, where transitions between BO states take place. We characterize electronic decoherence in terms of the electronic reduced density operator (eRDO) ρe(t). The latter is obtained by tracing out the nuclear DOF from the pure-state density operator |(t)(t)|of the complete electron-nuclear system. Next, we insert the exact factorization of Eq. (1), and we write the electronic factor as r-representation of an abstract (t,R)-dependent state of the electronic Hilbert space, φ(r|t,R)=r|φ(t,R). This leads to ρe(t)=∫dR|χ(R,t)|2Pe(t,R).(5) Here, Pe(t,R)= |φ(t,R)φ(t,R)|is the R-resolved conditional electronic RDO corresponding to the pure state |φ(t,R), while the eRDO ρe(t) in general describes a mixed state. Next, we employ the classical limit [30,36,37] for the nuclear wave function leading to |χ(R,t)|2≈δ(R−Rc t), 043075-2 ELECTRONIC DECOHERENCE ALONG A SINGLE … PHYSICAL REVIEW RESEARCH 7, 043075 (2025) where Rc trepresents the classical nuclear configuration at time t. As a consequence of this approximation, the eRDO [Eq. (5)] reduces to a pure-state RDO, namely, ρe(t)=Pe(t,Rc t), evaluated along the classical trajectory Rc t. Importantly, in contrast to Ehrenfest dynamics, a mean-field-like factorization, which removes electron-nuclear entanglement, is never made. The electronic conditional amplitude of the exact factorization retains electron-nuclear entanglement through its parametric dependence on the nuclear coordinates. Explicitly, the EOM followed by |φ(t,R)reads [18–26] i∂t|φ(t,R)=[HBO(R)+μUen(R,t)]|φ(t,R),(6a) where Uen(R,t)=Uen K[φ]+Uen Q[φ,χ]− φ(t,R)|Uen K[φ] |φ(t,R)is called electron-nuclear correlation operator with Uen K[φ]= ν [−i∇ν−Aν[φ](R,t)]2 2Mν (6b) and Uen Q[φ,χ]= ν pQ ν[φ,χ](R,t)·[−i∇ν−Aν[φ](R,t)] Mν , (6c) in which pQ ν[φ,χ](R,t)=[−i∇νχ(R,t)/χ(R,t)+Aν[φ](R,t)] (6d) is the so-called nuclear momentum function with the Berryconnection-like vector potential Aν[φ](R,t)=φ(t,R)|− i∇νφ(t,R). The classical trajectory Rc tof the nuclei enters the dynamics of the pure-state electronic factor in as much as the nuclear momentum function [Eq. (6d)] is taken to the classical limit pQ ν[φ,χ](R,t)→Pc t≡Mν˙ Rc t,where the classical forces acting on the nuclei are given by Mν¨ Rc ναt=∂Aνα (Rc t,t) ∂t−∂εRc t,t ∂Rc ναt + νβ ˙ Rc νβt(∂AναRc t,t ∂Rc νβt −∂AνβRc t,t ∂Rc ναt),(7) the electriclike [the first two terms in Eq. (7)] and magneticlike [the last term in the summation of Eq. (7)] contributions [21,30]. Here, Rc t= {Rc ναt}with νand αenumerating the nucleus in the system and their spatial components, and ε(R,t)=φ(t,R)|[HBO(R)−i∂t]|φ(t,R)is the scalar potential. A generic difficulty in the described approach is that the (nonunitary) propagation of the electronic conditional wave function requires knowledge of the firstand secondorder derivatives of φ(r|t,R) with respect to the nuclear coordinates. In fact, it is exactly these derivatives that lead to the coupling of trajectories in the CTMQC algorithms. If φ(r|t,R) was known exclusively along a single classical trajectory, these derivatives could not be evaluated. To tackle this problem, we make use of the fact that the correlated electron-nuclear problem has a small parameter, the electronic-over-nuclear mass ratio μ=me M,(8) where meis the electron mass and Mis a reference nuclear mass that can be taken to be the proton mass mpor the average nuclear mass of the particular system at hand. Mνis the mass of the nucleus νmeasured in units of the reference mass M. Since the proton is the lightest nucleus, μis upper bounded by μ6me/mp≈5.4×10−4. A major advantage of the EF is that the electronic and nuclear DOFs are governed by separate EOMs, allowing us to treat the two EOMs by different approximations. As detailed above, the nuclear EOM is treated in classical approximation by a single trajectory. The smallness of μis not used there. In fact, the limit μ→0 of the nuclear EOM leads to static nuclei [38], which is not a desirable starting point for our purpose. In contrast, in Eq. (6) for the electrons, the smallness of the dimensionless μ strongly suggests to treat μUen as a perturbation, while HBO defines the unperturbed dynamics [39]. The nuclear dynamics is affected by this perturbative treatment only in as much as the electronic wave function |φ(t,R)enters the calculation of the classical forces on the nuclei. If the electronic wave function is evaluated to a given perturbative order, the forces on the nuclei are evaluated consistently to the same order. We write |φ(t,R)= |φ(0)(t,R)+μ|φ(1)(t,R)+O(μ2).(9) Based on the knowledge of the zeroth-order state |φ(0)(t,R)=e−itHBO (R)|φ(0,R), it is straightforward to compute the dominant nonadiabatic correction, |φ(1)(t,R)=−i∫t 0 dte−i(t−t)HBO(R)|ξ(0)(t,R),(10a) |ξ(0)(t,R)=Uen K[φ(0)]+Uen Q[φ(0),χ] − φ(0)(t,R)|Uen K[φ(0)]|φ(0)(t,R) × |φ(0)(t,R).(10b) The pure-state electronic dynamics along a single classical trajectory can thus be explicitly evaluated under the above approximation described by Eqs. (9) and (10) with the involved nuclear momentum function taken to the classical limit depicted by Eq. (7). III. NONUNITARY PURE-STATE ELECTRONIC DYNAMICS A. Nonunitarity and nonadiabaticity along a classical trajectory Our main interests regarding the electronic evolution specific to a single nuclear trajectory Rc trefer to time-dependent pure-state processes, |φ(t,R)|t=0,Rc 0→|φ(t,R)|t>0,Rc t.(11) In principle, such a time-changing state |φ(t,Rc t)is assembled from the full knowledge of |φ(t,R)for arbitrary Runder the substitute of Rby Rc t. The corresponding trajectory-bound electronic evolution operator K(t) is therefore defined by ∣∣φt,Rc t〉=K(t)∣∣φ0,Rc 0〉.(12) We now analyze properties of such trajectory-bound evolution operator K(t), where by construction Rc tis taken as input. As the general solution |φ(t,R)is not available, we practically expand the evolution operator in line with 043075-3 MATISSE WEI-YUAN TU AND E. K. U. GROSS PHYSICAL REVIEW RESEARCH 7, 043075 (2025) Eq. (9), namely, K(t)=K(0)(t)+μK(1)(t)+O(μ2), where the order-moperation is defined by |φ(m)(t,Rc t)= K(m)(t)|φ(0,Rc 0). The knowledge of |φ(m)(t,R)for m=0,1 for arbitrary Ris available once the initial state |φ(0,R) is specified. By starting from various initial BO states |φ(0,R)= |ϕk(R), the later-time states obtained via Eq. (10) are re-denoted by |φ(m)(t,R) → |φ(m) k(t,R), for m=0,1 with χin Eq. (10b) re-denoted by χ0k. This fixes the left-hand side of Eq. (12) to |φ(t,Rc t)= {|φ(0) k(t,R)+μ|φ(1) k(t,R)+O(μ2)}R=Rc ttogether with |φ(0,Rc 0)= |ϕk(Rc 0)fixed on the right-hand side of Eq. (12). Inserting the complete set ∑k|ϕk(Rc t)ϕk(Rc t)| on the right-hand side of Eq. (12), the matrix elements [K(t)]lk ≡ ϕl(Rc t)|K(t)|ϕk(Rc t)appear. Using the knowledge of ϕl(R)|φ(m) k(t,R)for arbitrary R, these matrix elements can then be extracted upon setting R=Rc t, assuming that we also know ϕk(Rc t)|ϕk(Rc 0). These algebra then lead to [K(t)K†(t)−I]lk =−iμ∫t 0 dte−i(t−t)[εl(Rc t)−εk(Rc t)]  ν 1 Mν[p(0)Q νkRc t,t−p(0)Q νlRc t,t]·Dν lk Rc t +Lν lk Rc t−[Lν kl Rc t]∗.(13) Here, Dν lk (R)=ϕl(R)|(−i∇ν−Aν[ϕk](R))|ϕk(R)and Lν lk (R)=ϕl(R)|1 2(−i∇ν−Aν[ϕk](R))2|ϕk(R)are the derivative couplings of first and second order. Equation (13) is the main result of our analysis. The right-hand side only vanishes for l=kimplying that the trajectory-bound first-order nonadiabatic evolution process is by itself nonunitary. Note further that p(0)Q νk(R,t)≡ pQ ν[φ(0) k,χ0k] is in general complex. Its imaginary part given by −∇ν|χ0k(R,t)|/|χ0k(R,t)|is known as the nuclear quantum momentum, which has been taken to be the most essential factor for decoherence [20,22]. Without that imaginary part, p(0)Q νk(Rc t,t) taken along the classical trajectory, Rc tis just the classical momentum. Importantly, even with the classical momentum alone, the nonunitarity contributed by [p(0)Q νk(Rc t,t)−p(0)Q νl(Rc t,t)] ·Dν lk (Rc t) in Eq. (13) is still effective. This analysis thus clarifies that the nonunitarity persists even after taking the classical limit of the nuclei without any consideration of nuclear quantum momentum. We will show in the following example with explicit numerical solution to nonunitary pure-state electronic dynamics [Eqs. (9) and (10)] coupled to single-trajectory classical nuclear dynamics [Eq. (7)] that this nonunitarity without nuclear quantum momentum readily produces decoherence effects. B. Pure-state nonunitarity and decoherence 1. Defining coherence A general unitary propagation has two separate properties. On one hand, it preserves the norm of the propagated state. On the other hand, it keeps initially orthogonal states orthogonal. The nonunitary approach presented here conserves the norm [40]. Nevertheless, it leads to the de-orthogonalisation of initially orthogonal states. In order to explore physical impacts of the abstract mathematical nonunitarity on electronic coherence, we explicitly study two-state systems. Then, the electronic Hilbert space of interest is spanned by two BO states {|ϕ0(R),|ϕ1(R)}. In this case, the eRDO takes the form ρe(t)=∫dR{1  k=0 |χk(R,t)|2|ϕk(R)ϕk(R)| +[χ∗ 1(R,t)χ0(R,t)|ϕ0(R)ϕ1(R)|+H.c.]},(14) where {χk(R,t)}are the coefficients of the Born-Huang expansion, i.e., (r,R,t)=∑kr|ϕk(R)χk(R,t) [41]. Importantly, expanding the electronic factor of the EF in a BornHuang expansion as well, |φ(t,R)=∑kCk(t,R)|ϕk(R), the relevant quantity in the second line of Eq. (14) can be written as χ∗ 1(R,t)χ0(R,t)=|χ(R,t)|2C∗ 1(t,R)C0(t,R).(15) This will simplify the quantification of electronic coherence. In the context of molecular physics, coherence is usually defined by (t)≡dRχ∗ 1(R,t)χ0(R,t) [20,22,42]. It is complex-valued in general and we call |(t)|the coherence magnitude. Approximating the nuclear density by a single classical trajectory |χ(R,t)|2≈δ(R−Rc t) and employing Eq. (15), then yields (t)≈C∗ 1(t,Rc t)C0(t,Rc t). For disambiguity, we denote the coherence computed from the above EF-based single-trajectory approach by EF (t). To gain more insight, it is useful to compare it to another wellestablished single-trajectory approach, namely, the Ehrenfest dynamics. The Ehrenfest electronic state |φEh(t)is described by HBO(Rc t)|φEh(t)=i|˙ φEh(t). The coherence obtained under the Ehrenfest dynamics is then given by Eh(t)= [CEh 1(t)]∗CEh 0(t) with CEh k(t)=ϕk(Rc t)|φEh(t). 2. Messages from the exact solutions of the full problem To be concrete, we focus on a widely studied representative scenario, in which the electronic coherence |(t)|rises and falls (decoherence) along with the passage of the nuclear wave packet through an avoided crossing [20,43]. We take a simple modeling of the avoided crossing by HBO(R)= R(/2)σz+gσx. Here, Ris the nuclear coordinate and σx,σz are the Pauli matrices of the two-state electronic Hilbert space. The parameters gand determine the shape of the avoided crossing. The exact numerical solution to the TDSE of the full problem for this model is available. For the typical scenario of studying the dynamics of electronic coherence along passing the nuclear wave packet through the avoided crossing, the initial state is such that the nuclear density is a Gaussian with the electrons sitting on the upper BO surface [44]. To deepen our understanding of the mechanisms behind electronic decoherence, we now turn to exact numerical solutions of the full electron-nuclear dynamics across three qualitatively distinct nonadiabatic regimes, overviewed in Fig. 1. These regimes differ primarily in the extent of population transfer between the BO electronic states: Regime I corresponds to negligible transfer (<1%), Regime IIA to modest but non-negligible transfer (sub-50%), and Regime IIB 043075-4 ELECTRONIC DECOHERENCE ALONG A SINGLE … PHYSICAL REVIEW RESEARCH 7, 043075 (2025) FIG. 1. Distinct regimes exhibited by the avoided-crossing system used here. Upper panel: the regime of weak nonadiabatic transition (Regime I). The inter-BO-surface population transfer is almost negligible (red-solid/black dashed for upper/lower-surface population). The blue dashdot line in the upper right plot indicates unity population for eye guide. Lower panel: strong nonadiabatic transitions with moderate transfer (Regime IIA) and nearly 50% transfer (Regime IIB). Here, we demonstrate these regimes by setting the model parameters to M=20, =1.0, Pc 0=1.04, and Rc 0= −0.1 for Regime I and M=15,=3.0, Pc 0=15, and Rc 0=−1 for Regime IIA, while Regime IIB uses M=5 with other parameters being the same as Regime IIA. to strong transfer (approaching or exceeding 50%). The correlation between nonadiabaticity and nonunitarity, addressed analytically through Eq. (13) above, will be made more concrete in terms of decoherence later. In Fig. 2, we examine the dynamics of coherence and associated nuclear densities on the upper (|χ1(R,t)|2) and the lower (|χ0(R,t)|2) BO surfaces in these regimes. The first row of this figure plots the exact coherence magnitude |(t)|over time for each regime. While all three exhibit a pattern of rise and fall, the physical mechanisms underlying these temporal profiles actually differ and are more deeply revealed in the FIG. 2. Coherence magnitude |(t)|is shown by solid green lines on the first row. The second/third row takes snapshots of nuclear densities at the time indicated by the black/blue arrow on the time axis of the first row. The upper-surface densities |χ1(t)|2are red solid curves calibrated by the right vertical axis and the lower-surface densities |χ0(t)|2are black dashed lines calibrated by the left vertical axis. snapshots of |χ1/0(R,t)|2shown in the second and third rows, taken respectively at the instances with contrasted large and small magnitudes of coherence. In Regime I (left column), the observed rise and fall of coherence is strongly correlated with population dynamics: A small fraction of the upper population is temporarily transferred to the lower BO surface and then returns to the upper surface. So the product |χ1(R,t)χ0(R,t)|naturally follows this pattern of rise and fall. The second and third row plots show that the region in which the density |χ1(R,t)|2is nonzero overlaps very much with that for |χ0(R,t)|2at the moment where the coherence is high (middle row) as well as to the moment the coherence has decayed to nearly zero (the last row). In this weakly coupled regime, the apparent decay of coherence is a direct consequence of reversible population transfer. The relative phases between χ1(R,t) and χ0(R,t) as well as the time changing overlap between the two wave packets, as the most commonly attributed factors behind decoherence, are not essential in this regime. In contrast, Regime IIB (right column) showcases a more canonical case of electronic decoherence. Here, as the coherence temporarily becomes smaller, the nuclear wave packets on different surfaces also temporarily become spatially separated (see the changes from the second row at higher coherence to the last row with diminished coherence). This aligns with the commonly invoked physical picture of decoherence due to spatial separation and diminishing overlap between nuclear wave packets. Interestingly, Regime IIA (middle column) presents a more nuanced scenario. The coherence decay in this case is not associated with visible spatial separation of the wave packets. In fact, comparing the snapshots for Regime IIA, the spatial overlap between χ0(R,t) and χ1(R,t) is greater at the later time (low coherence) than at the earlier time (high coherence). Hence, the reduction of coherence here cannot be attributed to the loss of spatial overlap. Though the integrand in (t) remains locally nonzero, the overall integral vanishes due to cancellation across R. The well-established dephasing/phase jittering effects are consistent with this observation [37,45]. Altogether, the insights from Fig. 2establish that the observed coherence decay can result from multiple distinct mechanisms. These include reversible population transfer (Regime I), cancellation of contributions across R(Regime IIA) as well as nuclear wave packet separation (Regime IIB). These findings motivate the subsequent exploration of decoherence factors intrinsic to the single-trajectory dynamics of the electronic subsystem that are not directly identifiable from examining the BO-surface-resolved nuclear densities. This sets the stage for our main inquiry in the following subsection. 3. Single-trajectory factors behind the dynamics of coherence Having identified multiple underlying causes for coherence decay in exact quantum dynamics, we now examine to which extent such phenomenon can be captured using singletrajectory methods. Specifically, we compare the coherence dynamics obtained from three approaches: (1) the EF-based single-trajectory scheme developed in this work (EF (t)= C∗ 1(t,Rc t)C0(t,Rc t)), (2) the Ehrenfest method (Eh(t)= [CEh 1(t)]∗CEh 0(t)), and (3) the exact quantum dynamics used as 043075-5 MATISSE WEI-YUAN TU AND E. K. U. GROSS PHYSICAL REVIEW RESEARCH 7, 043075 (2025) FIG. 3. The first row shows trajectories from exact solution (cyan solid), EF-based single-trajectory method (green dashed), and the Ehrenfest method (red dotted). The exact coherence dynamics is replicated by cyan solid lines in the second and the third row. The Ehrenfest coherence with self-consistently obtained trajectories of the first row is shown as the red dashed lines, while that obtained by input exact trajectories is shown as the black dotted lines. The third row shows coherence from the EF-based method with self-consistently obtained trajectories (green dashed) and with exact trajectories (blue dotted). Despite slight deviations in these two trajectories in Regime IIB, the coherence from the EF-based method associated with these two trajectories agrees with each other. a reference benchmark (t)=dRχ∗ 1(R,t)χ0(R,t). Figure 3 summarizes these comparisons across the same three regimes discussed previously. The exact solution calculates the trajectory as the expectation value of the nuclear coordinate, namely, Rc t= dR|χ(R,t)|2R=dR|χ0(R,t)|2R+dR|χ1(R,t)|2R. The first row of Fig. 3demonstrates near-perfect agreement with results obtained from trajectories computed within each respective method. This indicates that the classical approximations on the nuclei from the single-trajectory methods (1) and (2) are appropriate. In the right column of Fig. 3, a slight deviation of the trajectory given by method (1) from the exact solution is observed. To see how such deviations from the exact classical trajectories affect coherence dynamics, in addition to computing coherence using methods (1) and (2) coupled to self-consistently obtained trajectories within those methods themselves, we also input the exact trajectories Rc tto the electronic propagations of these methods to obtain the electronic coherence (see captions for the second and the third rows of Fig. 3). We now discuss the subsequent coherence dynamics. The second and third rows of Fig. 3plot the results from the Ehrenfest and EF-based methods, respectively. Each plot includes the exact coherence (cyan solid line) for direct comparison. The Ehrenfest results in Regime I are close to the exact coherence. In Regimes IIA and IIB, however, the Ehrenfest coherence exhibits a clear failure: After reaching its peak, it decays only partially and then plateaus at a stable nonzero value, in contrast to the complete or near-complete decay seen in the exact solution. This discrepancy persists whether one uses the self-consistent Ehrenfest trajectory (the red-dashed line) or the exact nuclear trajectory as input (the black dotted line). This result reinforces the wellknown failure of Ehrenfest dynamics in capturing electronic decoherence [46–49]. The EF-based single-trajectory method performs differently. While it overestimates the peak magnitude of coherence in Regimes IIA and IIB, it captures the subsequent decay behavior correctly, namely, the coherence decays to nearly zero on a timescale comparable to that of the exact result. The deviation in the classical trajectory of the EF-based method from the exact solution in Regime IIB shows no visible effects on the coherence dynamics, the results in coherence computed from self-consistently obtained classical trajectories (green dashed lines) and that from taking the exact nuclear trajectory as input (blue dotted lines) agree with each other. These observations indicate that something intrinsic to the structure of the EF-based single-trajectory electronic dynamics—namely, its nonunitarity—accounts for the observed decoherence. This conclusion is reinforced by the final row of Fig. 3, which shows that in moving from Regime IIA to IIB, the increasing nonadiabatic transition is accompanied by increased rapidity of losing coherence (higher coherence magnitude is reached earlier in time and decays faster to zero in Regime IIB than Regime IIA). This echoes with the analysis in Sec. III A, which shows that nonadiabatic transitions are responsible for driving nonunitarity in the electronic evolution that is responsible for decoherence. That the unitary Ehrenfest dynamics fails to reproduce coherence decay, while the nonunitary EF-based method succeeds, strongly suggests that nonunitarity triggered by nonadiabatic transitions is the trajectory-native factor behind decoherence. We therefore have verified that the structure of the EF electronic dynamics along a single nuclear classical trajectory inherently encodes decoherence mechanisms, independent of ensemble averaging or environmental entanglement. This constitutes the central result of our investigation. IV. CONCLUSIONS In many mixed quantum-classical approaches, decoherence along a single nuclear trajectory is introduced by manually resetting the electronic states at each time step [6–17]. These decoherence corrections make the electronic evolution effectively nonunitary over each time step, in line with the general understanding that decoherence is associated with departures from unitarity [35]. However, the form and application of such corrections are typically heuristic, relying on intuitive models that utilize the finite widths of the wave packets to determine how coherence should be reduced [6–17]. In contrast, the nonunitary evolution investigated in this work emerges intrinsically from the structure of the EF formalism, even when the nuclear degrees of freedom are described by a single classical trajectory. No a priori resetting procedure is required, nor is there any reliance on the spatial widths of nuclear wave packets. Instead, nonadiabatic transitions directly induce nonunitarity in the dynamics of the electronic conditional state, which in turn leads to the decay of electronic coherence—capturing decoherence natively at the single-trajectory level. 043075-6 ELECTRONIC DECOHERENCE ALONG A SINGLE … PHYSICAL REVIEW RESEARCH 7, 043075 (2025) We expect future applications to periodic solids where the Born-Huang expansion is not feasible in practice, while our perturbative approach is readily applicable, especially when the electronic structure part is treated with the EF version of density functional theory [50,51]. Finally, the approach can be easily generalized to other composite systems since the EF framework is valid for arbitrary twoor multicomponent systems [3]. ACKNOWLEDGMENTS We thank Aaron Kelly, Jonathan Mannouch, and Sebastian de la Pena for useful discussions. M.W.-Y.T. acknowledges the hospitality of Angel Rubio at the Max Planck Institute for the Structure and Dynamics of Matter, Germany. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. ERC-2017-AdG-788890). E.K.U.G. acknowledges support as Mercator Fellow at the University Duisburg-Essen within SFB 1242 funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project No. 278162697. DATA AVAILABILITY The data that support the findings of this article are not publicly available. The data are available from the authors upon reasonable request. [1] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Cambridge University Press, Cambridge, 2002). [2] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information: 10th Anniversary Edition (Cambridge University Press, Cambridge, 2010). [3] X. Gonze, J. S. Zhou, and L. Reining, Variations on the “exact factorization” theme, Eur. Phys. J. B 91, 224 (2018). [4] J. L. Alonso, X. Andrade, P. Echenique, F. Falceto, D. PradaGracia, and A. Rubio, Efficient formalism for large-scale ab initio molecular dynamics based on time-dependent density functional theory, Phys. Rev. Lett. 101, 096403 (2008). [5] X. Andrade, A. Castro, D. Zueco, J. L. Alonso, P. Echenique, F. Falceto, and Á. Rubio, Modified Ehrenfest formalism for efficient large-scale ab initio molecular dynamics, J. Chem. Theory Comput. 5, 728 (2009). [6] B. J. Schwartz, E. R. Bittner, O. V. Prezhdo, and P. J. Rossky, Quantum decoherence and the isotope effect in condensed phase nonadiabatic molecular dynamics simulations, J. Chem. Phys. 104, 5942 (1996). [7] O. V. Prezhdo and P. J. Rossky, Mean-field molecular dynamics with surface hopping, J. Chem. Phys. 107, 825 (1997). [8] O. V. Prezhdo and P. J. Rossky, Evaluation of quantum transition rates from quantum-classical molecular dynamics simulations, J. Chem. Phys. 107, 5863 (1997). [9] M. D. Hack and D. G. Truhlar, A natural decay of mixing algorithm for non-Born-Oppenheimer trajectories, J. Chem. Phys. 114, 9305 (2001). [10] M. J. Bedard-Hearn, R. E. Larsen, and B. J. Schwartz, Meanfield dynamics with stochastic decoherence (MF-SD): A new algorithm for nonadiabatic mixed quantum/classical moleculardynamics simulations with nuclear-induced decoherence, J. Chem. Phys. 123, 234106 (2005). [11] A. W. Jasper and D. G. Truhlar, Electronic decoherence time for non-Born-Oppenheimer trajectories, J. Chem. Phys. 123, 064103 (2005). [12] G. Käb, Fewest switches adiabatic surface hopping as applied to vibrational energy relaxation, J. Phys. Chem. A 110, 3197 (2006). [13] G. Granucci, M. Persico, and A. Zoccante, Including quantum decoherence in surface hopping, J. Chem. Phys. 133, 134111 (2010). [14] J. E. Subotnik and N. Shenvi, A new approach to decoherence and momentum rescaling in the surface hopping algorithm, J. Chem. Phys. 134, 024105 (2011). [15] J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi, Understanding the surface hopping view of electronic transitions and decoherence, Annu. Rev. Phys. Chem. 67, 387 (2016). [16] M. P. Esch and B. G. Levine, State-pairwise decoherence times for nonadiabatic dynamics on more than two electronic states, J. Chem. Phys. 152, 234105 (2020). [17] Y. Shu and D. G. Truhlar, Decoherence and its role in electronically nonadiabatic dynamics, J. Chem. Theory Comput. 19, 380 (2023). [18] A. Abedi, N. T. Maitra, and E. K. U. Gross, Exact factorization of the time-dependent electron-nuclear wave function, Phys. Rev. Lett. 105, 123002 (2010). [19] Y. Suzuki, A. Abedi, N. T. Maitra, K. Yamashita, and E. K. U. Gross, Electronic Schrödinger equation with nonclassical nuclei, Phys. Rev. A 89, 040501 (2014). [20] S. K. Min, F. Agostini, and E. K. U. Gross, Coupledtrajectory quantum-classical approach to electronic decoherence in nonadiabatic processes, Phys. Rev. Lett. 115, 073001 (2015). [21] C. Li, R. Requist, and E. K. U. Gross, Energy, momentum, and angular momentum transfer between electrons and nuclei, Phys. Rev. Lett. 128, 113001 (2022). [22] E. V. Arribas and N. T. Maitra, Electronic coherences in molecules: The projected nuclear quantum momentum as a hidden agent, Phys. Rev. Lett. 133, 233201 (2024). [23] A. Abedi, N. T. Maitra, and E. K. U. Gross, Correlated electron-nuclear dynamics: Exact factorization of the molecular wavefunction, J. Chem. Phys. 137, 22A530 (2012). [24] S. K. Min, A. Abedi, K. S. Kim, and E. K. U. Gross, Is the molecular Berry phase an artifact of the Born-Oppenheimer approximation? Phys. Rev. Lett. 113, 263004 (2014). [25] R. Requist, F. Tandetzky, and E. K. U. Gross, Molecular geometric phase from the exact electron-nuclear factorization, Phys. Rev. A 93, 042108 (2016). [26] R. Requist, C. R. Proetto, and E. K. U. Gross, Asymptotic analysis of the Berry curvature in the E⊗eJahn-Teller model, Phys. Rev. A 96, 062503 (2017). [27] C. L. Lopreore and R. E. Wyatt, Quantum wave packet dynamics with trajectories, Phys. Rev. Lett. 82, 5190 (1999). [28] S. Garashchuk, V. Rassolov, and O. Prezhdo, Semiclassical bohmian dynamics, in Reviews in Computational Chemistry, edited by K. B. Lipkowitz (John Wiley & Sons, Inc., Hoboken, NJ, 2010), Vol. 27, pp. 287–368. 043075-7 MATISSE WEI-YUAN TU AND E. K. U. GROSS PHYSICAL REVIEW RESEARCH 7, 043075 (2025) [29] F. Talotta, F. Agostini, and G. Ciccotti, Quantum trajectories for the dynamics in the exact factorization framework: A proof-of-principle test, J. Phys. Chem. A 124, 6764 (2020). [30] F. Agostini, A. Abedi, and E. K. U. Gross, Classical nuclear motion coupled to electronic non-adiabatic transitions, J. Chem. Phys. 141, 214101 (2014). [31] S. K. Min, F. Agostini, I. Tavernelli, and E. K. U. Gross, Ab initio nonadiabatic dynamics with coupled trajectories: A rigorous approach to quantum (de)coherence, J. Phys. Chem. Lett. 8, 3048 (2017). [32] J.-K. Ha, I. S. Lee, and S. K. Min, Surface hopping dynamics beyond nonadiabatic couplings for quantum coherence, J. Phys. Chem. Lett. 9, 1097 (2018). [33] P. Vindel-Zandbergen, L. M. Ibele, J.-K. Ha, S. K. Min, B. F. E. Curchod, and N. T. Maitra, Study of the decoherence correction derived from the exact factorization approach for nonadiabatic dynamics, J. Chem. Theory Comput. 17, 3852 (2021). [34] J.-K. Ha and S. K. Min, Independent trajectory mixed quantumclassical approaches based on the exact factorization, J. Chem. Phys. 156, 174109 (2022). [35] W. H. Zurek, Decoherence, einselection, and the quantum origins of the classical, Rev. Mod. Phys. 75, 715 (2003). [36] F. A. Bornemann, P. Nettesheim, and C. Schütte, Quantumclassical molecular dynamics as an approximation to full quantum dynamics, J. Chem. Phys. 105, 1074 (1996). [37] G. A. Fiete and E. J. Heller, Semiclassical theory of coherence and decoherence, Phys. Rev. A 68, 022112 (2003). [38] F. G. Eich and F. Agostini, The adiabatic limit of the exact factorization of the electron-nuclear wave function, J. Chem. Phys. 145, 054110 (2016). [39] There are several possibilities to develop perturbation theories based on this small parameter, differing by the choice of unperturbed limit and the physical processes to be treated [38,52–58]. [40] This can be seen straightforwardly by examining ∂t(φ(t,R)|φ(t,R)) using Eq. (6) for the full electronic factor |φ(t,R). Alternatively, one can substitute the expansion of Eq. (9) to arbitrary order m>1 to ∂t(φ(t,R)|φ(t,R)) and yield ∑m n=0φ(n)(t,R)|φ(m−n)(t,R)=0. This result upon taking m=1 translates to φ(t,R)|φ(t,R)=1+O(μ2), which can also be verified explicitly through Eq. (10). Notably, these results hold regardless of how one treats the nuclear factor χ. [41] M. Born and K. Huang, Dynamical Theory of Crystal Lattices (Oxford University Press, Oxford, 1954). [42] E. Villaseco Arribas, N. T. Maitra, and F. Agostini, Nonadiabatic dynamics with classical trajectories: The problem of an initial coherent superposition of electronic states, J. Chem. Phys. 160, 054102 (2024). [43] A. Abedi, F. Agostini, Y. Suzuki, and E. K. U. Gross, Dynamical steps that bridge piecewise adiabatic shapes in the exact time-dependent potential energy surface, Phys. Rev. Lett. 110, 263001 (2013). [44] Since the actual electron mass meas well as the associated Bohr radius evoked in the atomic units used for general analysis is not really involved in models like the one exemplified here, for numerical calculations, it is more convenient to take the units of energy/time and coordinate/momentum directly from the defining parameters of the model. It is natural to take gas the energy unit. Noticeably, when =0, then the non-adiabatic coupling d(R)≡ ϕ0(R)|∂ϕ1(R)/∂Rvanishes identically for all Rs. In order to compare situations of various strengths of non-adiabatic couplings, one should allow to vary. Therefore, though g/has the physical dimension of length, it is not used for the unit of length. The situation =0 corresponds to that the nuclear particle is a free particle with the internal electronic states fixed to one of the two eigenstates of σx, independent of R. As a free particle, the initial gaussian width σconveniently serves as the length unit, which can still serve the purpose well for =0. Then, is measured in unit of g/σand the nuclear mass Mis in unit of σ−2/g(with ¯h=1). [45] M. Vacher, L. Steinberg, A. J. Jenkins, M. J. Bearpark, and M. A. Robb, Electron dynamics following photoionization: Decoherence due to the nuclear-wave-packet width, Phys. Rev. A 92, 040502 (2015). [46] J. C. Tully, Mixed quantum-classical dynamics, Faraday Discuss. 110, 407 (1998). [47] D. V. Makhov, W. J. Glover, T. J. Martinez, and D. V. Shalashilin, Ab initio multiple cloning algorithm for quantum nonadiabatic molecular dynamics, J. Chem. Phys. 141, 054110 (2014). [48] D. A. Fedorov and B. G. Levine, Nonadiabatic quantum molecular dynamics in dense manifolds of electronic states, J. Phys. Chem. Lett. 10, 4542 (2019). [49] M. P. Esch and B. G. Levine, An accurate, non-empirical method for incorporating decoherence into Ehrenfest dynamics, J. Chem. Phys. 155, 214101 (2021). [50] R. Requist and E. K. U. Gross, Exact factorization-based density functional theory of electrons and nuclei, Phys. Rev. Lett. 117, 193001 (2016). [51] R. Requist, C. R. Proetto, and E. K. U. Gross, Exact factorization-based density functional theory of electronphonon systems, Phys. Rev. B 99, 165136 (2019). [52] M. Born and R. Oppenheimer, Zur quantentheorie der molekeln, Ann. Phys. 389, 457 (1927). [53] G. A. Hagedorn, High order corrections to the time-dependent Born-Oppenheimer approximation I: Smooth potentials, Ann. Math. 124, 571 (1986). [54] L. A. Nafie, Velocity-gauge formalism in the theory of vibrational circular dichroism and infrared absorption, J. Chem. Phys. 96, 5687 (1992). [55] A. Scherrer, R. Vuilleumier, and D. Sebastiani, Nuclear velocity perturbation theory of vibrational circular dichroism, J. Chem. Theory Comput. 9, 5305 (2013). [56] A. Scherrer, F. Agostini, D. Sebastiani, E. K. U. Gross, and R. Vuilleumier, Nuclear velocity perturbation theory for vibrational circular dichroism: An approach based on the exact factorization of the electron-nuclear wave function, J. Chem. Phys. 143, 074106 (2015). [57] A. Schild, F. Agostini, and E. K. U. Gross, Electronic flux density beyond the Born-Oppenheimer approximation, J. Phys. Chem. A 120, 3316 (2016). [58] A. Scherrer, F. Agostini, D. Sebastiani, E. K. U. Gross, and R. Vuilleumier, On the mass of atoms in molecules: Beyond the Born-Oppenheimer approximation, Phys. Rev. X 7, 031035 (2017). 043075-8