An accurate and efficient Ehrenfest dynamics approach for calculating linear and nonlinear electronic spectra

Nonlinear optical spectroscopy is a powerful tool for probing the structure and dynamics of chemical systems,^1–3 and it typically provides significantly more information than its linear counterparts. For example, 2D electronic spectroscopy (2DES)^4,5 is frequently used to elucidate inter-chromophoric couplings, energy transfer and relaxation processes, and environmental effects in systems such as photosynthetic molecular aggregates with femtosecond time resolution.^6–9 However, interpreting these spectra in terms of the specific molecular structures and motions that give rise to the electronic states and their relaxation pathways often requires the assistance of simulations to disentangle spectral features such as short-time coherent oscillations and long-time population decay pathways. Accurately and efficiently simulating 2DES signals from atomistic simulations and uncovering how they arise from the underlying quantum dynamics of the nuclear and electronic states remains a significant challenge.

Simulating nonlinear optical spectra requires obtaining higher-order response functions, and hence, the number of methods that have successfully been applied to generate them is significantly more limited than for the less computationally demanding linear spectra. Exact methods such as the hierarchical equations of motion (HEOM)^10–13 and multiconfiguration time-dependent Hartree (MCTDH)^14,15 provide important insights and benchmarks for approximate methods. However, the desire to treat large condensed-phase systems containing many electronic states with a fully atomistic treatment of the nuclear motions and even on-the-fly evaluation of the electronic surfaces has spurred the recent development of approximate dynamics methods. When approximate methods are used, 2DES has been shown to provide a much stricter test of the quantum dynamics method employed than linear electronic spectroscopy. For example, when using Redfield theory,^16,17 it has been shown that in parameter regimes where the linear electronic spectrum can be quantitatively captured, the 2DES spectrum significantly deviates from the exact results,¹⁸ although this can be somewhat alleviated by freezing the low-frequency bath modes.^18,19

Trajectory-based methods provide a particularly appealing approach to simulate 2DES since they are typically compatible with a fully atomistic treatment of the nuclear motions and even on-the-fly evaluation of the electronic surfaces on which the nuclei evolve. Many recent applications of trajectory-based methods to 2DES have been based on semiclassical treatments of the Meyer–Miller–Stock–Thoss^20,21 mapping of the electronic degrees of freedom. In particular, the partially linearized density matrix (PLDM) approach²² and its more recent spin-PLDM variant^23,24 have been shown to produce accurate 2DES spectra in parameter regimes where perturbative methods break down.^25,26 Further work has shown how mapping-based semiclassical methods can be used to simulate 2DES beyond the perturbative treatment of the field–matter interaction.²⁷ The optimized mean trajectory approach has also been recently coupled with trajectories based on the Meyer–Miller Hamiltonian to compute two-dimensional vibrational-electronic spectra.^28,29 Other trajectory-based approaches include the numerical integration of the Schrödinger equation (NISE)^30–32 and the stochastic Liouville equation methods,^33,34 which involve an explicit treatment of the bath degrees of freedom but do not account for the back-reaction of the electronic degrees of freedom on the bath, leading to inaccuracies in linear and nonlinear spectra. To correct this, NISE has been combined with the fewest switches surface hopping approach, allowing it to incorporate the quantum back-reaction on the bath and thus improve its accuracy.³⁵ 

Here, we present an accurate method to simulate 2DES using Ehrenfest dynamics and contrast it with a previously reported method detailed in Ref. 36. By applying our method, we show that significantly more accurate 2DES spectra can be obtained from Ehrenfest dynamics if one uses an appropriate initialization of coherence and mixed states. In particular, we express all initial conditions as sums of pure states, which have important implications for the accuracy of the resulting dynamics and allow the simulation to be done unambiguously in either the wavefunction or density matrix formulations of Ehrenfest theory. Using a judicious splitting of density matrices into pure states, we are able to keep the computational cost under control by limiting the number of pure states that need to be time-propagated. We show that when formulated in this way, Ehrenfest theory can closely reproduce the HEOM result and that this result is not achieved when initialization is not performed from pure states.

