# Approximate but Accurate Quantum Dynamics from the Mori Formalism: I. Nonequilibrium Dynamics

## Abstract

We present a formalism that explicitly unifies the commonly used Nakajima-Zwanzig approach for reduced density matrix dynamics with the more versatile Mori theory in the context of nonequilibrium dynamics. Employing a Dyson-type expansion to circumvent the difficulty of projected dynamics, we obtain a self-consistent equation for the memory kernel which requires only knowledge of normally evolved auxiliary kernels. To illustrate the properties of the current approach, we focus on the spin-boson model and limit our attention to the use of a simple and inexpensive quasi-classical dynamics, given by the Ehrenfest method, for the calculation of the auxiliary kernels. For the first time, we provide a detailed analysis of the dependence of the properties of the memory kernels obtained via different projection operators, namely the thermal (Redfield-type) and population based (NIBA-type) projection operators. We further elucidate the conditions that lead to short-lived memory kernels and the regions of parameter space to which this program is best suited. Via a thorough analysis of the different closures available for the auxiliary kernels and the convergence properties of the self-consistently extracted memory kernel, we identify the mechanisms whereby the current approach leads to a significant improvement over the direct usage of standard semi- and quasi-classical dynamics.

## 1Introduction

The continued effort to develop accurate and efficient approaches for the calculation of the dynamics of many-body quantum systems has produced a rich variety of methods, ranging from the numerically exact to the approximate. While exact methods provide important benchmark results for model systems,[1] their computational cost makes them impractical for realistic multidimensional systems. Conversely, approximate methods, whether perturbative[9] or based on quasi-[12] or semi-classical[17] approaches, tend to scale more gracefully with system size and can address both model and realistic systems, albeit at the expense of general applicability and accuracy. For cases where one is interested in the dynamics of a limited number of degrees of freedom, the Nakajima-Zwanzig (NZ) equation[24] provides a useful starting point for a plethora of methods based on generalized quantum master equations (GQME).

The NZ equation, which may be derived via the projection operator technique,[26] dictates the evolution of the reduced density matrix (RDM) for the portion of the Hilbert space denoted as the system. The influence of the complementary subspace, referred to as the bath, on the RDM dynamics appears in the form of a memory term, full knowledge of which is tantamount to solving the original problem. The apparent simplicity of the NZ equation, however, belies the complexity of the memory term, which can be formidably difficult if not impossible to calculate exactly. Despite the seeming conservation of difficulty, different treatments of the memory kernel have led to manageable and often very successful approximate and numerically exact schemes.[10]

A major difficulty in the calculation of the memory kernel lies in the fact that its dynamical evolution involves the “projected” propagator, , where is the projector that defines the reduced dynamics. To circumvent the problem of projected dynamics, Shi and Geva[31] proposed a self-consistent expansion of the memory kernel, which requires the calculation of auxiliary kernels evolved with the *normal* rather than the *projected* propagator. From an exact perspective, this approach is useful only if the numerical effort necessary for the calculation of the auxiliary memory kernels is less than that required for the direct calculation of the system dynamics. Using the numerically exact quasi-adiabatic path integral (QUAPI) method,[2] Shi and Geva have shown that the memory kernel for the spin-boson model can decay up to 10 times faster than the system dynamics,[31] lending credence to the feasibility of the self-consistent approach. More recently, a similar scheme has been used by Rabani and co-workers within a path integral framework for the study of quantum transport problems.[34] Just as importantly, applications of the method have successfully used semi- and quasi-classical theories to calculate the auxiliary kernels.[32] In particular, Kelly, Markland, and coworkers have illustrated the impressive accuracy and robustness of this approach in both model and realistic problems.[39] These studies have led to two important conclusions: (i) the memory kernels are short-lived in a wide region of parameter space for canonical problems such as the spin-boson model and (ii) the self-consistent solution of the memory kernel using approximate dynamics can yield impressively more accurate results than *direct* simulation of the RDM dynamics using the very same approximate method.

Despite these important results, questions of general applicability still remain. For instance, the conditions that lead to short-lived memory kernels remain unknown. Further, it is still unclear how the breakdown of approximate methods (in unfavorable parameter regimes) affects the quality of the self-consistently extracted memory function. Perhaps most importantly, an understanding of the necessary and sufficient conditions from which the improvement in the accuracy of approximate dynamics from use in conjunction with the memory function formalism as observed in Refs. is lacking. Finally, the convergence properties of different versions of the auxiliary memory kernels arising from the alternative closures in the self-consistent expansion of the memory function have not yet been explored. We expect these convergence properties to differ when approximate methods are employed in the calculation of the auxiliary kernels.

The remarkable utility of the NZ equation notwithstanding, objects beyond single-time nonequilibrium dynamics are either cumbersome to obtain within this framework. Despite this difficulty, recent work has generalized the NZ equation to multi-time correlation functions.[42] In contrast, the more flexible Mori formulation permits direct extension to multiple-time and equilibrium correlation functions, which are essential in the treatment of linear[43] and non-linear spectroscopy,[44] and the calculation of chemical rate constants[45] and kinetic coefficients,[46] to name a few examples. For this reason, in paper I (this paper) of this series, we provide a unified Mori-type framework to approach single-time nonequilibrium correlation functions, and address several of the open questions listed above. In a second paper, we present the a similar framework to treat equilibrium correlation functions. It should be noted that a major advantage of the Mori formulation is that it can naturally address problems where no system-bath distinction exists, such as spin and fermion lattice models,[47] and quantum fluids.[50]

