Generalized quantum master equations can improve the accuracy of semiclassical predictions of multitime correlation functions

Multitime correlation functions, $C(t1,t2,…,tj)$⁠, are fundamental dynamical objects that characterize chemical and physical systems. When a system is subject to external stimuli, its response can be formulated as a multitime correlation that relates directly to experimental observables, as in multidimensional spectroscopy¹ and quantum control protocols.^2–4 Through this formalism, the energetic structure and dynamics of all manner of physical systems become available, from nanocrystals to extended solids and from gas phase molecules to solvated chromophores. Beyond materials science, multitime correlation functions also reveal the extent of quantum information scrambling in model systems,⁵ chemical reactions,^6,7 black holes,⁸ and many-body localized systems⁹ via out-of-time-ordered correlators. Hence, developing accurate and efficient theoretical and simulation frameworks to predict these objects is of paramount importance in chemistry, materials science, condensed matter physics, and astrophysics. Mathematically, while numerically exact solvers can predict the dynamics of few-state systems in contact with Gaussian heat baths,^10–17 going beyond effective harmonic environments and accounting for non-Gaussian, strong coupling to anharmonic motions remains a fundamental challenge. Such non-Gaussian interactions abound, emerging in charge transfer within structured environments,^18–21 optical spectroscopy of liquids,^22,23 interacting spin systems,^24,25 shot noise in measurements of current,^26,27 and radiation pressure.²⁸ Hence, a dynamical method should ideally capture faithful representations of fermion- and spin-nuclear couplings at the ab initio level. Path integrals and semiclassical methods^29,30 can, in principle, achieve this. However, in practice, these methods must balance accuracy and efficiency.³¹ Indeed, some techniques are accurate but prohibitively expensive, while others are largely inaccurate but (somewhat) affordable. It is thus crucial to develop reliable tools to reduce the computational cost and improve the accuracy of such schemes.

Over the past decade, the generalized quantum master equation (GQME) framework^32–34 has proven transformative in this regard for 1-time correlation functions: GQMEs can both drastically reduce the cost of path integral, semiclassical, and classical simulations—sometimes by multiple orders of magnitude—and also boost the accuracy of even starkly incorrect semiclassical results to within graphical accuracy of benchmark calculations. Such improvements have been shown for the 1-time nonequilibrium and equilibrium dynamics of Gaussian,^35–44 fermionic,^45–48 and atomisitic⁴⁹ models of charge and energy transfer in the condensed phase, and even for many-level photosynthetic complexes.^50,51 The GQME can even be constructed using mean-field theory (MFT) simulations,^{38,40,42,43,50,52} arguably the least accurate semiclassical method.⁵³ 

The GQME exactly rewrites a high-dimensional dynamics problem as a low-dimensional equation of motion of $C(t)$⁠, which encodes the correlation of the variables of interest. At the heart of the GQME is the choice of the projection operator,^54,55 which accounts for the initial preparation of the chemical system (i.e., equilibrium vs nonequilibrium) and selects a few variables to track dynamically. The projection operator also defines the memory kernel, $K(t)$⁠, which encapsulates the effect of all dynamical variables not being explicitly tracked and places an upper limit on how much simulation time is required to predict the dynamics of $C(t)$⁠. The efficiency gains arise from adequately choosing a projection operator that minimizes the memory kernel lifetime. What leads to the accuracy gains is less clear. While early work posited that the accuracy improvements arise from the short-time accuracy combined with short-lived memory kernels, later work⁴⁰ showed one could obtain efficiency gains even in memory kernels that are as long-lived as $C(t)$⁠. In fact, Ref. 40 proposed that the improved accuracy in 1-time nonequilibrium averages and correlation functions arises from the exact semiclassical sampling of bath operators at t = 0, which corrects the memory kernel through its self-consistent extraction from bath-augmented auxiliary kernels.

We have recently shown that one can formulate GQMEs that drastically reduce the cost of multitime simulations,⁵⁶ but is it also possible to use these to improve semiclassical multitime simulations? After all, recent advances in semiclassical theory have enabled the prediction of multitime correlation functions,^57–65 and the GQME could offer the ability to leverage this technology at a lower cost and greater accuracy. Yet, this is not a foregone conclusion. After all, our presumed source of improved accuracy, that is, the exact semiclassical sampling of bath correlations, is only strictly true at t = 0. This raises important and revealing questions:

1. Does a bath-sampling GQME based on approximate semiclassical dynamics allow us to improve accuracy across external temporal interactions?

2. How does additional bath sampling at later times impact the accuracy of the GQME predictions?

