Chinese Journal of Chemical Physics  2018, Vol. 31 Issue (4): 575-583

#### The article information

Yan-ying Liu, Ya-ming Yan, Meng Xu, Kai Song, Qiang Shi

Exact Generator and its High Order Expansions in Time-Convolutionless Generalized Master Equation: Applications to Spin-Boson Model and Exictation Energy Transfer

Chinese Journal of Chemical Physics, 2018, 31(4): 575-583

http://dx.doi.org/10.1063/1674-0068/31/cjcp1806146

### Article history

Accepted on: July 27, 2018
Exact Generator and its High Order Expansions in Time-Convolutionless Generalized Master Equation: Applications to Spin-Boson Model and Exictation Energy Transfer
Yan-ying Liu, Ya-ming Yan, Meng Xu, Kai Song, Qiang Shi
Dated: Received on June 20, 2018; Accepted on July 27, 2018
Beijing National Laboratory for Molecular Sciences, State Key Laboratory for Structural Chemistry of Unstable and Stable Species, CAS Research/Education Center for Excellence in Molecular Sciences, Institute of Chemistry, Chinese Academy of Sciences, Beijing 100190, China; University of Chinese Academy of Sciences, Beijing 100049, China
*Author to whom correspondence should be addressed. Qiang Shi, E-mail:qshi@iccas.ac.cn
Abstract: The time-convolutionless (TCL) quantum master equation provides a powerful tool to simulate reduced dynamics of a quantum system coupled to a bath. The key quantity in the TCL master equation is the so-called kernel or generator, which describes effects of the bath degrees of freedom. Since the exact TCL generators are usually hard to calculate analytically, most applications of the TCL generalized master equation have relied on approximate generators using second and fourth order perturbative expansions. By using the hierarchical equation of motion (HEOM) and extended HEOM methods, we present a new approach to calculating the exact TCL generator and its high order perturbative expansions. The new approach is applied to the spin-boson model with different sets of parameters, to investigate the convergence of the high order expansions of the TCL generator. We also discuss circumstances where the exact TCL generator becomes singular for the spin-boson model, and a model of excitation energy transfer in the Fenna-Matthews-Olson complex.
Key words: Time-convolutionless generalized master equation    Spin-boson model    Exact generator    High order expansion
Ⅰ. INTRODUCTION

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 [1-4]. Numerically exact approaches [5-11] are usually limited to some simple models, while approximate approaches [12-16] 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 [17-20].