The structure of this paper is as follows: In Section 2, we present the formalism for nonequilibrium correlation functions from the Mori perspective and show that, with the appropriate choice of projection operator, one recovers equations identical to those arising from the conventional NZ treatment. Section 2 introduces the spin-boson model and the proposes two types of projection operators for this model. Section 3 discusses different closures of the memory kernel based on different placements of the projection operator and on the use of time derivatives. To illustrate the arguments related to convergence, we implement the mean field Ehrenfest method, as first proposed in Ref. , to obtain the auxiliary memory kernels. We henceforth refer to the use of the Ehrenfest method coupled to the self-consistent extraction of the memory kernels as the GQME+MFT approach. In the Section 4, show the different properties of the memory kernels associated with the Redfield- and NIBA-type projectors, explore the performance of the GQME+MFT approach to SB models whose system-bath coupling is characterized by Ohmic and Debye spectral densities, and investigate the convergence properties of the distinct closures introduced in Section 3. In Section 5, we conclude.

## 2Mori Approach

For illustrative purposes, we focus on the spin-boson (SB) Hamiltonian, which is representative of typical open quantum systems, but note that the current approach is general and may be applied to any Hamiltonian.[54] In particular, a major advantage of the Mori formalism developed here over the traditional NZ approach is the ability to treat systems with no natural system-bath separation, such as spin-chains and lattice models[47] or liquids.[50] We reserve these applications for later work.

The SB Hamiltonian takes the form . It contains a system part consisting of two sites offset by an energy bias 2 and with off-diagonal coupling , which is assumed to be independent of the bath coordinates,

where corresponds to the Pauli matrix. The bath part of the Hamiltonian consists of independent harmonic oscillators,

where and are the mass weighted momenta and coordinates for the harmonic oscillator, respectively, and is the frequency of the mode. The coupling between the system and bath is assumed to be linear in the bath coordinates and diagonal and antisymmetric with respect to the system basis,

where is the coupling constant between the system and the oscillator, and depending on the definition of the model. The spectral density, , fully characterizes the system-bath interaction, and takes the form,

It is common to assume one of several forms for the spectral density. Two important cases describe Ohmic dissipation in condensed phase systems where is proportional to as . These are the standard Ohmic spectral density[11] characterized by an exponential cutoff, and the Debye spectral density characterized by a Lorentzian cutoff:

where the cutoff frequency determines the correlation time of the bath at sufficiently high temperatures.[55] The reorganization energy, , is a measure of the strength of the system-bath coupling. In the case of the Ohmic spectral density, the Kondo parameter, , is often used instead to gauge the coupling strength.

To assess the applicability of the formalism presented here, we compare our results to exact nonequilibrium population dynamics of the SB model,

where is the equilibrium density operator for the uncoupled bath, is the inverse thermal energy, and the initial condition corresponds to a Frank-Condon transition.

### 2.1Generalized Nakajima-Zwanzig-Mori Equation

Here, we deviate from the derivations commonly given for the NZ equation or the Mori equation of motion (EOM) for an operator and instead focus on a generalized EOM for the full propagator,

where is the projection operator that defines the subsystem whose dynamics we seek and is the complementary projection operator. This equation is general and can be employed within both the NZ and Mori approaches. For instance, applying Eq. (Equation 3) on an operator , yields the Mori EOM for that operator. Conversely, taking the Hermitian conjugate of Eq. (Equation 3), applying it on the initial density of the system and bath, , and acting the projection operator from the left followed by a trace over the bath degrees of freedom yields the NZ equation for the system’s RDM.

Essentially all work on the self-consistent expansion of the memory kernel has employed the thermal (Argyres-Kelley) projection operator ,[31] where is a bath operator with unit trace and corresponds to partial trace over the bath. To use the Mori formulation, it is convenient to rewrite the thermal projector, in the Heisenberg picture using Liouville notation,^{2}

Using the thermal type projector, applying Eq. (Equation 3) to , and closing on the left with , where is again a bath operator with unit trace that may be different from in the projection operator, yields the following set of EOMs for system observables,

where is a static rotation matrix, corresponds to nonequilibrium averages of populations and coherences with all possible factorizable initial conditions, and is the so-called the inhomogeneous term. The elements of the memory kernel are given by,

When , the inhomogeneous term disappears, . Often, is chosen such that , which correspond to the harmonic oscillator bath at equilibrium with the ground electronic state or with one of the two excited states, respectively. Initial conditions of the form , where is an arbitrary system operator and is taken from the set above, correspond either to a Frank-Condon excitation where the bath is in the electronic ground state also called the spectroscopic initial condition, or a charge transfer initial condition where the bath is in equilibrium with one of the excited states.

We henceforth refer to the thermal projector above with the additional restriction that as Redfield-type.[9] The reason for this name is that truncation of the memory kernel at second order in is equivalent to a second-order perturbative expansion of the memory kernel with respect to the system-bath coupling, which corresponds to Redfield theory. We further note that the choice for remains flexible.

An important alternative for the projection operator consists of restricting the set to the system populations , and choosing . When using this projection operator, a similar second-order truncation of the memory kernel with respect to leads to equations that are equivalent to the non-interacting blip approximation (NIBA), which is a second order expansion in the electronic coupling as opposed to the system-bath coupling.[11] Accordingly, we hereafter refer to this projector as NIBA-type. As in the case of the Redfield-type projector, the choice for determines whether the inhomogeneous term is zero or finite.