For a dynamical method to be well-defined, it is desirable for it to yield unique results based on its formulation. Previous work has demonstrated that the wavefunction (i.e., pure state) formulation of Ehrenfest offers such uniqueness.³⁷ Hence, here we employ a pure-state decomposition of the off-diagonal initial conditions required for the calculation of spectroscopic response functions. This allows the time-dependent self-consistent mean field approximation at the heart of the Ehrenfest method to evolve physically well-defined initial states, allowing one to construct the time-dependent statistical average from their evolution. In particular, we expand any given initial density matrix ρ₀ as a sum of pure states,

For an arbitrary initial density matrix ρ₀ = |a⟩⟨b|, the pure state decomposition can be accomplished as follows:

where

are pure states. We will refer to this approach as pure state Ehrenfest.

Electronic spectroscopy typically involves a ground state initial condition, |g⟩⟨g|, upon which a dipole operator μ is applied n times, with n determined by the type of response function being computed. In our implementation of pure state Ehrenfest, the initial electronic condition, ρ₀ = |g⟩⟨g| [see Eq. (17)], is maintained as a pair of vectors (|g⟩ and (⟨g|)^†) whose outer product yields the desired initial electronic density matrix, |g⟩⟨g|. The application of μ from the left (right) is then accomplished via a dot product between the μ matrix and the ket (bra) vector, resulting in a new pair of vectors. For example, applying the dipole operator from the right would be accomplished as follows:

where ⟨x| = ⟨ g|⋅ μ. The resultant |g⟩⟨x| can then be split into four pure states as described in Eq. (2). Since |g⟩⟨x| is still maintained as the vector pair |g⟩ and (⟨x|)^†, the pure state components ρ_aa, ρ_bb, $ρab+$⁠, and $ρab−$ are readily obtained by replacing |a⟩ and |b⟩ in Eq. (3) with |g⟩ and (⟨x|)^† = |x⟩, respectively. These pure states are then normalized, and the coefficients are saved so as to multiply them back in when computing the final response function. This ensures that the time propagation always begins from a density matrix whose norm equals unity, in line with the requirements of the density matrix formulation of Ehrenfest dynamics. Additionally, since |g⟩⟨x| is decomposed into four normalized pure states, we circumvent the explosion in terms that would occur were we to create a separate trajectory for each coherence formed when the dipole operator is applied, e.g.,

where j is the set of all possible states that result from a single excitation.

Each subsequent application of the dipole matrix is conducted on all four of the pure state components [ρ_aa(t), ρ_bb(t), $ρab+(t)$⁠, and $ρab−(t)$]. The application of the dipole matrix usually results in the formation of a new set of non-pure states. If we require more time propagation, each resulting non-pure state is further spawned into four pure state components using the procedure described above. Otherwise, we compute the response function by summing up all pure state components, taking care to multiply back the coefficients used to normalize them. This procedure spawns four final branches for the first-order response function and 64 for the third-order response.

We will contrast the pure state Ehrenfest method outlined above with the unmodified alternative, which we refer to as projected Ehrenfest, where all initial conditions |a⟩⟨b| are propagated directly without having been normalized or decomposed into pure states.

We consider a Frenkel exciton model of coupled chromophores where the full matter Hamiltonian is

Here, $Ĥs$ is the system Hamiltonian, $Ĥb$ is the bath Hamiltonian, and $Ĥsb$ is the system–bath coupling. The system Hamiltonian, $Ĥs$⁠, consists of individual chromophores, each with a site energy ϵ_m, that couple electronically via the transfer integrals J_mn,

The bath Hamiltonian, $Ĥb$⁠, consists of independent sets of phonon modes coupled to each chromophore. It is expressed in terms of the phonon mode frequencies ω and their mass-weighted momenta and coordinates $P̂$ and $Q̂$ as

Each chromophore site |m⟩⟨m| is linearly coupled to its local bath of phonon modes such that $Ĥsb$ takes the form

where c_mj is the coupling strength of the jth phonon mode attached to chromophore m. The characteristics of the baths and their effect on the chromophores are specified via the spectral density,

All chromophore sites are assumed to have identical spectral densities (J_m(ω) = J(ω)). Here, we use an Ohmic spectral density with a Lorentzian cutoff (Debye),

where ω_c is the bath cut-off frequency and $λ=(ℏπ)−1∫0∞dωJ(ω)/ω$ is the reorganization energy.

