Electronic decoherence processes in molecules and materials are usually thought and modeled via schemes for the system-bath evolution in which the bath is treated either implicitly or approximately. Here we present computations of the electronic decoherence dynamics of a model many-body molecular system described by the Su-Schrieffer-Heeger Hamiltonian with Hubbard electron-electron interactions using an exact method in which both electronic and nuclear degrees of freedom are taken into account explicitly and fully quantum mechanically. To represent the electron-nuclear Hamiltonian in matrix form and propagate the dynamics, the computations employ the Jordan-Wigner transformation for the fermionic creation/annihilation operators and the discrete variable representation for the nuclear operators. The simulations offer a standard for electronic decoherence that can be used to test approximations. They also provide a useful platform to answer fundamental questions about electronic decoherence that cannot be addressed through approximate or implicit schemes. Specifically, through simulations, we isolate basic mechanisms for electronic coherence loss and demonstrate that electronic decoherence is possible even for one-dimensional nuclear bath. Furthermore, we show that (i) decreasing the mass of the bath generally leads to faster electronic decoherence; (ii) electron-electron interactions strongly affect the electronic decoherence when the electron-nuclear dynamics is not pure-dephasing; (iii) classical bath models with initial conditions sampled from the Wigner distribution accurately capture the short-time electronic decoherence dynamics; (iv) model separable initial superpositions often used to understand decoherence after photoexcitation are only relevant in experiments that employ delta-like laser pulses to initiate the dynamics. These insights can be employed to interpret and properly model coherence phenomena in molecules.

## I. INTRODUCTION

Electronic decoherence is a basic feature of correlated electron-nuclear states^{1} and accompanies photoexcitation,^{2} passage through conical intersections,^{3} electron transfer,^{4} or any other dynamical process that creates superpositions of electronic diabatic states in molecules. Understanding electronic decoherence is central to the description of basic processes such as photosynthesis, vision, and electron transport,^{3–5} to the development of approximation schemes to the vibronic evolution of molecules,^{6,7} and to the isolation of superposition states with robust coherence properties that can subsequently be used in quantum technologies.^{8,9}

Developing deep insights into electronic decoherence requires detailed understanding of the dynamics that entangles the electrons with their nuclear environment. To see this, consider the coherence properties of electrons in a general entangled vibronic state

where the $$*E*_{n}$$ refers to electronic eigenstates (which are obtained by diagonalizing the Born-Oppenheimer electronic Hamiltonian at a particular nuclear geometry) and the $$*χ*_{n}$$ refer to the nuclear wavepackets associated with each $$*E*_{n}$$. The electronic reduced density matrix associated with such a state is given by

where the trace is over the nuclear coordinates. Importantly, note that the coherence between electronic states *n* and *m* is governed by the nuclear wavepacket overlap $Snm=\u27e8\chi m(t)\u2009\chi n(t)\u27e9$. The electron-nuclear entanglement that is generated as the *S*_{nm} decay (for $$*χ*_{m}$$, $$*χ*_{n}$$ ≠ 0) is reflected at the purely electronic level as a decay in the coherences between states *n* and *m*. Standard basis-independent measures of decoherence capture this precisely. For instance, the purity

decays with the overlaps between the environmental states *S*_{nm}. Detailed understanding of electronic decoherence can thus be developed by investigating the dynamics of the *S*_{nm}.^{1}

Because of the difficulty in following the vibronic evolution of molecules exactly, most of our insights into electronic decoherence have emerged from models where the nuclear bath is taken into account implicitly or approximately. In implicit approaches, decoherence effects are modeled via master equations^{10} where the bath is represented using a spectral density with adjustable parameters chosen to reproduce experimental findings. In turn, explicit approaches offer a detailed description of the nuclear dynamics albeit at an approximate level. Approximations that have been employed to model decoherence include surrogate Hamiltonians,^{11} path-integral techniques,^{12} frozen Gaussian approaches,^{13} semi-classical methods,^{14} and quantum-classical methods.^{2,15–20}

