Approximate but accurate quantum dynamics from the Mori formalism. II. Equilibrium time correlation functions

By encoding a system’s response to a weak perturbation,¹ equilibrium time correlation functions (ECFs) provide direct access to important dynamical quantities such as transport coefficients,^2,3 absorption spectra,⁴ and chemical rate constants.^3,5,6 It is therefore not surprising that the development of accurate and efficient approaches to the calculation of ECFs has been a focus of intense theoretical investigation, spawning a rich array of analytical and computational methods.

Attempts at calculating ECFs face the challenge of striking a balance between computational tractability and accuracy. On the more computationally expensive side lie numerically exact schemes, which are generally restricted to small, idealized models and tend to scale unfavorably, and at worst exponentially, with simulation time.^7–15 Approximate methods, on the other hand, are more scalable to realistic, multidimensional systems but suffer from limited accuracy. These include perturbative approaches,^16–19 analytic continuation,^20–25 quantum mode coupling theory,^26–30 and quasi-^31–35 and semi-classical^36–54 schemes. While each method has its virtues, all suffer from limited applicability, either due to the violation of parameter regime restrictions, the breakdown of uncontrolled approximations, or convergence problems. As a result, despite the availability of many methods for the calculation of ECFs, there remains a clear need for the development of widely applicable, accurate, and computationally efficient approaches.

A successful marriage of computational efficiency and accuracy may be enabled through a judicious use of the Mori formalism. Based on the projection operator technique,^55,56 the Mori approach provides a simple, low-dimensional equation of motion for the ECF, called a generalized quantum master equation (GQME). The reduced dimensionality of the GQME corresponds to that of the space spanned by the observables probed in the ECF, while the influence of the excluded degrees of freedom is encoded in the memory kernel. Calculation of the memory kernel, in turn, is beset by two difficulties: application of the projected propagator, $eiQLt$⁠, and the full dimensionality of the original problem. The former can be sidestepped using a Dyson-type expansion, which leads to a self-consistent equation for the memory kernel that requires only the calculation of projection-free auxiliary kernels.⁵⁷ The latter, however, continues to plague the calculation of the auxiliary kernels. Nevertheless, this approach has been vigorously pursued to obtain nonequilibrium averages in the context of impurity-type models where the auxiliary kernels have been obtained via numerically exact^57–64 or approximate methods.^65–70 These studies have shown that the short lifetime of the memory kernels can lead to dramatic increases in computational efficiency, regardless of the method used to compute the memory kernels.^57–70 Just as importantly, when approximate methods are used, the Mori approach has also provided impressive boosts in accuracy over the bare approximate dynamics.^65–70 Hence, the approach based on the GQME coupled to the self-consistent solution of the memory kernel shows great promise as a means of increasing the efficiency and, when appropriate, accuracy of dynamical methods.

Here we argue that the remarkable boosts in efficiency and accuracy afforded by the GQME approach can be extended to arbitrary systems and dynamical quantities beyond simple nonequilibrium averages. To show the viability of the approach, we specialize the Mori treatment to the symmetrized ECFs for the spin variables of the spin-boson (SB) model.^18,71 In the same spirit as the first paper of this series,⁶⁹ we calculate the auxiliary kernels necessary for the self-consistent solution of the memory kernel via the mean field theory (MFT) or Ehrenfest method and assess the potential benefits of this approach in terms of increases in efficiency and accuracy (we henceforth refer to this framework as the GQME+MFT approach). In this work, we also endeavor to elucidate the dependence of the GQME dynamics on the choice of closure and try to provide further evidence for the claim that the source of the improvement over bare semiclassical dynamics afforded by the GQME framework depends, at least partially, on the exact sampling of distinct initial conditions necessary in the calculation of the auxiliary kernels. It also bears remarking that the Mori approach can be easily generalized to multi-time correlation functions and is applicable to a wide variety of systems. In fact, a major advantage of the Mori formulation is that it can naturally address problems where the system-bath dichotomy is absent, such as spin and fermion lattice models,^72–74 and quantum fluids.^30,75–77

The paper is organized as follows. In Sec. II, we briefly introduce the SB model and the projection operator that allows for the investigation of symmetrized ECFs of the Pauli matrices. Sec. III A compares the dynamics obtained via the GQME+MFT approach to numerically exact results for the symmetrized spin-spin correlation function $Czz(t)=Re⟨σz(0)σz(t)⟩$⁠. Sec. III B explores the dependence of GQME results on the choice of closure and the accuracy of the initial conditions. Sec. IV is devoted to our concluding remarks.