To obtain linear and non-linear spectra, we treat the light–matter interaction perturbatively,⁵ i.e.,

where E(t) is the classical electric field and μ is the total dipole operator,

We consider a system of two coupled chromophores with ϵ₁ = 50 cm⁻¹, ϵ₂ = −50 cm⁻¹, and J₁₂ = 100 cm⁻¹, consistent with that used in previous studies.^18,38 In this system, there are four states, all of which are accessible: the ground state, two singly excited states, and one doubly excited state. The bath coordinates and momenta are initially sampled from the Wigner transform of the Boltzmann distribution on the ground state. We report results for two sets of bath parameters: a slow bath with a relaxation time $ωc−1=300$ fs and a fast bath with a relaxation time $ωc−1$ = 17.7 fs, both at k_BT = 208 cm⁻¹. Each bath was discretized using 300 phonon modes via the method employed in Ref. 39, which is designed to yield the exact reorganization energy regardless of the number of modes used. The reorganization energy λ was set to 50 cm⁻¹. All spectra were calculated with a transition dipole matrix, where μ₁/μ₂ = −5. The slow bath parameter regime was previously investigated via HEOM and Redfield theory and its frozen mode variant in Ref. 18.

Our HEOM and projected Ehrenfest calculations were conducted using the Python package PyRho.⁴⁰ For the HEOM calculations, converged results for the regimes studied here were obtained by employing the high-temperature approximation, i.e., using zero Matsubara frequencies (K = 0), and truncating the auxiliary density matrices at L = 15. Both the HEOM and projected Ehrenfest equations of motion were propagated using a fourth-order Runge–Kutta integrator, while the pure state Ehrenfest method was propagated via the split evolution method, allowing for the use of larger time steps. In the slow bath parameter regime, we used time steps of 10 fs for HEOM and 2 fs for both versions of Ehrenfest. Both versions of Ehrenfest were run at a time step of 2 fs to ensure consistency with projected Ehrenfest, which required the smaller time step to converge the fourth-order Runge–Kutta integrator. In the fast bath parameter regime, we used a time step of 2.5 fs for HEOM and 10 fs for pure state Ehrenfest. The linear spectra were generated using ∼20 000 trajectories of length 400 fs, while the nonlinear spectra used ∼20 000 trajectories of length 300 fs (in t₁ and t₃) for the slow bath regime and 200 fs for the fast bath regime.

Here, we compare the Ehrenfest results for both the linear and nonlinear optical spectra of the Frenkel exciton model outlined in Sec. III to the spectra obtained using the numerically exact HEOM approach. We demonstrate that, while in linear absorption spectra, the differences between the pure state and projected Ehrenfest are subtle, they become much more pronounced for nonlinear spectra.

Linear absorption spectra can be obtained from the Fourier transform of the first-order response function, χ(t),⁵ 

Applying the dipole operator at t = 0 results in optical coherences between the ground state and singly excited states, $ρ̃(0)=μ̂(0)ρ(0)$⁠. This initial condition can be represented either as an unmodified projected state (i.e., |e⟩⟨g|) or as a sum of four pure states as outlined in Eq. (1). Absorption spectra computed via pure state Ehrenfest, projected Ehrenfest, and the numerically exact HEOM for both the fast and slow bath parameter regimes are shown in Fig. 1. In both parameter regimes, the two peaks in the spectrum correspond to the two exciton states in the system. Their energies correspond approximately to the eigenenergies obtained upon diagonalization of the single-exciton section of the electronic Hamiltonian. In the slow bath regime [Fig. 1(a)], where Ehrenfest theory is expected to be accurate, both projected Ehrenfest and pure state Ehrenfest quantitatively reproduce the HEOM result. In the fast bath parameter regime [Fig. 1(b)], despite the qualitative agreement with the HEOM result, there are small inaccuracies for both Ehrenfest methods, with the pure state Ehrenfest method giving a more accurate result than the projected Ehrenfest method. These results can be explained by examining the dynamics arising from the coherence initial condition, ρ_S(0) = |0⟩⟨1|, shown in Figs. 2(a) and 2(b). |0⟩⟨1| is one of the coherence initial conditions (the others being ρ_S(0) = |0⟩⟨2| and their complex conjugates) that is propagated to give rise to the first-order response function, χ(t). In the slow bath regime, there is near quantitative agreement between both versions of Ehrenfest and the HEOM result, while in the fast bath regime, there are more pronounced differences, with the pure state Ehrenfest dynamics matching the HEOM result more closely and the projected Ehrenfest result being underdamped. These differences mirror and point to the source of inaccuracies in the absorption spectra for the fast bath parameter regime. The overall similarity in coherence dynamics across different methods (Fig. 2) also reveals why the linear spectra appear to be relatively insensitive to the initialization scheme used for Ehrenfest dynamics: the dynamics starting from a density matrix corresponding to a coherence with the ground state are not as sensitive to the initialization scheme. The relative insensitivity of the linear absorption spectra to methods that inaccurately treat coherence dynamics has previously been observed for Redfield theory, which was shown to yield accurate linear absorption spectra despite getting incorrect population dynamics,¹⁸ albeit with a different dynamical method rooted in perturbation theory.