Here we present accurate numerical simulations of electronic decoherence in a model molecular system using an exact method that takes into account both electrons and nuclei explicitly and fully quantum mechanically. As a model, we adopt the Su-Schrieffer-Heeger (SSH) Hamiltonian^{21} for *trans*-polyacetylene (PA) because it captures the essential vibronic phenomena of molecules. To quantify the effects of electron-electron interactions on the electronic decoherence,^{22} we augment this model with a Hubbard electron-electron interaction term.^{23} This Hamiltonian is represented in matrix form by the Jordan-Wigner transformation^{24–26} of the electrons and by the discrete variable representation^{27} (DVR) of the nuclei. The decoherence dynamics of such a system is propagated using the Crank-Nicolson method^{28} for a neutral SSH chain in Fock space with 4 electrons coupled to 2 vibrational modes. The advantage of this method is that it does not invoke any physical approximations, facilitating the interpretation of the results. Furthermore, because it solves the many-body problem exactly, it can access regimes that are challenging for other methods such as those encountered when the nuclear mass is small [where the Born-Oppenheimer approximation and even the very concept of a potential energy surface (PES) fail] or when the electron correlation is large. These simulations complement recent efforts to capture electronic decoherence dynamics in molecules using semiclassical approximations,^{2,15,16,18,29} the multiconfiguration time-dependent Hartree (MCTDH) method,^{30–32} and a recently proposed generalized theory for the time scale of the electronic decoherence in the condensed phase.^{33}

In addition to providing a standard for electronic decoherence dynamics in closed molecular systems, the exact modeling serves as a platform from which several basic questions about decoherence in molecules can be addressed. We focus this discussion around seven questions that explore the main requirements for the emergence of decoherence, its basic phenomenology, and the applicability of approximate schemes to capture decoherence. While the simulations pertain to a particular model system, the generic nature of the employed Hamiltonian permits interpreting the insights that result from this model in a broader sense.

Specifically, we investigate the basic mechanisms for the electronic coherence loss (Sec. III A) and the bath size requirements for its emergence (Sec. III B). We also investigate how timescales of electronic and nuclear decoherences compare (Sec. III C) and establish conditions for the accuracy of classical bath models (Sec. III D). In Secs. III E and III F, we study the effect of changing the nuclear mass and the initial-time preparation method on the electronic decoherence dynamics. Finally, in Sec. III G, we study the largely unexplored problem of how electronic interactions modulate the electronic decoherence.^{22} We summarize our main findings in Sec. IV.

## II. MODEL AND METHODS

### A. Model Hamiltonian

As an exemplifying model of a molecular system, we adopt the SSH tight-binding model for PA augmented with a Hubbard electron-electron interaction term, with Hamiltonian

Here,

are the electronic, nuclear, and electron-nuclear interaction components of the Hamiltonian, respectively. In this model, each site *n* represents a CH unit along a PA chain with *N* repeating units. In the electronic component, the fermion creation (annihilation) operators $\u0109n,\Delta \u2020(\u0109n,\Delta )$ create (or annihilate) an electron on site *n* with spin Δ and are subject to the usual fermionic anti-commutation rules. The electronic Hamiltonian describes the hopping process of the *π* electrons with the same spin between different sites with hopping strength *t*_{0} and a Hubbard electron-electron repulsion (*U*$\u2a7e$ 0) when two electrons occupy the same site. In the nuclear component, the *û*_{n} refers to the displacement of the *n*th CH unit from position *na*, where *a* is the lattice constant, $p^n$ is the momentum conjugate to *û*_{n}, *M* is the mass of the CH group, and *K* is the effective spring constant. The electron-nuclear interaction modulates the hopping integral as neighboring nuclei come closer together. In this paper, we use the standard SSH parameters for PA: *α* = 4.1 eV/Å, *K* = 21 eV/Å^{2}, *t*_{0} = 2.5 eV, *M* = 1349.14 eV fs^{2}/Å^{2}, and *a* = 1.22 Å. The quantity *U* is taken to be *U* = 0 except in Sec. III G, where the effect of the electronic correlations on the decoherence is investigated.

Below, we investigate the exact vibronic dynamics for a neutral SSH chain with 4 electrons. The end atoms of the molecule are taken to be clamped, resulting in two vibrational modes: a high frequency C–C stretching optical mode where the two middle nuclei move in the opposite direction with identical amplitude; and a lower frequency acoustic mode where the two middle nuclei move in the same direction with identical amplitude. Because the vibronic coupling is proportional to (*û*_{n} − *û*_{n−1}), the optical mode couples strongly to the electronic degrees of freedom, while the acoustic mode is weakly coupled.

### B. Matrix representation of the Hamiltonian