To address these questions, we first describe how one must formulate 1-time semiclassical GQMEs to improve efficiency and accuracy. We then introduce our extension of this method to 2-time simulations. By treating the simplest multitime case, we can proceed using fully continuous, rather than discretized, equations. For problems with additional intermediate-time interactions, the theory is trivially extendable, but by remaining at a single measurement, we preempt challenges of time and storage cost that add an additional layer of complexity, as in Ref. 56. We employ a spin-boson model with weak system–bath coupling and a fast decorrelating bath to compare our MFT and GQME results with an exact HEOM calculation. We refer the reader to  Appendix A for computational details.

The GQME shifts attention from the correlation function $C(t)$ to the memory kernel $K(t)$⁠. That is, as long as one knows $K(t)$⁠, it is easy to integrate Eq. (2) to obtain the dynamics of $C(t)$⁠. In general, this memory kernel has a shorter lifetime, $τK$⁠, than the correlation function’s equilibration time, t_eq, setting an upper limit on the amount of information required to characterize the evolution of $C(t)$⁠. That is, $τK<teq$⁠, where $K(t>τK)=0$⁠.

Alternatively, one can use the tensor transfer method,⁶⁶ which is formally equivalent. Thus, all relevant quantities can be extracted from reference dynamics for $C(t)$⁠.

For the GQME framework to improve the accuracy of 1-time semiclassical and quantum-classical dynamics, one must avoid fully decomposing the auxiliary kernels into $C(t)$ and its time derivatives. Exceptions arise when using methods whose approximate propagator, $eapp−iLt$⁠, commutes with the Liouvillian, $L$⁠, or that preserves the canonical distribution when computing equilibrium correlation functions. In such cases, the GQME framework cannot lead to any improvement in accuracy,⁴¹ but can still be used to improve computational efficiency. That is, the quality of the dynamics obtained is the same as that of the original semiclassical method used to parameterize it. However, since these exceptions do not arise for most methods in the semiclassical hierarchy^67–77 or that arise from the quantum-classical Liouville equation,^78–88 one may expect a suitable formulation of the semiclassical GQME to lead to improved dynamics.

In anticipation of the multitime derivation, and to contextualize the extent of GQME-based improvement that one may expect in the multitime response, we first illustrate how the semiclassical GQME can improve the accuracy of the MFT electronic dynamics in a parameter regime corresponding to a biased, underdamped spin connected to a fast decorrelating bath. Figure 1 compares the performance of the MFT against benchmark HEOM calculations and illustrates the improvement one obtains from the GQME sampling bath correlations during the MFT simulation. The GQME brings the approximate MFT dynamics to within the graphical accuracy of the exact HEOM. We have chosen this parameter regime as it places particular strain on MFT: similar to many other semiclassical methods, MFT violates detailed balance, a flaw that becomes most evident in biased systems, and the classical approximation for the bath deteriorates the most for fast, high-frequency baths. While a lower temperature would have rendered the classical treatment of the bath in MFT more questionable and heightened the disagreement between MFT and HEOM, anticipating a nearly prohibitive cost of converging the multitime HEOM benchmark, we have adopted an intermediate temperature that allows faster convergence with the number of Matsubara frequencies.

These data demonstrate why the argument claiming that accuracy improvement arises from the short-time accuracy of MFT (leading to memory kernels that are accurate over their short lifetimes) is incomplete. Indeed, the GQME kernel is more accurate than the MFT equivalent even before the first minimum in the dynamics at very early times.

Semiclassical multitime GQME: We now consider the case where the system experiences an impulsive external field at time t₁, which we choose to be the σ_x Pauli matrix. This corresponds to a photon-induced spin-flip of a qubit in quantum information protocols.² 

In general, time evolution causes the electronic system and bath to become entangled, making it impossible to write the full density matrix as a simple product of a system and bath component. Because this density matrix differs from the form in the projection operator—that is, multiplicatively separable with the bath always at equilibrium with the ground electronic state—introducing an interaction at intermediate times prevents the resulting multitime correlation function from being decomposable into products of 1-time correlation functions.⁹¹ In the mean-field method, the time-evolved density matrix of the entire system remains multiplicatively separable at all times. However, upon interaction with a photon, the system component of the density matrix collapses into a mixed state that cannot be separated into a simple product of a bath operator and a pure subsystem state.⁹² Thus, even with the simplifying approximation in the Ehrenfest method that neglects global correlation, the multitime equation does not simplify to the 1-time GQME.