For calculations of the dynamics within the reduced system dynamics framework, the NakajimaZwanzig (NZ) (or time-convolution, TC) master equation, an integro-differential equation with a superoperator called "memory kernel" derived from the full Hilbert space using the projection operator techniques, is often employed [21-24]. 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 spin-boson [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 system-bath coupling or a small coupling parameter in the system Hamiltonian [28-30]. It is noted that, numerical methods to calculate exact NZ memory kernel [31-33], 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 integro-differential equation to describe the bath effects, an alternative form of the generalized master equation is the Hashitsume-ShibataTakahashi or time-convolutionless (TCL) approach [35-37]. Also being formally exact, this approach describes the evolution of the RDO using a first order differential equation with time dependent coefficients [18, 36-38]. The key quality in the TCL generalized master equation is a super-operator 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 super-operator 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 non-perturbative 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 multi-dimensional integrals. It is then applied to the spin-boson 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 spin-boson 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 spin-boson 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 time-convolutionless master equation

We 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{eq-hamiltonian} \end{array}$ (1)

where $H_\textrm{S}$ is the system Hamiltonian, $H_\textrm{B}$ is the bath Hamiltonian, and $H_{\textrm{SB}}$ denotes the system-bath coupling. Evolution of the total system and bath density matrix is given by the quantum Liouville equation,

 $\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 ${\rho}_\textrm{T}$ is the density operator of the total system, ${\cal L} \rho_\textrm{T}(t)$=$[H_\textrm{T}, \rho_\textrm{T}(t)]$, and we set $\hbar$=1 throughout this work. The reduced system density operator $\rho_\textrm{S}(t)$ is defined as the partial trace over the bath degrees of freedom, $\rho_\textrm{S}(t)$=${\rm Tr}_\textrm{B} {\rho}_\textrm{T}(t)$.

To derive the generalized master equation, we separate the density operator into relevant and irrelevant parts by means of a projection operator $\mathcal{P}$. The residual projection operator ${\mathcal Q}$ is defined as ${\mathcal Q}$=1$-$${\mathcal P}. Both of {\mathcal P} and {\mathcal Q} obey the fundamental requirement of a projection, {\mathcal P}^2={\mathcal P}, {\mathcal Q}^2={\mathcal Q}. Then the TCL form of GQME [18, 36, 37] can be written as  \begin{array}{l} \label{eq-tclprojector1} \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} (t-s)}{\mathcal Q} {\cal L} {\mathcal P}\textrm{e}^{i {\cal L} (t-s)}\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} (t-t_0)} {\mathcal Q} \;\; \end{array} (5) and \rho_\textrm{T}(t_0) is the initial density operator of the total system. In this work, we choose a factorized initial condition \rho_\textrm{T}(0)=\rho_\textrm{S}(0)$$\otimes$$\rho_\textrm{B}(0) (we set t_0=0), and assume that {\mathcal P} \rho_\textrm{T}(0)=\rho_\textrm{T}(0), such that the inhomogeneous term {\mathcal I}(t) vanishes [44]. It is also assumed that \mathcal{PL}\mathcal{P}=0, which holds for the projection operator used later. Eq.(3) thus can be written as  \begin{array}{l} \label{eq-tclprojector2} \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 [56-58] for the total density operator, \rho_\textrm{T}$$\equiv$$|\rho_\textrm{T}\rangle\rangle, and define the product \langle \langle A | B \rangle \rangle$$\equiv$${\rm Tr}_\textrm{S} {\rm Tr}_\textrm{B} \{A^\dagger B\} where {\rm Tr}_\textrm{S} and {\rm Tr}_\textrm{B} are partial traces over the system and bath degrees of freedom, respectively. We denote an eigenstate of the quantum subsystem by | j \rangle , a Liouville space state of the quantum subsystem is then given by | jk \rangle \rangle=| j\rangle \langle k|. For simplicity, we use |j\rangle\rangle to denote |jj\rangle\rangle, and \langle\langle j|=\{|j\rangle\rangle\}^{\dagger }. The population of the quantum subsystem on state |j\rangle can be written as P_{j}(t)={\rm Tr}_\textrm{S}{\rm Tr}_{\textrm{B}}\{|j\rangle\langle j|\rho_{\textrm{T}}(t) \}=\langle\langle j|\rho_{\textrm{T}}(t)\rangle\rangle. We employ the following projection operator which has been used in many previous studies [56, 58-60]  \begin{array}{l} \label{projector} {\mathcal P} = \sum \limits_j | j \rho_j^\textrm{B} \rangle \rangle \langle \langle j | \;\; \end{array} (7) where | j\rho_j^\textrm{B}\rangle\rangle=| j \rangle \langle j |$$\otimes$$\rho_j^\textrm{B} , and the locally equilibrated (to be specified later) bath density operator \rho_j^\textrm{B} is taken to be associated with the system eigenstate |j\rangle . Thus in the reduced subsystem space, Eq.(6) can be written as  \begin{array}{l} \label{eq-tclpopulation1} \frac{\textrm{d}}{\textrm{d}t}P (t) = {\cal R}(t)P (t)\;\; \end{array} (8) where P(t) is a vector whose jth element is the population on state | j \rangle, and \mathcal{R}(t) is now a matrix whose elements are defined by {\cal R}_{jk}(t)=-i\langle\langle j|\mathcal{L}\theta(t) |k\rho_{k}^{B}\rangle\rangle. Considering the time evolution of {\rho}_\textrm{T}(t), we have  \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{eq-p} \end{array} (9) With the above projector operator \cal P and the relationship {\cal Q} {\rho}_\textrm{T}(0)=0, and in the reduced subsystem space, Eq.(9) can be written as  \begin{array}{l} P(t) = {\cal U}_\textrm{S}(t)P(0) \label{eq-evolution1} \end{array} (10) where the matrix element of \mathcal{U}_{\textrm{S}},  \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{eq-usdefinition} \end{array} (11) involves the propagator of the system only. Taking the time derivative operation of P(t) in Eq.(10) and making the inversion of {\mathcal U}_\textrm{S}(t), we get  \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{eq-tclpopulation2} \end{array} (12) provided the inverse does exist. By comparing with the TCL quantum master equation in Eq.(8), the TCL generator {\cal R}(t) is given by  \begin{array}{l} \label{eq-rt} {\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 {\cal U}_\textrm{S}(t) and the time inversion of {\cal U}_\textrm{S}(t). To get the concrete form of \dot{\cal U}_\textrm{S}(t), we take the derivative of Eq.(11) and obtain:  \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{eq-dus} \end{array} (14) In the above derivations, the condition {\cal P} {\cal L} {\cal P}=0 has been used. B. The spin-boson model and HEOM method For simplicity, we restrict ourselves to the spin-boson model for most part of this work. The spin-boson Hamiltonian, in the form of Eq.(1), describes a two-state 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 \epsilon and the coupling of the two states \Delta. The Hamiltonian of harmonic bath and the coupling to the system are  \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, {\sigma}_z and {\sigma}_x are the Pauli matrices. x_j and p_j are the mass-weighted coordinate and momentum of the jth bath mode with frequency {\omega}_j, respectively. The coupling coefficient between the system operator {\sigma}_z and the jth bath mode coordinate x_j is written as c_j. F=\displaystyle\sum_jc_jx_j is the collective bath coordinate. We set {\beta}=1/(k_\textrm{B}T) throughout this paper. The system-bath 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 C(t) is related to the spectral density via the fluctuation dissipation theorem [2, 62]:  \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{eq-corelation} \end{array} (19) where Z_\textrm{B}={\rm Tr} \textrm{e}^{-{\beta}H_\textrm{B}} is the partition function of the uncoupled harmonic bath. In this work, we will employ the Debye spectral density  \begin{array}{l} J(\omega) = \frac{\eta {\omega}_\textrm{c} \omega}{{\omega}^2+{{\omega}_\textrm{c}}^2}\;\; \end{array} (20) where \eta describes the coupling strength between system and bath, \omega_\textrm{c} is the cut-off frequency of the bath. Thus C(t) in Eq.(19) can be written into a sum of exponential decaying functions in time:  \begin{array}{l} C(t>0) = \sum\limits_k d_k\textrm{e}^{-{{\omega}}_kt} \label{eq-corelation1} \end{array} (21) where {{\omega}}_0={{\omega}_\textrm{c}} is the longitudinal relation constant, {{\omega}}_k=2k{\pi}/{\beta} is the Matsubara frequency and  \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, 63-67] or the stochastic Liouville equation approach [68-70]:  \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{eq-heom1} \end{array} (24) The subscript \textbf{n} denotes a set of index \{n_1, n_2, ..., n_k, ...\}, with the integer number n_k$$\geq$0 associated with the $k$th exponential terms in $C(t)$ of Eq.(21). The subscript ${\textbf{n}}_k^{\pm}$ differs from $\textbf{n}$ only by changing the specified $n_k$ to $n_k\pm$1. ${\rho}_{\textbf{0}}$ with $\textbf{0}$=$\{0, 0, ...\}$ corresponds to the RDO and the other ${\rho}_{\textbf{n}}s$ are the auxiliary density operators (ADOs).

