The ability to efficiently and accurately calculate equilibrium time correlation functions of many-body condensed phase quantum systems is one of the outstanding problems in theoretical chemistry. The Nakajima-Zwanzig-Mori formalism coupled to the self-consistent solution of the memory kernel has recently proven to be highly successful for the computation of nonequilibrium dynamical averages. Here, we extend this formalism to treat symmetrized equilibrium time correlation functions for the spin-boson model. Following the first paper in this series [A. Montoya-Castillo and D. R. Reichman, J. Chem. Phys. **144**, 184104 (2016)], we use a Dyson-type expansion of the projected propagator to obtain a self-consistent solution for the memory kernel that requires only the calculation of normally evolved auxiliary kernels. We employ the approximate mean-field Ehrenfest method to demonstrate the feasibility of this approach. Via comparison with numerically exact results for the correlation function $Czz(t)=Re\u27e8\sigma z(0)\sigma z(t)\u27e9$, we show that the current scheme affords remarkable boosts in accuracy and efficiency over bare Ehrenfest dynamics. We further explore the sensitivity of the resulting dynamics to the choice of kernel closures and the accuracy of the initial canonical density operator.

## I. INTRODUCTION

By encoding a system’s response to a weak perturbation,^{1} equilibrium time correlation functions (ECFs) provide direct access to important dynamical quantities such as transport coefficients,^{2,3} absorption spectra,^{4} 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.^{57} 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,^{69} 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\u27e8\sigma z(0)\sigma z(t)\u27e9$. 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.

## II. MORI APPROACH

As in Paper I,^{69} 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.^{76}

The SB Hamiltonian takes the form *H* = *H*_{S} + *H*_{B} + *H*_{SB}, where

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

The bath Hamiltonian consists of independent harmonic oscillators,

where *P*_{k}, *Q*_{k}, and $\omega 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 *k*th oscillator and $\alpha =\xb11$. 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^{18} with an exponential cutoff. Here, $\omega c$ is the cutoff frequency, which determines the correlation time for the bath at finite temperature. The Kondo parameter, $\xi $, describes the strength of the system-bath coupling and is proportional to the reorganization energy, $\lambda =\xi \omega c/2=\pi \u22121\u222b0\u221ed\omega J(\omega )/\omega $, 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.

### A. Mori-type GQME

As established in Paper I,^{69} 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\u2212P$ is the complementary projection operator.

