Approximate but Accurate Quantum Dynamics from the Mori Formalism: I. Nonequilibrium Dynamics
Abstract
We present a formalism that explicitly unifies the commonly used NakajimaZwanzig approach for reduced density matrix dynamics with the more versatile Mori theory in the context of nonequilibrium dynamics. Employing a Dysontype expansion to circumvent the difficulty of projected dynamics, we obtain a selfconsistent 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 spinboson model and limit our attention to the use of a simple and inexpensive quasiclassical 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 (Redfieldtype) and population based (NIBAtype) projection operators. We further elucidate the conditions that lead to shortlived 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 selfconsistently extracted memory kernel, we identify the mechanisms whereby the current approach leads to a significant improvement over the direct usage of standard semi and quasiclassical dynamics.
I Introduction
The continued effort to develop accurate and efficient approaches for the calculation of the dynamics of manybody 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,Makri (1992); Makri and Makarov (1995a, b); Meyer, Manthe, and Cederbaum (1990); Beck et al. (2000); Wang, Thoss, and Miller (2001); Thoss, Wang, and Miller (2001); Ishizaki and Fleming (2009) their computational cost makes them impractical for realistic multidimensional systems. Conversely, approximate methods, whether perturbativeBloch (1957); Redfield (1965); Leggett et al. (1987) or based on quasiGerber, Buch, and Ratner (1982); Stock (1995); Tully (1971, 1990, 1998) or semiclassicalMeyer and Miller (1979); Stock and Thoss (1997); Sun, Wang, and Miller (1998); Shi and Geva (2003a); Miller (2001a, b, 2009) 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 NakajimaZwanzig (NZ) equationNakajima (1958); Zwanzig (1960) 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,Grabert (1982); Fick and Sauermann (1990) 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.Redfield (1965); Bloch (1957); Leggett et al. (1987); Jang (2011); Tanimura, Suzuki, and Kubo (1989); Tanimura (2015)
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 GevaShi and Geva (2003b, 2004a); Zhang, Ka, and Geva (2006) proposed a selfconsistent 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 quasiadiabatic path integral (QUAPI) method,Makri and Makarov (1995a, b) Shi and Geva have shown that the memory kernel for the spinboson model can decay up to 10 times faster than the system dynamics,Shi and Geva (2003b) lending credence to the feasibility of the selfconsistent approach. More recently, a similar scheme has been used by Rabani and coworkers within a path integral framework for the study of quantum transport problems.Cohen and Rabani (2011); Cohen, Wilner, and Rabani (2013); Wilner et al. (2013, 2014a, 2014b) Just as importantly, applications of the method have successfully used semi and quasiclassical theories to calculate the auxiliary kernels.Shi and Geva (2004a) In particular, Kelly, Markland, and coworkers have illustrated the impressive accuracy and robustness of this approach in both model and realistic problems.Kelly and Markland (2013); Kelly, Brackbill, and Markland (2015); Pfalzgraff, Kelly, and Markland (2015) These studies have led to two important conclusions: (i) the memory kernels are shortlived in a wide region of parameter space for canonical problems such as the spinboson model and (ii) the selfconsistent 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 shortlived memory kernels remain unknown. Further, it is still unclear how the breakdown of approximate methods (in unfavorable parameter regimes) affects the quality of the selfconsistently 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. Shi and Geva, 2004a; Kelly and Markland, 2013; Kelly, Brackbill, and Markland, 2015 is lacking. Finally, the convergence properties of different versions of the auxiliary memory kernels arising from the alternative closures in the selfconsistent 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 singletime nonequilibrium dynamics are either cumbersome to obtain within this framework. Despite this difficulty, recent work has generalized the NZ equation to multitime correlation functions.Ivanov and Breuer (2015) In contrast, the more flexible Mori formulation permits direct extension to multipletime and equilibrium correlation functions, which are essential in the treatment of linearNitzan (2006) and nonlinear spectroscopy,Mukamel (1995) and the calculation of chemical rate constantsYamamoto (1960); Miller (2001a) and kinetic coefficients,Mahan (1990) to name a few examples. For this reason, in paper I (this paper) of this series, we provide a unified Moritype framework to approach singletime 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 systembath distinction exists, such as spin and fermion lattice models,Rasetti (1991); Sachdev (2011); Giamarchi (2004) and quantum fluids.Feenberg (1969); Pines and Nozieres (1966); Poulsen, Nyman, and Rossky (2005); Markland et al. (2011)
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. Sec. II introduces the spinboson model and the proposes two types of projection operators for this model. Sec. III 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. Kelly, Brackbill, and Markland, 2015, to obtain the auxiliary memory kernels. We henceforth refer to the use of the Ehrenfest method coupled to the selfconsistent extraction of the memory kernels as the GQME+MFT approach. In the Sec. IV, show the different properties of the memory kernels associated with the Redfield and NIBAtype projectors, explore the performance of the GQME+MFT approach to SB models whose systembath 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.
Ii Mori Approach
For illustrative purposes, we focus on the spinboson (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.Weiss (1992) 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 systembath separation, such as spinchains and lattice modelsRasetti (1991); Sachdev (2011); Giamarchi (2004) or liquids.Feenberg (1969); Pines and Nozieres (1966); Poulsen, Nyman, and Rossky (2005); Markland et al. (2011) 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 offdiagonal coupling , which is assumed to be independent of the bath coordinates,
(1) 
where corresponds to the Pauli matrix. The bath part of the Hamiltonian consists of independent harmonic oscillators,
(2) 
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,
(3) 
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 systembath interaction, and takes the form,
(4) 
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 densityLeggett et al. (1987) characterized by an exponential cutoff, and the Debye spectral density characterized by a Lorentzian cutoff:
(5)  
(6) 
where the cutoff frequency determines the correlation time of the bath at sufficiently high temperatures.Aslangul, Pottier, and SaintJames (1985) The reorganization energy, , is a measure of the strength of the systembath 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,
(7) 
where is the equilibrium density operator for the uncoupled bath, is the inverse thermal energy, and the initial condition corresponds to a FrankCondon transition.
ii.1 Generalized NakajimaZwanzigMori 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,
(8) 
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. (8) on an operator , yields the Mori EOM for that operator. Conversely, taking the Hermitian conjugate of Eq. (8), 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 selfconsistent expansion of the memory kernel has employed the thermal (ArgyresKelley) projection operator ,Shi and Geva (2003b, 2004a); Zhang, Ka, and Geva (2006); Cohen and Rabani (2011); Cohen, Wilner, and Rabani (2013); Kelly and Markland (2013); Kelly, Brackbill, and Markland (2015) 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}^{2}2For an introduction to this notation, see Chapter 2 in Ref. Mukamel, 1995 as , where contains all outer product states spanning the system. For the spinboson model, .
Using the thermal type projector, applying Eq. (8) 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,
(9) 
where is a static rotation matrix, corresponds to nonequilibrium averages of populations and coherences with all possible factorizable initial conditions, and is the socalled the inhomogeneous term. The elements of the memory kernel are given by,
(10) 
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 FrankCondon 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 Redfieldtype.Bloch (1957); Redfield (1965) The reason for this name is that truncation of the memory kernel at second order in is equivalent to a secondorder perturbative expansion of the memory kernel with respect to the systembath 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 secondorder truncation of the memory kernel with respect to leads to equations that are equivalent to the noninteracting blip approximation (NIBA), which is a second order expansion in the electronic coupling as opposed to the systembath coupling.Leggett et al. (1987); Sparpaglione and Mukamel (1988) Accordingly, we hereafter refer to this projector as NIBAtype. As in the case of the Redfieldtype projector, the choice for determines whether the inhomogeneous term is zero or finite.
Iii SelfConsistent Expansions for
As mentioned in Sec. I, the main difficulty associated with the memory kernel, Eq. (10), is the presence of the projected propagator, . To circumvent this problem, Shi and Geva proposed the use of the Dyson identity,
(11)  
(12) 
yielding a selfconsistent expansion of the memory kernel that only involves unprojected dynamics.Shi and Geva (2003b); Zhang, Ka, and Geva (2006) Not surprisingly, Eq. (8) can also be derived using the Dyson identity, Eqs. (11) and (12). Despite considering only timeindependent Hamiltonians in this work, extension of the formalism to timedependent Hamiltonians is simple and only requires the timeordered 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 of and type of auxiliary kernels employed.Zhang, Ka, and Geva (2006) 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 selfconsistent solution of the resulting integral equations.
iii.1 Bare expansions: Backward and Forward
The expansion of the memory kernel, Eq. (10), takes the following form,
(13) 
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
(14)  
(15) 
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 threemember expansions are equivalent to the twomember 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.Shi and Geva (2003b); Kelly and Markland (2013); Kelly, Brackbill, and Markland (2015)
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. (10) leads to the following expansion upon insertion of the Dyson identity,
(16) 
We note that has the form given in Eq. (14), and has the following form,
(17) 
It bears remarking that Eqs. (16) and (13) differ in the placement of under the integral and Eqs. (17) and (15) differ in whether or its Hermitian conjugate act on operators that require sampling only at or at finite times.
iii.2 Expansions using timederivatives
Because the auxiliary kernels given by Eqs. (14) and (15) or (17) require sampling of additional bath operators at and at finite times, convergence of these functions, at least within the context of semi and quasiclassical 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.Cohen and Rabani (2011); Cohen, Wilner, and Rabani (2013); Wilner et al. (2013, 2014b, 2014a); Kidon, Wilner, and Rabani (2015)
Here we focus on three types of auxiliary memory kernels that exploit different placements of the timederivative. 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,
(18)  
(19)  
(20) 
The second type focuses on replacing the Liouvillian operating on static operators with the time derivative, yielding the following expressions,
(21)  
(22)  
(23) 
The final type replaces all Liouvillians with time derivatives,
(24)  
(25)  
(26) 
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 , 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. (21) and (22) to solve for in Eq. 13.
Iv Results
The recent success achieved in using semi and quasiclassical schemes to calculate the auxiliary kernels required in the selfconsistent extraction of memory kernels underscores the importance of understanding the properties of this program in more detail. Consequently, we employ a simple quasiclassical method, namely Ehrenfest dynamics,Tully (1998); Grunwald, Kelly, and Kapral (2009) to obtain the auxiliary memory kernels and study the performance of the Redfield and NIBAtype 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.