As in Paper I,⁶⁹ we apply the Mori formalism for ECFs to the SB model, which is representative of typical condensed phase systems that exhibit nontrivial decoherence and dissipation patterns.^18,71 We note, however, that the Mori approach is general and may be applied to any Hamiltonian system, including those where the system-bath distinction is absent. This point is of crucial importance for the modeling of quantum ECFs such as those associated with absorption spectra in liquids.⁷⁶ 

The SB Hamiltonian takes the form H = H_S + H_B + H_SB, where

corresponds to the system part of the Hamiltonian. Here, $2ε$ corresponds the energy difference between the two sites, Δ, which is assumed to be static, represents the tunneling matrix element, and $σi$ corresponds to the i^th Pauli matrix.

The bath Hamiltonian consists of independent harmonic oscillators,

where P_k, Q_k, and $ωk$ are the mass-weighted momenta, coordinates, and frequency for the k^th harmonic oscillator, respectively. The system-bath coupling is assumed to be of the form,

where c_k is the coupling constant describing the strength of the interaction between the system and the kth oscillator and $α=±1$⁠. The system-bath interaction is fully characterized by the spectral density,

which encodes the frequency-resolved coupling between the system and the oscillators that compose the bath. The second line in Eq. (4) corresponds to the often used Ohmic form for the spectral density¹⁸ with an exponential cutoff. Here, $ωc$ is the cutoff frequency, which determines the correlation time for the bath at finite temperature. The Kondo parameter, $ξ$⁠, describes the strength of the system-bath coupling and is proportional to the reorganization energy, $λ=ξωc/2=π−1∫0∞dωJ(ω)/ω$⁠, which represents energy dissipated after the system makes a Frank-Condon transition. We also note that the present approach is neither limited to the SB model nor to the Ohmic form of the spectral density.

As established in Paper I,⁶⁹ the Mori (as well as the Nakajima-Zwanzig) approach stems from the generalized Nakajima-Zwanzig-Mori equation of motion for the propagator,

where $P$ is the projection operator, which consists of the dynamical operators whose correlations we seek, and $Q=1−P$ is the complementary projection operator.

For ECFs, the Mori-type projector traditionally consists of an operator B and its time derivative $B˙=iLB$⁠.^29,70,78,79 Instead of taking this approach, we employ a projector that spans the entire Hilbert space of the spin and orthonormalize the components of the projector to satisfy the idempotency requirement,

where $An∈{𝟏S,σx,σy,σz}$⁠, $ρ=e−βH/Tr[e−βH]$ is the canonical density operator, $β=[kBT]−1$ is the inverse thermal energy, and the transformation matrices $T$ and $R$ enforce the orthonormalization of the members of the projector (⁠$R(𝐀|𝐀)T=𝟏$⁠). The inner product is defined so as to recover the symmetrized form of the correlation function,

where $O$ is an operator and A and B are vectors in Liouville space. The form of the projector in Eq. (6) serves a dual purpose. It allows us to obtain any ECF of the spin variables, and it includes all potentially slow degrees of freedom associated with the subsystem that could result in a long-lived memory kernel, even if this does not preclude the possibility of slow contributions to the memory kernel coming from the bath.

Using the components of the projection operator to close Eq. (5) from both sides, we obtain the Mori-type GQME,

where $Cnm(t)=Tr[ρ{An†,Am(t)}]$ is the ECF of all members of A. The memory kernel takes the following form:

Given the difficulties associated with treating the dynamics required by the projected propagator in Eqs. (9) and (10), we follow Shi and Geva⁵⁷ and employ the Dyson decomposition,

to obtain the $Q$-forward (f) and $Q$-backward (b) self-consistent expansions of the memory kernel,

where the normally evolved auxiliary kernels

can be obtained via direct simulation. For more details about the derivation, we refer the reader to Ref. 69.

As discussed in the first paper of this series,⁶⁹ different closures may lead to different results when using approximate dynamics to calculate the auxiliary kernels. For that reason, we also explore the effect of the three additional closures that selectively replace the action of the Liouvillian with time-derivatives. For completeness, these alternative closures are reproduced below. In the first set, we replace the action of the Liouvillian acting on operators that require dynamic sampling with the time-derivative,

The second type replaces the action of the Liouvillian on static operators with the time derivative,

The final type replaces all Liouvillian operators with time derivatives,