## 3Self-Consistent Expansions for

As mentioned in Section 1, the main difficulty associated with the memory kernel, Eq. (Equation 5), is the presence of the projected propagator, . To circumvent this problem, Shi and Geva proposed the use of the Dyson identity,

yielding a self-consistent expansion of the memory kernel that only involves unprojected dynamics.[31] Not surprisingly, Eq. (Equation 3) can also be derived using the Dyson identity, Eqs. (Equation 6) and ( ?). Despite considering only time-independent Hamiltonians in this work, extension of the formalism to time-dependent Hamiltonians is simple and only requires the time-ordered form for the propagators in Eqs. (Equation 6) and ( ?).

It is important to remark that in the literature, different variants of the Dyson expansion have led to a menagerie of seemingly distinct expressions, which differ with respect to the number of and type of auxiliary kernels employed.[33] When these distinct expressions are evaluated via exact methods, all expansions yield equivalent results, up to numerical errors. However, when the auxiliary kernels are computed via approximate methods, different expansions can lead to memory kernels with different properties. In the following, we show that there is only a limited number of integral equation forms that can yield numerically distinct approximate memory kernels from the self-consistent solution of the resulting integral equations.

### 3.1Bare expansions: Backward and Forward

The expansion of the memory kernel, Eq. (Equation 5), takes the following form,

where the superscript refers to the placement of the in the projected propagator as “backward” with respect to the placement of the Liouvillian, i.e., . The auxiliary memory kernels take the forms

A seemingly distinct type of closure that is commonly used in the literature involves a *third* auxiliary memory kernel. In Appendix B, we show that the three-member expansions are equivalent to the two-member expansions given by Eqs. (Equation 7) and (Equation 9). We further note that the auxiliary kernels we obtain are equivalent to those used by others in the field.[31]

A second set of closures makes use of the identity . To indicate that the is to the right of the Liouvillian in the propagator, we have used the superscript (indicating a “forward” placement). Inserting the previous identity in Eq. (Equation 5) leads to the following expansion upon insertion of the Dyson identity,

We note that has the form given in Eq. (Equation 8), and has the following form,

It bears remarking that Eqs. (Equation 9) and (Equation 7) differ in the placement of under the integral and Eqs. (Equation 10) and ( ?) differ in whether or its Hermitian conjugate act on operators that require sampling only at or at finite times.

### 3.2Expansions using time-derivatives

Because the auxiliary kernels given by Eqs. (Equation 8) and ( ?) or (Equation 10) require sampling of additional bath operators at and at finite times, convergence of these functions, at least within the context of semi- and quasi-classical methods, necessitates the sampling of a larger number of bath realizations than for the direct simulation of , making the initial step of the calculation more expensive, even if trivially parallelizable. To avoid this added complexity and expense, the expressions for the auxiliary kernels can be rewritten as time derivatives of simpler correlation functions, including itself. Indeed, this observation has been made in recent work.[34]

Here we focus on three types of auxiliary memory kernels that exploit different placements of the time-derivative. The first type replaces the Liouvillian acting on operators that require dynamic sampling and leaves the Liouvillian acting on the static parts intact. In this scheme, remains unchanged. The other auxiliary memory kernels may be expressed as follows,

The second type focuses on replacing the Liouvillian operating on static operators with the time derivative, yielding the following expressions,

The final type replaces all Liouvillians with time derivatives,

It should be noted that Eqs. (Equation 11)–( ?) are exact identities in the context of exact quantum dynamics but *may yield different results when approximate quantum dynamics are employed.*

For clarity in the subsequent discussion, we henceforth refer to the different closures via abbreviations of the form , where , and where denotes the bare expansion and 1, 2, 3 the three types of expansion in the present section. For example, the closure refers to the use of Eqs. (Equation 12) and ( ?) to solve for in Equation 7.

## 4Results

The recent success achieved in using semi- and quasi-classical schemes to calculate the auxiliary kernels required in the self-consistent extraction of memory kernels underscores the importance of understanding the properties of this program in more detail. Consequently, we employ a simple quasi-classical method, namely Ehrenfest dynamics,[16] to obtain the auxiliary memory kernels and study the performance of the Redfield- and NIBA-type projectors and of the different closures available for the kernels. The procedural steps we follow can be summarized as follows:

Calculate the various auxiliary kernels via a dynamical method of choice. Here, we use the approximate Ehrefest approach.

Solve Eqs. (Equation 7) or (Equation 9) iteratively until the relative error becomes negligible. Our threshold is .

Numerically integrate Eq. (Equation 4) subject to the appropriate initial conditions. In the following we use a second-order Runge-Kutta algorithm.

The memory kernel decays to zero for a large sections of parameter space for the SB and other impurity-type models.^{3}

When approximate methods, such as semi- and quasi-classical schemes, are used to calculate the auxiliary kernels, the extracted memory function can accrue errors that are expected to grow with increasing simulation time. Hence, the decay of the memory kernel may not be accurately captured by these methods. Previous applications of the memory function approach have implemented a cutoff time for the memory kernel, after which all its components are set to zero. Sensitivity of the results to this cutoff time will be discussed in Section 4.3. For all results showing only one GQME+MFT curve, the cutoff time, , was chosen at the point where the extracted dynamics reached a plateau of stability. For a more thorough discussion of the Ehrenfest method and its implementation for obtaining the auxiliary kernels, see Appendices Appendix C and Appendix D.