The matrix representation of the SSH Hamiltonian is constructed by the tensor product of its electronic and nuclear components. To represent the fermionic creation and annihilation operators in matrix form, we adopt the Jordan-Wigner transformation.^{24–26} In this transformation, the fermionic annihilation operators at site *n* with spin Δ can be represented by a string matrix $ei\varphi n,\Delta $ times the corresponding spin-$12$ Pauli lowering matrix *σ*_{n,Δ} at the same site and with the same spin, i.e., $\u0109n,\Delta \u2250ei\varphi n,\Delta \sigma n,\Delta $. Similarly, the fermionic creation operator is represented as $\u0109n,\Delta \u2020\u2250\sigma n,\Delta \u2020e\u2212i\varphi n,\Delta $, where $\sigma n,\Delta \u2020$ is the spin-$12$ Pauli raising matrix at site *n* with spin Δ. Here, the phase matrix *ϕ*_{n,Δ} contains the sum over all the occupation matrices to the left of (*n*, Δ), i.e., $\varphi n,\Delta =\pi \u2211(k,\Delta \u2032)<(n,\Delta )\sigma k,\Delta \u2032\u2020\sigma k,\Delta \u2032$. In our notation, we assume that the spin up (+) is at the left of spin down (−) for each site.

The spin number operator $\sigma n,\Delta \u2020\sigma n,\Delta $ is idempotent. Furthermore, $\sigma n,\Delta \u2020\sigma n,\Delta $ with different (*n*, Δ) commute. Thus, by Taylor expanding each component of $ei\varphi n,\Delta =\u220f(k,\Delta \u2032)<(n,\Delta )ei\pi \sigma k,\Delta \u2032\u2020\sigma k,\Delta \u2032$, the fermionic annihilation operators can be expressed as

where *I*_{i} is the 2 × 2 identity matrix in the *i*th subspace. The fermionic creation operators at site *n* with spin Δ can be represented simply by replacing *σ*_{n,Δ} by $\sigma n,\Delta \u2020$ in Eq. (5). In this way, the resulting operators satisfy the desired fermionic anti-commutation rules (${\u0109n,\Delta ,\u0109n\u2032,\Delta \u2032\u2020}=\delta nn\u2032\delta \Delta \Delta \u2032,{\u0109n,\Delta ,\u0109n\u2032,\Delta \u2032}={\u0109n,\Delta \u2020,\u0109n\u2032,\Delta \u2032\u2020}=0$) because spin-$12$ Pauli matrices satisfy the following relations: ${\sigma n,\Delta ,\sigma n,\Delta \u2020}=1$, ${1\u22122\sigma n,\Delta \u2020\sigma n,\Delta ,\sigma n,\Delta \u2020}=0$, and ${1\u22122\sigma n,\Delta \u2020\sigma n,\Delta ,\sigma n,\Delta}=0$.

The Fock space in which the creation and annihilation operators are defined includes all possible electronic number states. To reduce the computational effort, we project the electronic Fock space to a Hilbert space with a fixed *n*_{e} number of electrons. This is possible since the SSH Hamiltonian commutes with the electron number operator $N^e=\u2211n=1N\u2211\Delta =\xb1\u0109n,\Delta \u2020\u0109n,\Delta $,

and thus the dynamics preserves *n*_{e}. The size of the net electronic basis is $22ne$ and can represent any many-body state in the electronic system with a fixed number of electrons *n*_{e}.

To represent the nuclear operators in terms of matrices, we employ the Discrete Variable Representative (DVR) as proposed in Ref. 27. For simplicity, we illustrate this method with one nuclear degree of freedom. For this degree of freedom, a basis consisting of grid points, {|*i*⟩}, is employed. The matrix elements of the kinetic energy operator in this basis are

where Δ*x* is the grid spacing. Correspondingly, the matrix elements of the position dependent function *V*(*û*) are

This idea can be extended to many degrees of freedom and is used to represent the nuclear component of the SSH Hamiltonian in matrix form.

DVR methods have been proved to be highly accurate to solve a variety of problems in molecular quantum dynamics and vibration-rotation spectroscopy.^{34} In this study with two nuclear degrees of freedom, the method provides convergence with respect to the grid spacing at Δ*x* = 0.02 Å for both degrees of freedom. These results cannot be achieved by the general grid basis method in which the second derivative in kinetic term is represented by a tridiagonal matrix because this method requires much smaller grid spacing leading to large memory needs. All results presented here have been tested for convergence in the grid spacing and the range of space considered in the simulation.

### C. Dynamical propagation

To propagate the dynamics, we employ the Crank-Nicolson scheme.^{28} In it, the evolution operator *Û*(*t* + Δ*t*, *t*) from time *t* to *t* + Δ*t* is computed by a first order Padé approximation^{35} as

In this way, the propagation from |Ψ(*t*)⟩ to |Ψ(*t* + Δ*t*)⟩ is transformed into the solution of the linear equation

where

This linear equation is solved using a biconjugate gradient stabilized iterative method,^{36} which is stable and faster than direct methods such as Gaussian elimination. Unless specified otherwise, results are reported with time step Δ*t* = 0.01 fs. For the model with standard SSH parameters prepared in separable superpositions states, this time step offers converged results (within 2% error) up to ∼500 fs, as verified by performing the dynamics with Δ*t* = 0.001 fs. For larger times, it offers qualitatively correct dynamics.