In the following, we refer to these closures using the contracted notation cfx or cbx, where $x∈{0,1,2,3}$⁠. In this notation, cf (cb) represents the $Q$-forward (⁠$Q$-backward) closures, which are obtained with the forward (backward) placement of the $Q$ in the projected propagator in the memory kernel in Eq. (9) (Eq. (10)). The number at the end of the notation denotes if and which Liouvillian in the expressions for the auxiliary kernels has been replaced with a numerical time derivative. In other words, cb2 refers to the closure where the memory kernel is extracted using Eq. (14) and the auxiliary kernels are given by Eqs. (21) and (22).

It bears repeating that all closures presented above are permitted by quantum mechanics and, when exact methods are used to calculate the auxiliary kernels, the memory kernels and GQME dynamics produced will be equivalent, regardless of the closure. However, as was shown in the nonequilibrium case,⁶⁹ the memory kernels and GQME dynamics obtained from auxiliary kernels calculated via approximate schemes can differ from each other, depending on the closure. Explicit expressions for the various closures in the context of the SB model can be found in  Appendix B.

As in the first paper of this series,⁶⁹ we employ the Ehrenfest method to calculate the auxiliary kernels necessary for the extraction of the memory kernels. The Ehrenfest approach is a simple quasi-classical scheme where the system (bath) evolves in the mean field of the bath (system). While most implementations of the Ehrenfest method (and other low-level quasi- and semi-classical theories such as the surface hopping^33,34 and linearized semiclassical initial value representation (LSC-IVR) schemes^38–40) use approximate forms for the fully correlated Boltzmann factor,^38,80–82 we employ a numerically exact representation which relies on the path integral approach developed in Ref. 83. The importance of the exact rendering of the canonical density operator will become evident in Sec. III A. Within this path integral framework, it is possible to produce expressions for the canonical distribution that increase in accuracy with the number of path integral slices, N, used. The number of path integral steps taken in the rendering of a given realization of $ρ$ is determined by the number of slices necessary for the convergence of static and dynamic data as a function of N. Details regarding the path integral approach to the canonical density are included in  Appendix A and those regarding the implementation of the Ehrenfest method are included in  Appendix B.

The protocol for the extraction of $K(t)$ can be summarized as follows. Using the approximate auxiliary kernels obtained via the Ehrenfest scheme, we solve Eqs. (13) and (14) iteratively until the relative error R.E. is negligible, i.e., $R.E.<10−10$⁠, where $R.E.=max[abs[Kn+1(t)−Kn(t)]]$ is the maximum absolute difference between two subsequent iterations of the memory kernel. With a converged solution for the memory kernel, we numerically integrate Eq. (8) using a second-order Runge-Kutta procedure, subject to the appropriate initial conditions, namely, $Cii(t=0)=1$⁠. To assess the viability of the GQME+MFT procedure in the equilibrium context, we compare the solution of Eq. (8) using the self-consistently extracted $K(t)$ to numerically exact results for $Czz(t)$ for the SB model with an Ohmic spectral density. We further specify that, in the same vein as previous implementations of the GQME coupled to the self-consistent solution of the memory kernel,^57–70 the finite lifetime of the memory kernel is explicitly included in the GQME evolution algorithm as a cutoff time, $τc$⁠, beyond which the memory kernel is set to zero. In practice, $τc$ is identified as any time point in the stability plateau corresponding to the range of time during which the memory kernel ceases to influence the GQME dynamics. When numerically exact methods are used to obtain the auxiliary kernels, this stability plateau is infinitely long.⁸⁴ Conversely, when approximate methods are used instead, the limited accuracy of the underlying method may cause the stability plateau to be finite or, in extreme cases, nonexistent. In all cases studied in this work, a well-defined stability plateau exists and all cutoff times for the memory kernels are specified (see Sec. III C for a more detailed exploration of how to determine the plateau of stability for difficult cases).

For nonequilibrium averages, the GQME+semiclassics approach has been shown to be uniformly beneficial in reducing the cost of the dynamics via short-lived memory kernels and generally advantageous in correcting approximate dynamics, especially for biased systems coupled weakly to the bath.^65–70 In the following, we show that similar benefits can be reaped in the equilibrium case.