### 4.1Projection Operators

In this section, we restrict our attention to the Ohmic spectral density, a model whose performance has already been studied extensively using the Redfield-type projection operator and closure scheme by Kelly, Markland, and co-workers.[40] The purpose of this section is mainly to provide an analysis of the NIBA-type memory kernels and show the viability of the GQME+MFT approach using both the Redfield- and NIBA-type projectors. Here and in Section 4.2, we show results only for the closure, and postpone the discussion of the closure dependence of the results to Section 4.3.

We first compare the different properties of the Redfield- and NIBA-type memory kernels for a realization of the biased () spin-boson model characterized by weak system-bath coupling (), low temperature (), and a moderately fast bath (). Figure 1 shows a representative set of components of the Redfield-type memory kernel, , in panels (a)-(d), and all components of the NIBA-type memory kernel in panels (e)-(f). Comparison of the two types of type memory kernel reveals the different timescales associated with their decay. Although seemingly noisy at longer times, the Redfield-type memory kernel has a short lifetime (), while the NIBA-type memory kernel decays much more slowly, .

Figure 2 illustrates the dynamics for the parameters used in Figure 1. Despite capturing the correct oscillation frequency and amplitude decay, the Ehrenfest dynamics (green dashes) fails to capture the long-time limit of the populations. Indeed, because of the assumption of a classical bath, the Ehrenfest method is known to violate detailed balance.[61] Instead, the dynamics resulting from both the Redfield- and NIBA-type GQMEs with the closure quantitatively agree (to within graphical accuracy) with the exact results, showing that either method presented here is viable for recovering highly accurate dynamics from approximate dynamics.

To understand the difference in the lifetimes of the Redfield- and NIBA-type memory kernels we recall that the Mori approach to Brownian motion,[62] which focuses on the properties of a massive particle suspended in a bath of lighter particles, relies on the separation of timescales for the dynamics of the heavy and light particles. This separation of timescales is made effective via the projection operator, which must be chosen such that it contains all slow variables associated with the massive particle. Appropriate inclusion of all slow variables in the projector ensures that the memory kernel decays on a shorter timescale than the system dynamics.[63] The Redfield-type projector spans the entire Hilbert space of the system, whereas the NIBA-type projector excludes projections onto coherences, where . While coherences often decay faster than populations, their decay is often slower than bath correlations, as long as the bath dimensionality is large. Hence, the slower time-scale associated with the decay of the coherences induces the slow decay of the NIBA-type memory kernels. This conclusion further suggests that the NIBA-type projector may be most useful in instances of fast system relaxation, such as strong system-bath coupling cases at high temperatures.

It has been suggested that the success of the memory function program, where the auxiliary kernels are calculated via semi- and quasi-classical schemes like Ehrenfest mean-field theory,[40] the momentum-jump solution to the quantum-classical Liouville equation,[39] and the linearized semiclassical initial value representation (LSC-IVR) scheme,[64] relies primarily on the confluence of two important factors. First, the memory kernels are short-lived in comparison to the desired system dynamics. Second, approaches based on semi-classical arguments are more accurate at short times. Considered in tandem, these factors imply that the present scheme can lead to highly accurate short-time memory kernels, thus avoiding problems associated with the long-time dynamics produced by these approximate methods by virtue of fast memory decay. However, the ability of the slowly decaying NIBA-type memory kernel to nearly recover the exact dynamics raises an important question: given the long lifetime of the NIBA-type kernel, how can it remain sufficiently accurate at long times to correct the long-time behavior of the bare quasi-classical dynamics? To answer this question, it is necessary to scrutinize the form of the auxiliary kernels. The NIBA-type auxiliary kernels for closure include two types of correlation functions, and given by Eqs. (Equation 19) and ( ?). Clearly, is the Ehrenfest version of , while involves a new type of correlator that requires the sampling of an additional bath operator, , at . At this point, two possible reasons for the improvement afforded by the NIBA-type approach seem likely. First, it may be that the Ehrenfest method describes the the dynamics of coherences, which are required as input in the auxiliary kernels, better than those of populations. Second, contains exact sampling of the bath operator at , which may encapsulate important information about the system-bath interaction. With the information above, however, it is difficult to decide on which hypothesis is more likely. We return to this discussion in Section 4.3. The previous questions notwithstanding, the ability of the long-lived NIBA-type memory kernel to produce dynamics that are comparably accurate to those obtained via the Redfield-type approach underscores the fact that *a rapidly decaying memory function is not required for the success of the GQME+MFT approach*.

Figure 3 provides a more thorough test of the performance of the NIBA-type projector. Here we compare the NIBA-type GQME+MFT dynamics to exact results for the cases addressed by Kelly *et al.*[40] in their recent work characterizing the performance of the Redfield-type GQME+MFT approach. For convenience, we include the results from the Redfield-type projector as well. We further remark that the focus of these cases is the performance of the present approach to biased systems coupled weakly () to a bath characterized by varying timescales, evident in the range of . As is clear from Figure 3, direct use of the Ehrenfest MFT method consistently leads to incorrect long-time values of the population difference. In agreement with the work of Kelly and co-workers,[40] the Redfield-type GQME+MFT method quantitatively corrects the dynamics in all considered cases. The NIBA-type approach generally provides clear improvement over direct use of MFT, but is slightly less accurate than the Redfield-type GQME+MFT. Regardless, the improvement of the dynamics produced by the NIBA-type projector is remarkable not just because of the fact that the memory function is long-lived, but also because such an approach is not tailored for the weak system-bath coupling limit as is the Redfield-type projector where the benchmark calculations of Figure 3 have been performed.