C. The exact TCL generator and its high order expansions

As shown in Eq.(13), to obtain the exact TCL generator, we first need to calculate the value of ${\cal U}_\textrm{S}(t)$, defined in Eq.(11). The matrix element of ${\cal U}_\textrm{S}(t)$ can be written as

 $\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{eq-us2} \end{array}$ (25)

where $\sigma_k(t)$=${\rm Tr}_\textrm{B} \{ \textrm{e}^{-i{\cal L}t} | k \rangle \langle k | \otimes \rho_k^\textrm{B} \}$ and the subscript $jj$ denotes the $j$th diagonal matrix element of the system RDO. It is noted that the initial state of bath to compute $[\sigma_k(t)]_{jj}$ is the relaxed equilibrium associated with the $k$th state of system. The following projection operator for the spin-boson model is employed,

 $\begin{array}{l}\label{Eq:sb-projction} \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 $\rho_{j}^{\textrm{B}}$ is defined as the locally equilibrated bath density operator associated with state $|j\rangle$, and $|1\rangle$ denotes donor state, $|2\rangle$ denotes acceptor state.

 $\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 $H^{(j)}$=$\pm({\epsilon}$+$\sum\limits_{\alpha}c_{\alpha}x_{\alpha})$+$H_\textrm{B}$ (+ for state $|1\rangle$, $-$ for state $|2\rangle$, respectively). Therefore, within the spin-boson model, we can calculate ${\cal U}_\textrm{S}(t)$ using the HEOM approach presented in the above subsection. By taking time inversion of ${\cal U}_\textrm{S}(t)$, we can get ${\cal U}^{-1}_\textrm{S}(t)$.