Before continuing to the ECFs, we discuss the behavior of different elements of the memory kernels (see Eqs. (7), (13), and (14)) for two different regimes. Panels (a) and (b) in Fig. 1 correspond to an unbiased (⁠$ε=0$⁠) system coupled to a moderately fast bath (⁠$ωc=2.5$⁠), and panels (c) and (d) to a strongly biased system (⁠$ε=−2$⁠) coupled to a fast bath (⁠$ωc=15$⁠). For both figures, the diagonal elements (panels (a) and (c)) are strongly peaked at t = 0 and quickly decay to approximately zero, except for what appears to be persistent, small amplitude oscillations. The off-diagonal elements for the unbiased realization of the model, shown in panel (b), display long-lived small amplitude oscillations. In contrast, the off-diagonal oscillations of the biased version in panel (d) display a strongly peaked behavior at t = 0 which subsequently decays to zero. We also note that all other elements of the memory kernel in the pre-rotation basis (i.e., before application of the rotation matrices $R$ and $T$ are necessarily zero, as ensured by the transformation matrix $Y$ in Eqs. (B1) and (B2) in  Appendix B. Figs. 2(a) and 3(d) present the GQME+MFT dynamics obtained from the memory kernels in panels (a)–(b) and (c)–(d), respectively. While it may be appealing to interpret the specific memory kernel elements in terms of distinct processes associated with different elements of the two-level system degrees of freedom, it is important to note that at the conclusion of the calculation the basis used to evolve the GQME must be rotated to the orthonormal basis where all kernel elements are non-zero. Finally, it is worth noting that although the memory kernels appear to contain long-lived oscillations both in the original and rotated bases, the GQME dynamics to which they give rise are insensitive to these oscillations, underlining the robustness of the GQME+MFT method.

We now turn to the dynamics obtained using the GQME+MFT approach. As Figs. 2 and 3 illustrate, the advantages afforded by the GQME formalism for nonequilibrium averages are transferable to equilibrium situations. For example, Fig. 2 illustrates that in the weak coupling regime shown in panels (a) and (b), the GQME approach is able to easily reproduce the already accurate dynamics produced via the Ehrenfest method. In these instances, the GQME+MFT scheme also affords boosts in efficiency, since the memory kernel cutoff necessary to recover the appropriate dynamics is $τc=1.0$ for panels (a) and (b). Naturally, the limited applicability of the Ehrenfest method can also lead to dynamics that significantly deviate from the exact results. Indeed, by virtue of its mean-field character, the Ehrenfest scheme is known to fail for cases where the system-bath coupling is large, as shown in panels (c) and (d). Here the Ehrenfest scheme leads to overly fast relaxation. In contrast, the GQME+MFT method is able to produce dynamics that are in notable agreement with the numerically exact results, albeit at a similar computational cost as the direct Ehrenfest calculation for panels (c) and (d), where $τc≥2.5$⁠. It bears remarking that the choice of $τc$ for systems characterized by strong system-bath coupling is somewhat more difficult. For panels (c) and (d), the choice of $τc$ indicates that the dynamics are insensitive for $τc≥2.5$ over the pertinent time scale for the ECF. If, instead, one was interested in going beyond the time shown, the cutoff times would have to be increased until a suitable plateau of stability is found for the desired time range. At worst, a plateau of stability may not exist. For these cases, the GQME+MFT approach is not appropriate. Finally, it is noteworthy that even for unbiased cases, where the Ehrenfest method already performs well for nonequilibrium dynamics, the GQME+semiclassics approach can offer marked improvements in the equilibrium case.

The approximations that underlie the Ehrenfest approach also lead to incorrect dynamics in other parameter regimes. For instance, the classical treatment of the bath implies that the Ehrenfest scheme is most accurate for slow baths where quantum effects are negligible.³² In addition, the mean-field character and classical treatment of the bath in the Ehrenfest framework results in the breaking of detailed balance, a problem that becomes most pronounced in the dynamics of biased systems.^85,86 Cases that lie beyond the region of validity of the Ehrenfest method are also those where the improvements afforded by the GQME+MFT approach are most dramatic, as Fig. 3 illustrates. Given the above considerations, it is not surprising that the Ehrenfest method can only qualitatively capture the relaxation of the ECF, leading to overly fast relaxation and, in the case of the strongly biased system in panel (d), leading to the wrong oscillation frequency. In contrast, the dynamics produced by the GQME+MFT approach are in remarkable agreement with the numerically exact results. Importantly, as the figure indicates, the GQME dynamics are robust with respect to the choice of cutoff time, $τc$⁠. Indeed, as the results demonstrate, it is possible to change $τc$ by an order of magnitude without affecting the resulting dynamics. Hence, as in the nonequilibrium case, the present method represents an important tool that has the potential to correct the semiclassical dynamics of systems characterized by nonzero bias or coupling to fast baths.

One question continues to emerge from the continued success of the GQME+semiclassics formalism: where do such improvements in accuracy over bare approximate dynamics come from? In the first paper of this series,⁶⁹ we conjectured that the source was, at least to some extent, the exact sampling of distinct initial conditions required for the calculation of the auxiliary kernels. As our added emphasis indicates, this improvement appears to rely on two different factors. The latter, in the nonequilibrium case, was identified as the “correction terms” arising from the Wigner transform of the product of the bath part of the system-bath interaction with the bath distribution, which took the form $ζ𝑛𝑜𝑛𝑒𝑞W=−α∑kckPktanh(βωk/2)/ωk$ (see Eq. (D6) in Appendix D of Ref. 69). Our treatment of the equilibrium problem produces the analogous bath operator given by Eq. (A19) in  Appendix B. The inclusion of the correlation functions that require the sampling of this operator in the auxiliary kernels (see Eqs. (B7)–(B9) and (B20)) is indeed necessary to obtain the improved dynamics shown in Figs. 2 and 3. The importance of the first factor, i.e., the exact nature of the sampling at t = 0, in the context of thermal equilibrium can only be confirmed by using representations of the canonical density that can be made arbitrarily accurate. For realizations of the SB model where the system and bath are weakly coupled and at high temperatures, the canonical distribution can be captured by a simple factorization approximation implemented earlier in the literature.³⁸ The two simplest approximations consist of rewriting the Boltzmann factor as simple products of system and bath operators, e.g., $e−βH≈e−βHSe−βHB$ and $e−βH≈e−β(HB+HSB)/2e−βHSe−β(HB+HSB)/2$⁠, which we refer to in our path integral notation as containing N = 0 and N = 1 path integral slices, respectively. For cases where the system-bath coupling is strong (large $ξ$⁠) or temperatures are low (large $β$⁠), the Boltzmann factor can no longer be captured by such crude approximations. In these cases, we adopt the path integral scheme developed in Ref. 83 to obtain highly accurate expressions for the canonical density, $ρ$⁠.

In Fig. 4, two sets of dynamics are presented. Panels (a) and (c) correspond to the direct Ehrenfest treatment of C_zz(t) for two different realizations of the SB model using an increasing number, N, of path integral slices. As is evident in panel (a), which corresponds to a case with strong system-bath coupling (⁠$ξ=2.55$⁠) and moderate to low temperatures (⁠$β=1.6$⁠), a large number of path integral steps (N = 6) are required to achieve convergence. Conversely, the high temperature (⁠$β=0.806$⁠) and weak coupling (⁠$ξ=0.2$⁠) case presented in panel (c) only requires two path integral steps (N = 2) to achieve convergence. More importantly, panels (a) and (c) illustrate that the accurate representation of the Boltzmann factor can significantly improve the dynamics produced via the bare Ehrenfest method. Panels (b) and (d) present the GQME+MFT dynamics obtained using the different approximations for $ρ$⁠. As the figure suggests, the accurate rendering of the canonical distribution is critical in obtaining improved accuracy via the GQME framework. Moreover, the sensitive dependence of the improvements afforded by the GQME formalism on the accurate rendering of $ρ$ lends credence to our conjecture⁶⁹ that the exact sampling of distinct initial conditions required by the auxiliary kernels is largely responsible for the improvements in accuracy afforded by the GQME framework.

As was demonstrated in the first paper of this series,⁶⁹ use of different closures can strongly influence the quality of the GQME+semiclassics dynamics for nonequilibrium averages. We demonstrate that the equilibrium problem is somewhat more sensitive to the choice of forward or backward closures than its nonequilibrium counterpart, but the difference between them is generally minor. In addition, the results in this section provide further numerical support for the analytical arguments furthered in Ref. 70, which states that closures that re-express the auxiliary kernels in terms of the original correlator and its time derivatives and numerically equivalent correlators must recover the original Ehrenfest dynamics (closures cf2, cf3, cb2, and cb3), while the other closures yield GQME dynamics that are distinct from the bare Ehrenfest result.

Initially, we focus on the differences between dynamics produced via the $Q$-forward and $Q$-backward closures. Taking the zeroth order closures cf0 and cb0, which do not replace the action of the Liouvillian with numerical derivatives, as representative of the differences between the forward and backward closures, it is clear that the forward closure performs slightly better than the backward closure. This manifests in the difference between the GQME results in panels (a)–(d) of Fig. 5. The reason for this discrepancy, however, is not well understood and will require further analysis beyond the scope of the present work. These differences notwithstanding, both the cf0 and cb0 closures clearly yield dynamics that are significantly more accurate than the bare quasiclassical results.

We now focus on the dynamics generated using the $Q$-forward closures. As is clear in panels (a) and (b) in Fig. 7, both the cf0 and cf1 closures yield dynamics that are in better agreement with the numerically exact results than the bare Ehrenfest dynamics. This is consistent with results obtained for the nonequilibrium problem.⁶⁹ In contrast, the cf2 and cf3 closures recover the direct Ehrenfest result, as shown in panels (c) and (d). Importantly, this agreement provides additional numerical support for the analytical arguments put forth in Ref. 70. Specifically, because the auxiliary kernels for closure cf3 (and cb3) can be written in terms of the original ECF in Eq. (8) and its time derivatives, the memory kernel obtained using this closure must recover the bare Ehrenfest result. Because the cf2 (and cb2) closure also reproduces the Ehrenfest dynamics, as was the case in the nonequilibrium problem,⁶⁹ it is clear that the Ehrenfest method permits replacing the action of the Liouvillian operator acting on a dynamically evolved operator with its time derivative. In turn, the fact that the cf0 and cf1 (and cb0 and cb1) closures yield GQME dynamics that are distinct from the Ehrenfest results is consistent with the observation that the Ehrenfest procedure is not ensemble-conserving. Nevertheless, the arguments in Ref. 70 do not explain the improvements afforded by the cf0 and cf1 closures. Instead, we again emphasize that the most likely reason for the improvement afforded by the cf0 and cf1 closures lies in the exact sampling of distinct bath operators alluded to in Sec. III A.

Again restricting our attention to the forward closures, the kernels produced via the various closures can shed light on the factors that lead to distinct GQME dynamics. Like the dynamics for the different closures shown in Fig. 7, the memory kernels shown in Fig. 6 suggest that the closures can be subdivided into two main sets which yield similar memory kernels: the first group consists of closures cf0 and cf1, while the second consists of cf2 and cf3. The main difference between the kernels produced by these two groups is clearest at intermediate times, where the cf2 and cf3 closures display a more strongly oscillatory behavior than the cf0 and cf1 closures. These differences, in fact, lead to the different dynamical behavior which either reproduces the bare Ehrenfest result (cf2 and cf3 closures) or yields GQME dynamics that are in close agreement with the numerically exact results (cf0 and cf1 closures). It bears repeating that we have only shown the four kernel elements which are analytically (and numerically, for closures cf0 and cf1) known to be nonzero. Nevertheless, the cf2 and cf3 generally lead to nonzero intermediate-time results for the other kernel elements (not shown). These non-zero contributions, while clearly undesirable if one wants to recover the numerically exact result, are essential for reproducing the bare Ehrenfest dynamics by means of the GQME.

A particularly appealing feature of the GQME approach is the promise of long-time dynamics for a comparatively small computational expense. As shown in Refs. 67 and 69, the GQME+MFT approach has proven to be capable of correctly recovering the long-time limit for nonequilibrium averages. The equilibrium formulation of the GQME+MFT approach is similarly capable of recovering accurate results in the long-time limit. The two panels of Fig. 8 show the dynamical behavior of the ECF for the parameters used in Fig. 7 for intermediate times (panel (a)) as well as over very long times (panel (b)). Panel (b) contains several important pieces of information. First, the violet region demarcates the range of time where exact dynamics are available, and the dashed line marks the correct long-time limit of the correlation function, $Re⟨σz(t→∞)σz(0)⟩=⟨σz⟩2≈0.75$⁠. This panel demonstrates that the closures that lead to improvements in accuracy for intermediate-time equilibrium dynamics and nonequilibrium averages, i.e., closures cb0. cb1, cf0, and cf1, also recover a significantly improved estimate of the long-time limit for the ECFs. It also bears remarking that the GQME dynamics are remarkably robust with respect to the choice of the cutoff time, $τc$⁠.

A final point of interest concerns the use of the approach presented here to extend the temporal range of standard semiclassical simulations. Due to the fact that the statistical error in quasi- and semi-classical methods increases with simulation time, converging quasi- and semi-classical dynamics on large time-scales can become prohibitively expensive. To ameliorate this problem, one can instead use the GQME+MFT dynamics obtained using the cf3 closure, which is known to reproduce the direct Ehrenfest result. In this way, it is possible to extend the Ehrenfest dynamics to nontrivially long times for a much smaller expense. For instance, a $10|Δ|−1$-long Ehrenfest calculation was sufficient to recover the Ehrenfest dynamics for $t>100|Δ|−1$ shown in Fig. 8.

In this paper, we have extended the GQME+MFT method introduced in the first paper of this series⁶⁹ to treat ECFs within the context of the Mori formalism. As Ref. 69 and the present work demonstrate, the Mori-based formulation furthered here provides a flexible framework to accurately study problems both in and out of equilibrium for SB-type models. We further emphasize that this approach is general and can be easily generalized to arbitrarily complex systems and dynamical quantities.

In the same spirit as earlier work on nonequilibrium dynamics,^57–70 we have demonstrated that equilibrium memory kernels for the SB model can be short-lived in comparison to the ECF lifetime. Indeed, it is possible to say that impurity-type systems at relatively high temperatures where the bath consists of a large number of modes with a broad distribution of energies have short-lived memory kernels, as long as the projection operator includes all system states. We have also shown that, as in the nonequilibrium case, the GQME+semiclassics approach is capable of yielding impressive improvements in computational efficiency and accuracy over direct application of bare quasi- and semi-classical methods. Consequently, at least for these types of models, we foresee that the GQME approach can become an invaluable tool for obtaining highly accurate nonequilibrium and equilibrium dynamics at a lower computational cost than would be necessary with conventional methods that go well beyond SB-type models, including the dynamics of liquid state systems with no system-bath distinction.

Via a systematic analysis of the dependence of the GQME dynamics on the kernel closures, we have further confirmed the analytical arguments posited in Ref. 70. Specifically, the results obtained from the cf3 and cb3 closures, which reproduce the direct Ehrenfest dynamics, confirm the proof presented in Ref. 70 which states that when the auxiliary kernels can be written in terms of the original ECF and its time derivatives, the GQME results will be equivalent to those arising from direct application of the dynamical method. The fact that the cf2 and cb2 closures also reproduce the Ehrenfest results is in agreement with the observation that the numerical time derivative of a correlation function is equivalent to the action of the Liouvillian on the dynamically sampled operator in the Ehrenfest method. We have also demonstrated that the cf0 closure yields the most accurate dynamics over a wide region of parameter space. Further, although the cb0 closure has been shown to produce slightly less accurate dynamics for some sets of parameters, it also affords a significant improvement over the bare mean-field dynamics, as do the cf1 and cb1 closures. The distinct results obtained from the cf0, cb0, cf1, and cb1 closures have also been shown to be consistent with the fact that the Ehrenfest method is not ensemble-conserving.⁷⁰ 

To conclude, we emphasize that the Mori formulation presented in this work and Ref. 69 is equally applicable to impurity-type models like the SB model and systems such as spin and fermion lattice models,^72–74 and quantum fluids.^30,75–77 In future papers, we explore the extension of the Mori framework to multiple-time correlation functions and systems coupled to arbitrary baths as well as problems where the system-bath distinction is absent.

The authors would like to thank Aaron Kelly and Tom Markland for helpful discussions. D.R.R. acknowledges support from NSF Grant No. CHE-1464802. A.M.C. thanks Hsing-Ta Chen for useful conversations.

After a partial Wigner transform with respect to the bath degrees of freedom, it is possible to write the system components of the canonical density operator as follows:⁸³ 

where $ρ=e−βH/Tr[e−βH]$⁠, $a,b∈{0,1}$ correspond to the two system states of the SB model, N_ab is a temperature dependent normalization constant, and $Ra,bW(𝐱,𝐩)$ is a bath operator of unit partial trace, i.e., $∫d𝐱d𝐩Ra,bW(𝐱,𝐩)=1$⁠, which can be interpreted as the bath distribution function corresponding to the subsystem state $|a⟩⟨b|$⁠. We henceforth drop the dependence of the bath distribution function on the bath coordinates and momenta, (x, p), for notational clarity.

To obtain expressions for N_ab and $Ra,bW$⁠, we follow the numerically exact framework presented in Ref. 83, which relies on the path integral expansion of the Boltzmann factor. This expansion employs the Trotter factorization, allowing us to rewrite the full Boltzmann factor as an N-membered product of basic path integral units,

Rigorously, this approximation becomes exact in the limit of $N→∞$⁠, but in practical calculations convergence is achieved with a small N.

By substituting Eq. (A2) into Eq. (A1), introducing resolutions of the identity $𝟏S=∑a|a⟩⟨a|$ and $𝟏B=∫d𝐪|𝐪⟩⟨𝐪|$ for the system and bath subspaces, respectively, and integrating over the bath degrees of freedom one can obtain expressions for N_ab and $Ra,bW$⁠. While the previous procedure formally eliminates the dependence on the bath degrees of freedom introduced with the resolutions of the identity, the same cannot be said of the system degrees of freedom. In fact, the number of paths grows exponentially with the number of path integral steps, and each path can be associated with a sequence of indices that indicate the electronic states visited along said path. To illustrate this point, consider the simpler case of applying the path integral formalism to the Boltzmann factor corresponding to the isolated two-level system,

where $kn∈{0,1}$ corresponds to the state used in the n^th path integral step of the expansion. For this decomposition, one may consider the path to consist of a string of numbers describing the identity of the states taken in the sequence, e.g., in an N = 3 expansion with the endpoints fixed, there are 4(=2^N−1) such sequences, {k₃,0,0,k₀}, {k₃,0,1,k₀}, {k₃,1,0,k₀}, and {k₃,1,1,k₀}. In the following, we refer to the sum over the set ${k1,…,kN−1}$ as the sum over paths.

Before providing expressions for N_ab and $Ra,bW$⁠, we provide a few definitions that render the notation simpler. We begin with the following basic definitions:

and the following path-independent quantities:

where A^(l) is a tridiagonal $N−1×N−1$ matrix whose diagonal and off-diagonal entries are equal to 2 and $−sech(2θl)$⁠, respectively.

A few quantities depend on the path taken in configuration space. These quantities, labeled by a tilde, take the following forms:

Using these definitions,

where $W∼a…b=S∼a,bexp[−∑lΛ∼(l)],$ and $Wa,b=∑pathsW∼a,b$⁠.

Finally, in  Appendix B, the Wigner transform of the product of $ρ$ and $VB=α∑lclxl$ will require the evaluation of the following operator:

We first introduce a notation based on static and dynamic matrices that facilitate the construction of the auxiliary kernels. Given the definition of the inner product in Eq. (7), one can easily verify that

where $Y←=RYR−1$⁠, $Y→=T−1YT$⁠, $Ynm=−2εznm$ is a static transformation matrix, $εijk$ is the Levi-Civita tensor (which is only defined with respect to the indices x, y, and z), $Ynm=0$ when n or m = 1, and

For convenience, we also introduce the following dynamical variables in the original (pre-rotation) basis:

where

To construct the auxiliary kernels, we need to rotate these dynamical matrices to the orthonormal basis of the projector,

In this notation, $C(t)=C(00)(t)$ and $⟨𝐕𝐁⟩=C(01)(t=0)=C(10)(t=0)$⁠.

Using these expressions, one can rewrite the auxiliary kernels thus,

In analogy to the definition of nonequilibrium averages,⁸⁷ one may express a correlation function within the Ehrenfest formalism as

where the superscript W indicates the partial Wigner transform of the operator with respect to the bath degrees of freedom, and N is the number of degrees of freedom over which the Wigner transform is being performed. The Wigner transform of an operator is defined as⁸⁸ 

In the following, we only perform the partial Wigner transform with respect to the bath degrees of freedom, as is required by the Ehrenfest procedure.

For equilibrium correlation functions of the form given by Eq. (B6), it is necessary to obtain an expression for the partially transformed canonical density operator $ρ=e−βH/Tr[e−βH]$⁠. Here we adopt the approach presented in Ref. 83, and express the partial Wigner transform of the canonical operator as

where $Sa$ is a pure system operator, and $BaW$ corresponds to a normalized bath density operator, i.e., $∫dΓBaW(𝐱,𝐩)=1$ (for detailed expressions, see  Appendix A). For convenience, we re-express $Sa$ in terms of a convenient basis, namely, $Sa=∑brabAb$⁠, where $Ab∈{𝟏S,σx,σy,σz}$⁠. Using the above notation, we rewrite the expressions for the dynamic matrices in Eq. (B6) as follows:

where

Noting that the Wigner transform of operator products can be obtained via the Moyal expansion⁸⁸ 

where

we can express the commutator and anticommutator of the bath density operator and the bath part of the interaction as

where

The expression for $ζBW$ can be found in Eq. (A19).

Substituting expressions (B18) into (B13) yields