### 4.2Debye Spectral Density

Due to its slower decay at large frequencies, the Debye spectral density is generally considered a more challenging case for trajectory-based dynamical methods.[65] Here we show that the conclusions drawn from the Ohmic case, namely that the GQME+MFT method can significantly improve the problematic MFT dynamics for weakly coupled, biased systems at low temperatures, are similarly applicable to the Debye case. A few differences are worth mentioning, chief among them that the highly oscillatory nature of the Redfield-type memory kernels for this spectral density generally means that a larger number of trajectories to achieve convergence are required.

Panels (a)-(d) of Figure 4 show the components for the Redfield-type memory kernel, while panels (e)-(f) show all components of the NIBA-type memory kernel. While the NIBA-type memory kernels do not show any obvious differences from their Ohmic counterparts, the Redfield-type memory kernel displays recurrent beating alongside overall decay. The presence of this much stronger oscillatory behavior requires that the dynamics of high frequency modes in the Ehrenfest procedure be treated more accurately than is necessary in the Ohmic case. We discuss this issue in more depth in the next section, where we explore the convergence properties of the difference closures.

Figure 5 presents some illustrative examples of the capability of the GQME+MFT approach to yield accurate dynamics for the biased SB model at low temperatures, over a wide range of . Panel (a), which corresponds to an *unbiased* case characterized by weak system-bath coupling at low temperature () and an intermediate bath frequency (), shows nearly perfect agreement between the Ehrenfest and exact dynamics. Both the Redfield- and NIBA-type GQME+MFT approaches are able to recover the remarkable agreement between the Ehrenfest and exact dynamics. Panels (b)-(d) correspond to biased cases, spanning a wide range of bath frequencies ( 0.25, 5.0) and temperatures (0.5, 50.0). As expected, the bare Ehrenfest method leads to incorrect long-time limits for all three biased cases. As in the Ohmic case, both the Redfield- and NIBA-type approaches yield results in almost quantitative agreement with the exact dynamics. Slight deviations are evident, as in panel (b), where the NIBA-type GQME slightly underestimates the long-time limit of the population difference. Perhaps the most difficult case for the current approach, panel (d), shows that the NIBA-type GQME+MFT treatment leads to overly damped oscillations at long times, whereas the Redfield-type approach yields results in near quantitative agreement with the exact dynamics. In short, the results in Figure 5 illustrate the robustness of the approach for weak-coupling cases over a wide range of bath frequencies and temperatures.

### 4.3Memory Kernel Closures and Dynamics

In Section 3, we introduced *eight* different closures of the memory kernels. These include two subsets consisting of the -forward and -backward closures, which are further subdivided into the bare expansion ( and ) and three expansions that use numerical derivatives of the simulated correlation functions, (– and –). While the resulting dynamics do not differ when the auxiliary kernels are calculated via exact methods, the same claim is not necessarily true when using approximate dynamics. Here, we continue to use the Ehrenfest method to illustrate the sensitivity of the results that occur across the spectrum of closures.

To inform the discussion on the properties of different closures of either the Redfield- or NIBA-type memory kernels, we first provide an overview of the underlying types of correlation functions that are employed in the calculation of the auxiliary kernels. These are summarized in Eqs. (Equation 19)–( ?) of Appendix B. For convenience, we reproduce these expressions, within the Ehfenfest approximation, below,

where and .

Inspection of Eqs. (Equation 14)–( ?) reveals that there are two main types of functions containing bath operators: those that require their sampling exclusively at [Eqs. (Equation 14), ( ?) and ( ?)], and those containing both statically and dynamically sampled bath operators [Eq. ( ?), ( ?), and ( ?)]. We recall the important fact that at , the Ehrenfest method is exact and the accuracy of the method diminishes with increasing simulation time.[66] What is not clear, however, is whether the accuracy associated with the *dynamical* sampling of a single system operator is the same as that of a product of system and bath operators, as is required in Eqs. ( ?), ( ?), and ( ?). Because the bath is treated classically, dynamical sampling of bath operators may indeed accrue larger errors. Since the classical approximation is most problematic for high frequency modes, this problem may be exacerbated by fast baths characterized by broad spectral densities, namely the Debye spectral density. Instead, when bath operators are sampled statically, as is the case in Eqs. (Equation 14), ( ?), and ( ?), the weighting of trajectories of the correlation functions is captured exactly. One may also distinguish the correlation functions above on the basis of sampling of *distinct* bath operators not normally included, whether explicitly or implicitly, in . Naively, one may suppose that Eqs. ( ?)–( ?) contain information distinct from that in contained in , but the Ehrenfest evolution algorithm requires sampling of , which contributes a dynamic component to the system’s bias energy . Consequently, only Eqs. ( ?) and ( ?), which sample , contain *distinct* information about the system which is not already included in the calculation of . Indeed, these correlation functions may contain additional information about the system-bath interaction via statically sampled bath operator, , that facilitate the improvement over the bare Ehrenfest dynamics afforded by the memory function approach.

