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.

Electronic decoherence is a basic feature of correlated electron-nuclear states1 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

|Ω(t)=n|En|χn(t),
(1)

where the En 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 En. The electronic reduced density matrix associated with such a state is given by

ρ^e(t)=TrN[|ΩΩ|]=nmχm(t)|χn(t)EnEm,
(2)

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=χm(t)χn(t). The electron-nuclear entanglement that is generated as the Snm 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

P(t)=Tr[ρ^e2(t)]=nm|Snm|2
(3)

decays with the overlaps between the environmental states Snm. Detailed understanding of electronic decoherence can thus be developed by investigating the dynamics of the Snm.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 equations10 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) Hamiltonian21 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 transformation24–26 of the electrons and by the discrete variable representation27 (DVR) of the nuclei. The decoherence dynamics of such a system is propagated using the Crank-Nicolson method28 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.

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

ĤSSH=Ĥe+ĤN+Ĥint.
(4)

Here,

Ĥe=t0n=1N1Δ=±(ĉn+1,Δĉn,Δ+ĉn,Δĉn+1,Δ)+Un=1Nĉn,+ĉn,+ĉn,ĉn,,ĤN=n=1Np^n22M+K2n=1N1(u^n+1u^n)2,Ĥint=αn=1N1Δ=±(ĉn+1,Δĉn,Δ+ĉn,Δĉn+1,Δ)(ûn+1ûn)

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 ĉn,Δ(ĉn,Δ) 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 t0 and a Hubbard electron-electron repulsion (U 0) when two electrons occupy the same site. In the nuclear component, the ûn refers to the displacement of the nth 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, t0 = 2.5 eV, M = 1349.14 eV fs22, 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.

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ϕn,Δ times the corresponding spin-12 Pauli lowering matrix σn at the same site and with the same spin, i.e., ĉn,Δeiϕn,Δσn,Δ. Similarly, the fermionic creation operator is represented as ĉn,Δσn,Δeiϕn,Δ, where σn,Δ 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., ϕn,Δ=π(k,Δ)<(n,Δ)σk,Δσk,Δ. In our notation, we assume that the spin up (+) is at the left of spin down (−) for each site.

The spin number operator σn,Δσn,Δ is idempotent. Furthermore, σn,Δσn,Δ with different (n, Δ) commute. Thus, by Taylor expanding each component of eiϕn,Δ=(k,Δ)<(n,Δ)eiπσk,Δσk,Δ, the fermionic annihilation operators can be expressed as

ĉ1+σ1+I1I2+IN+INĉ112σ1+σ1+σ1I2+IN+INĉ2+12σ1+σ1+12σ1σ1σ2+IN+INĉN12σ1+σ1+12σ1σ112σ2+σ2+12σN+σN+σN,
(5)

where Ii is the 2 × 2 identity matrix in the ith subspace. The fermionic creation operators at site n with spin Δ can be represented simply by replacing σn by σn,Δ in Eq. (5). In this way, the resulting operators satisfy the desired fermionic anti-commutation rules ({ĉn,Δ,ĉn,Δ}=δnnδΔΔ,{ĉn,Δ,ĉn,Δ}={ĉn,Δ,ĉn,Δ}=0) because spin-12 Pauli matrices satisfy the following relations: {σn,Δ,σn,Δ}=1, {12σn,Δσn,Δ,σn,Δ}=0, and {12σn,Δσn,Δ,σn,Δ}=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 ne number of electrons. This is possible since the SSH Hamiltonian commutes with the electron number operator N^e=n=1NΔ=±ĉn,Δĉn,Δ,

[ĤSSH,N^e]=0,
(6)

and thus the dynamics preserves ne. 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 ne.

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

i|T^|i=2(1)ii2M(Δx)2π2/3,i=i2(ii)2,ii,
(7)

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

i|V(û)|i=V(u(i))δii.
(8)

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.

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é approximation35 as

Û(t+Δt,t)=1iΔt2Ĥ(t+Δt2)1+iΔt2Ĥ(t+Δt2).
(9)

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

L^|Ψ(t+Δt)=|b,
(10)

where

L^=1+iΔt2Ĥ(t+Δt2),
(11)
|b=(1iΔt2Ĥ(t+Δt2))|Ψ(t).
(12)

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.

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.

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

|Ω(t=0)=12(|E0(u0)+|E1(u0))|χ(t=0),
(13)