## III. RESULTS AND DISCUSSION

The decoherence dynamics is investigated by following purity *P*(*t*) [Eq. (3)] which is a well-defined basis-independent measure of decoherence. The quantity *P* = 1 for a pure state, *P* < 1 for a mixed state, and *P* = 1/*M* for a maximally entangled state of *M* levels with equal populations. Note that a decay in purity directly signals coherence loss. This contrasts with electric polarization-based measures of coherence used in laser spectroscopy in which in addition to an overall decay that signals decoherence these measurements exhibit oscillations determined by Bohr transition frequencies in the system that are not linked to decoherence. Since the vibronic dynamics is solved exactly, the simulations take into account all possible potential energy surfaces (PESs) for the molecule and can access regimes where the Born-Oppenheimer picture or mixed quantum-classical (MQC) schemes are inadequate.

The following analysis is framed around seven questions on the electronic decoherence that are addressed in the light of the exact method.

### A. What are the main mechanisms for the electronic coherence loss?

There has been significant discussion in the literature about the main mechanisms for electronic decoherence in molecules.^{1,10,14–16,22,29–31,37} From Eq. (3), it is clear that the decoherence arises due to the nuclear wavepacket evolution in alternative diabatic PESs that lead to a decay in the nuclear wavepacket overlaps and thus to electron-nuclear entanglement. A recently proposed theory of the electronic decoherence shows that, in the short time, this can be divided into pure-dephasing dynamics, transitions between diabatic states, and their interference.^{33} The pure-dephasing dynamics refers to the electron-nuclear evolution that does not involve transitions into other electronic diabatic states.

To connect with these previous efforts and illustrate the main mechanisms for coherence loss in this model, we first consider the case in which the SSH chain is prepared in a separable tensor product of the form

where the electrons are initially in a superposition state of the ground and the first excited electronic state determined at the equilibrium nuclear coordinates *u*_{0} of the ground PES, while the initial nuclear state $$*χ*(*t* = 0)$$ is taken to be the ground vibrational state of the ground PES. As shown in Fig. 1 (left panel), these two PESs are well separated in energy and there are no avoided crossings, making the pure-dephasing model of the decoherence applicable. In this model, the vibronic evolution leads to an entangled vibronic state of the form

with initial condition $|\chi 0(t=0)=|\chi 1(t=0)=\u2009\chi (t=0)$. For this particular case, since $\chi 0(t)\u2009$ is stationary, $P=12+12|\u27e8\chi 0(t)|\chi 1(t)\u27e9|2=12+12|\u27e8\chi 0(0)|\chi 1(t)\u27e9|2=12+12|\u27e8\chi 1(0)|\chi 1(t)\u27e9|2=12+12|A(t)|2$, where $A(t)=\u27e8\chi 1(0)\u2009\chi 1(t)\u27e9$ is the autocorrelation function of the excited state nuclear wavepacket. Thus, *P*(*t*) is determined by *A*(*t*) provided that the dynamics is pure-dephasing. The autocorrelation function can be computed without propagating the quantum state as

where {*ϵ*_{n}} is the vibrational energy spectrum of the first excited PES and {*a*_{n} = ⟨*ϕ*_{n}$$*χ*_{1}(0)⟩} are the components of the initial state $$*χ*_{1}(*t* = 0)$$ projected along the vibrational eigenstates {|*ϕ*_{n}⟩} of the excited PES.

As shown in Fig. 2, the purity decay computed by the autocorrelation function is in quantitative agreement with the decoherence dynamics obtained via dynamical propagation with time step Δ*t* = 0.001 fs of the chain starting from the superposition state in Eq. (13). For this time step, the dynamics is essentially exact. Thus, a pure-dephasing picture is valid to illustrate the main mechanisms of the electronic coherence loss when starting from the state in Eq. (13).

We identify three distinct regions for the electronic decoherence. In the first 4 fs, the purity exhibits a Gaussian decay that arises due to the initial wavepacket motion on the excited state diabatic PES. Such initial dynamics leads to a decay in the nuclear overlap $|\u27e8\chi 0(0)|\chi 1(t)\u27e9|2$ as the snapshots of the nuclear wavepacket at *t* = 6 fs and 10 fs show. This initial decay of wavepacket overlap is the mechanism for the short-time decoherence that is captured by theories for decoherence timescales^{29,33,37} and that is expected to be dominant in condensed phase environments. Nevertheless, in this model, since the excited PES is bounded and of low dimensionality, the nuclear wavepacket eventually returns to its starting point leading to a recurrence in the purity corresponding to the snapshot at *t* = 32 fs. However, due to the anharmonicity of the first excited PES, the recurrences in the purity are never complete and this leads to an overall decay in the purity of the system. For longer times, the nuclear wavepacket is spread along the optical mode of the PES as the snapshot at *t* = 15000 fs illustrates, which leads to a stable decay in purity with small high frequency oscillations. Thus, this model recovers the initial Gaussian decay of purity that is expected to be dominant in the condensed phase.^{29–31,37} Furthermore, it captures recurrences that can be observed in small molecular systems and their decay in this case due to anharmonicities in the PES. At even longer times, beyond these two regions, fractional revivals of purity are observed (see Fig. 3), which is consistent with quantum revival theory.^{38} This long time behavior is beyond the applicability of the method, but can be estimated using Eq. (15). This fractional revival structure has also been encountered in experiments investigating the vibrational wavepacket evolution of Br_{2} molecules in the presence of a solid Ar environment via ultrafast pump-probe spectroscopy^{39} and in those studying the evolution of Rydberg electronic wavepackets in K atoms by photoionization measurements.^{40}

In this paper, we focus on the first two regions since they are expected to be the most relevant for molecules.^{15,30,31,41} Determining on general grounds how many recurrences are expected for a given molecule is a challenging problem. For the SSH model, mixed quantum-classical simulations in Ref. 15 indicate that as the number *N* of carbon atoms increases the number of visible recurrences decreases. In these simulations, for *N* = 50, only the initial Gaussian decay is observed. Furthermore, other examples in the literature suggest that for large molecules the initial Gaussian decay is often the dominant feature. This has been seen in Refs. 30 and 31 for water and paraxylene and in Ref. 33 for a displaced harmonic oscillator model in the high-temperature limit. However, there may be particular instances where several recurrences are observed even in large molecules, as claimed in Ref. 41.

### B. How large should a bath be in order for decoherence to emerge?

Another basic question in electronic decoherence is to determine the size of the bath required for decoherence to emerge. Investigations for a spin coupled to a bath of spins^{42} indicate that a bath as small as 20 spins is sufficient to generate a Gaussian decay in the overlap of the environmental states. The question is as follows: In a typical molecule, is decoherence only salient in the condensed phase where the electrons couple to a macroscopic number of bath degrees of freedom or are a few vibrational coordinates enough, and if so, how many?

As shown in Fig. 2, two vibrational degrees of freedom are enough to generate decoherence. However, as discussed below, in fact, just one vibrational coordinate is enough for decoherence to emerge since only the optical mode plays an important role during the dynamics.

To see this, consider Fig. 4 which shows the purity decay computed through the autocorrelation function along the optical mode on the first excited PES (1D model) and that on the full first excited PES (2D model). In this figure, the purity of the two models essentially coincides. Thus, we conclude that the coupling of the electrons to one vibrational degree of freedom (in this case, the optical high frequency mode) is sufficient to induce the electronic decoherence. While a macroscopic condensed phase is not required for the emergence of electronic coherence loss, a larger number of bath degrees of freedom can prevent the emergence of partial recurrences in the purity, as observed in Refs. 15, 30, and 31.

### C. How do the timescales for electronic and vibrational decoherence compare?

Conventional wisdom indicates that electronic decoherence is fast (∼10 fs), while vibrational decoherence is slow (∼10^{2}–10^{3} fs).^{43} However, in this model, the electronic decoherence rate is identical to the nuclear one, i.e., *P*_{e}(*t*) = *P*_{N}(*t*). This is a consequence of the Schmidt theorem^{9,44} (or the Carlson-Keller theorem^{45}), which indicates that for closed system-bath systems the purity of the system and the bath coincides (see Ref. 1 for a simple derivation of this fact). What is the origin of this apparent discrepancy?

To resolve this, consider a molecule in the presence of a solvent. For short times, the total purity decay of molecular vibrational degrees of freedom is just the product of the purity decay due to entanglement with electrons and with the solvent^{37}

where *P*_{N−e}(*t*) is the purity decay due to intra-molecular electron-nuclear coupling (*N* − *e*) and *P*_{N−s}(*t*) is that due to coupling to the solvent (*N* − *s*). Here, $\tau d(e)$ and $\tau d(s)$ are the corresponding decoherence times. The purity decay due to *N* − *e* interactions is usually faster than that due to *N* − *s* interactions (i.e., $\tau d(e)\u226a\tau d(s)$) as a consequence of the difference in timescales of electronic and solvent motion.

The origin of the apparent discrepancy is clarified by considering which bath generates the vibrational decoherence. Reference 43 considers decoherence of a vibrational state associated with a single electronic crude Born-Oppenheimer state caused by a solvent. In this case, there is no appreciable entanglement between the electrons and vibrations, i.e., *P*_{N−e} ≈ 1, and the purity decay due to solvent dominates, i.e., *P*_{N} ≈ *P*_{N−s}. Thus, in this scenario, the vibrational decoherence would be slower than the electronic decoherence.

By contrast, consider now the case in which the vibrations entangle both with electrons and the solvent, as would be the case for a molecule prepared in state Eq. (13) and immersed in a solvent. In this case, the purity decay due to the *N* − *e* interaction dominates as $\tau d(e)\u226a\tau d(s)$ and *P*_{N} ≈ *P*_{N−e} for short times. Thus, in this case, the vibrational decoherence is expected to have an initial timescale for coherence loss identical to one of the electronic decoherences as the Schmidt theorem indicates, followed by a slower decay due to the solvent.

### D. Are classical decoherence models accurate?

From a practical perspective, explicit decoherence modeling requires approximate description of the system-bath dynamics. One common approximation of practical importance in molecules is to model the nuclear degrees of freedom classically.^{2,15–20,46,47} In this case, decoherence is captured by propagating an ensemble of quantum-classical trajectories, each one evolving unitarily, with initial conditions sampled from an appropriate classical distribution meant to mimic the initial nuclear quantum state.^{15,16} The corresponding ensemble average of unitary quantum-classical evolutions mimics the nonunitary evolution of the density matrix of the system. This should be contrasted with true decoherence where a single-quantum system becomes entangled with environmental degrees of freedom and the unitary deterministic evolution of the system plus environment leads to a nonunitary evolution of the reduced density matrix of the system.

To what extent are classical bath models able to capture quantum decoherence processes in molecules?

To test this, we contrasted the exact decoherence dynamics of the SSH chain obtained in Sec. III A as prepared in state Eq. (13) with results obtained in a mixed quantum-classical (MQC) approximation where the nuclei move along a given fixed PES and the electrons instantaneously respond to the nuclear coordinates. In the MQC approximation, we choose the initial conditions for the nuclei by sampling from the Wigner distribution of the ground vibrational state of the ground PES. At time *t*, the electronic density matrix for the *i*th trajectory represented in the diabatic basis |*E*_{0}(*u*_{0})⟩, |*E*_{1}(*u*_{0})⟩ is

where *E*_{0}(*u*^{(i)}(*t*)) and *E*_{1}(*u*^{(i)}(*t*)) are the ground and first excited electronic energy at nuclear geometry *u*^{(i)}(*t*). To represent the electronic density operator in matrix form, as in Eq. (17), it is supposed that in the region where the nuclear wavepacket is distributed the electronic states are well approximated by diabatic states (with no dependence on the nuclear coordinates) obtained by diagonalizing the Born-Oppenheimer electronic Hamiltonian at the ground-state minimum energy geometry. The quantities *E*_{0} and *E*_{1} are determined by diagonalizing the full SSH Hamiltonian [Eq. (4)] at fixed nuclear positions *u*^{(i)}(*t*) encountered during the dynamics. The *u*^{(i)}(*t*) for each trajectory is determined by solving Newton’s equations of motion in a given potential *V*(*u*(*t*)). The electronic density matrix of the ensemble is taken to be the average over *N*_{traj} trajectories,

Using the average electronic density matrix, the purity is computed, as in Eq. (3). Figure 5 compares the electronic decoherence dynamics generated by moving the nuclei classically along the ground PES (*V* = *E*_{0}(*u*(*t*))), first excited PES (*V* = *E*_{1}(*u*(*t*))), the mean field PES between the two (*V* = (*E*_{0}(*u*(*t*)) + *E*_{1}(*u*(*t*)))/2), and a flat PES (*V* = 0) with the exact results. For short times (top panel), electronic decoherence dynamics generated by MQC along *any* potential essentially coincides with the exact method.

These simulations numerically validate a recent theoretical analysis^{37} which indicates that MQC methods correctly capture the initial purity decay when the initial state is sampled from the Wigner distribution, irrespective of the potential that is employed in the dynamics. Beyond short times, the simulations show that MQC schemes can capture some of the quantum recurrences, albeit the dynamics beyond short times depends sensitively on the potential *V* and severely overestimates the decoherence for this model.

Thus, MQC schemes with initial Wigner sampling for decoherence are expected to be accurate in the condensed phase, where the decoherence time is believed to be governed by the initial decay of purity.

### E. How does electronic decoherence vary with the mass of the nuclear bath?

It is challenging to intuitively predict what would be the effect of changing the mass of the nuclear bath on electronic decoherence timescales. Heavier nuclei move slower, and hence the decoherence is expected to be slower too. However, increasing the mass (decreasing the frequency) of the bath also makes the energy spectrum of the bath denser and this is expected to lead to faster decoherence. Which of these two processes is dominant? Furthermore, what happens as the nuclear mass becomes comparable to the mass of the electron and the system is beyond the regime of applicability of the Born-Oppenheimer approximation?

Insights into these problems were obtained by computing the purity dynamics for SSH chains initially prepared as in Eq. (13) with varying masses (0.01M, 0.05M, 0.25M, 0.5M, 1M, 2M, 3M; *M* = 1349.14 eV fs^{2}/Å^{2}); see Fig. 6. In different time windows, for 0.25M-3M, as mass decreases, we observe a faster decay of the initial purity decay, an earlier appearance of first recurrence peak, and a faster decay to the asymptotic purity behavior. All of these observations support our first hypothesis but contradict the second hypothesis. To differentiate the effects of these two possible contributions to the purity dynamics, we eliminate the effect of the speed of nuclear bath motion on the decoherence by studying how the decoherence changes versus the number of nuclear oscillation periods by defining a dimensionless time $\tau =\omega t2\pi $. Here, the frequency of the optical mode is used and $\omega =3KM$. The envelope of the purity decay versus *τ* for various masses is shown in the bottom panels of Fig. 6. It is clear that for 0.25M-3M the different decoherence versus *τ* plots overlap with one another. Thus, for these masses, we conclude that the effect of the energy spectrum of the bath on the decoherence is not important, even though the vibrational energy spectrum for heavier masses is, in fact, denser (see Fig. 7).

As shown in Fig. 6(b), for even smaller masses (0.01M and 0.05M), the purity decay is more complex and deviates from that observed for the other masses as the Born-Oppenheimer approximation starts to break down. However even in these cases, we observe that a lighter bath leads to faster short time decoherence.

### F. What is the effect of preparation on the decoherence dynamics?

In Sec. III A, III B, III D, and III E, we simulated the electronic decoherence dynamics by choosing a separable superposition state of the molecule as an initial state. This is a common strategy when defining decoherence time.^{2,15,16,29–31,48,49} However, experiments routinely use lasers to excite molecules and prepare the initial state. Are initial separable superpositions representative of the experimental situation? How are the decoherence dynamics affected by initial state preparation?

To understand this, we investigated the decoherence dynamics generated by photoexcitation with lasers of different durations. Specifically, we consider the SSH molecule in interaction with a laser in dipole approximation with radiation-matter interaction $\u0124RM=\u2212\mu \u2192^\u22c5\u03f5\u21920E(t)$. Here $\mu \u2192^$ is the dipole operator of the molecule, $\u03f5\u21920$ is the polarization of light, and *E*(*t*) is its electric field. The laser’s electric field is obtained from

where *A*(*t*) is the vector potential defined as

where *A*_{0} is the amplitude of the vector potential, $tw$ is the width of the Gaussian envelope centered at *t*_{0}, and *ω*_{0} is the central frequency of the laser pulse.

To study the effect of preparation on decoherence, a short $tw$ = 2 fs and a long $tw$ = 100 fs pulse are used to resonantly photoexcite the exact ground eigenstate of the SSH Hamiltonian. The laser frequency *ω*_{0} is chosen to be at resonance with the energy difference between the ground PES and first excited PES at the ground state equilibrium geometry. The frequency content of the short laser pulse (Δ*ω* ∼ 1 eV) can excite several vibronic transitions, while the long pulse (Δ*ω* ∼ 0.02 eV) can only excite a single (0-0) vibronic transition.

The decoherence dynamics induced by the two pulses is shown in Fig. 8. Before photoexcitation, the purity is close but not exactly 1 because the exact ground state of the molecule is entangled. The short pulse (top panel) creates a vibrational wavepacket in the excited state. In the limit in which the pulse is a delta kick and the Frank-Condon approximation is valid, this wavepacket would be identical to the one in the ground state. Since this pulse is short enough, the subsequent purity decay is reminiscent to that in Fig. 2. By contrast, the 100 fs laser pulse (bottom panel) performs a state-to-state photoexcitation and no recurrences are observed. The lower purity after photoexcitation reflects the entanglement properties of the final state.

As shown, the preparation mode has a strong influence on the decoherence dynamics. Importantly, model separable initial superposition states often used to understand decoherence after photoexcitation are only relevant in the experimental situation where a short delta-like laser pulse is employed to initiate the dynamics.

### G. How does the electron-electron interaction affect the electronic decoherence?

We now numerically address the largely unexplored connection between electronic interactions and decoherence. According to the formal analysis in Ref. 22, electron-electron interactions can influence the decoherence rate only when the decoherence dynamics is not pure-dephasing in nature, i.e., when [*Ĥ*_{e}, *Ĥ*_{int}] ≠ 0. Below we numerically investigate this claim and assess the quantitative effect of changing the degree of electron repulsion on the decoherence. For this, we vary *U* from 0 eV to 0.8 eV in Eq. (4) and investigate its effect on the purity when the system is prepared in (i) the superposition in Eq. (13) whose dynamics is well approximated by a pure-dephasing model and (ii) a superposition of the form

which involves electronic states with avoided crossings and therefore is not expected to be pure-dephasing in nature. Here, |*E*_{9}(*u*_{0})⟩ is the ninth excited electronic state determined at the equilibrium nuclear coordinates *u*_{0} of the ground PES (see Fig. 1).

Consider first the initial state in Eq. (13). The electronic decoherence dynamics with various electron-electron interaction strengths *U* is shown in Fig. 9(a). As shown, the decoherence dynamics observes minor changes when the electron-electron interaction strength *U* is increased in this case, which is in agreement with the analysis in Ref. 22. From the perspective of the PES, the reason for this behavior is that the electron-electron interaction does not introduce significant changes to the shape of the ground and first excited PES (Fig. 1).

By contrast, when the system is initially prepared in state Eq. (21), changing *U* has a strong effect on the purity dynamics. Specifically, as shown in Fig. 9(b) (bottom panel), the short-time purity dynamics is not affected by changing *U*. However, after 5 fs, *P*(*t*) changes strongly as *U* is varied because changing *U* changes the shape of the PESs and the magnitude of the avoided crossings involved. The fact that changing *U* strongly influences *P*(*t*) is in agreement with the formal analysis in Ref. 22. The numerical observation that for the first 5 fs *P*(*t*) is not affected by changing *U* indicates that this segment of the dynamics is approximately pure-dephasing in nature. In fact, an analysis of the population of the diabatic states (not shown) indicates that significant population transfer to other diabatic states (4 as labeled in Fig. 1) only occurs after 5 fs.

Thus, changing *U* can influence the decoherence dynamics by changing the shape of the PESs and by suppressing or enhancing avoided crossings. The numerical observations are consistent with the proposal in Ref. 22 that the electron correlation and the electronic decoherence only couple for non-pure-dephasing dynamics.

## IV. CONCLUSIONS

In this work, we have presented numerical simulations of electronic decoherence in a model vibronic SSH molecule with 4 electrons and 2 vibrations using an exact method that takes into account electrons and nuclei explicitly and fully quantum mechanically. The simulations serve as a standard of electronic decoherence in molecules and address several fundamental questions about electronic decoherence.

Specifically, we show the following:

Coupling to one anharmonic vibrational coordinate is sufficient for electronic decoherence to emerge. While in condensed phase environments no recurrences are expected, in small molecules the decoherence dynamics from initial separable superposition states exhibits recurrences in the purity and its decay.

While vibrational decoherence is usually considered to be slower than electronic decoherence, their early-time decoherence timescales can coincide, even in the condensed phase, when electrons and vibrations get entangled during the dynamics.

A class of mixed quantum-classical (MQC) methods with initial Wigner sampling capture the early-time decoherence correctly, irrespective of the potential employed in the propagation. This is in agreement with the theoretical developments in Ref. 37 and suggests that MQC can be adequate to capture electronic decoherence in the condensed phase.

Decreasing the nuclear mass generally leads to faster decoherence. During the initial Gaussian purity decay, this feature was observed even for masses for which the Born-Oppenheimer approximation is expected to fail.

The initial preparation has a strong influence on the decoherence dynamics. Model separable initial superpositions often used to understand decoherence after photoexcitation are only relevant in experiments that employ delta-like laser pulses to initiate the dynamics.

Electron-electron interactions can strongly affect the electronic decoherence when the electron-nuclear dynamics is not pure-dephasing. Electronic interactions influence the decoherence by changing the shape of the PES and modulating the strength of non-adiabatic effects in the dynamics. These numerical observations agree with and give insights into the formal argument in Ref. 22.

While the results pertain to a specific model system, the generic nature of the employed Hamiltonian permits interpreting them in a broader sense. As such, these insights can be used to interpret and properly model coherence phenomena in chemistry.

## ACKNOWLEDGMENTS

This material is based upon the work supported by the National Science Foundation under No. CHE-1553939. I.F. thanks Heiko Appel, Lorenzo Stella, and Angel Rubio for useful discussions.