A distinct source of error that may affect the accuracy of closures that implement numerical time-derivatives lies in the accuracy of the time-derivatives themselves. If the correlation functions calculated via the Ehrenfest procedure are sufficiently smooth and well converged and the time-step is sufficiently small, this error can be expected to be minimal. However, correlation functions containing bath operators that require sampling at finite times tend to be highly oscillatory, especially for fast baths, which may lead to less accurate results.

Armed with these considerations, it is possible to explore the differences associated with the different closures of the auxiliary kernels. Because the Redfield- and the NIBA-type projectors require different combinations of the aforementioned correlators, Eqs. (Equation 14)–( ?), as input for their auxiliary kernels, we discuss the behaviors of the different closures for the two projectors separately.

Focusing first on the Redfield-type kernel closures, we assess the effect of the -backward and -forward approaches by focusing first on the and closures. Inspection of Eqs. ( ?) and (Equation 10) reveals that the only difference between the -forward and -backward closures lies in the fact that contains the time-evolved bath operator , whereas requires sampling of the static bath operator, . Consistent with the above discussion, we may expect that the closure will lead to less accurate results than . As Figure 6 shows, the difference between the two closures is minimal. To understand the smallness of the difference between these two closures, it is sufficient to consider that, while each closure has a different form for , both closures share the same form for , which requires the sampling of bath operators both at and at finite times. Hence, any error associated with the explicit inclusion of time evolved bath operators would affect both closures, and , and any benefit derived from the exclusive sampling of bath operators is also maintained in the form for . Further, the similarity in performance of the and closures indicates that the error associated with the dynamical sampling of products of system and bath operators is often similar to the error associated with the exclusive sampling of system operators. Because the difference between the results of the -forward and -backward closures is small, we henceforth exclusively address the differences among the -backward closures, , , , and .

Consideration of the remaining three -backward closures, , , and , likewise requires close scrutiny of the types of correlation functions that are used in each. First we note that, while the form of in the closure avoids the sampling of bath operators at finite time, still samples bath operators at finite times. Instead, the closure completely avoids the sampling of bath operators at , while still benefiting from sampling of static bath operators for both auxiliary kernels. In contrast, the and closures explicitly avoid sampling of static bath operators, other than the density operator for the bath. As is evident from panels (a) and (b) in Figure 7, the closure performs as well or better than the closure. However, because *the closure leads to equations that are easier to converge and auxiliary kernels that are easier to calculate*, it should be preferred over the closure.

Inspection of panels (c) and (d) of the same figure shows that the and essentially recover the bare Ehrenfest behavior. Thus, we have demonstrated the remarkable fact that different closures for the memory function, all of which are exact when implemented with exact input, can yield markedly different results when combined with approximate dynamical input. The analytical proof that, for example, the (and ) closure must yield correlators that are identical to the use of the bare input dynamics is provided in the companion paper.[67] Further, the fact that the closure also recovers the bare Ehrenfest dynamics clearly indicates that the first of the two criteria specified in Ref. is satisfied by the Ehrenfest method. Indeed, the numerical data show that within the Ehrenfest approach the action of the Liouvillian acting on a dynamically sampled operator is equivalent to the numerical time-derivative of the analogous correlation function. In addition, given the violation of the second criterion (), it is not surprising that the and closures yield GQME dynamics that are *distinct* from the the results of direct application of the Ehrenfest method. However, the violation of the second criterion does not explain the reason for marked improvement in the dynamics afforded by the and closures. These results also lend additional credence to the claim that the success of the GQME+MFT approach does *not* rely on the short-time accuracy of Ehrenfest dynamics. Further, they illustrate that improvement within the memory formalism over the bare quasi-classical theory depends sensitively on the correlation functions calculated as input for the auxiliary memory kernels. In Section 4.1 we suggested that the Ehrenfest method might capture the dynamics of coherences more accurately than that of the populations, and that the self-consistent extraction of the memory kernel would include corrections to the population dynamics afforded by the ostensibly more accurate coherence dynamics. However, the recovery of Ehrenfest dynamics by closures and , which also use the dynamics of coherences, implies that this cannot be the root cause of the improvement of dynamics within the memory function approach. Instead, the more likely explanation is that the sampling of static bath operators in Eqs. ( ?) and ( ?) contributes important information about the system-bath interaction, which leads to far greater accuracy in the extracted memory kernels themselves.

The differences in the dynamics resulting from the different closures are also evident in the extracted memory kernels. Direct comparison of the Redfield-type memory kernels for cases characterized by the Debye spectral density fails to reveal much, since their highly oscillatory behavior obfuscates subtle differences among the memory kernels obtained from different closures. However, the fast decay of the Ohmic spectral density, which results in quickly decaying memory kernels with minor oscillations, makes discerning qualitative and quantitative differences between the extracted kernels possible. Figure 8 compares a representative set of Redfield-type kernel elements arising from different closures. As is clear from the figure, all closures agree within numerical and sampling error in their values. However, the and closures display a stronger oscillatory behavior, in contrast to the and closures, which correctly recover accurate dynamics (see Figure 2). We further note that the difference in behavior is greatest at intermediate times.