Two-dimensional electronic spectra offer a much stricter test than absorption spectroscopy of a method’s accuracy because they require one to correctly capture both the population and coherence dynamics for a wider range of initial conditions.

2DES spectra can be obtained from a Fourier transform over the rephasing (R_rp) and non-rephasing (R_nr) terms of third-order response function under the rotating wave approximation,⁵ 

where

with Φ₁(t₃, t₂, t₁) and Φ₄(t₃, t₂, t₁) yielding the stimulated emission (SE) contributions, Φ₂(t₃, t₂, t₁) and Φ₅(t₃, t₂, t₁) yielding the ground state bleaching (GSB) contributions, and Φ₃(t₃, t₂, t₁) and Φ₆(t₃, t₂, t₁) yielding the excited state absorption (ESA) contributions. These are computed as

Here, $G(t)ρ=e−iHmatt/ℏρeiHmatt/ℏ$ is the Liouville-space propagator. t_i∈{1,2,3} are time intervals between successive light–matter interactions, and

A procedure for obtaining the third-order response function from Ehrenfest dynamics was previously outlined in Ref. 36, which employs ensemble averages over the ground state wavefunction, |g⟩. In this scheme, the dipole operator is applied via a product with individual dipole matrix elements, time propagation is done only over select population states, and the effect of coherences is modeled via an average over two wavefunctions formed from the ket and bra. However, the use of specific dipole matrix elements as well as time-propagation over select population states makes it difficult to generalize that method to arbitrary systems with multiple quantum states. Here, we present a generalizable method that exploits the linearity of density matrices to compute third-order response functions as follows:

1. Apply the first dipole operator to ρ₀ at t = 0 and split the resultant density matrix into four pure states as described in Sec. II. Via the Condon approximation, the bath degrees of freedom are untouched by the dipole operator. Each resultant pure state inherits its own copy of the bath. Propagate states and their corresponding baths independently through t₁.

2. Apply the second dipole operator at t = t₁, this time to all four propagated states. This operation usually results in a new set of non-pure states, which are subsequently split into four of their respective pure states, bringing the total number of branches to 16. Each branch adopts a copy of bath coordinates and momenta from its parent branch, ensuring continuity from t = 0. Propagate the pure states and baths independently through the waiting time, t₂.

3. Apply the third dipole operator at t = t₁ + t₂ to the 16 propagated states. This splits them again into a total of 64 branches, each a pure state. Propagate all branches through t₃.

4. Consolidate all 64 branches at t = t₁ + t₂ + t₃, and apply the final dipole operator to obtain the third-order response function.

Multiple instances of this procedure are run in order to average over an ensemble of initial bath states. From the above, it is clear that the pure state Ehrenfest method requires more sampling than projected Ehrenfest, which does not involve any splitting of initial conditions into pure states. As we will see below, this splitting into pure states is required to obtain accurate 2DES spectra. The procedure can also be trivially parallelized since all trajectories in the scheme run independently. Additionally, as demonstrated in Sec. II, we decompose each initial density matrix |a⟩⟨b| as a sum of only four pure states in the same Liouville space.