For ECFs, the Mori-type projector traditionally consists of an operator *B* and its time derivative $B\u02d9=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\u2208{\U0001d7cfS,\sigma x,\sigma y,\sigma z}$, $\rho =e\u2212\beta H/Tr[e\u2212\beta H]$ is the canonical density operator, $\beta =[kBT]\u22121$ is the inverse thermal energy, and the transformation matrices $T$ and $R$ enforce the orthonormalization of the members of the projector ($R(\mathbf{A}|\mathbf{A})T=\U0001d7cf$). 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[\rho {An\u2020,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^{57} 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,^{69} 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\u2208{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, *cb*2 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,^{69} 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.

## III. RESULTS

As in the first paper of this series,^{69} 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 $\rho $ 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\u221210$, where $R.E.=max[abs[Kn+1(t)\u2212Kn(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, $\tau c$, beyond which the memory kernel is set to zero. In practice, $\tau 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.^{84} 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).

### A. Improvements in efficiency and accuracy

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 ($\epsilon =0$) system coupled to a moderately fast bath ($\omega c=2.5$), and panels (c) and (d) to a strongly biased system ($\epsilon =\u22122$) coupled to a fast bath ($\omega 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 $\tau 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 $\tau c\u22652.5$. It bears remarking that the choice of $\tau c$ for systems characterized by strong system-bath coupling is somewhat more difficult. For panels (c) and (d), the choice of $\tau c$ indicates that the dynamics are insensitive for $\tau c\u22652.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.^{32} 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, $\tau c$. Indeed, as the results demonstrate, it is possible to change $\tau 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,^{69} 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 $\zeta \mathit{n}\mathit{o}\mathit{n}\mathit{e}\mathit{q}W=\u2212\alpha \u2211kckPktanh(\beta \omega k/2)/\omega 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.^{38} The two simplest approximations consist of rewriting the Boltzmann factor as simple products of system and bath operators, e.g., $e\u2212\beta H\u2248e\u2212\beta HSe\u2212\beta HB$ and $e\u2212\beta H\u2248e\u2212\beta (HB+HSB)/2e\u2212\beta HSe\u2212\beta (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 $\xi $) or temperatures are low (large $\beta $), 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, $\rho $.

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 ($\xi =2.55$) and moderate to low temperatures ($\beta =1.6$), a large number of path integral steps (*N* = 6) are required to achieve convergence. Conversely, the high temperature ($\beta =0.806$) and weak coupling ($\xi =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 $\rho $. 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 $\rho $ lends credence to our conjecture^{69} 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.

### B. Memory kernel closures

As was demonstrated in the first paper of this series,^{69} 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 *cf*2, *cf*3, *cb*2, and *cb*3), 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 *cf*0 and *cb*0, 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 *cf*0 and *cb*0 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 *cf*0 and *cf*1 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.^{69} In contrast, the *cf*2 and *cf*3 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 *cf*3 (and *cb*3) 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 *cf*2 (and *cb*2) closure also reproduces the Ehrenfest dynamics, as was the case in the nonequilibrium problem,^{69} 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 *cf*0 and *cf*1 (and *cb*0 and *cb*1) 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 *cf*0 and *cf*1 closures. Instead, we again emphasize that the most likely reason for the improvement afforded by the *cf*0 and *cf*1 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 *cf*0 and *cf*1, while the second consists of *cf*2 and *cf*3. The main difference between the kernels produced by these two groups is clearest at intermediate times, where the *cf*2 and *cf*3 closures display a more strongly oscillatory behavior than the *cf*0 and *cf*1 closures. These differences, in fact, lead to the different dynamical behavior which either reproduces the bare Ehrenfest result (*cf*2 and *cf*3 closures) or yields GQME dynamics that are in close agreement with the numerically exact results (*cf*0 and *cf*1 closures). It bears repeating that we have only shown the four kernel elements which are analytically (and numerically, for closures *cf*0 and *cf*1) known to be nonzero. Nevertheless, the *cf*2 and *cf*3 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.

### C. Long-time dynamics

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\u27e8\sigma z(t\u2192\u221e)\sigma z(0)\u27e9=\u27e8\sigma z\u27e92\u22480.75$. This panel demonstrates that the closures that lead to improvements in accuracy for intermediate-time equilibrium dynamics and nonequilibrium averages, i.e., closures *cb*0. *cb*1, *cf*0, and *cf*1, 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, $\tau 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 *cf*3 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|\Delta |\u22121$-long Ehrenfest calculation was sufficient to recover the Ehrenfest dynamics for $t>100|\Delta |\u22121$ shown in Fig. 8.

## IV. CONCLUSIONS

In this paper, we have extended the GQME+MFT method introduced in the first paper of this series^{69} 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 *cf*3 and *cb*3 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 *cf*2 and *cb*2 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 *cf*0 closure yields the most accurate dynamics over a wide region of parameter space. Further, although the *cb*0 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 *cf*1 and *cb*1 closures. The distinct results obtained from the *cf*0, *cb*0, *cf*1, and *cb*1 closures have also been shown to be consistent with the fact that the Ehrenfest method is not ensemble-conserving.^{70}

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.

## ACKNOWLEDGMENTS

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.

### APPENDIX A: PATH INTEGRAL TREATMENT OF THE WIGNER TRANSFORMED CANONICAL DENSITY OPERATOR

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:^{83}

where $\rho =e\u2212\beta H/Tr[e\u2212\beta H]$, $a,b\u2208{0,1}$ correspond to the two system states of the SB model, *N*_{ab} is a temperature dependent normalization constant, and $Ra,bW(\mathbf{x},\mathbf{p})$ is a bath operator of unit partial trace, i.e., $\u222bd\mathbf{x}d\mathbf{p}Ra,bW(\mathbf{x},\mathbf{p})=1$, which can be interpreted as the bath distribution function corresponding to the subsystem state $|a\u27e9\u27e8b|$. 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\u2192\u221e$, but in practical calculations convergence is achieved with a small *N*.

By substituting Eq. (A2) into Eq. (A1), introducing resolutions of the identity $\U0001d7cfS=\u2211a|a\u27e9\u27e8a|$ and $\U0001d7cfB=\u222bd\mathbf{q}|\mathbf{q}\u27e9\u27e8\mathbf{q}|$ 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\u2208{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*_{3},0,0,*k*_{0}}, {*k*_{3},0,1,*k*_{0}}, {*k*_{3},1,0,*k*_{0}}, and {*k*_{3},1,1,*k*_{0}}. In the following, we refer to the sum over the set ${k1,\u2026,kN\u22121}$ 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\u22121\xd7N\u22121$ matrix whose diagonal and off-diagonal entries are equal to 2 and $\u2212sech(2\theta 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\u223ca\u2026b=S\u223ca,bexp[\u2212\u2211l\Lambda \u223c(l)],$ and $Wa,b=\u2211pathsW\u223ca,b$.

Finally, in Appendix B, the Wigner transform of the product of $\rho $ and $VB=\alpha \u2211lclxl$ will require the evaluation of the following operator:

### APPENDIX B: EXPRESSIONS FOR THE AUXILIARY KERNELS

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\u2190=RYR\u22121$, $Y\u2192=T\u22121YT$, $Ynm=\u22122\epsilon znm$ is a static transformation matrix, $\epsilon 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 $\u27e8\mathbf{V}\mathbf{B}\u27e9=C(01)(t=0)=C(10)(t=0)$.

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

#### 1. Quasi-classical treatment of $Cnm(jk)(t)$

In analogy to the definition of nonequilibrium averages,^{87} 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^{88}

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 $\rho =e\u2212\beta H/Tr[e\u2212\beta 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., $\u222bd\Gamma BaW(\mathbf{x},\mathbf{p})=1$ (for detailed expressions, see Appendix A). For convenience, we re-express $Sa$ in terms of a convenient basis, namely, $Sa=\u2211brabAb$, where $Ab\u2208{\U0001d7cfS,\sigma x,\sigma y,\sigma 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^{88}

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 $\zeta BW$ can be found in Eq. (A19).

Importantly, the system operators *X*_{pn} and *Y*_{pn} can be written as linear combinations of the subsystem outer products, i.e., ${|j\u27e9\u27e8k|}$, which makes direct evaluation of these functions straightforward via the Ehrenfest method. We direct the reader to Ref. 69 for instructions on how to properly treat the system initial conditions for these functions.

#### 2. Ehrenfest method

In the Ehrenfest method,^{31,32,35,87} the system (bath) evolves in the mean field of the bath (system). Additionally, the dynamics of the bath are assumed to be accurately captured by classical mechanics. Naturally, the mean field and classical bath approximations imply that the Ehrenfest method is not applicable to cases where the system and bath are strongly coupled (large $\lambda $) and baths characterized by high frequencies (large $\omega c$) or low temperatures, respectively.

Classical treatment of the bath requires that the continuous spectral density be discretized into individual oscillators. We use the approach outlined in Ref. 49 where the spectral density is decomposed into *f* oscillators. The frequency of the *k*th oscillator takes the form,

and the coupling constant,

Evolution of the system and bath within the Ehrenfest scheme is given by the following equations of motion. For the system, the Liouville equation under the influence of a modified Hamiltonian determines the evolution of the system’s wavefunction,

where

is the modified system Hamiltonian and $\lambda cl(t)=\alpha \u2211kckQk(t)$ is the classical fluctuation that modulates the system bias. Here, $\rho S(0)=|\psi S(0)\u27e9\u27e8\psi S(0)|$ is the initial density matrix for the system, which must be initialized in a pure state (for a detailed discussion, see Ref. 69).

The evolution of the bath variables is given by Hamilton’s equations subject to Hamiltonian of the free bath modified by a quantum force arising from the time-dependent populations in the system,

where

and $\sigma \xafz(t)=TrS[\rho S(t)\sigma z].$

A second-order Runge-Kutta scheme was used to calculate $Cnm(ij)(t)$. During individual time steps, $\sigma \xafz(t)$ is kept constant for the evolution of the bath, while $\lambda cl(t)$ is kept constant during the evolution of the system. Over a half time step, the equations for the classical variables take the forms,

and

where

Approximately 10^{6} trajectories were necessary to achieve convergence.

## REFERENCES

It should be noted that there are problems where the memory kernels never decay to zero. A clear example of this is the case of a SB model where the bath consists of one oscillator, where the persistent quantum coherence results in a memory kernel that never decays. For problems where decoherence and dissipation are expected to be significant, a non-decaying memory kernel is a sign that the projection operator has to be chosen differently. Clearly, the GQME formalism cannot provide any advantages with respect to computational efficiency for cases where the memory kernel never decays.