Numerically integrate Eq. (9) subject to the appropriate initial conditions. In the following we use a secondorder RungeKutta algorithm.
The memory kernel decays to zero for a large sections of parameter space for the SB and other impuritytype models.^{3}^{3}3We are currently exploring cases for which the memory kernel decays to zero very slowly or decays to a finite constant. We refer to the timescale 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 quasiclassical 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.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 C and D.
iv.1 Projection Operators
In this section, we restrict our attention to the Ohmic spectral density, a model whose performance has already been studied extensively using the Redfieldtype projection operator and closure scheme by Kelly, Markland, and coworkers.Kelly, Brackbill, and Markland (2015) The purpose of this section is mainly to provide an analysis of the NIBAtype memory kernels and show the viability of the GQME+MFT approach using both the Redfield and NIBAtype projectors. Here and in Sec. IV.2, we show results only for the closure, and postpone the discussion of the closure dependence of the results to Sec. IV.3.
We first compare the different properties of the Redfield and NIBAtype memory kernels for a realization of the biased () spinboson model characterized by weak systembath coupling (), low temperature (), and a moderately fast bath (). Fig. 1 shows a representative set of components of the Redfieldtype memory kernel, , in panels (a)(d), and all components of the NIBAtype 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 Redfieldtype memory kernel has a short lifetime (), while the NIBAtype memory kernel decays much more slowly, .
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 longtime limit of the populations. Indeed, because of the assumption of a classical bath, the Ehrenfest method is known to violate detailed balance.Parandekar and Tully (2006) Instead, the dynamics resulting from both the Redfield and NIBAtype 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 NIBAtype memory kernels we recall that the Mori approach to Brownian motion,Mori (1965) 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.Reichman and Charbonneau (2005) The Redfieldtype projector spans the entire Hilbert space of the system, whereas the NIBAtype 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 timescale associated with the decay of the coherences induces the slow decay of the NIBAtype memory kernels. This conclusion further suggests that the NIBAtype projector may be most useful in instances of fast system relaxation, such as strong systembath 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 quasiclassical schemes like Ehrenfest meanfield theory,Kelly, Brackbill, and Markland (2015) the momentumjump solution to the quantumclassical Liouville equation,Kelly and Markland (2013) and the linearized semiclassical initial value representation (LSCIVR) scheme,Shi and Geva (2004b) relies primarily on the confluence of two important factors. First, the memory kernels are shortlived in comparison to the desired system dynamics. Second, approaches based on semiclassical arguments are more accurate at short times. Considered in tandem, these factors imply that the present scheme can lead to highly accurate shorttime memory kernels, thus avoiding problems associated with the longtime dynamics produced by these approximate methods by virtue of fast memory decay. However, the ability of the slowly decaying NIBAtype memory kernel to nearly recover the exact dynamics raises an important question: given the long lifetime of the NIBAtype kernel, how can it remain sufficiently accurate at long times to correct the longtime behavior of the bare quasiclassical dynamics? To answer this question, it is necessary to scrutinize the form of the auxiliary kernels. The NIBAtype auxiliary kernels for closure include two types of correlation functions, and given by Eqs. (54) and (57). 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 NIBAtype 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 systembath 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.3. The previous questions notwithstanding, the ability of the longlived NIBAtype memory kernel to produce dynamics that are comparably accurate to those obtained via the Redfieldtype 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 NIBAtype projector. Here we compare the NIBAtype GQME+MFT dynamics to exact results for the cases addressed by Kelly et al.Kelly, Brackbill, and Markland (2015) in their recent work characterizing the performance of the Redfieldtype GQME+MFT approach. For convenience, we include the results from the Redfieldtype 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 Fig. 3, direct use of the Ehrenfest MFT method consistently leads to incorrect longtime values of the population difference. In agreement with the work of Kelly and coworkers,Kelly, Brackbill, and Markland (2015) the Redfieldtype GQME+MFT method quantitatively corrects the dynamics in all considered cases. The NIBAtype approach generally provides clear improvement over direct use of MFT, but is slightly less accurate than the Redfieldtype GQME+MFT. Regardless, the improvement of the dynamics produced by the NIBAtype projector is remarkable not just because of the fact that the memory function is longlived, but also because such an approach is not tailored for the weak systembath coupling limit as is the Redfieldtype projector where the benchmark calculations of Fig. 3 have been performed.
iv.2 Debye Spectral Density
Due to its slower decay at large frequencies, the Debye spectral density is generally considered a more challenging case for trajectorybased dynamical methods.Berkelbach, Reichman, and Markland (2012) 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 Redfieldtype memory kernels for this spectral density generally means that a larger number of trajectories to achieve convergence are required.
Panels (a)(d) of Fig. 4 show the components for the Redfieldtype memory kernel, while panels (e)(f) show all components of the NIBAtype memory kernel. While the NIBAtype memory kernels do not show any obvious differences from their Ohmic counterparts, the Redfieldtype 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.
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 . Panel (a), which corresponds to an unbiased case characterized by weak systembath coupling at low temperature () and an intermediate bath frequency (), shows nearly perfect agreement between the Ehrenfest and exact dynamics. Both the Redfield and NIBAtype 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 longtime limits for all three biased cases. As in the Ohmic case, both the Redfield and NIBAtype approaches yield results in almost quantitative agreement with the exact dynamics. Slight deviations are evident, as in panel (b), where the NIBAtype GQME slightly underestimates the longtime limit of the population difference. Perhaps the most difficult case for the current approach, panel (d), shows that the NIBAtype GQME+MFT treatment leads to overly damped oscillations at long times, whereas the Redfieldtype 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 weakcoupling cases over a wide range of bath frequencies and temperatures.
iv.3 Memory Kernel Closures and Dynamics
In Sec. III, 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 NIBAtype 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. (54)–(59) of Appendix B. For convenience, we reproduce these expressions, within the Ehfenfest approximation, below,
(27)  
(28)  
(29)  
(30)  
(31)  
(32) 
where and .
Inspection of Eqs. (27)–(32) reveals that there are two main types of functions containing bath operators: those that require their sampling exclusively at [Eqs. (27), (29) and (30)], and those containing both statically and dynamically sampled bath operators [Eq. (28), (31), and (32)]. We recall the important fact that at , the Ehrenfest method is exact and the accuracy of the method diminishes with increasing simulation time.Golosov and Reichman (2001) 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 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. (28)–(32) 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. (30) and (32), 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 systembath 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 timederivatives lies in the accuracy of the timederivatives themselves. If the correlation functions calculated via the Ehrenfest procedure are sufficiently smooth and well converged and the timestep 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 NIBAtype 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 Redfieldtype kernel closures, we assess the effect of the backward and forward approaches by focusing first on the and closures. Inspection of Eqs. (15) and (17) reveals that the only difference between the forward and backward closures lies in the fact that contains the timeevolved 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 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 , 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 Fig. 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.Kelly et al. (2016) Further, the fact that the closure also recovers the bare Ehrenfest dynamics clearly indicates that the first of the two criteria specified in Ref. Kelly et al., 2016 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 timederivative 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 shorttime accuracy of Ehrenfest dynamics. Further, they illustrate that improvement within the memory formalism over the bare quasiclassical theory depends sensitively on the correlation functions calculated as input for the auxiliary memory kernels. In Sec. IV.1 we suggested that the Ehrenfest method might capture the dynamics of coherences more accurately than that of the populations, and that the selfconsistent 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. (29) and (31) contributes important information about the systembath 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 Redfieldtype 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 Redfieldtype 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 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 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” Wignertransformed 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. (78). Results for this approximation are shown in Fig. 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 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 Redfieldtype closures, the NIBAtype kernels only contain two types of correlation functions, and . As mentioned in the discussion of the Redfieldtype closures, contains the same information as Ehrenfest version of , whereas contains exact information about the systembath 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.Kelly et al. (2016) Indeed, as Fig. 10 shows, the and closures are able to quantitatively correct the Ehrenfest dynamics. In contrast to the Redfield case, the and NIBAtype 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 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 NIBAtype 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 shortlived 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 timestep 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 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.
V Conclusions
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 singletime nonequilibrium populations and coherences as well as more complicated dynamical objects, such as multitime, equilibrium and nonequilibrium correlation functions. We have shown that use of the Redfieldtype projector recovers the conventional NZ treatment previously used by Shi and GevaShi and Geva (2003b, 2004b); Zhang, Ka, and Geva (2006) and Kelly, Markland, and coworkersKelly and Markland (2013); Kelly, Brackbill, and Markland (2015); Pfalzgraff, Kelly, and Markland (2015) in the context of the SB model, and Rabani and coworkersCohen and Rabani (2011); Cohen, Wilner, and Rabani (2013); Wilner et al. (2013, 2014b, 2014a) for more general models.
While previous applications of the GQME+semiclassical approachShi and Geva (2004b); Kelly and Markland (2013); Kelly, Brackbill, and Markland (2015) have been limited to the Redfieldtype projector and have focused on the improvement over the bare semiclassical 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 NIBAtype projector, do not result in an inaccurate description of the GQME dynamics. This demonstrates that the success of the GQME+semiclassical approach is not a function of the shorttime 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 semiclassical 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,Kelly et al. (2016) 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 Moribased formulation furthered in this work provides a flexible framework to accurately study problems that go beyond the scope of nonequilibrium dynamics for SBtype models. For instance, the Mori formalism can easily address equilibrium and multitime correlation functions in systems coupled to harmonic and anharmonic baths, as well as problems where the systembath distinction is absent. Work in this latter direction will be pursued in future papers.
Vi Acknowledgements
The authors would like to thank Aaron Kelly and Tom Markland for helpful discussions. D.R.R. acknowledges support from NSF No. CHE1464802. A.M.C. thanks HsingTa Chen and Guy Cohen for useful conversations.
Appendix A FourierLaplace Analysis of Closures
We begin by introducing the FourierLaplace transform of a timedependent function ,
(33) 
Its first and second timederivatives take the following form,
(34)  
(35) 
In their original paper, Shi and GevaShi and Geva (2003b) derived the following identity for the Redfieldtype (thermal) projector,
(38)  
(39) 
where is the Liouvillian corresponding to the systembath 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,
(40)  
(41) 
Replacing these expressions for the projected propagator in Eq. (10) followed by use of the Dyson decomposition leads to the following three rather than twomembered expansions,
(42)  
(43) 
where the second auxiliary kernel also contains the projected propagator,
(44)  
(45) 
In FourierLaplace space, the previous expressions for the memory kernel become,
(46)  
(47) 
In the timedomain, the second auxiliary kernels may be expanded as follows,
(48)  
(49) 
where . Transforming Eqs. (48) and (49), we may solve for
(50)  
(51) 
Substitution of Eqs. (50) and (51) into Eqs. (46) and (47) for the forward and backward closures yields,
(52)  
(53) 
which are clearly equivalent to Eqs. (36) and (37), implying that the threemembered closures are equivalent to the twomembered closures.
Appendix B Expressions for Auxiliary Memory Kernels
Here we provide explicit expressions for the components of the memory kernel using the Redfield and NIBAtype projection operators. Before going further, however, we introduce for notational clarity the following correlation functions.
(54)  
(55)  
(56)  
(57)  
(58)  
(59) 
where and indicate symmetrized (anticommutator) or antisymmetrized (commutator) bath products.
Using the Redfieldtype projector,
(60) 
where , , and , the elements of the memory kernels take the following forms,
(61)  
(62)  
(63) 
where and .
Using the NIBAtype projector,
(64) 
where , and , the elements of the memory kernels take the following forms,
(65)  
(66)  
(67) 
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 NIBAtype, should be . Instead, has contributions of first and second order in . Indeed, the proper NIBAtype projector has the following form,
(68) 
where .
Appendix C Initial Conditions in the Ehrenfest method
In an open quantum system where a subsystem interacts weakly with a heat bath, the Ehrenfest methodGerber, Buch, and Ratner (1982); Stock (1995); Tully (1998); Grunwald, Kelly, and Kapral (2009) 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 wavefunctionTully (1998) and density matrix formulations.Grunwald, Kelly, and Kapral (2009) Because of its clarity, the derivation based on the density matrix and the quantumclassical 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 spinboson model
(69) 
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. (69) into sums of correlation functions with proper initial conditions,
(70)  
(71)  
(72) 
where .
Returning to the wavefunctionbased derivation of the Ehrenfest approach requires that any subsystem initial condition correspond to a pure state. Referring again to Eq. (69), we may rewrite the trace over the subsystem in the eigenbasis of the subsystem’s initial condition,
(73) 
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,