An alternative method for obtaining accurate Ehrenfest dynamics from non-pure state density matrices is outlined in our previous work in the Appendixes of Ref. 41, and it involves stochastically sampling an auxiliary wavefunction. For example, the initial condition |0⟩⟨1| would be obtained from

where ϕ is a parameter that is sampled over the range 0–2π. In this stochastic scheme, sampling is required both over values of ϕ and initial positions and momenta of the bath. In the context of a third-order response function, ϕ would require fresh sampling at each light–matter interaction (i.e., at times t = 0, t = t₁, and t = t₁ + t₂) because the auxiliary wavefunction would otherwise return incorrect dynamics. In contrast, the pure state Ehrenfest method detailed above only requires the usual sampling over initial bath conditions, which can be done once at t = 0 since the baths are continuous through t = t₁ + t₂ + t₃.

Figure 3 shows the 2DES spectra obtained from both the pure state and projected forms of Ehrenfest and how they compare to the HEOM spectra. At t₂ = 0, we see two partially resolved diagonal peaks centered at energies −130 and 140 cm⁻¹, which roughly match the eigenenergies of single-exciton states. We also observe two cross-peaks that indicate inter-chromophoric coupling, with their positive amplitudes, suggesting that they arise from SE and GSB contributions involving two different one-exciton states. The cross-peak in the upper diagonal region is weaker at all times because the ESA contributions largely cancel out the SE + GSB contributions. Spectral diffusion is observed at later times (t₂ = 200, 600 fs), and in the HEOM spectra, the lower-energy exciton state evolves to have a higher population than the higher-energy exciton state, consistent with the detailed balance expectation.

The response functions used to obtain the Ehrenfest spectra were truncated (windowed using a step function) as outlined in the supplementary material, Sec. I. From Fig. 3, we observe that pure state Ehrenfest recovers almost quantitative agreement with the HEOM result, with just a slight deterioration as the waiting time t₂ increases. In contrast, projected Ehrenfest yields much worse agreement with the HEOM spectrum. While most of the spectral features from the HEOM spectrum are present at t₂ = 0, the projected Ehrenfest spectrum becomes progressively worse as t₂ increases and bears little resemblance to the HEOM result by the time t₂ = 600 fs. Additionally, it is well known that the Ehrenfest method is unable to capture detailed balance correctly.⁴² Indeed, the t₂ = 600 fs panel in Fig. 3 shows that for pure state Ehrenfest, the intensity of the higher-energy diagonal peak is larger than would be normally expected in comparison to the exact HEOM spectrum, demonstrating that pure state Ehrenfest inherits the breakdown of detailed balance observed in the projected Ehrenfest formulation.

Although one could imagine that the poor performance of the projected Ehrenfest method could be attributed to underconverged sampling of initial conditions and trajectories, in the supplementary material, Sec. V, we demonstrate that increasing the sampling does not remove the incorrect oscillatory features observed in the projected Ehrenfest 2DES spectrum. Instead, we look for explanations of the inaccuracies by examining the contributions to the response function that make up the 2DES spectrum. Figures 4 and 5 illustrate this using plots of Φ₂, with the rest of the contributions to the response function being shown in the supplementary material, Figs. 2–13. Here, we observe that, as expected, pure state Ehrenfest accurately reproduces the HEOM result, while projected Ehrenfest fails to do so and that the discrepancy is much more striking in the imaginary part of the response function. These discrepancies can be further traced to the single-time dynamics obtained from the coherence initial conditions between excited states. Since the Frenkel exciton model is in the global ground state ρ₀ = |0⟩⟨0| at t = 0, such coherence states can only be obtained when the dipole operator is applied multiple times, as is the case when computing the third-order response. As such, the dynamics of coherence initial conditions between excited states (e.g., |1⟩⟨2|) are relevant to nonlinear spectra but not linear spectra. Figure 6 shows both the population and coherence dynamics for |1⟩⟨2|. In both the slow and fast bath parameter regimes, we see from the first two rows that the population dynamics obtained from pure state Ehrenfest more closely matches the HEOM result, especially at longer times. In contrast, projected Ehrenfest result matches the HEOM result until ∼200 fs, after which there are notable deviations. This poor performance of the projected Ehrenfest is also observed to a lesser extent in the coherence dynamics of the |1⟩⟨2| initial state (bottom two rows of Fig. 6). We emphasize that the coherence dynamics presented here are distinct from those involved in the absorption spectra since they do not involve the ground state. The worsening performance of projected Ehrenfest with time in reproducing the dynamics arising from |1⟩⟨2| also explains why the respective 2DES spectrum degrades progressively as the waiting time t₂ increases (Fig. 3). As such, the differences in the performance between projected and pure state Ehrenfest can be attributed to inaccuracies in the projected Ehrenfest dynamics when populations, and coherences to a lesser extent, arise from an initial condition of coherences between two excited states.

