The article information
 Yanying Liu, Yaming Yan, Meng Xu, Kai Song, Qiang Shi
 柳燕莺, 闫亚明, 许猛, 宋凯, 史强
 Exact Generator and its High Order Expansions in TimeConvolutionless Generalized Master Equation: Applications to SpinBoson Model and Exictation Energy Transfer^{†}
 时间无卷积广义主方程中核函数的精确解和高阶微扰展开:在自旋玻色模型和激发态能量转移中的应用
 Chinese Journal of Chemical Physics, 2018, 31(4): 575583
 化学物理学报, 2018, 31(4): 575583
 http://dx.doi.org/10.1063/16740068/31/cjcp1806146

Article history
 Received on: June 20, 2018
 Accepted on: July 27, 2018
Quantum dynamics plays a significant role in many chemical and physical processes of condensed phases, while the development of accurate and efficient schemes to simulate such processes remains an important challenge in theoretical chemistry [14]. Numerically exact approaches [511] are usually limited to some simple models, while approximate approaches [1216] generally have limited range of applicability. Therefore, searching for methods that are accurate, efficient, and general still stays in the forefront of quantum dynamics in condensed phases. A popular idea to overcome this issue is to split the total system into two parts: the quantum system which we are the most interested in, and the environment or bath that exchanges energy and particles with the system. The quantum system is described specifically using the reduced density operator (RDO) and its interaction with the bath is included in the equation of motion of the RDO [1720].
For calculations of the dynamics within the reduced system dynamics framework, the NakajimaZwanzig (NZ) (or timeconvolution, TC) master equation, an integrodifferential equation with a superoperator called "memory kernel" derived from the full Hilbert space using the projection operator techniques, is often employed [2124]. With the NZ generalized master equation, the complexity of treating the bath effects is reduced to the calculation of the memory kernel which completely determines the dynamics of the system [18]. Although formally exact memory kernels can be written using the projection operators, they are hard to compute explicitly except for several simple models. For example, no analytical exact kernels are available for the commonly used spinboson [2, 25] and Anderson impurity [26, 27] models. Therefore a large number of studies using NZ master equation have been based on perturbation expansions with respect to the strength of the systembath coupling or a small coupling parameter in the system Hamiltonian [2830]. It is noted that, numerical methods to calculate exact NZ memory kernel [3133], are proposed in many recent studies. In a recent work, we have also investigated systematically the convergence of high order expansions of the NZ memory kernel [34].
Instead of using an integrodifferential equation to describe the bath effects, an alternative form of the generalized master equation is the HashitsumeShibataTakahashi or timeconvolutionless (TCL) approach [3537]. Also being formally exact, this approach describes the evolution of the RDO using a first order differential equation with time dependent coefficients [18, 3638]. The key quality in the TCL generalized master equation is a superoperator referred to as the TCL kernel or generator, which eliminates the integration over the system history and is local in time [18]. However, the exact calculation of the TCL generator is intractable for it relies on performing the inversion of a superoperator in the full Hilbert space [18]. Closed expressions for TCL master equations have been obtained for several analytically solvable models such as a harmonic oscillator coupled bilinearly to a harmonic bath [2, 19, 39, 40], the JaynesCummings model [41, 42], and the resonant energy level model in charge transport [43]. Although there are recent attempts to obtain the nonperturbative TCL generators numerically [44, 45], many applications of the TCL generalized master equations rely on perturbative schemes, especially the second order perturbation approximations due to their simplicity. A major problem of the second order TCL master equations is that their validity is not guaranteed in the case of strong coupling. The analytical fourth order perturbation approaches [18, 46], however, are usually rather cumbersome and their extensions to higher orders are difficult. An attempt to formulate a recursive approach to calculating high orders of the TCL generator has been proposed recently [47], but has not been applied to realistic calculations of reduced quantum dynamics.
It should be noted that, when solved exactly, the TCL and NZ master equations give the same results [48]. While comparing with the NZ form of the generalized quantum master equation (GQME), the TCL form of the GQME usually has wider applications in several problems such as in calculating spectroscopic signals [49, 50]. Especially, the second order TCL GQME becomes exact for the pure dephasing model where a two level system with zero interstate coupling is coupled linearly to a harmonic bath [51, 52], while the second order NZ GQME does not. On the computational cost to obtain the memory kernels and generators in the two forms of GQME, it is found that both kernels/generators usually decay on similar time scale and thus, on this point, neither formulations is more advantageous than the other.
In this work, we adopt a population propagator approach based on recent studies [42, 44, 45], and the hierarchical equation of motion (HEOM) method [10, 53, 54] to calculate the exact TCL generator. Then, similar to the case of the high order expansions of the NZ memory kernel [34], a new approach to calculating the high order expansions of the TCL generator is proposed based on the previously developed extended HEOM method to calculate perturbation expansion of open system dynamics [55]. The new approach utilizes the results from the extended HEOM, but does not in volve the cumbersome multidimensional integrals. It is then applied to the spinboson model to investigate the convergence of the high order TCL generators. We also show that, in several examples, the exact TCL generators may become singular, which may limit the application of the TCL generalized master equations in certain circumstances.
The remainder of this paper is organized as follows. In Section Ⅱ.A and Section Ⅱ.B, we briefly review the TCL generalized master equation, the spinboson model, and the HEOM method. Methods to derive the exact TCL generator and its high order expansions are given in Section Ⅱ.C. In Section Ⅲ, we apply the proposed method to the spinboson model and present the numerical results for the exact TCL generators and their high order expansions. Examples demonstrating the singularity of the TCL generator are also presented. Conclusions and discussions are made in Section Ⅳ.
Ⅱ. THEORY A. The timeconvolutionless master equationWe consider a general total Hamiltonian describing a system coupled to a bath
$ \begin{array}{l} H = H_\textrm{S} + H_\textrm{B} + H_{\textrm{SB}} \;\; \label{eqhamiltonian} \end{array} $  (1) 
where
$ \begin{array}{l} \frac{\textrm{d}}{\textrm{d}t}\rho_\textrm{T}(t) = i[H, \rho_\textrm{T}(t)] \\ = i{\cal L} \rho_\textrm{T}(t) \end{array} $  (2) 
where
To derive the generalized master equation, we separate the density operator into relevant and irrelevant parts by means of a projection operator
$ \begin{array}{l} \label{eqtclprojector1} \frac{\textrm{d}}{\textrm{d}t}{\mathcal P} \rho_\textrm{T}(t) = i\mathcal{PL}\theta(t)\mathcal{P}\rho_\textrm{T}(t) +{\cal I}(t){\mathcal Q} \rho_\textrm{T}(t_0)\;\; \end{array} $  (3) 
where
$ \begin{array}{l} {\theta}(t) = {\left[1+i \int_{t_0}^t \textrm{d}s \textrm{e}^{i {\mathcal Q} {\cal L} (ts)}{\mathcal Q} {\cal L} {\mathcal P}\textrm{e}^{i {\cal L} (ts)}\right]}^{1} \;\; \end{array} $  (4) 
$ \begin{array}{l} {\cal I}(t) = i{\mathcal P} {\cal L} \theta (t) \textrm{e}^{i {\mathcal Q} {\cal L} (tt_0)} {\mathcal Q} \;\; \end{array} $  (5) 
and
$ \begin{array}{l} \label{eqtclprojector2} \frac{\textrm{d}}{\textrm{d}t}{\mathcal P}\rho_\textrm{T}(t) = i\mathcal{PL}\theta(t)\mathcal{P}\rho_\textrm{T}(t)\;\; \end{array} $  (6) 
In the following parts of this work, we adopt the Liouville space notation [5658] for the total density operator,
We employ the following projection operator which has been used in many previous studies [56, 5860]
$ \begin{array}{l} \label{projector} {\mathcal P} = \sum \limits_j  j \rho_j^\textrm{B} \rangle \rangle \langle \langle j  \;\; \end{array} $  (7) 
where
$ \begin{array}{l} \label{eqtclpopulation1} \frac{\textrm{d}}{\textrm{d}t}P (t) = {\cal R}(t)P (t)\;\; \end{array} $  (8) 
where
$ \begin{array}{l} {\mathcal P}{\rho}_\textrm{T} (t) = {\mathcal P}\textrm{e}^{i{\cal L} t} {\rho}_\textrm{T}(0) \\ = {\mathcal P}\textrm{e}^{i{\cal L} t} ({\mathcal P}+{\mathcal Q}) {\rho}_\textrm{T}(0)\ \label{eqp} \end{array} $  (9) 
With the above projector operator
$ \begin{array}{l} P(t) = {\cal U}_\textrm{S}(t)P(0) \label{eqevolution1} \end{array} $  (10) 
where the matrix element of
$ \begin{array}{l} {\cal U}_{\textrm{S};jk}(t) = \langle\langle j \textrm{e}^{i{\cal L} t} k\rho_k^{\textrm{B}}\rangle\rangle \label{equsdefinition} \end{array} $  (11) 
involves the propagator of the system only. Taking the time derivative operation of
$ \begin{array}{l} \dot{P}(t) = \dot{{\cal U}}_\textrm{S}(t)P(0) \\ = \dot{{\cal U}}_\textrm{S}(t) {\cal U}_\textrm{S}^{1}(t)P(t)\ \label{eqtclpopulation2} \end{array} $  (12) 
provided the inverse does exist. By comparing with the TCL quantum master equation in Eq.(8), the TCL generator
$ \begin{array}{l} \label{eqrt} {\cal R}(t) = \dot{{\cal U}}_\textrm{S}(t){\cal U}_\textrm{S}^{1}(t)\ \end{array} $  (13) 
Thus, in order to obtain the TCL generator, one has to compute the time derivation of
$ \begin{array}{l} \dot{{\cal U}}_{\textrm{S};jk}(t) = i\langle\langle j \textrm{e}^{i{\cal L} t}\mathcal{L} k\rho_k^{\textrm{B}}\rangle\rangle \label{eqdus} \end{array} $  (14) 
In the above derivations, the condition
For simplicity, we restrict ourselves to the spinboson model for most part of this work. The spinboson Hamiltonian, in the form of Eq.(1), describes a twostate system. The Hamiltonian of the system
$ \begin{array}{l} H_\textrm{S} = {\epsilon} {\sigma}_z + {\Delta} {\sigma}_x \ \end{array} $  (15) 
is characterized by the energy bias
$ \begin{array}{l} H_\textrm{B} = \sum\limits_j \left(\frac{p_j^2}{2} + \frac{1}{2} {\omega}_j^2 x_j^2 \right) \end{array} $  (16) 
$ \begin{array}{l} H_{\textrm{SB}} =  \sum\limits_{j} c_jx_j \otimes {\sigma}_z\equiv F \otimes {\sigma}_z\ \end{array} $  (17) 
Here,
The systembath interaction is described by the spectral density, which is defined as: [1, 61]
$ \begin{array}{l} J(\omega) = \frac{\pi}{2} \sum\limits_j \frac{c_j^2}{{\omega}_j} \delta (\omega {\omega_j}) \end{array} $  (18) 
The correlation function
$ \begin{array}{l} C(t>0) = \frac{1}{Z_\textrm{B}} {\rm Tr}[\textrm{e}^{{\beta}H_\textrm{B}}\textrm{e}^{iH_\textrm{B}t}F\textrm{e}^{iH_\textrm{B}t}F] \\ = \frac{1}{\pi} \int_{\infty}^{\infty} \textrm{d}{\omega}J ({\omega})\frac{\textrm{e}^{i{\omega}t}}{1\textrm{e}^{{\beta}{\omega}}}\ \label{eqcorelation} \end{array} $  (19) 
where
$ \begin{array}{l} J(\omega) = \frac{\eta {\omega}_\textrm{c} \omega}{{\omega}^2+{{\omega}_\textrm{c}}^2}\;\; \end{array} $  (20) 
where
$ \begin{array}{l} C(t>0) = \sum\limits_k d_k\textrm{e}^{{{\omega}}_kt} \label{eqcorelation1} \end{array} $  (21) 
where
$ \begin{array}{l} d_0 = \frac{1}{2}{\eta}{{\omega}_\textrm{c}}\left[\cot\left(\frac{\beta{\omega}_\textrm{c}}{2}\right)i\right] \end{array} $  (22) 
$ \begin{array}{l} d_k = \frac {4k\pi\eta{\omega}_\textrm{c}}{{(2k\pi)}^2 {(\beta{\omega}_\textrm{c})}^2}\ , \, {\rm for}\, \, k\geq1 \end{array} $  (23) 
Then the HEOM can be derived using the path integral technique, [2, 10, 53, 6367] or the stochastic Liouville equation approach [6870]:
$ \begin{array}{l} \frac{\partial}{{\partial}t} {\rho}_{\textbf{n}}(t) =  \left(i{\cal L} + \sum\limits_k n_k {{\omega}}_k\right) {\rho}_{\textbf{n}}(t)  \\ i\left[ {\sigma}_z, \sum\limits_k {\rho}_{{\textbf{n}}_k^+}(t)\right]  \\ i \sum\limits_kn_k \left(d_k {\sigma}_z {\rho}_{{\textbf{n}}_ k^}(t) d_k^* {\rho}_{{\textbf{n}}_k^}(t) {\sigma}_z \right)\ \label{eqheom1} \end{array} $  (24) 
The subscript
As shown in Eq.(13), to obtain the exact TCL generator, we first need to calculate the value of
$ \begin{array}{l} {[{\cal U}_\textrm{S}(t)]}_{jk} = \langle \langle j  \textrm{e}^{i {\cal L}t}  k \rho_k^\textrm{B} \rangle \rangle \\ = {\rm Tr}_\textrm{S}{\rm Tr}_\textrm{B} \{ j \rangle \langle j  \textrm{e}^{i{\cal L}t}  k \rangle \langle k \otimes\rho_k^\textrm{B} \} \\ = [\sigma_k(t)]_{jj} \;\; \label{equs2} \end{array} $  (25) 
where
$ \begin{array}{l}\label{Eq:sbprojction} \mathcal{P} = 1\rho_1^{\textrm{B}}\rangle\rangle\langle\langle 1 + 2\rho_2^{\textrm{B}}\rangle\rangle\langle\langle 2 \;\; \end{array} $  (26) 
where
$ \begin{array}{l} {\rho}_j^\textrm{B} = \frac{\textrm{e}^{{\beta}H^{(j)}}}{{\rm Tr}_\textrm{B}\{\textrm{e}^{{\beta}H^{(j)}}\}} \end{array} $  (27) 
with
To calculate the derivative of
$ \begin{array}{l} H_\textrm{S} = H_0 + H_1 \end{array} $  (28) 
where
$ \begin{array}{l} {[\dot{\cal U}_\textrm{S}(t)]}_{jk} = {(1)}^{j+k+1}{\Delta}{[{\sigma}_z{\sigma}_ {y, k}(t)]}_{jj} \;\; \label{equs3} \end{array} $  (29) 
where
In the rest part of this section, we will deduce the high order expansions of
$ \begin{array}{l} {\cal R}(t) = \sum\limits_{N=2}^{\infty}{\Delta}^N{\cal R}^{(N)}(t) \end{array} $  (30) 
$ \begin{array}{l} \label{equs5} {\cal U}_\textrm{S}(t) = {\mathbf I}+\sum\limits_{N=2}^{\infty}{\Delta}^N {\cal U}_\textrm{S}^{(N)}(t) \;\; \end{array} $  (31) 
$ {{\mathcal{ \dot U}}_{\rm{S}}}(t) = \sum\limits_{N = 2}^\infty {{\Delta ^N}} {\mathcal{\dot U}}_{\rm{S}}^{(N)}(t) $  (32) 
By exploiting the identity
$ \begin{array}{l} {\cal U}_\textrm{S}^{1}(t) = \sum\limits_{n=0}^{\infty}{(1)}^n {\left[\sum\limits_{N=2}^{\infty}{\Delta}^N{\cal U}_\textrm{S}^{(N)}(t) \right]}^n\ \end{array} $  (33) 
Thus we can obtain from Eq.(13),
$ \begin{array}{l} \sum\limits_{N=2}^{\infty}{\Delta}^N{\cal R}^{(N)}(t) = \left\{\sum\limits_{N=2}^{\infty}{\Delta}^N\dot{\cal U}_ \textrm{S}^{(N)}(t)\right\} \Bigg\{\sum\limits_{n=0}^{\infty}{(1)}^n\cdot \\ \left[\sum\limits_ {N=2}^{\infty}{\Delta}^N{\cal U}_\textrm{S}^{(N)}(t) \right]^n\Bigg\} \end{array} $  (34) 
By expanding this equation and comparing coefficients for each order of
$\begin{array}{l} {\cal R}^{(2)}(t) = \dot{\cal U}_\textrm{S}^{(2)}(t) \;\; \end{array} $  (35) 
$ \begin{array}{l} {\cal R}^{(4)}(t) = \dot{\cal U}_\textrm{S}^{(4)}(t)  \dot{\cal U}_\textrm{S}^{(2)}(t) {\cal U}_\textrm{S}^{(2)}(t) \;\; \end{array} $  (36) 
$ \begin{array}{l} {\cal R}^{(6)}(t) = \dot{\cal U}_\textrm{S}^{(6)}(t)  \dot{\cal U}_\textrm{S}^{(4)}(t) {\cal U}_\textrm{S}^{(2)}(t) + \\ \dot{\cal U}_\textrm{S}^{(2)}(t) {\cal U}_\textrm{S}^{(2)}(t) {\cal U}_\textrm{S}^{(2)}(t)  \dot{\cal U}_\textrm{S}^{(2)}(t) \ {\cal U}_\textrm{S}^{(4)}(t) \;\; \quad\\ \cdots \end{array} $  (37) 
The odd order terms in the expansion of the exact generator
$ \begin{array}{l} \label{eqrecursion} {\cal R}^{(2n)}(t) = \dot{\cal U}_\textrm{S}^{(2n)}(t) \sum\limits_{m=1}^{n1} [{\cal R}^{(2m)}(t) {\cal U}_\textrm{S}^{[2(nm)]}(t)]\ \end{array} $  (38) 
In order to obtain the expansions in Eq.(31) and Eq.(32), the extended HEOM method [xu17] is employed. More specifically, the perturbation expansion of
$ \begin{array}{l} {\sigma}_{y, k}(t) = \sum\limits_{N=0}^{\infty} \frac{1}{N!} {\sigma}_{y, k}^{(N)}(t){\Delta}^N \;\; \label{eqsigma} \end{array} $  (39) 
The
$ \begin{array}{l} \frac{\partial}{{\partial}t}{\rho}_{\textbf{n}}^{(N)}(t) = i{\cal L}_0{\rho}_{\textbf{n}}^{(N)}(t)iN[{\sigma}_x, {\rho}_ {\textbf{n}}^{(N1)}(t)] \\ \sum\limits_k(n_k{{\omega}}_k) {\rho}_{\textbf{n}}^{(N)}(t) i\left[{\sigma}_z, \sum\limits_k{\rho}_{{\textbf{n}}_k^+}^{(N)}(t)\right]  \\ i\sum\limits_kn_k\left[c_k{\sigma}_z{\rho}_{{\textbf{n}}_k^}^ {(N)}(t) d_k^*{\rho}_{{\textbf{n}}_k^}^{(N)}(t){\sigma}_z\right] \;\; \label{eqheom2} \end{array} $  (40) 
where
For all numerical calculations in this section, the initial state of the total system is assumed to be equilibrated on state
There are four numerical examples in this section. The first three are based on the spinboson model with different parameters, including the systembath coupling strength
In this subsection, we investigate the convergence of the high order perturbative expansions of the TCL generator presented in the above Section Ⅱ.C. In the first example, FIG. 1 shows the time evolution of population calculated by the HEOM method as well as the TCL generalized master equation with the exact generator. The parameters in this example are
The high order expansions up to the 12th order of the TCL generator
In the second example, the following parameters are used:
The high order expansions of the TCL generator for the second example are shown in FIG. 4(a). In contrast to the case shown in FIG. 2(a), the amplitude of the high order terms does not decrease as the perturbation order increases. Thus it is expected that the perturbation expansion at finite order will soon breakdown when the propagation time increases so that the generator expansion becomes hard to converge. This is exactly the case as shown in FIG. 4(b).
By comparing the high order expansions shown in FIG. 2 with that in FIG. 4, a good indication of the possible convergence is to compare the relative amplitudes of the second order perturbation term
To explore the convergence of the high order TCL generators in different parameter regimes, we also calculate the critical values of interstate coupling
$ \begin{array}{l} S_{11}^{(n)}(t/\pi=2.5)  S_{11}^{(10)}(t/\pi=2.5)\leq0.001 \;\; \end{array} $  (41) 
for the perturbation order
In contrast to the memory kernel in the NZ generalized master equation, the formally exact generator of the TCL generalized master equation in Eq.(3) involves the inverse of a Liouville space superoperator. And as discussed in Section Ⅱ, we assume that the inverse of
In this subsection, we show such an example by setting
This result can be understood using the following theoretical analysis. For the symmetric spinboson model, the
$ \begin{array}{l} \begin{array}{l} {\cal U}_{\rm{S}}(t) = \left( {\begin{array}{*{20}{c}} {a(t)}&{1  a(t)}\\ {1  a(t)}&{a(t)} \end{array}} \right)\end{array}\end{array} $  (42) 
$ \begin{array}{l} \begin{array}{l} {\cal R}(t) = \left( {\begin{array}{*{20}{c}} {b(t)}&{  b(t)}\\ {  b(t)}&{b(t)} \end{array}} \right) \end{array} \label{eqsingularity} \end{array} $  (43) 
When
Because of the existence of the singularities, it is reasonable to predict the divergence of the high order expansion of the TCL generator in this example. The results shown in FIG. 7 (a) and (b) verify this prediction. It can be seen from FIG. 7(a) that amplitudes of high order terms increase rapidly with the perturbation order. And the results in FIG. 7(b) show that the summations of perturbative terms diverge from the exact one after only a very short time.
It is interesting to investigate whether similar problem of singular generators also exists in more general model systems beyond the simple symmetric spinboson model. We apply the method mentioned in Section Ⅱ.C to the problem of excitation energy transfer in the FMO complex. The model Hamiltonian for this problem can be found in previous work [74], and the parameters used in this simulation are obtained from Refs.[74, 75]. FIG. 8(a) shows the time evolution of population when the initial state is assumed to be equilibrated on site 1 of the FMO complex (It is noted that most simulations in the literature use an unrelaxed initial state for the bath degrees of freedom). And FIG. 8(b) shows the corresponding exact generator
It should also be noted here that, the above behaviors of singular TCL generator are investigated using a generalized master equation of populations, using the projection operator defined in Eq.(7). The condition of singular generator will also be different if other types of projection operators are employed.
Ⅳ. CONCLUSIONS AND DISCUSSIONSIn this work, we have combined the population propagator formalism and the HEOM method to obtain the exact TCL generator, and proposed a new approach to calculate the high order expansions of the generator with the recently developed extended HEOM method to obtain high order expansions of the open system dynamics. A recursive relation is derived with which one can obtain the high order expansions of the TCL generator without the need to calculate the cumbersome and timeconsuming multidimensional integrals.
By using the spinboson model as an example, we have investigated the convergence of the high order expansion of the TCL generator. It is shown that the convergence of the high order generators depends on the system parameter, and the ratio of the fourth to second order terms can serve as a good indicator of the convergence. A potential problem of the exact TCL generator is that, it may become singular in certain parameter regimes, as has been demonstrated in the analytically solvable JaynesCummings model [42] and the resonant energy model [43] for charge transport. We have also shown that the problem of singular exact generator also exists in the case of spinboson model, and a commonly used model for excitation energy transfer in the FMO complex. Of course, the condition for singular generator also depends on the choice of projection operator, as well as the parameters. The current study is based on the TCL generalized master equation using system populations, and this topic may be further explored for other different forms of TCL GQME.
Ⅴ. ACKNOWLEDGEMENTSThis work is supported by the National Natural Science Foundation of China (No.21673246), and the Strategic Priority Research Program of the Chinese Academy of Sciences (No.XDB12020300).
[1]  A. Nitzan, Chemical Dynamics in Condensed Phases:Relaxation, Transfer and Reactions in Condensed Colecular Systems. New York: Oxford University Press (2006). 
[2]  U. Weiss, Quantum Dissipative Systems, 4th Edn. New Jersey: World Scientific, (2012). 
[3]  B. J. Berne, G. Ciccotti, and D. F. Coker, Eds., Classical and Quantum Dynamics in Condesed Phase Simulations, New Jersey: World Scientific, (1998). 
[4]  V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, 3rd Edn. Weinheim: WileyVCH, (2011). 
[5]  D. E. Makarov, and N. Makri, Chem. Phys. Lett. 221 , 482 (1994). DOI:10.1016/00092614(94)002754 
[6]  N. Makri, and D. E. Makarov, J. Chem. Phys. 102 , 4600 (1995a). DOI:10.1063/1.469508 
[7]  N. Makri, and D. E. Makarov, J. Chem. Phys. 102 , 4611 (1995b). DOI:10.1063/1.469509 
[8]  H. Wang, M. Thoss, and W. H. Miller, J. Chem. Phys. 115 , 2979 (2001). DOI:10.1063/1.1385561 
[9]  H. Wang, and M. Thoss, J. Chem. Phys. 119 , 1289 (2003). DOI:10.1063/1.1580111 
[10]  Y. Tanimura, and R. Kubo, J. Phys. Soc. Jpn. 58 , 101 (1989). DOI:10.1143/JPSJ.58.101 
[11]  A. Ishizaki, and G. R. Fleming, J. Chem. Phys. 130 , 234111 (2009). DOI:10.1063/1.3155372 
[12]  J. C. Tully, Faraday Discuss. 110 , 407 (1998). DOI:10.1039/a801824c 
[13]  J. C. Tully, J. Chem. Phys. 137 , 22A301 (2012). DOI:10.1063/1.4757762 
[14]  M. Thoss, H. B. Wang, and W. H. Miller, J. Chem. Phys. 115 , 2991 (2001). DOI:10.1063/1.1385562 
[15]  T. C. Berkelbach, D. R. Reichman, and T. E. Markland, J. Chem. Phys. 136 , 034113 (2012). DOI:10.1063/1.3671372 
[16]  A. MontoyaCastillo, T. C. Berkelbach, and D. R. Reichman, J. Chem. Phys. 143 , 194108 (2015). DOI:10.1063/1.4935443 
[17]  W. T. Pollard, A. K. Felts, and R. A. Friesner, Adv. Chem. Phys. 93 , 77 (1996). 
[18]  H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press on Demand, (2002). 
[19]  Y. Yan, and R. Xu, Annu. Rev. Phys. Chem. 56 , 187 (2005). DOI:10.1146/annurev.physchem.55.091602.094425 
[20]  K. Blum, Density Matrix Theory and Applications, Springer Science & Business Media, (2013). 
[21]  S. Nakajima, Prog. Theor. Phys. 20 , 948 (1958). DOI:10.1143/PTP.20.948 
[22]  R. Zwanzig, J. Chem. Phys. 33 , 1338 (1960). DOI:10.1063/1.1731409 
[23]  R. Zwanzig, Physica 30 , 1109 (1964). DOI:10.1016/00318914(64)901028 
[24]  H. Mori, Prog. Theor. Phys. 33 , 423 (1965). DOI:10.1143/PTP.33.423 
[25]  A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59 , 1 (1987). DOI:10.1103/RevModPhys.59.1 
[26]  G. D. Mahan, ManyParticle Physics, New York: Springer Science & Business Media, (2013). 
[27]  E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83 , 349 (2011). DOI:10.1103/RevModPhys.83.349 
[28]  A. Royer, Phys. Lett. A 315 , 335 (2003). DOI:10.1016/S03759601(03)010521 
[29]  D. R. Reichman, F. L. H. Brown, and P. Neu, Phys. Rev. E 55 , 2328 (1997). DOI:10.1103/PhysRevE.55.2328 
[30]  M. Aihara, H. M. Sevian, and J. L. Skinner, Phys. Rev. A 41 , 6596 (1990). DOI:10.1103/PhysRevA.41.6596 
[31]  Q. Shi, and E. Geva, J. Chem. Phys. 119 , 12063 (2003). DOI:10.1063/1.1624830 
[32]  E. Y. Wilner, H. Wang, M. Thoss, and E. Rabani, Phys. Rev. B 89 , 205129 (2014). DOI:10.1103/PhysRevB.89.205129 
[33]  J. Cerrillo, and J. Cao, Phys. Rev. Lett. 112 , 110401 (2014). DOI:10.1103/PhysRevLett.112.110401 
[34]  M. Xu, Y. Yan, Y. Liu, and Q. Shi, J. Chem. Phys. 148 , 164101 (2018). DOI:10.1063/1.5022761 
[35]  M. Tokuyama, and H. Mori, Prog. Theor. Phys. 56 , 1073 (1976). DOI:10.1143/PTP.56.1073 
[36]  N. Hashitsumae, F. Shibata, and M. Shingū, J. Stat. Phys. 17 , 155 (1977). DOI:10.1007/BF01040099 
[37]  F. Shibata, Y. Takahashi, and N. Hashitsume, J. Stat. Phys. 17 , 171 (1977). DOI:10.1007/BF01040100 
[38]  S. Chaturvedi, and F. Shibata, Z. Phys. B 35 , 297 (1979). DOI:10.1007/BF01319852 
[39]  H. Grabert, P. Schramm, and G. L. Ingold, Phys. Rep. 168 , 115 (1988). DOI:10.1016/03701573(88)900233 
[40]  B. Garraway, Phys. Rev. A 55 , 2290 (1997). DOI:10.1103/PhysRevA.55.2290 
[41]  M. Murao, and F. Shibata, J. Phys. Soc. Jpn. 64 , 2394 (1995). DOI:10.1143/JPSJ.64.2394 
[42]  A. Smirne, H. P. Breuer, J. Piilo, and B. Vacchini, Phys. Rev. A 82 , 062114 (2010). DOI:10.1103/PhysRevA.82.062114 
[43]  W. M. Zhang, P. Y. Lo, H. N. Xiong, M. W. Y. Tu, and F. Nori, Phys. Rev. Lett. 109 , 170402 (2012). DOI:10.1103/PhysRevLett.109.170402 
[44]  G. Nan, Q. Shi, and Z. Shuai, J. Chem. Phys. 130 , 134106 (2009). DOI:10.1063/1.3108521 
[45]  L. Kidon, E. Y. Wilner, and E. Rabani, J. Chem. Phys. 143 , 234110 (2015). DOI:10.1063/1.4937396 
[46]  S. Jang, J. Cao, and R. J. Silbey, J. Chem. Phys. 116 , 2705 (2002). DOI:10.1063/1.1445105 
[47]  G. Gasbarri, and L. Ferialdi, Phys. Rev. A 97 , 022114 (2018). DOI:10.1103/PhysRevA.97.022114 
[48]  D. Chruściński, and A. Kossakowski, Phys. Rev. Lett. 104 , 070406 (2010). DOI:10.1103/PhysRevLett.104.070406 
[49]  D. Egorova, M. Thoss, W. Domcke, and H. Wang, J. Chem. Phys. 119 , 2761 (2003). DOI:10.1063/1.1587121 
[50]  M. Schröder, U. Kleinekathöfer, and M. Schreiber, J. Chem. Phys. 124 , 084903 (2006). DOI:10.1063/1.2171188 
[51]  S. Mukamel, Principles of Nonlinear Optical Spectroscopy. New York: Oxford (1995). 
[52]  R. Doll, D. Zueco, M. Wubs, S. Kohler, and P. Hänggi, Chem. Phys. 347 , 243 (2008). DOI:10.1016/j.chemphys.2007.09.003 
[53]  A. Ishizaki, and Y. Tanimura, J. Phys. Soc. Jpn. 74 , 3131 (2005). DOI:10.1143/JPSJ.74.3131 
[54]  Q. Shi, L. P. Chen, G. J. Nan, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 130 , 084105 (2009). DOI:10.1063/1.3077918 
[55]  M. Xu, L. Song, K. Song, and Q. Shi, J. Chem. Phys. 146 , 064102 (2017). DOI:10.1063/1.4974926 
[56]  M. Sparpaglione, and S. Mukamel, J. Chem. Phys. 88 , 3263 (1988). DOI:10.1063/1.453922 
[57]  M. Sparpaglione, and S. Mukamel, J. Chem. Phys. 88 , 4300 (1988). DOI:10.1063/1.453789 
[58]  H. T. Chen, T. C. Berkelbach, and D. R. Reichman, J. Chem. Phys. 144 , 154106 (2016). DOI:10.1063/1.4946809 
[59]  M. G. Mavros, and T. Van Voorhis, J. Chem. Phys. 141 , 054112 (2014). DOI:10.1063/1.4891669 
[60]  R. Zwanzig, Phys. Rev. 124 , 983 (1961). DOI:10.1103/PhysRev.124.983 
[61]  S. Chakravarty, and A. J. Leggett, Phys. Rev. Lett. 52 , 5 (1984). DOI:10.1103/PhysRevLett.52.5 
[62]  Q. Shi, L. P. Chen, G. J. Nan, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 130 , 164518 (2009). DOI:10.1063/1.3125003 
[63]  Y. Tanimura, Phys. Rev. A 41 , 6676 (1990). DOI:10.1103/PhysRevA.41.6676 
[64]  Y. Tanimura, J. Phys. Soc. Jpn. 75 , 082001 (2006). DOI:10.1143/JPSJ.75.082001 
[65]  R. Xu, and Y. Yan, Phys. Rev. E 75 , 031107 (2007). DOI:10.1103/PhysRevE.75.031107 
[66]  Y. Tanimura, and S. Mukamel, J. Phys. Soc. Jpn. 63 , 66 (1994). DOI:10.1143/JPSJ.63.66 
[67]  A. Ishizaki, and Y. Tanimura, Chem. Phys. 347 , 185 (2008). DOI:10.1016/j.chemphys.2007.10.037 
[68]  Y. A. Yan, F. Yang, Y. Liu, and J. S. Shao, Chem. Phys. Lett. 395 , 216 (2004). DOI:10.1016/j.cplett.2004.07.036 
[69]  Y. Zhou, Y. A. Yan, and J. S. Shao, Europhys. Lett. 72 , 334 (2005). DOI:10.1209/epl/i2005102624 
[70]  Y. Zhou, and J. S. Shao, J. Chem. Phys. 128 , 034106 (2008). DOI:10.1063/1.2818095 
[71]  S. Jang, J. Chem. Phys. 131 , 164101 (2009). DOI:10.1063/1.3247899 
[72]  D. P. McCutcheon, and A. Nazir, New. J. Phys. 12 , 113042 (2010). DOI:10.1088/13672630/12/11/113042 
[73]  H. T. Chang, P. P. Zhang, and Y. C. Cheng, J. Chem. Phys. 139 , 224112 (2013). DOI:10.1063/1.4840795 
[74]  M. Cho, H. M. Vaswani, T. Brixner, J. Stenger, and G. R. Fleming, J. Phys. Chem. B 109 , 10542 (2005). DOI:10.1021/jp050788d 
[75]  S. I. Vulto, M. A. de Baat, R. J. Louwe, H. P. Permentier, T. Neef, M. Miller, H. van Amerongen, and T. J. Aartsma, J. Phys. Chem. B 102 , 9577 (1998). 