The previous discussion suggests that the main factor leading to highly accurate memory kernels is the *exact* sampling of specific static bath operators. A corollary question arises: do all correlation functions with statically sampled bath operators lead to a similar improvement? After all, Eq. ( ?) samples at , but its use in closure does not lead to any improvement over the bare Ehrenfest dynamics. This suggests either that the exact sampling of exclusively static bath operators adds an important correction to the auxiliary kernels, or that correlation functions that sample are not as important as those that sample . The idea that is neither sampled explicitly nor implicitly (via the Ehrenfest evolution protocol) in the calculation of provides some support to the latter claim. It is also fair to ask whether one can similarly benefit from “improperly” Wigner-transformed bath operator products sampled at . This question becomes particularly important when a functional form for the density operator of the bath is either not available or challenging to obtain. To see the importance of properly including the terms in the Wigner transformation, we take an approximate form for the Wigner transform of the product (defined in Appendix D). Our approximation truncates the Moyal expansion for the Wigner transform of a product of operators at zeroth order in , neglecting the second term (containing ) on the right side of Eq. (Equation 29). Results for this approximation are shown in Figure 9. As is evident in the figure, the benefits in closures and that originally led to the quantitative agreement between the GQME+MFT and exact dynamics are eliminated. Instead, the final result, while not unphysical, is not better than the standard Ehrenfest result. Since the and closures do not contain this neglected term, the results in panels (c) and (d) of Figure 9 are the same as those in panels (c) and (d) of Figure 7. This result underscores the importance of proper sampling of *all* contributions arising from the Wigner transform of operator products.

In comparison to the Redfield-type closures, the NIBA-type kernels only contain two types of correlation functions, and . As mentioned in the discussion of the Redfield-type closures, contains the same information as Ehrenfest version of , whereas contains *exact* information about the system-bath interaction at . For this reason, we expect the and closures to yield significantly better dynamics than the and closures, which should simply recover the Ehrenfest dynamics.[67] Indeed, as Figure 10 shows, the and closures are able to quantitatively correct the Ehrenfest dynamics. In contrast to the Redfield case, the and NIBA-type auxiliary kernels require the sampling of the same bath operators, which explains the lack of difference in the behaviors of the two closures.

Of paramount importance to the success of the memory approach is the finite lifetime of the memory kernel. For the GQME+MFT implementation used here, we have chosen a cutoff time for the memory kernel, , which lies in a stability plateau alluded to in Section 4. The range of the stability plateau can depend sensitively on the regime of parameter space explored. Figure 11 shows the dependence of the GQME+MFT dynamics for the Redfield- and NIBA-type projectors on the specific value of used. Panels (a) and (b), corresponding to a fast bath, high temperature, biased case, show the greatest sensitivity of the GQME dynamics to the exact cutoff time, . In contrast, the results in panels (c) and (d), which correspond to lower temperatures, are more stable. Despite the slight sensitivity of the results on the choice of for the examples shown, the GQME dynamics are clearly robust. We remark, however, that there are regions of parameter space for which this stability plateau is short-lived or nonexistent. In those cases, the present approach is clearly not appropriate.

Finally, we address the convergence properties of the GQME dynamics with respect to the time-step of the memory kernel and the GQME evolution algorithm. As Figure 12 shows, decreasing the time step used in the extraction of the memory kernel can greatly alter the accuracy of the GQME dynamics. The closures most sensitive to the time step are and closures. The sensitivity of the closure to the time step is not difficult to understand, since for large the second numerical time derivative becomes noisy, leading to the growing oscillations in panel (d). Because closure shown in panel (a) contains the correlation functions and , which require sampling of bath operators at and at finite times and are highly oscillatory functions, a smaller time step is required to achieve sufficiently accurate memory kernels. It is important to note that, as stated before, the and closures, when fully converged with respect to bath realizations and time step, yield identical results that greatly improve the Ehrenfest results, while closures and are only capable of recovering the Ehrenfest dynamics.

## 5Conclusions

In this paper we have developed a method to obtain the nonequilibrium population and coherence dynamics based on the Mori formalism. Our approach is general and, depending on the choice of projector, can treat arbitrary single-time nonequilibrium populations and coherences as well as more complicated dynamical objects, such as multi-time, equilibrium and nonequilibrium correlation functions. We have shown that use of the Redfield-type projector recovers the conventional NZ treatment previously used by Shi and Geva[31] and Kelly, Markland, and coworkers[39] in the context of the SB model, and Rabani and coworkers[34] for more general models.

While previous applications of the GQME+semi-classical approach[64] have been limited to the Redfield-type projector and have focused on the improvement over the bare semi-classical dynamics that the memory function formalism can afford, we have systematically explored the sensitivity of the results to the choice of projector and the type of closure employed. In doing so, we find two important facts. First, slowly decaying memory kernels, often observed when using the NIBA-type projector, do not result in an inaccurate description of the GQME dynamics. This demonstrates that the success of the GQME+semi-classical approach is *not* a function of the short-time accuracy of the approximate method used to calculate the auxiliary kernels. Second, we identify the types of closures that consistently lead to improvements over the bare semi-classical dynamics (, , , and ) and attribute this improvement, in part, to the sampling of static bath operators, , which do not appear in the evaluation of the approximate bare populations. Just as importantly, we identify the types of closures that recover the Ehrenfest dynamics (, , , and ). Our findings also provide numerical confirmation of the analytical proof included in the companion paper,[67] which indicates that use of the and closures *can only* recover the level of dynamics used to calculate the auxiliary kernels.