We now turn to our simulation results. First, we compute the benchmark 2-time correlation function (Fig. 2) using the numerically exact HEOM over a 10 × 10 Δ⁻¹ grid of points with a 0.01 Δ⁻¹ resolution. For the GQME, we do not need to directly simulate the entire grid, which is one of the central benefits of the method. Instead, we demarcate the first 2 × 2 Δ⁻¹ of the grid with black boxes and first simulate this region with MFT. We chose this value because our 1-time kernel lifetime is 0.7 Δ⁻¹ (see Fig. 1), and our previous work⁵⁶ suggests the 2-time kernel goes to zero at least as fast as the 1-time kernel. A 2 Δ⁻¹ range on both axes should be more than sufficient to fully parameterize the 2-time kernel, and thus generate the GQME dynamics for the full 10 × 10 Δ⁻¹ grid in Fig. 2. In this and all the following figures, the 16 panels represent the matrix elements of the displayed quantity, with the i = j = 0 element in the top left.

To quantify the improvement that the multitime semiclassical GQME can afford, we hone in on the difference between the full propagator and the 2-time correlation function constructed using only 1-time quantities, that is, the second term of Eq. (11). Figure 3 (left) shows the benchmark of this quantity obtained with HEOM. For these exact data, we employ a continuous-time algorithm (see  Appendix A of Ref. 56) to extract the 2-time kernel, shown in Fig. 3 (right) from the 2-time term. In this regime, the 2-time kernel has only two significant elements that are confined to small triangles of $∼0.5$ Δ⁻¹ base. This confirms our prediction that 2 Δ⁻¹ is an adequate simulation time for the MFT trajectories, to which we now turn.

Figure 4 presents both MFT and semiclassical GQME equivalents to the panels of Fig. 3. We employ the semiclassical GQME-improved version of the 1-time quantities presented in Fig. 1, which are effectively identical to the HEOM in this case. This enables us to isolate all errors in the 2-time term from the 2-time kernel. The top row of Fig. 4 shows that MFT is completely inaccurate, even at the shortest times. In particular, the 2-time term of Fig. 4(a) has no clear structure and is almost an order of magnitude too large. This means that, in this regime, discarding the 2-time term completely and just using 1-time improvements would be preferable to the result of the 2-time calculation. We construct the associated 2-time kernel in Fig. 4(b) from derivatives of the 2-time term and, therefore, obtain an object that is also of the wrong form and magnitude.

Interestingly, if one rotates only the first $QL$ term with $LSB$ when obtaining Eq. (13), that is, performing bath sampling only at t = 0, then the improvement is lost. Previous work⁴⁰ found that when replacing the remaining $Q$ with $1−P$⁠, the difference between the two terms was slow to converge as a function of the number of trajectories. Hence, we ascribe this observation to the increased difficulty of converging the small difference between two large terms subject to statistical noise from finite sampling.

Currently, the MFT multitime protocol is computationally expensive: our 2-time MFT simulations average over 10 000 trajectories, which has the same cost as performing 100 000 000 1-time trajectories. To obtain the 2-time analog of the 100 000 that satisfactorily converges the 1-time data in Fig. 1 would take four orders of magnitude longer. While converging the data over the entire time domain is feasible for the spin-boson model as semiclassical calculations are trivially parallelizable, performing this many simulations with an atomistic bath would become difficult as force calculations constitute the greatest expense in molecular dynamics. For this reason, we demonstrate that it is possible to obtain improved accuracy and efficiency through our semiclassical GQME even without fully converging the long-term data. Indeed, while at this level of convergence the 2-time kernel of Fig. 4(d) is still noisy compared to the HEOM result in Fig. 3 (right), the GQME 2-time term in Fig. 4(c) obtains semiquantitative agreement in form and magnitude. In addition, because it is a trajectory-based method, earlier time data is more converged, and we only need convergence in the region used in the memory kernel construction, that is, the smaller boxes in Fig. 4 (left). To be quantitative about the degree to which the GQME kernel is underconverged, in Fig. 5, we show the difference between GQME and HEOM 2-time terms on these two timescales.

The fact that there is little to no error in the 0.5 Δ⁻¹ region of Fig. 5 (right) directly represents an improvement in accuracy: while the 2-time term is not maximal in this region, one can see from the (1, 3) and (2, 0) elements of Fig. 3 (left) that it is non-zero and not trivial to obtain. The larger errors at later times are due to an accumulation of error from the underconverged tails of $Kσx(t2,t1)$⁠, consistent with the observation that semiclassical methods require increasingly larger numbers of trajectories to converge correlation functions evaluated at longer times. Since the size of the 2-time term itself is the same as the error, careful comparison shows this is mostly due to slightly shifted frequencies (relative positions of regions of positive and negative signal) rather than their intensities. We can therefore have some confidence that—since the magnitudes are correct, the frequencies are largely accurate, and the GQME 2-time kernel of Fig. 4(c) clearly goes to zero—it is reasonable to evolve the GQME to the full 10 × 10 Δ⁻¹ grid of Fig. 2. Using just the data in the 0.5 × 0.5 Δ⁻¹ region and the 1-time improvement (which is of insignificant cost in comparison), we obtain the GQME 2-time propagator and display it in Fig. 6.