where the electrons are initially in a superposition state of the ground and the first excited electronic state determined at the equilibrium nuclear coordinates u0 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

|Ω(t)=12(|E0(u0)|χ0(t)+|E1(u0)|χ1(t)),
(14)

with initial condition |χ0(t=0)=|χ1(t=0)=χ(t=0). For this particular case, since χ0(t) is stationary, P=12+12|χ0(t)|χ1(t)|2=12+12|χ0(0)|χ1(t)|2=12+12|χ1(0)|χ1(t)|2=12+12|A(t)|2, where A(t)=χ1(0)χ1(t) 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

|A(t)|2=|χ1(t=0)|χ1(t)|2=n,m=0|an|2|am|2ei(ϵnϵm)t/,
(15)

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

FIG. 1.

Slice of the potential energy surfaces of the SSH Hamiltonian along u2 = −u3 with electron-electron interactions U = 0 eV (left) and U = 0.8 eV (right). Here, u0 is the ground-state equilibrium geometry. The numbers 0-9 (in red) are used to label the first 10 electronic states. For definitiveness, the figure only shows those states with a zero net spin along the z direction.

FIG. 1.

Slice of the potential energy surfaces of the SSH Hamiltonian along u2 = −u3 with electron-electron interactions U = 0 eV (left) and U = 0.8 eV (right). Here, u0 is the ground-state equilibrium geometry. The numbers 0-9 (in red) are used to label the first 10 electronic states. For definitiveness, the figure only shows those states with a zero net spin along the z direction.

Close modal

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).

FIG. 2.

Decoherence dynamics for a neutral SSH chain with 4 electrons. The plots show the purity dynamics for a chain initially prepared in the superposition state with equal coefficients between the ground and first excited state in Eq. (13) using a time step Δt = 0.001 fs. The decay and recurrences in the purity reflect the vibrational dynamics in the excited state anharmonic PES. Snapshots of the nuclear probability density in the excited diabatic state are shown in the upper panels. The purity dynamics assuming a pure-dephasing model is shown in red.

FIG. 2.

Decoherence dynamics for a neutral SSH chain with 4 electrons. The plots show the purity dynamics for a chain initially prepared in the superposition state with equal coefficients between the ground and first excited state in Eq. (13) using a time step Δt = 0.001 fs. The decay and recurrences in the purity reflect the vibrational dynamics in the excited state anharmonic PES. Snapshots of the nuclear probability density in the excited diabatic state are shown in the upper panels. The purity dynamics assuming a pure-dephasing model is shown in red.

Close modal

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 |χ0(0)|χ1(t)|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 timescales29,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 Br2 molecules in the presence of a solid Ar environment via ultrafast pump-probe spectroscopy39 and in those studying the evolution of Rydberg electronic wavepackets in K atoms by photoionization measurements.40 

FIG. 3.

Fractional revivals in the decoherence dynamics of a SSH chain for the initial state in Fig. 2. The decoherence was computed through the autocorrelation function as P=12+12|A(t)|2.

FIG. 3.

Fractional revivals in the decoherence dynamics of a SSH chain for the initial state in Fig. 2. The decoherence was computed through the autocorrelation function as P=12+12|A(t)|2.

Close modal

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.

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 spins42 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.

FIG. 4.

The purity decay along the u2 = −u3 direction (optical mode) on the first excited PES (1D model, black) and along the full first excited PES (2D model, red). The red line cannot be seen as it coincides exactly with the 1D model black line. This implies that the nuclear evolution along the optical mode dominates the decoherence dynamics.

FIG. 4.

The purity decay along the u2 = −u3 direction (optical mode) on the first excited PES (1D model, black) and along the full first excited PES (2D model, red). The red line cannot be seen as it coincides exactly with the 1D model black line. This implies that the nuclear evolution along the optical mode dominates the decoherence dynamics.

Close modal

Conventional wisdom indicates that electronic decoherence is fast (∼10 fs), while vibrational decoherence is slow (∼102–103 fs).43 However, in this model, the electronic decoherence rate is identical to the nuclear one, i.e., Pe(t) = PN(t). This is a consequence of the Schmidt theorem9,44 (or the Carlson-Keller theorem45), 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 solvent37 

PN(t)=PNe(t)PNs(t)=expt2τd(e)2expt2τd(s)2,
(16)