To calculate the derivative of ${\cal U}_\textrm{S}(t)$, we define

 $\begin{array}{l} H_\textrm{S} = H_0 + H_1 \end{array}$ (28)

where $H_0$=$\epsilon {\sigma}_z$ and $H_1$=$\Delta {\sigma}_x$. The relation ${\mathcal Q}{\cal L}{\mathcal P}$=${\cal L}_1 {\mathcal P}$ holds for the above defined projection operator in Eq.(26). $\dot{\cal U}_{\textrm{S};jk}(t)$ defined in Eq.(14) can then be calculated using ${\cal L}_1 | k\rangle\rangle$=${(-1)}^ki{\Delta}{\sigma}_y$,

 $\begin{array}{l} {[\dot{\cal U}_\textrm{S}(t)]}_{jk} = {(-1)}^{j+k+1}{\Delta}{[{\sigma}_z{\sigma}_ {y, k}(t)]}_{jj} \;\; \label{eq-us3} \end{array}$ (29)

where ${\sigma}_{y, t}(t)$=${\rm Tr}_\textrm{B}\{\textrm{e}^{-i{\cal L}t}{\sigma} _y{\rho}_k^\textrm{B}\}$. Similar to ${\cal U}_\textrm{S}(t), \, \dot{\cal U}_\textrm{S}(t)$ can also be computed using the HEOM method. Therefore, by calculating $\dot{\cal U}_\textrm{S}(t)$ and ${\cal U}^{-1}_\textrm{S}(t)$, it is easy to get the value of the exact TCL generator with Eq.(13).

In the rest part of this section, we will deduce the high order expansions of ${\cal R}(t)$. By using the projection operator in Eq.(7), the coupling between the two states $\Delta$ is used as the expansion parameter. We first write the perturbative expansions of ${\cal R}(t)$, ${\cal U}_\textrm{S}(t)$, and $\dot{\cal U}_\textrm{S}(t)$, with respect to the inter-state coupling ${\Delta}$:

 $\begin{array}{l} {\cal R}(t) = \sum\limits_{N=2}^{\infty}{\Delta}^N{\cal R}^{(N)}(t) \end{array}$ (30)
 $\begin{array}{l} \label{eq-us5} {\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 ${(1+x)}^{-1}$=$\displaystyle\sum_{n=0}^{\infty} {(-1)}^nx^n$ (under the assumption that $| x |$$<1), we can invert Eq.(31) as  \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 \Delta results in an infinite series of equations:  \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 {\cal R}(t) are all zero. With the above equations, the following relations can be obtained for the even order terms:  \begin{array}{l} \label{eq-recursion} {\cal R}^{(2n)}(t) = \dot{\cal U}_\textrm{S}^{(2n)}(t)- \sum\limits_{m=1}^{n-1} [{\cal R}^{(2m)}(t) {\cal U}_\textrm{S}^{[2(n-m)]}(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 {\cal U}_\textrm{S}(t) can be calculated directly using the method presented in Ref.[55] with proper initial conditions. Expansions of \dot{\cal U}_\textrm{S}(t) with respect to \Delta can be calculated by expanding {\sigma}_{y, k}(t) defined in Eq.(29) into Taylor series:  \begin{array}{l} {\sigma}_{y, k}(t) = \sum\limits_{N=0}^{\infty} \frac{1}{N!} {\sigma}_{y, k}^{(N)}(t){\Delta}^N \;\; \label{eq-sigma} \end{array} (39) The {\sigma}_{y, k}^{(N)}(t) terms can be obtained using the extended HEOM method [55] with the initial condition {\rho}_\textrm{T}={\sigma}_y$$\otimes$${\rho}_k^\textrm{B}:  \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}}^{(N-1)}(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{eq-heom2} \end{array} (40) where {\cal L}_0 {\rho}=[{\epsilon}{\sigma}_z, {\rho}]. Ⅲ. RESULTS For all numerical calculations in this section, the initial state of the total system is assumed to be equilibrated on state | 1\rangle, i.e., {\rho}_\textrm{T}(0)=| 1 \rangle \langle 1|$$\otimes$${\rho}_1^\textrm{B}. Thus the inhomogeneous term in Eq.(3) vanishes. The HEOM approach is used to calculate {\cal U}_\textrm{S}(t) with Eq.(25) and the benchmark numerical exact population dynamics. Population dynamics of TCL quantum master equation is calculated from the differential equation in Eq.(8), which is solved by the fourth order Runge-Kutta method. The high order expansion terms of the TCL generator {\cal R}(t) are computed by the extended HEOM approach with Eqs. (29), (40), and (38). There are four numerical examples in this section. The first three are based on the spin-boson model with different parameters, including the system-bath coupling strength \eta, the cut-off frequency {\omega}_\textrm{c}, the intra-state coupling \Delta, the inverse temperature \beta, and the energetic bias \epsilon, while the last one is for the excitation energy transfer in the Fenna-Matthews-Olson (FMO) complex. The first two examples describe circumstances where the high order expansions of the TCL generator converge and diverge, respectively. The possible singularity of the exact generator is shown in the third example for the spin-boson model and the fourth example for the FMO complex. A. Convergence of high order expansions of the generator 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 {\beta}=0.5, {\omega}_\textrm{c}=5, {\eta}=5, {\Delta}=1, and {\epsilon}=0. It is noted that the population dynamics resulting from the two different approaches is exactly same, thus validating the method to obtain the exact generator. Population dynamics using the second order TCL generalized master equation is also shown for comparison. In this example, results from the second order TCL approximation also agree well with the exact ones.  FIG. 1 The time dependent population on state |1\rangle in the spin-boson model, obtained from the HEOM method, the exact TCL generalized master equation, and the TCL equation with the second order approximate generator, respectively. The parameters are {\beta}=0.5, {\omega}_\textrm{c}=5, {\eta}=5, {\Delta}=1, and {\epsilon}=0. The high order expansions up to the 12th order of the TCL generator {\cal R}(t) are shown in FIG. 2(a), and the fine details are shown in the inset. The corresponding summations of the nth order generators are given in FIG. 2(b). FIG. 2(a) shows that the amplitudes of high order terms decrease quickly when the perturbation order increases, making the expansion of the generator converge easily, as shown in FIG. 2(b).  FIG. 2 Panel (a) shows the second and high order expansion terms of the TCL generator and the fine details are presented in the inset. Panel (b) shows the corresponding summations to certain orders. The parameters are {\beta}=0.5, {\omega}_\textrm{c}=5, {\eta}=5, {\Delta}=1, and {\epsilon}=0. In the second example, the following parameters are used: {\beta}=1, {\omega}_\textrm{c}=1, {\eta}=2, {\Delta}=1, and {\epsilon}=0, which correspond to the case of weak system-bath coupling and a slow bath. FIG. 3 shows the exact population dynamics using the HEOM method, and the TCL generalized master equations with the exact and second order generators. It can be seen that, after a short time, results from the second order TCL generalized master equation quickly deviate from the exact ones.  FIG. 3 The time dependent population on state |1\rangle in the spin-boson model, obtained from the HEOM method, the exact TCL generalized master equation, and the TCL equation with the second order approximate generator, respectively. The parameters are {\beta}=1, {\omega}_\textrm{c}=1, {\eta}=2, {\Delta}=1, and {\epsilon}=0. 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).  FIG. 4 Panel (a) shows the second and high order expansion terms of the TCL generator. Panel (b) shows the corresponding summations to certain orders. The parameters are {\beta}=1, {\omega}_\textrm{c}=1, {\eta}=2, {\Delta}=1, and {\epsilon}=0. 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 {\cal R}(2) with the fourth order one {\cal R}(4). This has been discussed previously for the NZ form of the generalized master equations or polaron-transformed GQMEs [34, 58, 71-73]. The general trend is also applicable for the TCL generators. To explore the convergence of the high order TCL generators in different parameter regimes, we also calculate the critical values of interstate coupling {\Delta}_\textrm{c} as a function of the system-bath coupling strength \eta. The high order expansions converge when \Delta$$ <$$\Delta_\textrm{c}$. The convergence criterion is set to

 $\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 $n$ up to 28, where $S_{11}^{(n)}(t/{\pi})$ represents the first matrix element of the summation of the $n$th order generators at time $t$. The other parameters we used are ${\beta}$=0.5, ${\omega}_\textrm{c}$=5, and ${\epsilon}$=0. The results are shown in FIG. 5. It can be seen that, with the increase of the system-bath coupling strength $\eta$, $\Delta_\textrm{c}$ also increases.

 FIG. 5 The critical interstate coupling constant $\Delta_\textrm{c}$ for converged expansion of the TCL generator, as a function of the system-bath coupling strength $\eta$. The other parameters are ${\beta}$=0.5, ${\omega}_\textrm{c}$=5, and ${\epsilon}$=0.
B. Singularity of the exact generator ${\cal R}(t)$

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 ${\cal U}_\textrm{S}(t)$ does exist but there are some circumstances in which ${\cal U}_\textrm{S}(t)$ is not a full rank matrix, resulting in the singularity points of the generator at certain time $t$.

In this subsection, we show such an example by setting ${\omega}_\textrm{c}$=1, ${\eta}$=1, and the other parameters are the same as those in the second example. We can see from FIG. 6(b) that with this set of parameters, ${\cal R}(t)$ has two singularities. Comparing FIG. 6 (a) with (b), it is noted that when $P_{11}$ equals to $P_{22}$, the singularities appear.

 FIG. 6 Panel (a) shows the time dependent population on state $|1\rangle$ and $|2\rangle$ in the spin-boson model. Panel (b) shows the corresponding numerical result of the exact generator. The parameters are: ${\beta}$=1, ${\omega}_\textrm{c}$=1, ${\eta}$=1, ${\Delta}$=1, and ${\epsilon}$=0.

This result can be understood using the following theoretical analysis. For the symmetric spin-boson model, the ${\cal U}_\textrm{S}(t)$ and ${\cal R}(t)$ satisfy the following relations:

 $\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{eq-singularity} \end{array}$ (43)

When $a(t)$=1/2, the inverse of ${\cal U}_\textrm{S}(t)$ does not exist. However, in this case, by using the above Eq.(43), the TCL generalized master equation in Eq.(8) predicts that $\dot P_i(t)$=0 can not hold except in the equilibrium state. So for the symmetric spin-boson model, whenever the population curves cross, there will be singularity in the exact generator.

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.

 FIG. 7 Panel (a) shows the second and high order (4th to 12th) expansion terms of the TCL generator. Panel (b) shows the corresponding summation terms and the fine details are presented in the inset. The parameters are: ${\beta}$=1, ${\omega}_\textrm{c}$=1, ${\eta}$=1, ${\Delta}$=1, and ${\epsilon}$=0.

It is interesting to investigate whether similar problem of singular generators also exists in more general model systems beyond the simple symmetric spin-boson 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 $\mathcal{R}_{11}$. It is found that there are also several singularity points in the TCL generator. So that although the exact TCL generalized master equation is very appealing in the simulation of open system dynamics, the issue of singularity in the exact TCL generator may exist in many problems. In contrast, the memory kernel in the NZ generalized master equation is usually well behaved.

 FIG. 8 Panel (a) shows the time dependent population on sites 1 to 7 of the FMO complex. The initial state is equalibrated on site 1. Panel (b) shows the corresponding numerical result of the 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 DISCUSSIONS

In 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 time-consuming multi-dimensional integrals.

By using the spin-boson 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 Jaynes-Cummings 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 spin-boson 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.

Ⅴ. ACKNOWLEDGEMENTS

This 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).

Reference
 [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: Wiley-VCH, (2011). [5] D. E. Makarov, and N. Makri, Chem. Phys. Lett. 221 , 482 (1994). DOI:10.1016/0009-2614(94)00275-4 [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. Montoya-Castillo, 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/0031-8914(64)90102-8 [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, Many-Particle 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/S0375-9601(03)01052-1 [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/0370-1573(88)90023-3 [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/i2005-10262-4 [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/1367-2630/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).