Pump–probe spectra, generated by integrating over ω₁, are shown in Fig. 7. The results here mirror those of the 2DES spectrum, where the pure state Ehrenfest method accurately reproduces the HEOM result, albeit with slightly higher peaks at 200 cm⁻¹ for t₂ ≥ 200 fs. In contrast, while the projected Ehrenfest method yields accurate results at t₂ = 0, the accuracy of the pump–probe spectra degrades at t₂ ≥ 200 fs, with the peak around 200 cm⁻¹ losing definition.

Next, we examine how well pure state Ehrenfest performs in a more challenging fast bath regime. Figure 8 shows 2DES spectra for the fast bath parameters computed via HEOM and pure state Ehrenfest. The range of values for t₂ (0–150 fs) is more limited here than in the slow bath regime because we observed virtually no change in the 2DES spectrum beyond t₂ = 200 fs (see the supplementary material Fig. 14). As with the slow bath parameter regime, the pure state Ehrenfest response functions were truncated before being Fourier transformed to produce the 2DES spectra (see the supplementary material, Sec. I). In this more challenging regime, pure state Ehrenfest does not perform as well as it did for the slow bath parameters. This is not surprising, as it aligns with the well-known deficiencies of the Ehrenfest method, which struggles to capture the correct dynamics in the limit of fast nuclei. Nonetheless, it is able to qualitatively capture most of the features present in the HEOM spectrum. This result is corroborated by the pump–probe spectra shown in Fig. 9.

Here, we have introduced the pure state Ehrenfest method to obtain linear and nonlinear spectra. We have shown that exploiting the linearity of the density matrix to express all initial conditions as sums of pure states is essential to obtaining accurate coherence dynamics. Our approach to obtaining pure states following the application of the dipole operator ensures that the number of branches involved in the time propagation does not depend on the number of quantum states, making our approach suitable for systems with large numbers of excited states. While our approach provides modest improvements over projected Ehrenfest in computing linear spectra, it is able to dramatically improve the accuracy with which nonlinear spectra (2DES and pump–probe) can be obtained, especially at longer waiting times. The physical origin of this improvement arises from the ability of our pure state Ehrenfest approach to much more accurately capture the dynamics when the initial condition is a coherence between excited states. These conditions do not occur when calculating the linear spectra (where all the initial density matrices are coherences with the ground state) but play a vital role in obtaining the higher-order response functions needed to capture nonlinear spectroscopy. We have shown that pure state Ehrenfest gives excellent agreement with the exact linear and 2DES spectra in the slow bath regime where the Ehrenfest method is expected to be accurate and have also demonstrated that it still provides good results in a faster bath regime and preserves the qualitative features observed in the exact spectra. These results suggest that the pure state Ehrenfest approach provides a tractable and accurate approach to calculate nonlinear spectra for multichromophoric systems in a wide range of chemical environments.

See the supplementary material for details of the time cutoffs involved in the construction of the 2DES spectra, plots of individual components of the third-order response function for both the slow and fast bath parameter regimes, and a discussion of the 2DES spectra for the fast bath parameter regime at longer waiting times.

This work was funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences (Grant No. DE-SC0020203 to T.E.M.). A.O.A. acknowledges support from the Edward Curtis Franklin Award and the Stanford Diversifying Academia, Recruiting Excellence Fellowship. A.M.C. acknowledges the start-up funds from the University of Colorado, Boulder.

The authors have no conflicts to disclose.

Austin O. Atsango: Conceptualization (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Andrés Montoya-Castillo: Conceptualization (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). Thomas E. Markland: Conceptualization (equal); Funding acquisition (lead); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

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