Approximate but accurate quantum dynamics from the Mori formalism: I. Nonequilibrium dynamics

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–8 their computational cost makes them impractical for realistic multidimensional systems. Conversely, approximate methods, whether perturbative^9–11 or based on quasi-^12–16 or semi-classical^17–23 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,25 provides a useful starting point for a plethora of methods based on generalized quantum master equations (GQMEs).

The NZ equation, which may be derived via the projection operator technique,^26,27 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.^{9–11,28–30}

A major difficulty in the calculation of the memory kernel lies in the fact that its dynamical evolution involves the “projected” propagator, $ei(1−P)Lt$⁠, where $P$ is the projector that defines the reduced dynamics. To circumvent the problem of projected dynamics, Shi and Geva^31–33 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 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,3 Shi and Geva have shown that the memory kernel for the spin-boson (SB) model can decay up to 10 times faster than the system dynamics,³¹ 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–38 Just as importantly, applications of the method have successfully used semi- and quasi-classical theories to calculate the auxiliary kernels.³² In particular, Kelly, Markland, and coworkers have illustrated the impressive accuracy and robustness of this approach in both model and realistic problems.^39–41 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 SB 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 that lead to improvements in accuracy of approximate dynamics via the memory function formalism, as observed in Refs. 32 and 39–41, is lacking. Finally, the convergence properties of different versions of the auxiliary 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 cumbersome to obtain within this framework. Despite this difficulty, recent work has generalized the NZ equation to multi-time correlation functions.⁴² 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⁴³ and non-linear spectroscopies,⁴⁴ and the calculation of chemical rate constants^21,45 and kinetic coefficients,⁴⁶ to name a few examples. For this reason, in this paper (the first in a 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 forthcoming paper (the second in this series),⁶⁹ we present 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–49 and quantum fluids.^50–53 

The structure of this paper is as follows: In Sec. II, 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 II also introduces the SB model and proposes two types of projection operators for this model. Section III discusses different closures of the memory kernel based on different placements of the projection operator $Q=1−P$ and the use of time derivatives. To illustrate the arguments related to convergence, we implement the mean field Ehrenfest method, as first proposed in Ref. 40, to obtain the auxiliary 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 Sec. IV, we 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 Sec. III. In Sec. V, we conclude.

For illustrative purposes, we focus on the 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.⁵⁴ 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–49 and liquids.^50–53 We reserve these applications for later work.

The SB Hamiltonian takes the form H = H_S + H_B + V. 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 σ_i corresponds to the ith Pauli matrix. The bath part of the Hamiltonian consists of independent harmonic oscillators,

where P_k and Q_k are the mass weighted momenta and coordinates for the kth harmonic oscillator, respectively, and ω_k is the frequency of the kth 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 c_k is the coupling constant between the system and the kth oscillator, and α = ± 1, depending on the definition of the model. The spectral density, J(ω), 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 J(ω) is proportional to ω as ω → 0. These are the standard Ohmic spectral density¹¹ characterized by an exponential cutoff, and the Debye spectral density characterized by a Lorentzian cutoff

Here the cutoff frequency ω_c determines the correlation time of the bath at sufficiently high temperatures.⁵⁵ The reorganization energy, $λ=π−1∫0∞dωJ(ω)/ω$⁠, is a measure of the strength of the system-bath coupling. In the case of the Ohmic spectral density, the Kondo parameter, η = πξ/2 = πλ/ω_c, 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 ρ_B = e^−βH_B/Tr_B[e^−βH_B] is the equilibrium density operator for the uncoupled bath, β = [k_BT]⁻¹ is the inverse thermal energy, and the initial condition |1〉〈1|ρ_B corresponds to a Frank-Condon transition.

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

where $L=[H, …],$$P$ is the projection operator that defines the subsystem whose dynamics we seek, and $Q=1−P$ is the complementary projection operator. This equation is general and can be employed within both the NZ and Mori approaches. For instance, applying Eq. (8) on an operator $Aˆ$⁠, yields the Mori EOM for that operator. Conversely, taking the Hermitian conjugate of Eq. (8), applying it on the initial density matrix of the system and bath, ρ₀, and acting the projection operator $P$ from the left followed by a trace over the bath degrees of freedom yields the NZ equation for the system’s RDM.

To date, essentially all work on the self-consistent expansion of the memory kernel has employed the thermal (Argyres-Kelley) projection operator $P=RBTrB[…]$⁠,^{31–36,39,40} where $RB$ is a bath operator with unit trace and Tr_B[…] 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,⁵⁶ as $P=∑i|Ai〉〉〈〈RBiAi|$⁠, where {A_i} contains all outer product states spanning the system. For the SB model, A_i ∈ {|1〉〈1|, |2〉〈1|, |1〉〈2|, |2〉〈2|}.

Using the thermal type projector, applying Eq. (8) to |A_k〉〉, and closing on the left with $〈〈RBj′Aj|$⁠, where $RBj′$ is again a bath operator with unit trace that may be different from $RBj$ in the projection operator, yields the following EOM for system observables:

where $Xjk=〈〈RBjAj|L|Ak〉〉$ is a static rotation matrix, $Cjk(t)=〈〈RBj′Aj|eiLt|Ak〉〉=Tr[(RBj′)†Aj†Ak(t)]$ corresponds to nonequilibrium averages of populations and coherences with all possible factorizable initial conditions, and $Ijk(t)=i〈〈RBj′Aj|QeiLQtL|Ak〉〉$ is the so-called the inhomogeneous term. The elements of the memory kernel are given by

When $RBj′=RBj$⁠, the inhomogeneous term disappears, $I(t)=0$⁠. Often, $RBj=RBj/TrB[RBj]$ is chosen such that $RBj∈{e−βHB,e−β(HB±α∑kckQˆk)}$⁠, 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 $ρ(0)=ρS(0)RBi$⁠, where ρ_S(0) is an arbitrary system operator and $RBi$ 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 R_Bi = exp[ − βH_B] as Redfield-type.^9,10 The reason for this name is that truncation of the memory kernel at second order in $Q$ is equivalent to a second-order perturbative expansion of the memory kernel with respect to the system-bath coupling, which corresponds to Redfield theory.^9,10 The choice for $RB,i′$⁠, however, remains flexible.

An important alternative for the projection operator consists of restricting the set {A_i} to the system populations B_i = |i〉〈i|, and choosing $RBi = exp[−β(HB + (−1)i+1α∑kckQˆk)]$⁠. When using this projection operator, a similar second-order truncation of the memory kernel with respect to $Q$ 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,57 Accordingly, we hereafter refer to this projector as NIBA-type. As in the case of the Redfield-type projector, the choice for $RB,i′$ determines whether the inhomogeneous term is zero or finite.

As mentioned in Sec. I, the main difficulty associated with the memory kernel, Eq. (10), is the presence of the projected propagator, $eiLQt$⁠. 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,33 Not surprisingly, Eq. (8) can also be derived using the Dyson identity, Eqs. (11) and (12). 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. (11) and (12).

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 and type of auxiliary kernels employed.³³ 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 are only a limited number of kernel expansions and ways of expressing the auxiliary kernels that can yield numerically distinct approximate memory kernels from the self-consistent solution of the resulting integral equations.

Substitution of Eq. (11) into Eq. (10), where $A=iL$ and $B=−iPL$⁠, yields the following integral equation for the memory kernel:

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

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

A second set of closures makes use of the identity $eiQLtQ=QeiLQt$⁠. To indicate that the $Q$ is to the right of the Liouvillian $L$ in the propagator, we have used the superscript f (indicating a “forward” placement). Inserting the previous identity in Eq. (10) and expanding the projected propagator using the Dyson identity given by Eq. (12), with $A=iL$ and $B=−iLP$⁠, yields

We note that $K(1)(t)$ has the form given in Eq. (14), and $K(3f)(t)$ has the following form:

It bears remarking that Eqs. (16) and (13) differ in the placement of $K(t)$ under the integral and Eqs. (17) and (15) differ in whether $QL$ or its Hermitian conjugate act on operators that require sampling only at t = 0 or at finite times. We further note that the auxiliary kernels given by Eqs. (14), (15), and (17) no longer require the use of projected dynamics and can be simulated directly.

Because the auxiliary kernels given by Eqs. (14), (15), and (17) require sampling of additional bath operators at t = 0 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 $C(t)$⁠, 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 $C(t)$ itself. Indeed, this observation has been made in recent work.^34–38,58

Here we focus on three types of auxiliary 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, $K(3b)(t)$ remains unchanged. The other auxiliary 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. (18)–(26) 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 c(xy), where x ∈ {f(orward), b(ackward)}, and y ∈ {0, 1, 2, 3} where 0 denotes the bare expansion and 1, 2, 3 the three types of expansion in the present section. For example, the cb2 closure refers to the use of Eqs. (21) and (22) to solve for $K(t)$ in Eq. (13).

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,59 to obtain the auxiliary 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 thus:

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

2. Solve Eqs. (13) or (16) iteratively until the relative error becomes negligible. We define the relative error, R.E., as the maximum absolute difference between two subsequent iterations of the memory kernel in the self-consistent scheme, i.e., $R.E.=max[abs[Kn+1(t)−Kn(t)]]$⁠. Our threshold is 10⁻¹⁰.

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

The memory kernel decays to zero for large sections of parameter space for the SB and other impurity-type models.⁶⁰ We refer to the time scale that determines this decay as the memory lifetime. As mentioned in the Introduction, the computational efficiency of the memory function approach depends sensitively on this lifetime.

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 Sec. IV C. For all results showing only one GQME+MFT curve, the cutoff time, τ_c, was chosen at the point where the extracted $C(t)$ dynamics reached a plateau of stability. Importantly, our implementation of the Ehrenfest method employs Wigner-transformed initial conditions and is exact at t = 0, and, as is characteristic of all methods stemming from the semiclassical hierarchy, decreases in accuracy for increasing simulation time.⁶¹ For a more thorough discussion of the Ehrenfest method and its implementation for obtaining the auxiliary kernels, see  Appendices C and  D.

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 cb0 closure scheme by Kelly, Brackbill, and Markland.⁴⁰ 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 Sec. IV B, we show results only for the cb1 closure and postpone the discussion of the closure dependence of the results to Sec. IV C.

We first compare the different properties of the Redfield- and NIBA-type memory kernels for a realization of the biased (ε = 1) SB model characterized by weak system-bath coupling (ξ = 0.1), low temperature (β = 5.0), and a moderately fast bath (ω_c = 2.5). Fig. 1 shows a representative set of components of the Redfield-type memory kernel, $Kx2(t)$⁠, in panels (a)-(d), and all components of the NIBA-type memory kernel $K(t)$ in panels (e) and (f). Comparison of the two types of memory kernel reveals the different time scales associated with their decay. Although seemingly noisy at longer times, the Redfield-type memory kernel has a short lifetime (τ_c ∼ 2), while the NIBA-type memory kernel decays much more slowly, (τ_c > 15).

Fig. 2 illustrates the dynamics for the parameters used in Fig. 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.⁶² Instead, the dynamics resulting from both the Redfield- and NIBA-type GQMEs with the cb1 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,⁶³ which focuses on the properties of a massive particle suspended in a bath of lighter particles, relies on the separation of time scales for the dynamics of the heavy and light particles. This separation of time scales 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 time scale than the system dynamics.⁶⁴ For the SB and other impurity-type models, the Redfield-type projector spans the entire Hilbert space of the system, whereas the NIBA-type projector excludes projections onto coherences, |i〉〈j| where i≠j. 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,⁴⁰ the momentum-jump solution to the quantum-classical Liouville equation,³⁹ and the linearized semiclassical initial value representation (LSC-IVR) scheme,⁶⁵ 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 cb1 include two types of correlation functions, $qnm(00)(t)$ and $qnm(10a)(t)$ given by Eqs. (B1) and (B4). Clearly, $qnm(00)(t)$ is the Ehrenfest version of $C(t)$⁠, while $qnm(10a)(t)$ involves a new type of correlator that requires the sampling of an additional bath operator, $ζW=−α∑jcjPjtanh(βωj/2)/ωj$⁠, at t = 0. 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 dynamics of coherences, which are required as input in the auxiliary kernels, better than those of populations. Second, $qnm(10a)(t)$ contains exact sampling of the bath operator ζ^W at t = 0, 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 Sec. IV C. 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.

Fig. 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.⁴⁰ 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. The focus of these cases is the performance of the present approach to biased systems coupled weakly (λ = 2ξ/ω_c < 1) to a bath characterized by varying time scales, evident in the range of ω_c. As is clear from Fig. 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,⁴⁰ the Redfield-type GQME+MFT method quantitatively corrects the dynamics in all cases considered. 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 Fig. 3 have been performed.

Due to its slower decay at large frequencies, the Debye spectral density is generally considered a more challenging case for trajectory-based dynamical methods.⁶⁶ 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 is required to achieve convergence.

Panels (a)-(d) of Fig. 4 show the components $Kx2(t)$ for the Redfield-type memory kernel and panels (e) and (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 Subsection IV C, where we explore the convergence properties of the difference closures.

Fig. 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 ω_c. Panel (a), which corresponds to an unbiased case characterized by weak system-bath coupling at low temperature (β = 50.0) and an intermediate bath frequency (ω_c = 1.0), 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 (ω_c = 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 Fig. 5 illustrate the robustness of the approach for weak-coupling cases over a wide range of bath frequencies and temperatures.

In Sec. III, we introduced eight different closures for the memory kernels. These include two subsets consisting of the $Q$-forward and $Q$-backward closures, which are further subdivided into the bare expansion (cf0 and cb0) and three expansions that use numerical time derivatives of the simulated correlation functions, (cf1–cf3 and cb1–cb3). 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. (B1)–(B6) of  Appendix B. For convenience, we reproduce these expressions, within the Ehrenfest approximation, below,

where $VBW=∑kckxk$ and $ζW=−α∑jcjPjtanh(βωj/2)/ωj$⁠.

Inspection of Eqs. (27)–(32) reveals that there are two main types of functions containing bath operators: those that require their sampling exclusively at t = 0 [Eqs. (27), (29), and (30)], and those containing both statically and dynamically sampled bath operators [Eqs. (28), (31), and (32)]. We recall the important fact that at t = 0, the Ehrenfest method is exact and the accuracy of the method diminishes with increasing simulation time.⁶¹ 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. (28), (31), and (32). 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. (27), (29), and (30), the t = 0 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 $C(t)$⁠. Naively, one may suppose that Eqs. (28)–(32) contain information distinct from that contained in $C(t)$⁠, but the Ehrenfest evolution algorithm requires sampling of $VBW(t)$⁠, which contributes a dynamic component to the system’s bias energy $ε↦ε+VBW(t)$⁠. Consequently, only Eqs. (30) and (32), which sample ζ^W, contain distinct information about the system which is not already included in the calculation of $C(t)$⁠. Indeed, these correlation functions may contain additional information about the system-bath interaction via the statically sampled bath operator, ζ^W, that facilitates 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. (27)–(32), 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 $Q$-backward and $Q$-forward approaches by focusing first on the cb0 and cf0 closures. Inspection of Eqs. (15) and (17) reveals that the only difference between the $Q$-forward and $Q$-backward closures lies in the fact that $K(3f)(t)$ contains the time-evolved bath operator $VBW(t)$⁠, whereas $K(3b)(t)$ requires sampling of the static bath operators, $VBW(0)$ and ζ^W(0) (see also Eqs. (B8)–(B10)). Consistent with the above discussion, we may expect that the cf0 closure will lead to less accurate results than cb0. As Fig. 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 $K(3)(t)$⁠, both closures share the same form for $K(1)(t)$⁠, which requires the sampling of bath operators both at t = 0 and at finite times. Hence, any error associated with the explicit inclusion of time evolved bath operators would affect both closures, cb0 and cf0, and any benefit derived from the exclusive sampling of t = 0 bath operators is also maintained in the form for $K(1)(t)$⁠. Further, the similarity in performance of the cb0 and cf0 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 $Q$-forward and $Q$-backward closures is small, we henceforth exclusively address the differences among the $Q$-backward closures, cb0, cb1, cb2, and cb3.

Consideration of the remaining three $Q$-backward closures, cb1, cb2, and cb3, likewise requires close scrutiny of the types of correlation functions that are used in each. First we note that, while the form of $K(3b)(t)$ in the cb0 closure avoids the sampling of bath operators at finite time, $K(1)(t)$ still samples bath operators at finite times. Instead, the cb1 closure completely avoids the sampling of bath operators at t≠0, while still benefiting from sampling of static bath operators for both auxiliary kernels. In contrast, the cb2 and cb3 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 Fig. 7, the cb1 closure performs as well or better than the cb0 closure. However, because the cb1 closure leads to equations that are easier to converge and auxiliary kernels that are easier to calculate, it should be preferred over the cb0 closure.

Inspection of panels (c) and (d) of the same figure shows that the cb2 and cb3 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 cb3 (and cf3) closure must yield correlators that are identical to the use of the bare input dynamics is provided in Ref. 67. The slight discrepancy between the cb3 and the Ehrenfest results may be attributed to the finite precision of the first and second numerical time derivatives of $C(t)$⁠. Further, the fact that the cb2 closure also recovers the bare Ehrenfest dynamics clearly indicates that the first of the two criteria specified in Ref. 67, i.e.,

where ρ corresponds to an initial, normalized density matrix and A and B are arbitrary operators, 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 (⁠$[LEh,(eiLt)Eh]≠0$⁠), it is not surprising that the cb0 and cb1 closures yield GQME dynamics that are distinct from the results of direct application of the Ehrenfest method. However, the violation of the second criterion does not provide a reason for the marked improvement in the dynamics afforded by the cb0 and cb1 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 kernels. In Sec. IV A 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 cb2 and cb3, 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. (29) and (31) 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. Fig. 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 t = 0 values. However, the cb2 and cb3 closures display a stronger oscillatory behavior, in contrast to the cb0 and cb1 closures, which correctly recover accurate dynamics (see Fig. 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. (28) samples $VBW(0)$ at t = 0, but its use in closure cb3 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 $VBW(0)$ are not as important as those that sample ζ^W. The idea that ζ^W is neither sampled explicitly nor implicitly (via the Ehrenfest evolution protocol) in the calculation of $C(t)$ 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 t = 0. 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 ρ_BV_B (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 ζ^W) on the right side of Eq. (D5). Results for this approximation are shown in Fig. 9. As is evident in the figure, the benefits in closures cb0 and cb1 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 cb2 and cb3 closures do not contain this neglected term, the results in panels (c) and (d) of Fig. 9 are the same as those in panels (c) and (d) of Fig. 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 (see Eqs. (B12)–(B14)) only contain two types of correlation functions, q⁽⁰⁰⁾(t) and q^(10a)(t). As mentioned in the discussion of the Redfield-type closures, q⁽⁰⁰⁾(t) contains the same information as the Ehrenfest version of $C(t)$⁠, whereas q^(10a)(t) contains exact information about the system-bath interaction at t = 0. For this reason, we expect the c0b and c1b closures to yield significantly better dynamics than the cb2 and cb3 closures, which should simply recover the Ehrenfest dynamics.⁶⁷ Indeed, as Fig. 10 shows, the cb0 and cb1 closures are able to quantitatively correct the Ehrenfest dynamics. In contrast to the Redfield case, the cb0 and cb1 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, τ_c, which lies in a stability plateau alluded to in Sec. IV. The range of the stability plateau can depend sensitively on the regime of parameter space explored. Fig. 11 shows the dependence of the GQME+MFT dynamics for the Redfield- and NIBA-type projectors on the specific value of τ_c 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, τ_c. 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 τ_c 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 Fig. 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 cb0 and cb3 closures. The sensitivity of the cb3 closure to the time step is not difficult to understand, since for large dt the second numerical time derivative becomes noisy, leading to the growing oscillations in panel (d). Because closure cb0 shown in panel (a) contains the correlation functions q^(11s)(t) and q^(11a)(t), which require sampling of bath operators at t = 0 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 cb0 and cb1 closures, when fully converged with respect to bath realizations and time step, yield identical results that greatly improve the Ehrenfest results, while closures cb2 and cb3 are only capable of recovering the Ehrenfest dynamics.

In this paper we have developed a method to obtain 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,33,65 and Kelly, Markland, and coworkers^39–41 in the context of the SB model, and Rabani and coworkers^34–38 for more general models.

While previous applications of the GQME+semi-classical approach^39,40,65 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 (cb0, cb1, cf0, and cf1) and attribute this improvement, in part, to the sampling of static bath operators, ζ^W, 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 (cb2, cb3, cf2, and cf3). Our findings also provide numerical confirmation of the analytical proof included in Ref. 67, which indicates that use of the cb3 and cf3 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.

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 and Guy Cohen for useful conversations.

We begin by introducing the Fourier-Laplace transform of a time-dependent function f(t),

Its first and second time-derivatives take the following forms:

The kernel expansions in Eqs. (13) and (16) can be rewritten as follows:

In their original paper, Shi and Geva³¹ derived the following identity for the Redfield-type (thermal) projector:

where $Lsb$ is the Liouvillian corresponding to the system-bath interaction for the SB model $V=ασz∑kckxk$⁠. The second line is a simple extension of the derivation provided by Shi and Geva. These identities allow for the exact rewriting of the “$Q$-surrounded” projected propagator,

Replacing these expressions for the projected propagator in Eq. (10) 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 x ∈ {b, f}. Transforming Eqs. (A16) and (A17), we may solve for $K(2x)(ω)$

Substitution of Eqs. (A18) and (A19) into Eqs. (A14) and (A15) for the $Q$-forward and $Q$-backward closures yields

which are clearly equivalent to Eqs. (A4) and (A5), implying that the three-membered closures are equivalent to the two-membered closures.

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 s and a indicate symmetrized (anticommutator) or antisymmetrized (commutator) bath products.

Using the Redfield-type projector,

where A_i ∈ {|0〉〈0|, |1〉〈0|, |1〉〈0|, |1〉〈1|}, i ∈ {1, 2, 3, 4}, and ρ_B = e^−βH_B/Tr_B[e^−βH_B], the elements of the memory kernels take the following forms:

where X_n = 2(δ_n3 − δ_n2) and Y_n = 2(δ_n1 − δ_n4).

Using the NIBA-type projector,

where B_i ∈ {|0〉〈0|, |1〉〈1|}, i ∈ {1, 2} and ρ_B = e^−βH_B/Tr_B[e^−βH_B], the elements of the memory kernels take the following forms:

We employ a notation where some indices are squared since the q(t) functions are labelled using the indices corresponding to the A_j operators, which are related to the B_j operators in the following way: B₁ ↦ A₁ = A_1² and B₂ ↦ A₄ = A_2².

For the projector above to be truly of NIBA-type, $K(1)(t)$ should be $O(Δ2)$⁠. Instead, $K(1)(t)$ has contributions of first and second order in Δ. Indeed, the proper NIBA-type projector has the following form:

where $ρB(i)=e−β(HB−(−1)iVB)/TrB[e−β(HB−(−1)iVB)]$⁠.

In an open quantum system where a subsystem interacts weakly with a heat bath, the Ehrenfest method^12,13,16,59 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¹⁶ and density matrix formulations.⁵⁹ 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 SB model:

where i ∈ {x, y}, and the bath initial condition has unit trace, Tr_B[R_B] = 1. This corresponds to a nonequilibrium initial condition where the system is initially in a superposition of coherences. Immediately it is clear that ρ_S(0) = σ_i has zero norm. To remedy this, one may take advantage of the linearity of the problem and rewrite Eq. (C1) as sums of correlation functions with proper initial conditions,

where 1_S = |1〉〈1| + |2〉〈2|.

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. (C1), 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.

Fig. 13 shows the calculation of the nonequilibrium population dynamics given different system initial conditions, ρ_S(0) ∈ {σ_x, σ_y, |1〉〈1|}. For the initial conditions corresponding to the Pauli matrices, we implement three different decompositions given by Eqs. (C2)–(C4), labeled A, B, and C, respectively. We also include the results for two different decompositions for the initial condition |1〉〈1|, labeled A and B. 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 B underscores the problems associated with the lack of uniqueness in the density matrix based approach.

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

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

where the superscript W denotes the Wigner transform of the operator and Γ ≈ {Q_j,  P_j} is the set of all classical variables.

For the auxiliary kernels of the SB model, the following Wigner transforms are necessary: