Generalized quantum master equations (GQMEs) are an important tool in modeling chemical and physical processes. For a large number of problems, it has been shown that exact and approximate quantum dynamics methods can be made dramatically more efficient, and in the latter case more accurate, by proceeding via the GQME formalism. However, there are many situations where utilizing the GQME approach with an approximate method has been observed to return the same dynamics as using that method directly. Here, for systems both in and out of equilibrium, we provide a more detailed understanding of the conditions under which using an approximate method can yield benefits when combined with the GQME formalism. In particular, we demonstrate the necessary manipulations, which are satisfied by exact quantum dynamics, that are required to recast the memory kernel in a form that can be analytically shown to yield the same result as a direct application of the dynamics regardless of the approximation used. By considering the connections between these forms of the kernel, we derive the conditions that approximate methods must satisfy if they are to offer different results when used in conjunction with the GQME formalism. These analytical results thus provide new insights as to when proceeding via the GQME approach can be used to improve the accuracy of simulations.

Generalized quantum master equations (GQMEs) provide a formal framework to describe the time evolution of observables and correlation functions in complex many-body systems based on the projection operator method.1–3 The generalized master equation formalism has found extensive use both in allowing efficient and accurate calculations of material properties, including diffusion constants of liquids,4–8 density fluctuations in glasses,9–12 and structural relaxation in polymers.13,14 In addition, it has also been heavily exploited as an analysis tool to uncover the inherent time scales in complex chemical systems,15 as a dimensionality reduction technique in the development of coarse-grained molecular models,16 and more broadly in areas such as meteorological and financial time-series analysis and optimal prediction methods.17–20 The central quantity in the GQME formalism is the memory kernel, which encodes the effect of the projected dynamical degrees of freedom on the observable. However, the standard expressions for the memory kernel contain projected dynamical quantities that are impractical to simulate directly. This has led to the development of a number of ways to approximate the memory kernel which allows it to be recast in terms of unprojected dynamical quantities,21,22 assumed functional forms,23 or a given (perturbative or Markovian) limit.24–29 

Just over a decade ago, Shi and Geva derived a formally exact representation for the memory kernel of the Nakajima-Zwanzig GQME that requires only projection-free input.30 This representation opened the door to using either numerically exact or approximate methods to simulate the memory kernel. For exact treatments, whose computational cost increases severely with propagation time, exploiting the rapid decay of the memory kernel has been shown to allow for significant gains in the efficiency of simulating charge and energy transport in the condensed phase.30–40 

When approximate methods, such as those arising from the quantum-classical and semiclassical hierarchies,41–45 are used to calculate the memory kernel, significant improvements in accuracy have been observed when compared to their direct application.31,32,46–49 However, such improvements sensitively depend on how one calculates the projection-free partial kernels that are used to construct the memory kernel.49 These observations naturally raise questions as to why this is and when proceeding via the projection operator formalism will be advantageous.

Here we show, for both equilibrium and nonequilibrium systems, the conditions under which proceeding via the GQME formalism yields results that are guaranteed to be identical to the original dynamics used in the projection free input and suggest how this limitation can be overcome. To achieve this, we show how the memory kernel governing the evolution of equilibrium and nonequilibrium systems can be exactly recast in terms of unprojected correlation functions, which can be straightforwardly simulated using either exact or approximate methods. By analyzing these expressions, we derive the necessary requirements for an approximate dynamics to yield the same results when used directly and as an approximation to the memory kernel in the GQME approach. These results thus provide insights into when GQME methods might allow for improvement in accuracy or efficiency in equilibrium and nonequilibrium situations in a diverse set of systems.

To begin we consider the Nakajima-Zwanzig GQME,1,2 which provides a particularly straightforward route for describing the reduced dynamics of systems out of equilibrium. In this scheme the total system is decomposed into two parts: the system, which consists of all the degrees of freedom of interest in the problem, and the bath, which comprises the remaining degrees of freedom. For equilibrium systems the Mori approach is typically the preferred formulation.3 However, the Mori formalism, which renders distinction between system and bath unnecessary, is general and can be used for both equilibrium and nonequilibrium problems. Indeed, as observed in Ref. 49, both the Mori and Nakajima-Zwanzig approaches can be written in a unified formalism stemming from the projected equation of motion for the propagator,

(1)

where the Liouville operator is L=1ħ[Hˆ,], Hˆ is the Hamiltonian operator, P is the projection operator, and Q=1P is the complementary projection operator.

Central to the Mori approach is the appropriate choice of projection operator, P. In the following, we assume that it takes the form

(2)

where A contains a subset of observables which are of particular interest. The definition of the inner product is chosen such that (A|A) = 1, thus satisfying the idempotency condition, P2=P. Consequently, the complementary projector, Q, is by construction orthogonal to the subspace defined by P.

Using the Dyson operator identity

(3)

to expand the second term in the second equality of Eq. (1) yields

(4)

Restricting our attention to correlation functions defined by the inner product of the elements constituting P, the equation of motion for the propagator in Eq. (4) yields the following Mori-type GQME for the correlation function C(t)=(A|A(t)),

(5)

and the memory kernel, K(t), takes the form

(6)

Direct evaluation of the memory kernel in Eq. (6) is problematic, as it requires the action of the projected propagator, eiQLt.

To circumvent the difficulty of the projected propagator, the Dyson identity, Eq. (3), can be used to obtain a self-consistent expansion of the memory kernel

(7)

where the partial (auxiliary) kernels

(8)
(9)

no longer require the use of projected dynamics. The labels for the partial kernels of 1 and 3 are chosen so as to be consistent with earlier work.30,31,46–49 In principle, the memory kernel K(t) can be obtained by generating K1(t) and K3(t) from simulation and solving Eq. (7) numerically. Depending on the choice of projection operator and definition of the inner product, one can specialize this result to equilibrium correlation functions or nonequilibrium population dynamics.

In the equilibrium case, the Kubo-transformed correlation function is obtained by defining the inner product as

(10)

where O is a general superoperator in Liouville space (e.g., eiLt), ρeq=Z1eβHˆ is the canonical density operator, Z=Tr[eβHˆ] is the partition function, β = 1/kBT is the inverse of the thermal energy, and χAA = Tr[ρA(0)A(0)]. Commonly, the elements of the vector A are chosen to consist of a dynamical variable, a, which is a function of the coordinates and momenta of the system, and its time derivative, ȧ=iLa. For example, in the case of diffusion a could be chosen to be the position or velocity of some or all of the particles, and for infrared spectroscopy as the system dipole moment.

In nonequilibrium cases, the inner product may be defined as49 

(11)

where the set {An} spans a limited subspace of the total Hilbert space of the system, and RB is an operator that belongs to the complementary space, which is conventionally denoted as the bath. The normalization condition on RB requires that the trace over the bath degrees of freedom yields unity, TrB[RB]=1. The elements of A can, for example, be chosen such that {An} spans all outer products of the system states and RB=eβHB/TrB[eβHB] corresponds to the canonical distribution for the bath degrees of freedom. Such a choice for the projector provides access to nonequilibrium population and coherence dynamics of the system and has been used in the context of the spin-boson model,30,31,39,46,47,49 and in a wide class of quantum impurity models to access site occupation dynamics in the presence of one or more noninteracting fermionic or bosonic baths.34–38,40

The expressions in Eqs. (7)(9) have been shown to improve the accuracy of the dynamics produced by a number of approximate methods31,46–49 when used in the GQME formalism. This is shown in Fig. 1 for the spin boson problem, where the three different approximate semiclassical methods yield almost quantitatively exact results when used as an approximation to the memory kernel (solid lines) but fail markedly when used directly (dashed lines).

FIG. 1.

Evolution of the subsystem population difference versus time for the nonequilibrium relaxation of a spin-boson system initially prepared in the excited state. Numerically exact QUAPI55 (black dots), Ehrenfest mean field (blue), PBME53 (green), FBTS54 (red). Solid lines represent results using the GQME kernels in Eqs. (8) and (9), and dotted lines use those given in Eqs. (14) and (13). The dashed lines are direct dynamics. The simulation procedures and definition of the parameters are as used in Refs. 47 and 49.

FIG. 1.

Evolution of the subsystem population difference versus time for the nonequilibrium relaxation of a spin-boson system initially prepared in the excited state. Numerically exact QUAPI55 (black dots), Ehrenfest mean field (blue), PBME53 (green), FBTS54 (red). Solid lines represent results using the GQME kernels in Eqs. (8) and (9), and dotted lines use those given in Eqs. (14) and (13). The dashed lines are direct dynamics. The simulation procedures and definition of the parameters are as used in Refs. 47 and 49.

Close modal

Using manipulations that are exactly satisfied in quantum mechanics, but not necessarily by approximate methods, here we show how the memory kernel can be written in a form that returns a result that is identical to that obtained using the approximate method directly (irrespective of the approximate method employed). By analyzing the steps that are necessary to derive this expression, we are thus able to show the classes of approximate dynamics methods that will be guaranteed to always yield the same result when used directly or via the memory kernel formalism, no matter which expression for the memory kernel is used.

One begins by expanding the complementary projection operator Q and applying the Liouville operators.49 Carrying out the former operation for K3 (Eq. (9)) yields

(12)

The Liouville operator, L, can then be applied backwards on the static part or forwards to generate the time derivative. Doing the latter allows the partial kernel to be written purely as a function of C(t) and its time derivatives,

(13)

Proceeding similarly for K1 gives

(14)

where in the second and last equalities, the braces denote the anticommutator.

Due to the convolution-based structure of the equations relating the full memory kernel to the partial kernels, it is particularly convenient to consider its Fourier-Laplace representation. The Fourier-Laplace transform of C(t) is defined as

(15)

and its n-th time-derivative, C(n)(t)dndtnC(t), is

(16)

Application of the Fourier-Laplace transform to Eq. (7) gives

(17)

and to Eqs. (14) and (13) yields

(18)
(19)

where Ω(ω)=iω+Ċ(t=0). Using the Fourier-Laplace transforms of the memory kernels in Eqs. (18) and (19) in Eq. (17) gives

(20)

where we have introduced the superscript direct to highlight that Cdirect(ω) is the Fourier-Laplace transformed correlation function obtained from the direct dynamics. The corresponding GQME expression for the correlation function, CGQME(ω), can be obtained by Fourier-Laplace transforming and rearranging the equation of motion for the correlation function in Eq. (5),

(21)

where we have used C(t=0)=1, which follows from the idempotency property of the projector. Inserting the expression for the memory kernel in terms of the direct correlation function from Eq. (20) in this expression and rearranging gives CGQME(ω)=Cdirect(ω). Hence, when one uses expressions for the partial kernels that depend only on the original correlation function, C(t), and its time-derivatives, one is certain to recover exactly the same result as would be obtained from a direct application of the dynamics used to calculate the partial kernels.

This result therefore proves that no accuracy benefit can be obtained through the GQME for any approximate method using Eqs. (13) and (14).56 More explicitly, using manipulations that are exactly satisfied by quantum mechanics to recast the partial kernels in terms of the original correlation function and its time derivatives removes any potential benefit for gains in accuracy by proceeding via the GQME formalism. This is shown in Fig. 1 where, when Eqs. (13) and (14) are simulated using Ehrenfest mean field theory (blue dotted line), the same result is obtained as a direct application of mean field theory (blue dashed line).

It should be noted that all the expressions given for the partial kernels in the preceding equations (i.e., Eqs. (8) and (9) as well as Eqs. (14) and (13)) are guaranteed to give the same result as a direct application when exact quantum dynamics is used to generate them. However, if the partial kernels are shorter-lived than the correlation function, one can still obtain significant efficiency gains when using numerically exact methods, which scale poorly with simulation time.

The fact that Eqs. (13) and (14) cannot be used to obtain an increase in accuracy for any approximate method begs the question: What are the conditions that an approximate method must satisfy to guarantee that the same result will be obtained by using Eqs. (8) and (9) as Eqs. (14) and (13)? Identifying such conditions allows the identification of the classes of approximate methods which cannot gain accuracy benefits by combination with the GQME formalism. Indeed, we can immediately see that approximate methods only need to satisfy two such conditions. The first is that

(22)

where L indicates a Liouville operator corresponding to an approximate dynamics. Explicitly, this condition requires that the correlation function of A and the operator resulting from the action of the Liouvillian acting on A is equivalent to the time-derivative of the original correlation function. While the equality in Eq. (22) is trivially maintained for numerically exact methods, the same is not necessarily true for approximate methods. The second condition requires that the approximate Liouville operator commutes with the propagator,

(23)

These conditions are all that are required of an approximate method to perform the manipulations in Eq. (14) and likewise to obtain Eq. (13) from Eq. (9). Consequently, the GQME approach cannot yield improvements in accuracy when the approximate method used to calculate the partial kernels satisfies the conditions given by Eqs. (22) and (23).

In the case of equilibrium correlation functions, where the Liouville operator commutes with the canonical density matrix, Eq. (23) is equivalent to time-translational invariance. For example, dynamics methods such as centroid molecular dynamics (CMD),52 ring polymer molecular dynamics (RPMD),50,51 and purely classical mechanics satisfy these properties and thus are guaranteed not to obtain increased accuracy when used to approximate equilibrium correlation functions via the GQME approach. However, for nonequilibrium correlation functions, and for methods that do not satisfy these conditions, one can expect the dynamics resulting from the GQME to differ from the direct application of the approximate method, as has been shown in previous work where significant gains in both efficiency and accuracy have been achieved.31,32,46–49

In this paper, we have shown that when the memory kernel is formally cast in terms of the original correlation function and its time derivatives, the GQME is guaranteed to return the same result as that obtained from a direct simulation, regardless of the method used to construct the kernel. Further, we have used this to define the criteria that, if satisfied by an approximate dynamics method, guarantee that the method will produce GQME dynamics that are identical to that given by its direct dynamics, independent of the way the partial kernels are expressed. In these cases, where no such benefit in accuracy can be obtained, one may still expect the GQME scheme to yield improved efficiency if the memory kernels are short-lived. Conversely, violation of the criteria in Eqs. (22) and (23) implies that the resulting GQME dynamics will in general be different from direct simulation of the correlation functions. Previous studies31,46–49 and Fig. 1 have indeed confirmed that the GQME approach is capable of yielding significant gains both in efficiency and accuracy, when one does not invoke the formal recasting of the partial kernels in Eqs. (13) and (14), using a wide range of approximate theories. These insights thus outline a path for future applications of the GQME formalism by allowing one to assess the benefits that it may afford.

We greatly thank David Reichman for numerous insightful discussions. We also thank Will Pfalzgraff for helpful comments and a thorough reading of this manuscript, Mariana Rossi for helpful discussions and numerical verification of some of these analytic results, and Hsing-Ta Chen for useful conversations. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award No. DE-SC0014437. T.E.M. also acknowledges support from a Cottrell Scholarship from the Research Corporation for Science Advancement and an Alfred P. Sloan Research fellowship. A.K. and L.W. acknowledge postdoctoral fellowships from the Stanford Center for Molecular Analysis and Design.

1.
S.
Nakajima
,
Prog. Theor. Phys.
20
,
948
(
1958
).
2.
R.
Zwanzig
,
J. Chem. Phys.
33
,
1338
(
1960
).
3.
H.
Mori
,
Prog. Theor. Phys.
33
,
423
(
1965
).
4.
P. C.
Martin
and
S.
Yip
,
Phys. Rev.
170
,
151
(
1968
).
5.
G.
Harp
and
B.
Berne
,
Phys. Rev. A
2
,
975
(
1970
).
6.
D.
Levesque
and
L.
Verlet
,
Phys. Rev. A
2
,
2514
(
1970
).
7.
J.-P.
Boon
and
S. A.
Rice
,
J. Chem. Phys.
47
,
2480
(
1967
).
8.
T.-W.
Nee
and
R.
Zwanzig
,
J. Chem. Phys.
52
,
6353
(
1970
).
9.
E.
Rabani
and
D. R.
Reichman
,
J. Chem. Phys.
120
,
1458
(
2004
).
10.
E.
Rabani
and
D. R.
Reichman
,
Annu. Rev. Phys. Chem.
56
,
157
(
2005
).
11.
T. E.
Markland
,
J. A.
Morrone
,
B. J.
Berne
,
K.
Miyazaki
,
E.
Rabani
, and
D. R.
Reichman
,
Nat. Phys.
7
,
134
(
2011
).
12.
T. E.
Markland
,
J. A.
Morrone
,
K.
Miyazaki
,
B. J.
Berne
,
D. R.
Reichman
, and
E.
Rabani
,
J. Chem. Phys.
136
,
074511
(
2012
).
13.
B. U.
Felderhof
,
J. M.
Deutch
, and
U. M.
Titulaer
,
J. Chem. Phys.
63
,
740
(
1975
).
14.
K. S.
Schweizer
,
J. Chem. Phys.
91
,
5802
(
1989
).
15.
D.
Chandler
,
J. Chem. Phys.
68
,
2959
(
1978
).
16.
C.
Hijón
,
P.
Español
,
E.
Vanden-Eijnden
, and
R.
Delgado-Buscalioni
,
Faraday Discuss.
144
,
301
(
2010
).
17.
A. J.
Chorin
,
O. H.
Hald
, and
R.
Kupferman
,
Proc. Natl. Acad. Sci. U. S. A.
97
,
2968
(
2000
).
18.
D. T.
Schmitt
and
M.
Schulz
,
Phys. Rev. E
73
,
056204
(
2006
).
19.
M.
Niemann
,
T.
Laubrich
,
E.
Olbrich
, and
H.
Kantz
,
Phys. Rev. E
77
,
011117
(
2008
).
20.
D.
Venturi
and
G. E.
Karniadakis
,
Proc. R. Soc. A
470
,
20130754
(
2014
).
21.
B. J.
Berne
,
J. P.
Boon
, and
S. A.
Rice
,
J. Chem. Phys.
45
,
1086
(
1966
).
22.
B. J.
Berne
,
M. E.
Tuckerman
,
J. E.
Straub
, and
A. L. R.
Bug
,
J. Chem. Phys.
93
,
5084
(
1990
).
23.
X.
Chen
and
R. J.
Silbey
,
J. Chem. Phys.
132
,
204503
(
2010
).
25.
A. G. G.
Redfield
,
Adv. Magn. Reson.
1
,
1
(
1965
).
26.
27.
A. J.
Leggett
,
S.
Chakravarty
,
A.
Dorsey
,
M.
Fisher
,
A.
Garg
, and
W.
Zwerger
,
Rev. Mod. Phys.
59
,
1
(
1987
).
28.
U.
Weiss
,
Quantum Dissipative Systems
(
World Scientific
,
1992
).
29.
Y. C.
Cheng
and
R. J.
Silbey
,
J. Phys. Chem. B
109
,
21399
(
2005
).
30.
Q.
Shi
and
E.
Geva
,
J. Chem. Phys.
119
,
12063
(
2003
).
31.
Q.
Shi
and
E.
Geva
,
J. Chem. Phys.
120
,
10647
(
2004
).
32.
Q.
Shi
and
E.
Geva
,
J. Chem. Phys.
121
,
3393
(
2004
).
33.
M.-L.
Zhang
,
B. J.
Ka
, and
E.
Geva
,
J. Chem. Phys.
125
,
044106
(
2006
).
34.
G.
Cohen
and
E.
Rabani
,
Phys. Rev. B
84
,
075150
(
2011
).
35.
G.
Cohen
,
E. Y.
Wilner
, and
E.
Rabani
,
New J. Phys.
15
,
73018
(
2013
).
36.
G.
Cohen
,
E.
Gull
,
D. R.
Reichman
,
A. J.
Millis
, and
E.
Rabani
,
Phys. Rev. B
87
,
195108
(
2013
).
37.
E. Y.
Wilner
,
H.
Wang
,
G.
Cohen
,
M.
Thoss
, and
E.
Rabani
,
Phys. Rev. B
88
,
045137
(
2013
).
38.
E. Y.
Wilner
,
H.
Wang
,
M.
Thoss
, and
E.
Rabani
,
Phys. Rev. B
89
,
205129
(
2014
).
39.
E. Y.
Wilner
,
H.
Wang
,
M.
Thoss
, and
E.
Rabani
,
Phys. Rev. B
92
,
195143
(
2015
).
40.
L.
Kidon
,
E. Y.
Wilner
, and
E.
Rabani
,
J. Phys. Chem.
143
,
234110
(
2015
).
41.
42.
J. C.
Tully
,
J. Chem. Phys.
55
,
562
(
1971
).
43.
W. H.
Miller
,
J. Phys. Chem. A
105
,
2942
(
2001
).
44.
W. H.
Miller
,
J. Chem. Phys.
125
,
132305
(
2006
).
45.
R.
Kapral
,
J. Phys.: Condens. Matter
27
,
073201
(
2015
).
46.
A.
Kelly
and
T. E.
Markland
,
J. Chem. Phys.
139
,
014104
(
2013
).
47.
A.
Kelly
,
N.
Brackbill
, and
T. E.
Markland
,
J. Chem. Phys.
142
,
094110
(
2015
).
48.
W. C.
Pfalzgraff
,
A.
Kelly
, and
T. E.
Markland
,
J. Phys. Chem. Lett.
6
,
4743
(
2015
).
49.
A.
Montoya-Castillo
and
D. R.
Reichman
,
J. Chem. Phys.
144
,
184104
(
2016
).
50.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
121
,
3368
(
2004
).
51.
S.
Habershon
,
D. E.
Manolopoulos
,
T. E.
Markland
, and
T. F.
Miller
III
,
Annu. Rev. Phys. Chem.
64
,
387
(
2013
).
52.
S.
Jang
and
G. A.
Voth
,
J. Chem. Phys.
111
,
2371
(
1999
).
53.
H.
Kim
,
A.
Nassimi
, and
R.
Kapral
,
J. Chem. Phys.
129
,
084102
(
2008
).
54.
C.-Y.
Hsieh
and
R.
Kapral
,
J. Chem. Phys.
137
,
22A507
(
2012
).
55.
D. E.
Makarov
and
N.
Makri
,
Chem. Phys. Lett.
221
,
482
(
1994
).
56.

We note that a second decomposition of the memory kernel can be obtained from Eq. (6) through the use of the identity eiQLtQ=QeiLQt.49 A similar proof for this type of closure is straightforward and follows the same steps outlined above.