where PNe(t) is the purity decay due to intra-molecular electron-nuclear coupling (Ne) and PNs(t) is that due to coupling to the solvent (Ns). Here, τd(e) and τd(s) are the corresponding decoherence times. The purity decay due to Ne interactions is usually faster than that due to Ns interactions (i.e., τd(e)τ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., PNe ≈ 1, and the purity decay due to solvent dominates, i.e., PNPNs. 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 Ne interaction dominates as τd(e)τd(s) and PNPNe 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.

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 ith trajectory represented in the diabatic basis |E0(u0)⟩, |E1(u0)⟩ is

ρe(i)(t)12      12ei0t(E0(u(i)(t))E1(u(i)(t)))dt12ei0t(E1(u(i)(t))E0(u(i)(t)))dt      12,
(17)

where E0(u(i)(t)) and E1(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 E0 and E1 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 Ntraj trajectories,

ρe¯(t)=1Ntrajiρe(i)(t).
(18)

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 = E0(u(t))), first excited PES (V = E1(u(t))), the mean field PES between the two (V = (E0(u(t)) + E1(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.

FIG. 5.

Mixed quantum-classical description of electronic decoherence for the initial state Eq. (13). The plots show the purity of an ensemble of MQC trajectories initially sampled from the ground Wigner distribution of the nuclei. The nuclei evolve in the ground (red) and first excited (blue) PES, their mean (yellow) PES, and a flat potential (magenta). Note that the initial purity decay is independent of the potential on which the classical bath propagates and that it is a good approximation to the exact decoherence for short times. For longer times, the MQC scheme overestimates the overall decoherence time when recurrences are present.

FIG. 5.

Mixed quantum-classical description of electronic decoherence for the initial state Eq. (13). The plots show the purity of an ensemble of MQC trajectories initially sampled from the ground Wigner distribution of the nuclei. The nuclei evolve in the ground (red) and first excited (blue) PES, their mean (yellow) PES, and a flat potential (magenta). Note that the initial purity decay is independent of the potential on which the classical bath propagates and that it is a good approximation to the exact decoherence for short times. For longer times, the MQC scheme overestimates the overall decoherence time when recurrences are present.

Close modal

These simulations numerically validate a recent theoretical analysis37 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.

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 fs22); 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 τ=ωt2π. Here, the frequency of the optical mode is used and ω=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).

FIG. 6.

Dependence of the decoherence dynamics on the nuclear mass in units of M = 1349.14 eV fs22. Bottom panels plot the envelope of the purity decay [defined by the positions of the peaks in P(t)] versus dimensionless time τ=t2π3KM for all masses.

FIG. 6.

Dependence of the decoherence dynamics on the nuclear mass in units of M = 1349.14 eV fs22. Bottom panels plot the envelope of the purity decay [defined by the positions of the peaks in P(t)] versus dimensionless time τ=t2π3KM for all masses.

Close modal
FIG. 7.

Components {an = ⟨ϕnχ1(0)⟩} of the initial nuclear state χ1(0) along the vibrational eigenstates {|ϕn⟩} of the first excited PES. Note that the vibrational energy spectrum becomes more dense as the mass is increased.

FIG. 7.

Components {an = ⟨ϕnχ1(0)⟩} of the initial nuclear state χ1(0) along the vibrational eigenstates {|ϕn⟩} of the first excited PES. Note that the vibrational energy spectrum becomes more dense as the mass is increased.

Close modal

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.

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 ĤRM=μ^ϵ0E(t). Here μ^ is the dipole operator of the molecule, ϵ0 is the polarization of light, and E(t) is its electric field. The laser’s electric field is obtained from

E(t)=dA(t)dt,
(19)

where A(t) is the vector potential defined as

A(t)=A0e(tt0)2tw2sin(ω0t),
(20)

where A0 is the amplitude of the vector potential, tw is the width of the Gaussian envelope centered at t0, 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.

FIG. 8.

Effect of preparation on the decoherence dynamics. The plots show purity during and after laser photoexcitation with a 2 fs laser pulse (top panel) and a 100 fs laser pulse (bottom panel) resonantly tuned to the HOMO → LUMO transition. Excitation via a short pulse leads to decoherence dynamics which is reminiscent to that in Fig. 2. Excitation via a long pulse is dramatically different creating an almost-stationary entangled vibronic state. Insets: Long-time purity decay.

FIG. 8.

Effect of preparation on the decoherence dynamics. The plots show purity during and after laser photoexcitation with a 2 fs laser pulse (top panel) and a 100 fs laser pulse (bottom panel) resonantly tuned to the HOMO → LUMO transition. Excitation via a short pulse leads to decoherence dynamics which is reminiscent to that in Fig. 2. Excitation via a long pulse is dramatically different creating an almost-stationary entangled vibronic state. Insets: Long-time purity decay.

Close modal

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.

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

|Ω(t=0)=12(|E0(u0)+|E9(u0))|χ(t=0),
(21)

which involves electronic states with avoided crossings and therefore is not expected to be pure-dephasing in nature. Here, |E9(u0)⟩ is the ninth excited electronic state determined at the equilibrium nuclear coordinates u0 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).

FIG. 9.

Exact decoherence dynamics for the neutral SSH chain with 4 electrons under different electron-electron interaction strengths U. The chain is initially prepared in a separable superposition state between (a) the ground and first excited electronic states [Eq. (13)] and (b) the ground and the ninth excited state [Eq. (21)]. The PESs for U = 0 and U = 0.8 eV are shown in Fig. 1. Note that changing U only affects the electronic decoherence in case (b) where the dynamics is not pure-dephasing.

FIG. 9.

Exact decoherence dynamics for the neutral SSH chain with 4 electrons under different electron-electron interaction strengths U. The chain is initially prepared in a separable superposition state between (a) the ground and first excited electronic states [Eq. (13)] and (b) the ground and the ninth excited state [Eq. (21)]. The PESs for U = 0 and U = 0.8 eV are shown in Fig. 1. Note that changing U only affects the electronic decoherence in case (b) where the dynamics is not pure-dephasing.

Close modal

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.

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:

  1. 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.

  2. 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.

  3. 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.

  4. 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.

  5. 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.

  6. 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.

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.

1.
A. F.
Izmaylov
and
I.
Franco
,
J. Chem. Theory Comput.
13
,
20
(
2017
).
2.
I.
Franco
,
M.
Shapiro
, and
P.
Brumer
,
J. Chem. Phys.
128
,
244905
(
2008
).
3.
B. G.
Levine
and
T. J.
Martinez
,
Annu. Rev. Phys. Chem.
58
,
613
(
2007
).
4.
Y. C.
Cheng
and
G. R.
Fleming
,
Annu. Rev. Phys. Chem.
60
,
241
(
2009
).
5.
S. H.
Choi
,
C.
Risko
,
M. C. R.
Delgado
,
B.
Kim
,
J. L.
Brédas
, and
C. D.
Frisbie
,
J. Am. Chem. Soc.
132
,
4358
(
2010
).
7.
J. E.
Subotnik
and
N.
Shenvi
,
J. Chem. Phys.
134
,
244114
(
2011
).
8.
M.
Shapiro
and
P.
Brumer
,
Quantum Control of Molecular Processes
(
Wiley-VCH
,
2012
).
9.
M. A.
Nielsen
and
I. L.
Chuang
,
Quantum Computation and Quantum Information
(
Cambridge University Press
,
2011
).
10.
L. A.
Pachon
and
P.
Brumer
,
Phys. Chem. Chem. Phys.
14
,
10094
(
2012
).
11.
G.
Katz
,
M. A.
Ratner
, and
R.
Kosloff
,
Phys. Rev. Lett.
98
,
203006
(
2007
).
12.
P.
Huo
and
D. F.
Coker
,
J. Phys. Chem. Lett.
2
,
825
(
2011
).
13.
M. J.
Bedard-Hearn
,
R. E.
Larsen
, and
B. J.
Schwartz
,
J. Chem. Phys.
123
,
234106
(
2005
).
14.
G. A.
Fiete
and
E. J.
Heller
,
Phys. Rev. A
68
,
022112
(
2003
).
15.
I.
Franco
and
P.
Brumer
,
J. Chem. Phys.
136
,
144501
(
2012
).
16.
I.
Franco
,
A.
Rubio
, and
P.
Brumer
,
New J. Phys.
15
,
043004
(
2013
).
17.
S.
Shim
,
P.
Rebentrost
,
S.
Valleau
, and
A.
Aspuru-Guzik
,
Biophys. J.
102
,
649
(
2012
).
18.
M. I.
Mallus
,
M.
Schallwig
, and
U.
Kleinekathöfer
,
J. Phys. Chem. B
121
,
6471
(
2017
).
19.
A. V.
Akimov
and
O. V.
Prezhdo
,
J. Chem. Theory Comput.
10
,
789
(
2014
).
21.
A. J.
Heeger
,
S.
Kivelson
,
J. R.
Schrieffer
, and
W. P.
Su
,
Rev. Mod. Phys.
60
,
781
(
1988
).
22.
A.
Kar
,
L.
Chen
, and
I.
Franco
,
J. Phys. Chem. Lett.
7
,
1616
(
2016
).
23.
J.
Hubbard
,
Proc. R. Soc. A
276
,
238
(
1963
).
24.
P.
Jordan
and
E.
Wigner
,
Z. Phys.
47
,
631
(
1928
).
25.
A.
Altland
and
B. D.
Simons
,
Condensed Matter Field Theory
(
Cambridge University Press
,
2010
).
26.
J.
Flick
, “
Exact nonadiabatic many-body dynamics: Electron-phonon coupling in photoelectron spectroscopy and light-matter interactions in quantum electrodynamical density-functional theory
,” Ph.D. thesis,
Humboldt-Universität zu Berlin
,
2016
.
27.
D. T.
Colbert
and
W. H.
Miller
,
J. Chem. Phys.
96
,
1982
(
1992
).
28.
J.
Crank
and
P.
Nicolson
,
Math. Proc. Cambridge Philos. Soc.
43
,
50
(
1947
).
29.
O. V.
Prezhdo
and
P. J.
Rossky
,
Phys. Rev. Lett.
81
,
5294
(
1998
).
30.
C.
Arnold
,
O.
Vendrell
, and
R.
Santra
,
Phys. Rev. A
95
,
033425
(
2017
).
31.
M.
Vacher
,
M. J.
Bearpark
,
M. A.
Robb
, and
J. P.
Malhado
,
Phys. Rev. Lett.
118
,
083001
(
2017
).
32.
J.
Schulze
,
M. F.
Shibl
,
M. J.
Al-Marri
, and
O.
Kühn
,
J. Chem. Phys.
144
,
185101
(
2016
).
33.
B.
Gu
and
I.
Franco
,
J. Phys. Chem. Lett.
9
,
773
(
2018
).
34.
J. C.
Light
and
T.
Carrington
, “
Discrete-variable representations and their utilization
,” in
Advances in Chemical Physics
(
John Wiley & Sons, Inc.
,
2007
), pp.
263
310
.
35.
G. A.
Baker
and
P.
Graves-Morris
,
Padé Approximants
(
Cambridge University Press
,
1996
).
36.
H. A.
van der Vorst
,
SIAM J. Sci. Comput.
13
,
631
(
1992
).
37.
B.
Gu
and
I.
Franco
,
J. Phys. Chem. Lett.
8
,
4289
(
2017
).
39.
M.
Gühr
,
H.
Ibrahim
, and
N.
Schwentner
,
Phys. Chem. Chem. Phys.
6
,
5353
(
2004
).
40.
J. A.
Yeazell
,
M.
Mallalieu
, and
C. R.
Stroud
,
Phys. Rev. Lett.
64
,
2007
(
1990
).
41.
P. M.
Felker
and
A. H.
Zewail
,
Phys. Rev. Lett.
53
,
501
(
1984
).
42.
M.
Schlosshauer
,
Phys. Rev. A
72
,
012109
(
2005
).
43.
A.
Nitzan
,
Chemical Dynamics in Condensed Phases: Relaxation, Transfer, and Reactions in Condensed Molecular Systems
(
OUP Oxford
,
2013
).
44.
E.
Schmidt
,
Math. Ann.
63
,
433
(
1907
).
45.
B. C.
Carlson
and
J. M.
Keller
,
Phys. Rev.
121
,
659
(
1961
).
46.
J. C.
Tully
,
Faraday Discuss.
110
,
407
(
1998
).
47.
R.
Kapral
and
G.
Ciccotti
,
J. Chem. Phys.
110
,
8919
(
1999
).
48.
J. P.
Paz
,
S.
Habib
, and
W. H.
Zurek
,
Phys. Rev. D
47
,
488
(
1993
).
49.
V. S.
Batista
and
P.
Brumer
,
Phys. Rev. Lett.
89
,
143201
(
2002
).