Finally, we remark that the Mori-based formulation furthered in this work provides a flexible framework to accurately study problems that go beyond the scope of nonequilibrium dynamics for SB-type models. For instance, the Mori formalism can easily address equilibrium and multi-time correlation functions in systems coupled to harmonic *and* anharmonic baths, as well as problems where the system-bath distinction is absent. Work in this latter direction will be pursued in future papers.

## 6Acknowledgements

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

## AFourier-Laplace Analysis of Closures

We begin by introducing the Fourier-Laplace transform of a time-dependent function ,

Its first and second time-derivatives take the following form,

The kernel expansions in Eqs. (Equation 7) and (Equation 9) can be rewritten as follows,

In their original paper, Shi and Geva[31] derived the following identity for the Redfield-type (thermal) projector,

where is the Liouvillian corresponding to the system-bath interaction for the SB model . The second line is a simple extension of the derivation provided by Shi and Geva. These identities allow for the exact rewriting of the “-surrounded” projected propagator,

Replacing these expressions for the projected propagator in Eq. (Equation 5) followed by use of the Dyson decomposition leads to the following *three*- rather than *two*-membered expansions,

where the second auxiliary kernel also contains the projected propagator,

In Fourier-Laplace space, the previous expressions for the memory kernel become,

In the time-domain, the second auxiliary kernels may be expanded as follows,

where . Transforming Eqs. (Equation 17) and ( ?), we may solve for

Substitution of Eqs. (Equation 18) and ( ?) into Eqs. (Equation 16) and ( ?) for the -forward and -backward closures yields,

which are clearly equivalent to Eqs. (Equation 15) and ( ?), implying that the three-membered closures are equivalent to the two-membered closures.

## BExpressions for Auxiliary Memory Kernels

Here we provide explicit expressions for the components of the memory kernel using the Redfield- and NIBA-type projection operators. Before going further, however, we introduce for notational clarity the following correlation functions.

where and indicate symmetrized (anticommutator) or antisymmetrized (commutator) bath products.

Using the Redfield-type projector,

where , , and , the elements of the memory kernels take the following forms,

where and .

Using the NIBA-type projector,

where , and , the elements of the memory kernels take the following forms,

We employ a notation where some indices are squared since the functions are labelled using the indices corresponding to the operators, which are related to the operators in the following way: and .

For the projector above to be truly of NIBA-type, should be . Instead, has contributions of first and second order in . Indeed, the proper NIBA-type projector has the following form,

where .

## CInitial Conditions in the Ehrenfest method

In an open quantum system where a subsystem interacts weakly with a heat bath, the Ehrenfest method[12] treats the subsystem quantum mechanically and the bath classically. The validity of this approximation relies on two important assumptions: correlations between the system and bath are negligible, and the characteristic energy of the bath is smaller than the other energy scales in the problem, justifying the use of classical mechanics for the evolution of the bath.

The Ehrenfest method has been derived from complementary wavefunction[16] and density matrix formulations.[59] Because of its clarity, the derivation based on the density matrix and the quantum-classical Liouville equation has garnered much attention in the last decade. The density matrix formulation only requires that the subsystem and bath density matrices have norms equal to unity. However, the lack of restriction of the subsystem density matrix to pure states results in ambiguities in its implementation. To illustrate the source of the ambiguity, we focus on the following correlation function for the spin-boson model

where , and the bath initial condition has unit trace, . This corresponds to a nonequilibrium initial condition where the system is initially in a superposition of coherences. Immediately it is clear that has zero norm. To remedy this, one may take advantage of the linearity of the problem and rewrite Eq. (Equation 25) into sums of correlation functions with proper initial conditions,

where .

Returning to the wavefunction-based derivation of the Ehrenfest approach requires that any subsystem initial condition correspond to a pure state. Referring again to Eq. (Equation 25), we may rewrite the trace over the subsystem in the eigenbasis of the subsystem’s initial condition,

Although there are other ways of choosing a pure state so as to evaluate the above correlation function, the wavefunction formulation avoids the ambiguity fostered by the density matrix derivation.

Figure 13 shows the calculation of the nonequilibrium population dynamics given different system initial conditions, . For the initial conditions corresponding to the Pauli matrices, we implement three different decompositions given by Eqs. (Equation 26), ( ?), and ( ?), labeled , , and , respectively. We also include the results for two different decompositions for the initial condition , labeled and . As is clear from panels (a) and (b), all density matrix based decompositions agree in their short- and long-time limits, but disagree in their descriptions of intermediate-time behavior. More importantly, all density matrix decompositions disagree with the wavefunction-based result. Panel (c) shows the ability of the density matrix-based approach to recover the wavefunction-based result when the initial condition corresponds to a population rather than a coherence, but decomposition underscores the problems associated with the lack of uniqueness in the density matrix based approach.

## DEhrenfest Method: Correlation Functions

Unlike previous implementations of the Ehrenfest method, we are interested in two time correlation functions rather than nonequilibrium single quantity dynamics,

where () is a generic system (bath) operator. Under the quasi-classical approximation of Wigner dynamics, we may rewrite the above correlation function as follows,

where the superscript denotes the Wigner transform of the operator and is the set of all classical variables.

For the auxiliary memory kernels of the spin-boson model, the following Wigner transforms are necessary,

and, using the Moyal bracket for products of operators,[68]