The agreement between the HEOM and semiclassical GQME results for the multitime correlation function $Cσx(t2,t1)$ in Figs. 6 and 2 is excellent, even at later times. Compare, for example, the small red dots in the white “wells” of the (0, 0) element of both figures; these features represent, for example, the extent to which a quantum control sequence could be used to manipulate the coherence properties of this spin system. In contrast, Fig. 7 shows the “full MFT” multitime correlation function. To generate this figure, we do not perform the 2-time MFT simulation over the 10 × 10 Δ⁻¹ grid, as converging it with the same resolution becomes prohibitively expensive. Instead, we employ both 1-time and 2-time GQME formulas from our previous multitime GQME work,⁵⁶ which exclusively improve the efficiency but not the accuracy of semiclassics. Features such as these small dots are entirely missing in the 2-time MFT result, and the qualitative structure of many elements is incorrect [e.g., the (2, 3) element has a pattern of wells where it should be diagonal stripes], showing how critical it is to sample bath information when employing a GQME with MFT reference dynamics. Also noteworthy is that our parameter regime corresponds to weak electronic-nuclear coupling (small reorganization energy), which leads our 2-time memory kernel to offer a comparatively small correction. However, many experimental systems display intermediate to strong electronic-nuclear coupling, such as light-harvesting complexes,^93,94 and small polaron-forming materials such as polymers^95,96 and transition metal oxides.^97,98 In such cases, a 2-time kernel that is five times too large would lead to catastrophically incorrect 2-time propagators even with exact 1-time dynamics.

Our weak-coupling parameter regime is one where time convolutionless second-order perturbation (TCL2) theory^99–101 is expected to perform well for 1-time correlation functions, that is, where $max[2λ/βωc2,2λ/πωc]<1$⁠.¹⁰² There is, however, an additional Markovianity requirement for the accuracy of multitime TCL2. In particular, multitime TCL2 is consistent with the quantum regression theorem,¹⁰³ so it cannot capture non-Markovian correlation¹⁰⁴ across intermediate-time interactions.^105,106 Since our semiclassical GQME correctly predicts correlation across intermediate-time interactions via the multitime kernel, our approach offers a significant advance over standard multitime TCL2 and provides a means to improve the accuracy of multitime semiclassical techniques.

In conclusion, we showed a representative example of when the MFT method cannot quantitatively or even qualitatively capture the 2-time correction to a quantum multitime correlation function reporting on light-induced interactions in the spin-boson model. We then demonstrated that our semiclassical multitime GQME dramatically improved the accuracy of the 2-time kernel by employing classical sampling of bath interactions obtained from the same qualitatively and quantitatively incorrect MFT dynamics. Computationally, the semiclassical GQME gave a modest 10/0.7 ≃ 14 times efficiency gain over the 1-time window and impressive accuracy gains that brought the semiclassical MFT into agreement with the HEOM benchmark, consistent with previous work. For the 2-time problem, the 10/0.5 = 20 times improvement was squared, meaning that to produce the figures in this paper, in which the dynamics have not yet even reached equilibrium, our multitime semiclassical GQME reduced the cost of the simulation by 400 times. Just as impressively, extending the multitime simulation to equilibrium incurs only a trivial computational cost compared with the effort to do this with either HEOM or semiclassical theory alone. Our illustrative example also recapitulates⁵⁶ a general principle that since the higher-order memory kernels have lifetimes comparable to (and even along one axis, shorter than) the 1-time kernel, the multitime GQME obtains higher savings as the number of time indices in the correlation function increases. In upcoming work, we will demonstrate that a discrete-time analog of Eq. (4) can be united with the method introduced in Ref. 56 to affordably treat 3-time correlation functions. This offers a promising approach to utilizing semiclassical dynamics^107,108 to efficiently obtain accurate 2D spectra for complex systems.

This work was partially supported by an Early Career Award in CPIMS program in the Chemical Sciences, Geosciences, and Biosciences Division of the Office of Basic Energy Sciences of the U.S. Department of Energy under Award No. DE-SC0024154. We also acknowledge start-up funds from the University of Colorado Boulder.