Generalized Quantum Master Equation: A Tutorial Review and Recent Advances^†

Ⅰ.   INTRODUCTION

The dynamics in many important chemical processes, such as excitation energy, charge, and quantum coherence transfer in photosynthesis, photovoltaics, organic light-emitting diodes (OLEDs) as well as reactive photochemistry, can be described in terms of dynamics in open quantum systems [1, 2]. Open quantum system refers to a quantum system embedded in a quantum environment. The system is usually chosen to be the interested or experimentally accessible degrees of freedom (DOF) such as the electronic states in light-harvesting complexes or a tagged vibrational mode in a condensed-phase matter. The environment or the bath represents the rest of DOF that couple with the system and thus affect the system dynamics.

For example, in nonadiabatic dynamics, the electronic DOF are viewed as the system and the nuclear DOF serve as the thermal bath, and the dynamics of the electronic reduced density matrix (RDM) is the target reduced dynamics, which directly reports the population and coherence transfer. Quantum dynamics in condensed-phase systems can be used to theoretically describe the electronic energy and charge transfer pathways along with the reaction mechanisms in complex systems. Fundamental understanding of open quantum systems allows us to predict and control quantum systems in various applications ranging from quantum optics [3], quantum information and quantum computing [4, 5], quantum thermodynamics [6-9], quantum biology [10, 11], and hybrid quantum device [12].

However, it is usually not possible to solve the equation of motion for the overall system including the system and the bath due to the exponentially growing complexity with increasing DOF. Thus, it is useful and practical to formulate the problem to solve for the dynamics of the interested part only, of course, with the influence of the surroundings at minimum possible numerical cost. Projection operator technique, in particular, is designed for filtering out the irrelevant part from the interested information. Applying the projection operators to the equation of motion (i.e., quantum Liouville equation) for open quantum systems leads to the generalized quantum master equation (GQME) that describes the reduced system dynamics formally exactly with the influence from the bath absorbed in the so-called memory kernel.

The main goal of this work is to provide a tutorial review and along the way present and summarize recent advances on GQME methods. It is important to note that throughout the literature spanning many disciplines, GQME sometimes goes by other names such as the non-Markovian master equation, quantum dissipative theory [13], or sometimes simply referred to as one of the master equation [14]. Here, we restrict our introduction and discussion to focus on the non-Markovian GQME and will not cover the details of Markovian quantum master equation (QME) theory. Discussions regarding the relationship between the non-Markovian GQME and the Markovian QME can be found in Refs.[14-17].

The remainder of this review is organized as follows. In Section Ⅱ, we briefly sketch a historical background of GQME followed by the current status. In Section Ⅲ, we provide the derivation of the two forms of GQME including time-convoluted Nakajima-Zwanzig (NZ) GQME and time-convolutionless (TCL) GQME using the projection operators technique. In Section Ⅳ, we focus on introducing and discussing various schemes for evaluating the memory kernel or the TCL generator, as well as the inhomogeneous term in the GQME. In Section Ⅴ, we proceed with summarizing the expressions of the molecular projection-free inputs in terms of system-dependent time correlation function (TCF) and survey the common quantum dynamics numerical methods that can be used to evaluate them. In Section Ⅵ, we present several topical discussions on emerging trends and technology that can be used to accelerate and improve the utility of GQME. In Section Ⅶ, concluding remarks and future outlook are given. Appendix A provides the mathematical tools necessary for following the derivations presented in this review. Appendix B shows the proof for deriving the simplification of the memory kernel in the NZ-GQME with system-bath Hamiltonian.

Ⅱ.   BACKGROUND OF GQME

A   Historical development of GQME

Quantum master equations (QME) are a family of methods for obtaining the reduced dynamics in open quantum systems under the Markovian assumption that the system is memoryless on its previous dynamics history. The first QME was originally proposed by Pauli, which was based on the Born approximation for which there exists a separation of timescales between the system dynamics and the relaxation of the environment in addition to the Markovian assumption [18]. Later, more practical and popular forms of QME were proposed such as the Bloch-Redfield master equation [19, 20] and the Lindblad master equation [21, 22].

The generalization of the QME to the non-Markovian case is referred to as the GQME. A non-Markovian system has a memory of its previous dynamics. The first attempt for deriving GQME ansatz dated back to Van Hove's work in the late 1950s, where Fourier transform was applied directly to the Hamiltonian that unfortunately led to a rather complex expression involving contour integrals [23, 24]. In 1958 [25] and 1960 [26], Nakajima and Zwanzig independently derived the alternative version of GQME using projection operators, which is now known as Nakajima-Zwanzig GQME (NZ-GQME). The NZ-GQME is derived by employing the projection operators to the Liouville operator to solve the quantum Liouville equation via Laplace transform. With slightly different notation and use of determinants, Montroll derived a version that was equivalent to the NZ-GQME in 1960 [27, 28]. The introduction of the projection operator technique made clear various advantages of working in the density matrix formalism, and a timely summary was published [29]. In 1961, yet another formulation of GQME was produced by Prigonine and Résibois, referred to as the PR-GQME [30, 31]. Along with this development, in 1962 Swenson revisited, assessed, and simplified the rather complex mathematical expression and singular diagonals assumptions used in the original Van Hove's derivation, and yielded the so-called VHS-GQME [32, 33]. In 1963, the density matrix formalism for the VHS-GQME was derived by Peterson and Quay [34]. In 1964, Weinstock [35] also derived a version of GQME that use expansion in term of density, in contrast to other earlier versions that expand in term of interaction (i.e., the system-bath coupling) or Swenson's VHS-GQME that expand in terms of the two-body scattering matrix [33].

All the above formulations of GQME depend on the assumption that there is no correlation of the initial state (which is justified for most problems), nevertheless a more general one that could account for the arbitrary initial condition was also interesting and was presented by Weinstock in 1965 [36]. The PR-, Montroll-, and NZ-GQME are exact and equivalent to each other; although the VHS-GQME is also exact, it is not identical to the other three versions. Discussions on the equivalence between the PR-GQME and VHS-GQME can be found in Ref.[37], and the comparison and equivalence between the PR-, Montroll-, and NZ-GQME were also summarized by Zwanzig [28]. Out of the four early versions of GQMEs, the NZ-GQME formulation received the most attention and popularity owing to its simple form and ease of application, which we will discuss and elaborate further throughout this review. The early 1960s has witnessed the ripening of the fundamental formulation of GQME that prepared the ground for various applications toward solving the irreversibility problem in open quantum systems.

The first modeling for electron transport phenomena using GQME was proposed in 1961 [38]. The numerical demonstration for the aforementioned models under various parametric cases was also achieved not long after [39, 40]. The first review on the GQME to treat irreversible system was given by Chester [41] and soon followed by Prigogine et al. [42]. Soon after these exciting developments, the inconvenience of solving the integro-differential became apparent. To this end, in 1967, Fulinski derived the first time-convolutionless GQME (TCL-GQME) formalism, i.e., without explicit memory term or integral over time [43]. The original proposal was initially criticized [44] but then rigorously proved that the TCL-GQME is also formally exact as the other time-convoluted GQMEs [45]. Motivated by the time-convolutionless form, several investigators furiously developed and formulated different forms of TCL-GQME [46-51]. Following what has been done for NZ-GQME, the perturbation expansion approach was also reported for the TCL-GQME [43, 45, 48-51], and discussed in detail in Ref.[52]. The use of time-ordered operator to simplify the equations became popular [28, 53, 54]. Along with the advent of various numerical methods and increase in computation power, the problem of increasing complexity started to be addressed [55, 56].

B   GQME today

Nowadays GQME has been employed to investigate the quantum dynamics problems in complex molecular systems [57-59], light-harvesting complexes and photosynthetic systems [60-68], molecular junction [69], systems that involve intersystem crossing [70] and conical intersection [71, 72] also in the presence of electromagnetic field [73]. In quantum optics and information, GQME is capable of addressing linear and nonlinear spectroscopies [74-76], quantum computing's qubit dynamics [77], and more.

The two popular forms of GQME are the time-convoluted NZ-GQME and the TCL-GQME. Both of them can be derived using projection operators and are formally exact (i.e., without any assumptions). The two forms of GQME consist of three parts, i.e., the time-local generator of the system, the inhomogeneous term, and the memory kernel in NZ-GQME or the TCL generator in TCL-GQME. These parts narrow down to calculating the evolution operator (or matrix) and the time correlation functions that are used as projection-free inputs with molecular information. The numerical calculation can in turn be performed by various numerically exact or approximate methods for quantum dynamics as will be elaborated in Section Ⅴ.

Recently, in 2003, Shi and Geva proposed a scheme to transform the NZ-GQME memory kernel into the solution of two coupled Volterra equations with projection-free inputs consisting of various system-dependent bath correlation functions that can be evaluated numerically exactly or using approximate dynamical methods [78]. Later in 2006, Zhang, Ka, and Geva further simplified the NZ-GQME memory kernel into the solution of a single Volterra equation that can be solved iteratively [79]. In 2009, Nan, Shi, and Shuai simulated non-perturbative time-convolutionless GQME using path integral approach. In 2015, Pfalzgraff, Kelly, and Markland combined direct Ehrenfest mean-field dynamics with atomistic environment with the SG scheme of NZ-GQME to study aqueous ferrous-ferric electron transfer [80]. In 2015, Kidon, Wilner, and Rabani demonstrated the exact calculation of the TCL-GQME [81]. In 2016$-$2017, Montoya-Castillo and Reichman proposed alternative schemes where the projection operator acts before or after the projected propagator [82, 83]. In 2019, Mulvihill, Geva, and co-workers proposed the modified GQME which is formulated with the more general perturbation type of Hamiltonian, which does not rely on the partition of the system and bath [68, 84-86]. In 2019, Shi and co-workers simulated charge carrier transport in organic molecular crystals using GQME.

FIG. 1 depicts the common partition of the overall open quantum system into the system, the bath, and the system-bath (SB) coupling parts, whose Hamiltonian is given by the system-bath type:

Figure  1.  Illustrations of open quantum system where a system is coupled to the bath through system-bath coupling.

DownLoad: Full-Size Img PowerPoint

where $\hat H_ {\rm{S}}$, $\hat H_ {\rm{B}}$, $\hat H_{ {\rm{BS}}}$ are the system, the bath, and the system-bath coupling Hamiltonians, respectively. In literature, the bath is also named the environment, the reservoir, the noise, and so on depending on the problem being investigated. In the projection operator language, the system and the bath are also called the relevant and the irrelevant parts, respectively. More generally, we can combine the system and the bath Hamiltonians to be the unperturbed zero Hamiltonian, i.e., $\hat H_0$ = $\hat H_ {\rm{S}}$+$\hat H_ {\rm{B}}$, and the SB coupling is treated as the perturbation, $\hat H_ {\rm{I}}$ = $\hat H_{ {\rm{BS}}}$ such that the overall Hamiltonian is written as the perturbation type:

This type of Hamiltonian can also be defined in the following way: the zero Hamiltonian contains the diagonal part and the perturbation is the off-diagonal part, which is used in constructing the recently proposed modified GQME (M-GQME).

Ⅲ.   THEORY OF GQME

A   Quantum Liouville equation (QLE)

The quantum Liouville equation describes the dynamics of a closed quantum system in terms of the total density operator $\hat{\rho}$:

where $\mathcal{L}$ is the quantum Liouville superoperator defined as

In Eq.(3) above, $\hat{\rho}(t)$ represents the overall density operator at time $t$, and $[\hat{A},\hat{B}]$ is the quantum commutator between arbitrary quantum operators $\hat A$ and $\hat{B}$. Expressed in the Hilbert space $\{ |n\rangle\}$, the density matrix element is $\rho_{jk}(t)$ = $\langle j | \hat \rho(t) | k\rangle$ such that

and the elements of the Liouville superoperator is given in a tetradic notation [87, 88]:

Thus, the equation of motion for the density matrix elements is given by

For time-independent Hamiltonian $\hat H$, the evolution of the state vector (wavefunction) is described by the time-evolution operator $\hat{U}(t,t_0)$ = ${\rm{e}}^{-i\hat{H}(t-t_0) / \hbar}$ that is unitary (i.e., $\hat{U}^{-1}$ = $\hat{U}^{\dagger}$ and $\hat{U}\hat{U}^{\dagger}$ = $\hat{U}^{\dagger}U$ = $\hat{\mathit{\boldsymbol{{\rm{{I}}}}}}$, with $\hat{ {\rm{{I}}}}$ being the identity operator). The evolution operator $\hat{U}(t,t_0)$ propagates the ket state vector on the right side as in $|{\Psi}(t) \rangle$ = $\hat U(t,t_0) | {\Psi}(t_0)\rangle$, while the operator $\hat{U}^{\dagger}(t,t_0)$ propagates the bra state vector on the left side as in $\langle {\Psi}(t) |$ = $\langle {\Psi}(t_0)| \hat U^\dagger (t,t_0)$.

The time-dependent density matrix can be expressed in terms of the $\hat{U}(t,t_0)$ and $\hat{U}^\dagger(t,t_0)$:

or, it is convenient to define the time-evolution superoperator $\hat{\mathcal{U}}(t,t_0)$ for more general time-dependent Liouville as below

where $\hat{\mathcal{U}}(t,t_0)$ and its adjoint operator $\hat{\mathcal{U}}^{\dagger}(t,t_0)$ expressed in time-ordered exponential are given by

However, solving the quantum Liouville equation is usually impossible or not necessary for open quantum systems, since we are only interested in the reduced system dynamics. In what follows, we will often revisit QLE in Eq.(3) as the starting point for deriving various forms of the GQME. Being able to work with the system dynamics using the reduced density matrix (RDM) (i.e., $\hat \sigma$ = $\mathrm{Tr}_{ {\rm{B}}}\left\{\hat \rho \right\}$, where $\mathrm{Tr}_{ {\rm{B}}}\{\bullet\}$ is the trace over the bath DOF), and to deal conveniently with mixed states are some of the advantages of working in the density operator formalism of quantum mechanics. To this end, we will introduce the concept of projection operator (strictly speaking superoperator) in the next section.

B   Projection operator

The core of GQME techniques relies on the strategy that conceived the whole system could be partitioned into the relevant part (i.e., the system) and the irrelevant part (i.e., the bath), and that we care about the time-evolution of the relevant part. The projection operator is the key mathematical tool to realize the aforementioned strategy [25, 26, 28]. We can also think of this projection operator as a dimensionality reduction operation that allows us to deal efficiently only with what is essential.

Define a projection operator $\mathcal{P}$ that extracts the relevant information such that it satisfies $\mathcal{P}^2$ = $\mathcal{P}$ (idempotence relation), and its complimentary projection operator $\mathcal{Q}$ = 1$-\mathcal{P}$. The following properties are thus satisfied:

Moreover, the projection operators $\mathcal P$ and $\mathcal Q$ commute with time differentiation and time integral, e.g., $\mathcal{P} {\frac{\mathrm d}{\mathrm dt}}$ = ${\frac{\mathrm d}{\mathrm dt} \mathcal{P}}$ and $\int\mathrm d t \mathcal P$ = $\mathcal P$$\int\mathrm d t$.

For example, the projection operator that extracts the diagonal elements is given by

whose elements are

the projection operator $\mathcal P$ is expressed in the tetradic matrix representation. It can be shown that this projection operator satisfies $\mathcal{P}^2$ = $\mathcal{P}$ and Eq.(13).

Another example of projection operator is what actually used in the NZ-GQME schemes for system-bath Hamiltonian (Eq.(1)):

where the bath equilibrium density operator

with inverse temperature $\beta$ = $1/k_ {\rm{B}}T$ satisfies $\mathrm{Tr}_ {\rm{B}}\{ \hat \rho^{\mathrm{eq}}_ {\rm{B}}\}$ = 1. Eq.(17) satisfies $\mathcal{P}^2$ = $\mathcal{P}$ and Eq.(13) as well.

It is important to note that Eqs.(14) and (17) are not used in deriving NZ-GQME and TCL-GQME forms. In fact, the definition of the projection operator is not required in the derivations of NZ-GQME or TCL-GQME, which hold for any projection operator. However, to evaluate the GQME, one needs to adopt a scheme (Section IV), where the projection operators have to be specified. We will introduce the projection operators that are actually used later, but let us keep it general for now.

C   Derivations of two popular GQME forms

In this section, we take advantage of the projection operator technique to derive GQME in both NZ-GQME and TCL-GQME forms. In FIG. 2 we present the road map of how GQME can be derived, formulated, simplified using different schemes, and evaluated numerically, as will be introduced throughout this tutorial review. Specifically, we will start with the QLE as the starting point for deriving two forms of GQME with the help of the projection operator technique. Furthermore, different schemes and numerical evaluation will also be introduced and discussed in Section Ⅳ and Section Ⅴ, respectively.

Figure  2.  An illustrative road map on how generalized quantum master equation (GQME) can be derived, formulated, simplified using different schemes, and evaluated numerically.

DownLoad: Full-Size Img PowerPoint

1   Nakajima-Zwanzig GQME (NZ-GQME)

NZ-GQME is time-convoluted, time-nonlocal, and the chronological ordering prescription (COP) formulation of GQME [54]. These different names come from the way it can be derived, either using the Fano-Zwanzig-Mori projection operator [28, 89, 90] as will be shown in this section, or the Kubo's ordering prescription as will be mentioned in the comparison between NZ and TCL forms [54, 75, 91]. More recently, the derivation based on dynamical map approach has also been developed [92, 93].

To derive the NZ-GQME, we start by applying the projection operators $\mathcal{P}$ and $\mathcal{Q}$ on the quantum Liouville equation (Eq.(3)) and obtain the following two coupled equations

To solve the inhomogeneous first-order differential equations in Eqs.(19) and (20), we make use of the Laplace transform $\mathbb{L}\{\hat{\rho}(t)\}(s)$ = $\tilde{\mathbb \rho}(s)$ (see Appendix A) and apply projection operators $\mathcal{P}$+$\mathcal{Q}$ such that

To solve for the relevant part $\mathcal{P}\tilde{\rho}(s)$, we take Laplace transform on both sides of the QLE in Eq.(3), apply Eq.(A2), in our case, $s\tilde{\rho}(s)$$-$$\hat{\rho}(0)$ = $-i\mathcal{L}\tilde{\rho}(s)/\hbar$, and obtain

Taking Laplace transform on Eqs.(19) and (20), we have

Solving for Eq.(24) by arranging the term with $\mathcal{Q}\tilde{\rho}(s)$ on the left hand side, we have

where $\mathcal{Q}\tilde{\rho}(s)$ can be obtained in the form like $\tilde{\rho}(s)$ in Eq.(22):

Plugging in the above result for $\mathcal{Q}\tilde{\rho}(s)$ into the Eq.(23) leads to

Finally, perform the inverse Laplace transform on Eq.(27) and note that the product of two functions of the complex variable $s$ now turns into a convolution in time, which allows us to arrive at the NZ-GQME [28]:

where $\mathcal{I}^{\mathrm{NZ}}(t)$ and $\mathcal{K}^{\mathrm{NZ}}(t,\tau)$ are the inhomogeneous term and the memory kernel superoperator, respectively:

Note that in the above derivation, we use the idempotence property of the projection operators (i.e., $\mathcal{P}^2$ = $\mathcal{P}, \mathcal{Q}^2$ = $\mathcal{Q}$). Specifically, when defining $\mathcal{K}^{\mathrm{NZ}}(t,\tau)$, it is conventional throughout the literature to introduce an extra $\mathcal{P}$ for convenience. Also, recall that $\mathcal{L}$ is always a $t$-dependent Liouville superoperator, thus here the memory kernel $\mathcal{K}^{\mathrm{NZ}}(t,\tau)$ is not only a function of $\tau$ but also a function of $t$. In the form given by NZ-GQME in Eq.(28), the reduced density matrix (RDM) defined as

which reports the system dynamics under the influence of the bath.

We now discuss the composition and characteristics of Eq.(28). The first term is the time-local generator and describes the time evolution of the reduced density matrix as if the system is isolated from the bath, i.e., if $\mathcal{QLP = PLQ}$ = 0, Eq.(28) reduces to QLE for a closed system ${\frac{\mathrm d}{\mathrm dt} \hat{\sigma}(t)}$ = $- {\frac{i}{\hbar}\mathcal{P}\mathcal{L} \hat{\sigma}(t)}$. The second term is the inhomogeneous term that accounts for the effect of initial state preparation, $\hat \rho(0)$. This term vanishes if $\hat \rho(0)$ is invariant to projection (which is often the case), i.e., if $\mathcal P\hat{\rho}(0)$ = $\hat \rho(0)$, then $\mathcal Q \hat \rho(0)$ = 0 and so $\mathcal I^{\mathrm{NZ}}(t)$ = 0. The last term is the time convolution of the memory kernel $\mathcal{K}^{\mathrm{NZ}}(t,\tau)$ and the RDM $\hat \sigma(t)$ at earlier time, which reflects the non-Markovian effect where the RDM depends on its previous dynamics under the influence of the bath.

The memory kernel is at the heart of evaluating the GQME and is a central subject of this review that will often be revisited, reformulated, and discussed. Before proceeding with such an elaborate discussion of several schemes to evaluate the memory kernel, we will first introduce the time-convolutionless GQME in next section.

2   Time-convolutionless GQME (TCL-GQME)

The time-convolutionless GQME is also called time-local [1], memoryless, or partial ordering prescriptions (POP) formulation for the GQME [54]. In this section, we present that one can derive formally exact time-convolutionless form of GQME that is equivalent to the NZ-GQME in Eq.(28). The derivation for obtaining the TCL-GQME can be done through several approaches [43, 45, 46, 48-51].

Here we show the derivation of TCL-GQME using the partial ordering prescription. We start from the QLE in Eq.(3) and apply the projection operator $\mathcal{P}$ and $\mathcal{Q}$ to obtain two coupled linear differential equations in Eq.(19) and Eq.(20). The solution of Eq.(20), a first-order inhomogeneous ordinary differential equation, is given by

This time we do not plug this integrated result back to Eq.(19). Instead, we plug in the formal solution of the QLE (see Eq.(3)) as below [50, 90]

into Eq.(32), and obtain

To solve the above equation, we now define the damping operator $\Sigma(t)$ after the damping theory developed by Hashitsumae et al. [49, 50] as follows

We also define the forward propagator, $\mathcal{G}_{\mathcal Q}(t,t_0)$, for general time-dependent Liouville in terms of positive time-ordered expansion (Eq.(11))

and the backward propagator, $G(t, t_0)$, in terms of negative time-ordered expansion (Eq.(12))

Now, the damping operator can be simplified as

Next, we rewrite Eq.(34) using the damping operator and the forward/backward propagators and obtain

Collecting the term containing $\mathcal{Q}\hat{\rho} (t)$ on the left

and expressing $\mathcal{Q} \hat{\rho}(t)$ in terms of $\mathcal{P} \hat{\rho}(t)$ and $\mathcal{Q} \hat{\rho}(0)$, we have

Plugging in the above explicit form of $\mathcal{Q} \hat{\rho}(t)$ back to Eq.(19) yields the TCL-GQME formula as follows:

where the TCL generator or time-local generator $\mathcal{K}^{\mathrm{TCL}}(t)$ and the inhomogeneous term $\mathcal{I}^{\mathrm{TCL}}(t)$ are defined as

Here, we introduce an extra $\mathcal{P}$ to the end of $\mathcal{K}^{\mathrm{TCL}}(t)$ using the idempotence property.

3   Relationship between NZ-GQME and TCL-GQME

We discuss the similarities and differences between the memory kernel $\mathcal{K}^{\mathrm{NZ}} (t,\tau)$ in NZ-GQME and the time-local generator $\mathcal{K}^{\mathrm{TCL}} (t)$ in TCL-GQME, as well as when the inhomogeneous term vanishes.

First, we stress again that the two forms of GQME (as shown in Eq.(28) and Eq.(42)) are both formally exact and equivalent. However, the main difference between the two stems from the so-called memory effect that either enters into the GQME explicitly in the form of time-convolution or without such explicit dependence as in the time-convolutionless form. Efforts have been made to establish the connection between the two forms through approximate solutions based on expansions [1, 43, 45], numerical methods [92], and recently through generalized dephasing models [16].

Second, the memory kernel and the inhomogeneous term in NZ-GQME (Eq.(28)) can also be expressed in terms of the forward propagator $\mathcal{G}_{\mathcal Q}(t,\tau)$ defined in Eq.(36), which is used in the derivation of TCL-GQME:

Here, the evolution for the forward propagator $\mathcal{G}_{\mathcal Q}(t,0)$ is given by

with the initial condition $\mathcal{G}_{\mathcal Q}(0,0)$ = $\hat{ {\rm{I}}}$, but it is not straightforward to solve numerically exactly since it involves the complementary projection operator $\mathcal Q$.

Third, the inhomogeneous term is non-convoluted in both GQME forms and will vanish when the initial state of the system and the bath is uncorrelated. For example, if we adopt the projection operator defined in Eq.(17) and assume the initial density matrix can be factorized as follows

then

So, the inhomogeneous terms of NZ-GQME in Eq.(29) and of TCL-GQME in Eq.(44) vanish. Under this circumstance the time-evolution of the reduced density matrix is simplified to

Although the choice of the initial state as in Eq.(48) is widely applicable, there are cases where the inhomogeneous term does not vanish, and we will briefly discuss at the end of Sec. IV.

Lastly, it is noted that the derivations presented here for both NZ-GQME and TCL-GQME are not unique. There exist various approaches to derive the same results, for instance, through the path integral formalism [94], Mori equation of motion [82, 83], and diagrammatic expansion [92, 95].

Ⅳ.   MEMORY KERNEL AND INHOMOGENEOUS TERM

In this section, we focus on describing several schemes to solve the key terms in NZ-GQME---the memory kernel $\mathcal{K}^{\mathrm{NZ}}(t,\tau)$ and the inhomogeneous term $\mathcal{I}^{\mathrm{NZ}}(t)$. As we indicated earlier, directly solving any of the two superoperators either analytically or numerically is a challenge, notably because of the explicit participation of the complementary projection operator $\mathcal{Q}$ in the propagator ${\rm{e}}^{-i\mathcal{Q}\mathcal{L}\tau /\hbar}$. Therefore, it is desirable to formulate these superoperators in terms of projection-free inputs. To this end, the projection operator like Eq.(17) and the Hamiltonian type like the system-bath type (Eq.(1)) and the perturbation type (Eq.(2)) have to be specified first. We will begin the discussion with schemes for NZ-GQME memory kernel expressed in system-bath Hamiltonian and then followed by the alternative NZ-GQME called the modified GQME (M-GQME) that is formulated in perturbation Hamiltonian.

A   Memory kernel for system-bath Hamiltonian

Now, we adopt the system-bath Hamiltonian $\hat H$ = $\hat H_ {\rm{S}}$+$\hat H_ {\rm{B}}$+$\hat H_{ {\rm{BS}}}$ (Eq.(1)) and describe several schemes that provide formally exact memory kernel in NZ-GQME. The overall Liouville superoperator consists of three terms including the system, the bath, and the system-bath coupling such that

where $\mathcal L_ {\rm{S}}(\bullet)$ = $\left[\hat H_ {\rm{S}},\bullet\right]$, $\mathcal L_ {\rm{B}}(\bullet)$ = $\left[\hat H_ {\rm{B}},\bullet\right]$, and $\mathcal L_{ {\rm{BS}}}(\bullet)$ = $\left[\hat H_{ {\rm{BS}}},\bullet\right]$. We use the projection operator $\mathcal P(\bullet)$ = $\hat \rho^{\mathrm{eq}}_ {\rm{B}} \otimes \mathrm{Tr}_ {\rm{B}}(\bullet)$ (Eq.(17)) and assume that the initial state of the overall density matrix can be expressed in the direct product of equilibrium bath part and the system part given by $\hat{\rho}(0)$ = $\hat \rho^{\mathrm{eq}}_ {\rm{B}} \otimes \hat{\sigma}(0)$ (Eq.(48)). Moreover, we assume the following system operator vanishes

which can be made without loss of generality, because if $\mathrm{Tr}_ {\rm{B}}\left[\hat \rho_ {\rm{B}}^{\mathrm{eq}} \hat H_{ {\rm{BS}}} \right]$$\neq$0, one can always redefine the system and the system-bath coupling such that $\hat H_ {\rm{S}}'$ = $\hat H_ {\rm{S}}$+$\mathrm{Tr}_ {\rm{B}}\left[\hat \rho_ {\rm{B}}^{\mathrm{eq}} \hat H_{ {\rm{BS}}} \right]$ and $\hat H_{ {\rm{BS}}}'$ = $\hat H_{ {\rm{BS}}} - \mathrm{Tr}_ {\rm{B}}\left[\hat \rho_ {\rm{B}}^{\mathrm{eq}} \hat H_{ {\rm{BS}}} \right]$ which would guarantee that $\mathrm{Tr}_ {\rm{B}}\left[\hat \rho_ {\rm{B}}^{\mathrm{eq}} \hat H_{ {\rm{BS}}}' \right]$ = 0.

Now we apply the projection operator in Eq.(17) to the NZ-GQME (Eq.(28)) and perform a partial trace over the bath DOF on each term. The LHS of Eq.(28) is given by

The first term on RHS of Eq.(28) is the isolated system dynamics:

where we used that $[\hat H_{ {\rm{B}}}, \hat \rho_ {\rm{B}}^{\mathrm{eq}}]$ = 0 and $\mathrm{Tr}_ {\rm{B}} \left\{ [\hat H_{ {\rm{BS}}}, \hat \rho_ {\rm{B}}^{\mathrm{eq}}]\right\}$ = 0.

The second term on RHS of Eq.(28) is the memory kernel term that can be simplified to be (see proof in Appendix B)

The third term on RHS of Eq.(28) is the inhomogeneous term, which vanishes using the current choice of projection operator and initial state (see Eqs. (49) and (50)).

Combining the above terms leads to the NZ-GQME for system-bath Hamiltonian, which is the equation of motion for the RDM of an open quantum system:

where the formally exact memory kernel (a superoperator of system DOF) is given by

Here, $- {\frac{i}{\hbar}} \mathcal L_ {\rm{S}}\hat\sigma(t)$ = $- \frac{i}{\hbar}\left[\hat H_ {\rm{S}}, \hat \sigma(t) \right]$ corresponds to the bath-free and $- \int_0^t\mathrm d\tau \mathcal K(\tau) \hat\sigma$($t$$-$$\tau)$ corresponds to the bath-induced contributions to the system's reduced dynamics.

Note that the NZ memory kernel here is expressed as a simplified notation but still formally exact: $\mathcal{K}^{\mathrm{NZ}}(t,\tau)$$\longrightarrow$$\mathcal{K}(\tau)$ with only the relative lag time $\tau$ as argument shown explicitly. Moreover, the memory kernel in Eq.(58) can be written in the following equivalent forms [86]:

To facilitate discussions of the evaluation schemes, we denote the memory kernel representations in Eqs.(58)–(64) as the version 1$-$7 of the NZ-GQME memory kernel for system-bath Hamiltonian $\mathcal K(t)$. These expressions of memory kernel are equivalent and will generate the same reduced dynamics if evaluated with exact methods but the result may differ when approximate methods are used. In the following sections, we will introduce the Shi-Geva, Zhang-Ka-Geva, and Montoya-Castillo-Reichman schemes for calculating the NZ memory kernel for the system-bath Hamiltonian.

B   Shi-Geva (SG) scheme

The Shi-Geva (SG) scheme [78] provides a pathway to circumvent the dependency on operator $\mathcal{Q}$ in the propagator of the memory kernel of NZ-GQME. To this end, the Shi-Geva scheme starts from the memory kernel version 7 in Eq.(64), which relies on the following relations

whose proof is included in Appendix A of Ref.[78] and then apply Dyson decomposition (Eq.(A11) with $\mathcal{A}$ = $\mathcal{L}$ and $\mathcal{B}$ = $\mathcal L$$-$$\mathcal L_{ {\rm{BS}}} \mathcal P$, i.e.

Applying the above Dyson decomposition to Eq.(64), one can arrive at a Volterra equation of the second kind [78, 96]:

where the {auxiliary kernels} $\mathcal{K}_{1}(t)$ and $\mathcal{K}_{2}(t)$ are defined as

Here, $\mathcal{K}_{1}(t)$ is projection-free but $\mathcal{K}_{2}(t)$ is projection-dependent, which can be reformulated to projection-free by applying the Dyson decomposition again in Eq.(67) to $\mathcal{K}_2(t)$:

where the projection-free auxiliary kernel $\mathcal{K}_{3} (t)$ is given by

Thereby, the memory kernel $\mathcal{K}(t)$ can be calculated by obtaining the projection-free auxiliary kernels $\mathcal{K}_1(t)$ and $\mathcal{K}_3(t)$ as inputs, then solving the Volterra equation Eq.(71) for $\mathcal{K}_2(t)$, and finally solving the Volterra equation Eq.(68) for $\mathcal{K}(t)$ [78, 79, 84]. The projection-free inputs can be formulated in terms of time correlation functions to be obtained from numerical calculations using various dynamics methods, which will be discussed later.

C   Zhang-Ka-Geva (ZKG) scheme

The Zhang-Ka-Geva (ZKG) scheme is a further extension of the Shi-Geva scheme. In the ZKG scheme, the Dyson decomposition in Eq.(A11) is performed directly to the memory kernel of version 1 in Eq.(58) to yield a single Volterra equation as below

where the projection-free $\Phi(t)$ is defined as

It is noted that the term $\dot{\Phi}(t)$ = ${\frac{\mathrm d}{\mathrm dt}\Phi(t)}$ has similar form to the $\mathcal{K}_1(t)$ (cf. Eq.(69)) but with the last Liouvillian given by $\mathcal{L}$ instead of $\mathcal{L}_{ {\rm{BS}}}$,

This scheme is an improvement over the previously introduced SG scheme in that only one Volterra equation is to be evaluated taking in projection free-inputs $\Phi(t)$ and $\dot{\Phi}(t)$. The projection-free inputs can be formulated in terms of time correlation functions, which will be discussed later.

D   Montoya-Castillo-Reichman (MCR) scheme

The Montoya-Castillo-Reichman scheme [82, 83] can be derived in two ways starting from version 2 or version 5 of the memory kernel expressions. Depending on whether the complementary projection operator $\mathcal{Q}$ is applied before (Eq.(59)) or after (Eq.(62)) the projected propagator ${\rm{e}}^{-i\mathcal{Q}\mathcal{L}t/\hbar}$, the corresponding backward (MCR-b) and forward (MCR-f) memory kernels are obtained. Applying the Dyson identity in Eq.(A11) to Eq.(59) where the projection is applied first, we arrive at the MCR-b memory kernel:

and applying the Dyson identity in Eq.(146) with $\mathcal A$ = $\mathcal{LQ}$ and $\mathcal{B}$ = $\mathcal{L}$ to Eq.(62) where the projection is applied second, we arrive at the MCR-f memory kernel:

Here, $\mathcal{K}_1(t)$, $\mathcal{K}_3(t)$, and $\Phi(t)$ are the same as what are introduced earlier in Eq.(69), Eq.(72), and Eq.(74). Within the above MCR-b/f scheme, the projection-free auxiliary kernels are used as inputs to solve the single Volterra equation, and the MCR-b scheme where projection acts first, was shown to give more accurate result than the MCR-f version [97].

E   Modified GQME (M-GQME)

The modified GQME (M-GQME) scheme was proposed by Mulvihill, Geva, and co-workers in 2019 [68, 84-86], which provides an alternative way to describe the quantum dynamics for system with overall Hamiltonian expressed in a more general perturbation type (Eq.(2)). In quantum dynamics for condensed phase, a typical case is the study of nonadiabatic dynamics, where the electronic and the nuclear DOF are interrelated and defining a unique single bath is sometimes not possible. To this end, consider an $F$-electronic-state Hamiltonian

where $\hat{H}_{j}$ = $\hat{H}_{jj}$ = $\hat{{\bf{P}}}^2/2$+$\hat{V}_{j}(\hat{{\bf{R}}})$ is the nuclear Hamiltonian of the electronic state $|j\rangle$, $\hat{{\bf{P}}}$ = $(\hat{P}_1, \ldots ,\hat{P}_N)$ and $\hat{{\bf{R}}}$ = $(\hat{R}_1, \ldots ,\hat{R}_N)$ are the mass-weighted momenta and coordinates for the $N$ nuclear DOF, and $\hat{V}_{j}(\hat{{\bf{R}}})$ is the potential energy surface of the electronic state $|j\rangle$. In addition, $\hat{V}_{jk} (\hat{{\bf{R}}})$ is the coupling between state $|j\rangle$ and $|k\rangle$, which becomes independent of the nuclear configuration in the Condon approximation, i.e., $\hat V_{jk}$ = $V_{jk}$ = Const., whereas when the coupling depends on nuclear configuration, we have the general non-Condon case.

The initial state of the overall system is assumed to be in the factorized form

where $\hat \rho_ {\rm{N}}(t)$ = $\mathrm{Tr}_ {\rm{e}}\{ \hat \rho(t) \}$ and $\hat \sigma (t)$ = $\mathrm{Tr}_ {\rm{N}}\{ \hat \rho(t)\}$ are the nuclear and the electronic reduced density matrices, respectively, which are obtained by taking the partial trace over the electronic and the nuclear Hilbert spaces, $\mathrm{Tr}_ {\rm{e}}\{ \bullet \}$ and $\mathrm{Tr}_ {\rm{N}}\{ \bullet \}$, respectively.

In the M-GQME scheme, the projection operator $\mathcal{P}$ is defined as

Here, $\hat{\rho}_{ {\rm{N}}}^{\mathrm{ref}}$ is an arbitrary reference nuclear density operator, which satisfies $\mathrm{Tr}_{ {\rm{N}}} \{\hat{\rho}_{ {\rm{N}}}^{\mathrm{ref}} \}$ = 1. This also means that the initial state is not restricted to the thermal equilibrated bath (cf. Eq.(18)). Here, we choose

which leads to a vanishing inhomogeneous term, since $\mathcal{Q} \hat \rho_ {\rm{N}}(0)$ = 0 such that the inhomogeneous term

It is noteworthy that by choosing a specific reference, the resulting memory kernel will be explicitly dependent on the given initial condition, and altering the initial condition will lead to different dynamics.

In the non-Condon case, applying the projection operator (Eq.(80)) and Dyson decomposition (Eq.(A11)) to NZ-GQME (Eq.(28)), we have M-GQME as follows

where $\langle\mathcal{L}\rangle^{0}_{ {\rm{N}}}$ is an electronic superoperator as a result of taking the overall Liouvillian averaged over the initial state of the nuclear DOF, i.e.

and the memory kernel superoperator is given by

In the Condon case, the electronic couplings are constants, i.e., $\hat V_{jk}$ = $V_{jk}$ = Const., where we denote the diagonal part as the zero Hamiltonian and the off-diagonal part as the interaction Hamiltonian as follows

The Liouville operators in the Condon case are given by $\mathcal{L}_{\mathrm{zero}}(\bullet)$ = $[\hat{H}_{\mathrm{zero}},\bullet]$ and $\mathcal{L}_{\mathrm{int}}(\bullet)$ = $[\hat{H}_{\mathrm{int}},\bullet]$, where the $\mathcal{L}_{\mathrm{int}}$ is purely electronic superoperator since the off-diagonal elements are constant. In this case, the memory kernel becomes

Applying Dyson decomposition in Eq.(A11) to memory kernel for the non-Condon case in Eq.(85) and to memory kernel in the Condon case in Eq.(88), we obtain their corresponding Volterra equations [84, 97]:

The projection-free inputs are then defined by [97]

The projection free-inputs can be formulated to time-correlation functions, then they can be evaluated by various numerical methods. End-to-end formulation, technical discussions, and applications of the M-GQME can be found in Refs.[68, 85, 86, 97].

F   Inhomogeneous term

It is commonly assumed that in the initial state, the relevant and irrelevant parts are uncorrelated so that the overall density matrix can be factorized, such as Eq.(48) for system-bath NZ-GQME and Eq.(79) for M-GQME. As a consequence, the inhomogeneous term in the GQMEs vanishes. Situations for which the bath starts in thermal equilibrium while the system is created by some excitation exemplify this case.

However, the uncorrelated initial state is not always the case. For example, if the initial state is $\hat \rho(0)$ = ${\rm{e}}^{-\beta(\hat H_ {\rm{S}}+\hat H_ {\rm{B}}+\hat H_{ {\rm{BS}}})}/\mathrm{Tr}_ {\rm{B}}\{ {\rm{e}}^{-\beta(\hat H_ {\rm{S}}+\hat H_ {\rm{B}} + \hat H_{ {\rm{BS}}})}\}$, then it cannot be written in separated form. If the system and the bath are initially correlated, the inhomogeneous term persists, which makes the situation rather complicated. The first concern is that not every arbitrary inhomogeneous term could lead to physical dynamics of the whole system [98]. From more practical aspect, one is also concerned with how long the correlation lasts and how strong the correlation is. For a relatively weak and short correlation, one could address this issue through the perturbation expansion for the memory kernel and TCL generator. Illustrating for the perturbation expansion of the inhomogeneous term (Eq.(46)) up to the second order, we have [1, 99, 100]

A more general framework to account for the various initial conditions under the weak-coupling limit can be found in the work of de Vega and Alonso [101]. A general description applied to strong coupling regime has also been developed [56]. Furthermore, comprehensive discussions on the influence of correlated initial condition toward the validity of the consequent dynamics [4, 14, 101, 102], how to account for the bath effect [103], optimizing the initial state formulation [104], and preserving the exactness of the resulting inhomogeneous term, are actively being discussed. This issue is still under ongoing investigations and is interesting notably because in many instances the complete positivity relationship might be violated and leads to unphysical dynamics [4, 98]. This raises a question regarding the validity of the GQME dynamics when no well-defined boundary exists [105, 106]. Extension of this question includes dealing with other master equations that focus on decoherence, such as the Hu-Paz-Zhang master equation [107] in relation to investigating the nature of the inhomogeneous term [108].

Similar to the memory kernel, the inhomogeneous term could also be obtained from numerical methods by calculating the reduced density matrix, evolution operator, or the time correlation functions. In the next section, we summarize the numerical quantum dynamical methods applicable to calculating the above-mentioned quantities of our interest that enter the GQMEs.

G   Other approaches

The schemes presented above are merely the tip of the iceberg of how one can further develop and combine GQME techniques. A comprehensive summary that explores various road maps to solve for the memory kernel can be found in Ref.[97]. Beyond the schemes that are pertinent to the chemical physics community, there are many more schemes applicable for evaluating the memory kernel that is interested by a wide community such as quantum information [4, 98], quantum thermodynamics [6-9], and more. For example, there is a scheme that connects the NZ-GQME and TCL-GQME to rewrite the memory kernel in terms of the reduced propagator and effective Liouvillian [95]. Another scheme resembling the spirit of the M-GQME for electron-phonon has also been developed taking advantage of Fock space formulation [109]. There is a recent development using a completely different approach by applying diagrammatic expansion on the time-convolutionless master equation [110].

We would like to comment on the widely used perturbation expansion of the memory kernel, which is valid when the system-bath coupling is weak. Starting with Hamiltonian of the perturbation type in Eq.(2) and the projection operator in Eq.(17), we note for $n$ = 0, 1, 2, $\ldots$ that

Then the perturbation expansion in terms of the interaction Hamiltonian $\hat{H}_I$ up to the second oder expansion gives [1]

For the TCL generator as in Eq.(43), assuming that the inverse operator $\left(1-\Sigma(t)\right)^{-1}$ exists, we can expand the term $\left(1-\Sigma(t)\right)^{-1}$ = $\sum\limits_{n = 0}^{\infty} \left(\Sigma(t)\right)^n$ such that [1, 43, 45, 49-52]

Therefore, we obtain TCL generator up to the second order that is given by

with $\Sigma_1 (t)$ = $\int_0^t \mathrm d \tau \mathcal{Q}\mathcal{L}(\tau)\mathcal{P}$, which gives the same expression as in the second-order NZ-GQME memory kernel. A series of studies by many workers have shown that the second-order TCL-GQME scheme gives better numerical performance compared to the second-order NZ-GQME [1, 54, 81, 100, 111]. It is worth noting that efforts have been made to address how the TCL-GQME can be evaluated non-perturbatively [94, 112, 113]. Throughout the literature, it is common to refer to results obtained from the truncation procedure as a {partial resummation} of the perturbation expansion. It is generally known that expression obtained through the perturbation theory might not always be convergent [113-115]. To this end, different resummation techniques are typically employed to obtain better convergent results. Two of the typically employed resummation techniques are the Landau-Zener resummation [116] and the P$\acute{\mathrm{a}}$de resummation [117, 118]. The accuracy of the resummation technique could be a concern [113, 115]. More recently new schemes have been developed to perform better resummation of memory kernels [114].

Ⅴ.   PROJECTION-FREE INPUTS: SYSTEM-DEPENDENT BATH CORRELATION FUNCTIONS

In this section, we summarize the projection-free inputs for calculating the memory kernel in the NZ-GQME and M-GQME. These projection-free inputs carry the dynamical information of realistic molecular or condensed-phase system, which may be evaluated with numerically exact or approximate quantum dynamical methods. If these inputs are evaluated exactly, the resulting dynamics from NZ-GQME or M-GQME should be equal, whereas if the inputs are obtained approximately, the GQME reduced dynamics may differ.

For an $F$-electronic-state system, the memory kernel superoperator $\mathcal K(t)$ is an $F^2$$\times$$F^2$ matrix defined in Liouville space with basis $\{|j,k\rangle\rangle$ = $|j\rangle \langle k|; \;j,k = 1,\ldots,F\}$, and so are the other projection-free inputs. We define the elements of the memory kernel as

Here, $\mathrm{Tr}\{\bullet\}$ = $\mathrm{Tr}_ {\rm{S}}\{\mathrm{Tr}_ {\rm{B}}(\bullet)\}$ = $\mathrm{Tr}_ {\rm{e}}\{\mathrm{Tr}_ {\rm{N}}(\bullet)\}$ denotes the trace over all DOF.

A   Projection-free inputs for NZ-GQME

Projection-free inputs for the memory kernel of NZ-GQME for the system-bath Hamltonian in SG, ZKG, MCR schemes include $\mathcal{K}_{1}(t)$, $\mathcal{K}_{3}(t)$, $\Phi(t)$, and $\dot{\Phi}(t)$ that are defined in Eqs.(69), (72), (74), (75) and listed again as follows

Assume the system-bath coupling can be factorized as below

where $\hat \Lambda$ is a bath operator (so Eq.(53) implies $\langle \hat \Lambda \rangle_{\mathrm{eq}}$ = $\mathrm{Tr}_ {\rm{B}}\{ \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \Lambda \}$ = 0) and $\hat S$ = $S(\hat x)$ is a system operator that satisfies $S(\hat x) |x\rangle$ = $S_x | x\rangle$. Let us first look at $\mathcal{K}_1(t)$ that will be required in SG, MCR-b/f schemes. The elements of $\mathcal{K}_1(t)$ are given by

where we define the system-dependent bath correlation function (SDBCF) [78] with $\langle \bullet \rangle$ = $\mathrm{Tr}\left\{ \hat \rho_ {\rm{B}}^{\mathrm{eq}} \bullet \right\}$ and for any bath operator $\hat B$, we define

so, SDBCF in general satisfies

As will see later, all that one needs to provide to GQME as projection-free inputs can be formulated into SDBCFs, with nuclear operators $\hat A,\hat B$ including $\hat \Lambda$, $\hat I_B$ (bath identity operator).

Second, we look at $\mathcal{K}_3(t)$ that will be required in the SG and MCR-f schemes. Similarly, we can expand the quantum commutator due to $\mathcal L_{ {\rm{BS}}}$ and express each term using SDBCF.

Third, we look at $\Phi(t)$ that will be required in the ZKG and MCR-b schemes.

Fourth, we look at $\dot{\Phi}(t)$ that will be required in the ZKG scheme, which should be evaluated explicitly rather than taking numerical derivative of $\Phi(t)$. This will be expanded into many terms due to the existence of the overall Hamiltonian

B   Projection-free inputs for M-GQME

Now, we discuss the projection-free inputs for the memory kernel of M-GQME for the perturbation type of Hamltonian. They include $\mathcal F(t)$ and $\dot{\mathcal F}(t)$ for non-Condon case as well as $\mathcal F_1(t)$ and $\mathcal F_2(t)$ for Condon case that are defined in Eqs.(92–95):

First, we look at the $\mathcal{F}(t)$ for the non-Condon case.

Second, we look at the $\dot{\mathcal{F}}(t)$ for the non-Condon case.

Next, we look at the $\mathcal{F}_{1}(t)$ and $\mathcal{F}_{2}(t)$ for the Condon case.

In fact, the projection-free SDBCFs required for M-GQME can be formulated in the following form:

where the "time 0" nuclear operator $\hat \Omega$ is $\hat \rho_ {\rm{N}}(0)$ in $\mathcal F(t)$ and $\mathcal{F}_2(t)$, $\left(\hat H_u \hat \rho_ {\rm{N}}(0) - \hat \rho_ {\rm{N}}(0) \hat H_{v}\right)$ in $\mathcal{F}_1(t)$, as well as $\left(\hat H_u \hat \rho_ {\rm{N}}(0) - \hat \rho_ {\rm{N}}(0) \hat H_{v}\right)$, $\hat V_{au} \hat \rho_ {\rm{N}}(0)$, or $\hat \rho_ {\rm{N}}(0)\hat V_{vb}$ for $\dot{\mathcal F}(t)$ [68]. The "time $t$" nuclear operator $\hat B$ in $\mathcal F(t)$ and $\dot{\mathcal F}(t)$ is given by

where only the first option applies to $\mathcal{F}_1(t)$ and $\mathcal{F}_2(t)$.

C   Quantum dynamical methods

The projection-free inputs for obtaining the memory kernel boil down to the calculation of SDBCFs, which are generically expressed by Eq.(106) for system-bath formulation of NZ-GQME and by Eq.(114) for the M-GQME. The SDBCFs contain all the molecular information or the feature of the specific physical system under study that GQME needs as the inputs, therefore it is important to have accurate SDBCFs so as to obtain the accurate reduced system dynamics such as RDM $\hat \sigma(t)$. The quantum dynamical methods are referred to as the numerical techniques that either provide numerically exact evaluation of the SDBCFs or offer approximate estimation of them. Although numerically exactness is desirable since its combination with NZ-GQME or M-GQME will lead to exact reduced dynamics of the system, it is not always possible to perform numerically exact calculations for most of chemical systems with realistic all-atom Hamiltonians. Therefore, the applications of GQME have been divided into model Hamiltonians or all-atom realistic Hamiltonians, where the all-atom systems are mostly treatable by the approximate dynamical approaches.

To this end, we first introduce numerically exact approaches that could give rise to exact SDBCFs, which will be used as projection-free inputs for GQME. The numerically exact approaches include the path-integral formalism, wave-function formalism, and dynamical map formalism. First, within the path-integral formalism, there are methods such as the quasi-adiabatic propagator path integral (QuAPI) [119-121], the quantum-classical path integral (QCPI) [121], and the iterative blip-summed path integral [122]. One can also directly obtain the memory kernel from real-time path integral Monte Carlo [123, 124], though it is generally only possible for some special cases. Another prominent technique of the path-integral formalism is the hierarchical equations of motion (HEOM) [125-127], and the recently proposed generalized HEOM to account for arbitrary correlation function was also recently developed [128]. Second, under the hood of the wavefunction formalism, there are methods such as the time evolving density matrix using orthogonal polynomial algorithm (TEDOPA) [129], wave-packet dynamics [130], time-dependent density matrix renormalization group (t-DMRG) [131, 132], hierarchical stochastic Schrödinger equations (HSSE) [133], multiconfiguration time-dependent Hartree (MCTDH) [134-138] and the multilayer MCTDH (ML-MCTDH) [139, 140]. Third, the approach based on the dynamical map is typically realized through the recently popular technique, the tensor-train (TT) decomposition [141], and other methods developed based on it [142-144]. In addition, Cerrillo and Cao proposed the transfer tensor method (TTM) that generically fits the dynamical map using the initial reduced system dynamics as input, which circumvents the calculation of GQME [145, 146].

However, the numerically exact dynamical methods have restrictions on the Hamiltonian forms and scalability with increasing DOF, which make it difficult if not impossible to apply to realistic all-atom systems. Approximate numerical methods, on the other hand, provide more tractable ways to obtain the projection-free inputs for GQME calculation, such as mixed quantum-classical approaches including the mean-field (MF) Ehrenfest dynamics [68, 79, 80, 85, 147, 148], the mixed quantum-classical Liouville (MQCL) [149], as well as semiclassical approaches including the linearized semiclassical (LSC) dynamics and symmetrical quasiclassical (SQC) dynamics based on the Meyer-Miller mapping Hamiltonian [85]. Combining the assorted approximate methods with the GQME formalism is a strategy commonly employed, which has been shown to improve accuracy and numerical efficiency in contrast to directly obtaining the system dynamics from the approximate methods without GQME [68, 148, 150, 151].

Here, we take the mean-field Ehrenfest dynamics as an example to show how the SDBCF could be obtained. The Ehrenfest dynamics treats the electronic (system) DOF quantum-mechanically and the nuclear (bath) DOF classically, and there is a feedback loop between the two in the mean-field way [152, 153]. Traditionally, the quantum part in Ehrenfest method is described by pure-state wavefunctions,

and its equation of motion is the time-dependent Schrödinger equation

where the electronic Hamiltonian is the potential energy part and electronic couplings of the overall Hamiltonian $\hat{h}({\bf R})$ = $\sum\limits_{j,k}^{F}V_{jk}\left({{\bf{R}}}\right)|j\rangle\langle k|$, as in the overall Hamiltonian $\hat H$ = ${{\bf P}}^2/2$+$\hat h({\bf R})$, and thus the time evolution of the wavefunction coefficients $c_j(t)$ is given by

The nuclear DOF are propagated classically on the instantaneous mean-field PES that are defined as

which depends on the wavefunction coefficients $c_j(t)$ of the system.

Alternatively, for representing mixed states, we express the quantum system by the reduced density matrix whose equation of motion is given by quantum Liouville equation

Solving the above equation of motion leads to electronic propagation:

where the time evolution operator $\hat U(t+\delta t,t)$ in the diabatic basis can be expressed in terms of adiabatic energies through the transformation matrix ${\bf T}$ that diagonalizes $\hat h({\bf R}_t)$ such that ${\bf T}^\dagger \hat h({\bf R}_t) {\bf T}$ = $\mathrm{diag}(E_1,\ldots, E_F)$:

The mean-field potential is given by

which dictates the nuclear equations of motion

The initial state of the overall system is assumed to be given by the overall density operator $\hat \rho(0)$:

The initial sampling for the nuclear DOF is given by the Wigner transformed nuclear density as follows

whereas there is no initial sampling for the electronic DOF since the initial values of the electronic reduced density matrix $\sigma_{jk}(0)$ or the wavefunction coefficients $c_j(0)$ are known. In the pure-state case, $\sigma_{jk}$ = $c_j c_k^*$ and both versions of Ehrenfest dynamics (wavefunction and density matrix) are equivalent.

In the mean-field dynamics, the SDBCF defined in Eq.(106) or Eq.(114) can be evaluated in the following forms:

where the system starts out with $\hat \sigma(0)$ = $| u\rangle \langle v|$ with only $\sigma_{uv}(0)$ = 1.

Before we close this section, it is noted that in the weak coupling limit, the SDBCF becomes free-bath correlation function (FBCF) given by

which is only a property of the bath and will be used as input for GQME.

D   Model Hamiltonians

GQME has been used to learn about many open quantum systems described mostly by harmonic models such as the spin-boson model [68, 78, 79, 82-85], Caldeira-Leggett model [154-157], single-impurity Anderson model [123, 124], Frenkel excitonic model for light-harvesting Fenna-Matthwes-Olson (FMO) complex [61, 68, 158-160], one-dimensional Holstein model of organic molecule crystals [161]. Beyond harmonic models, Pfalzgraff, Kelly, and Markland combined direct Ehrenfest mean-field dynamics with atomistic environment with the SG scheme of NZ-GQME to study aqueous ferrous-ferric electron transfer [80].

For condensed phase systems, mapping all-atom anharmonic Hamiltonian to model harmonic Hamiltonians is typically performed for the GQME calculation [61, 65, 158-160, 162-166]. In particular, spin-boson model among others has been serving as the test ground for benchmarking purposes. Spin-boson model is a two-level system with the electronic Hilbert space defined in terms of ground state denoted as $|0\rangle$ and the excited state denoted as $|1\rangle$, and the electronic system bilinearly couples with a total of $N$ bath modes. The spin-boson Hamiltonian expressed in terms of the system-bath type (Eq.(1)) is given by

Here, $\hat{\sigma}_x$ = $|0\rangle \langle 1|$+$|1\rangle \langle 0|$ and $\hat{\sigma}_z$ = $|0\rangle \langle 0|$$-$$|1\rangle \langle 1|$ are the Pauli matrices; $\Gamma$ is the electronic coupling, $\Delta E$ = $-2\varepsilon$ is the reaction energy, $\{ \hat{R}_j, \hat{P}_j,\omega_j\}$ = $\{ \hat{R}_j, \hat{P}_j,\omega_j | j$ = $1,...,N \}$ are the mass-weighted coordinates, momenta and frequencies associated with the $N$ nuclear normal modes, respectively, $\{c_j\}$ are the system-bath coupling coefficients between the electronic DOF and the nuclear DOF; the nuclear operator $\hat \Lambda$ = $\sum\limits_{j = 1}^N c_j\hat R_j$ is half of the energy gap between the $|0\rangle$ and $|1\rangle$ states. Alternatively, the spin-boson Hamiltonian can be cast into the perturbation type (Eq.(2)) that will be used with M-GQME as below:

For two-state systems, the bath mode frequencies and the system-bath coupling coefficients, $\{\omega_j, c_j\}$ are straightforwardly determined by the spectral density, which is defined by

The spectral density can be assumed to be Ohmic, Debye, a sum of Lorentzian modes [56], a sum of Debye and oscillatory Debye modes [167], a sum of damped harmonic modes [168], or other empirical forms [163, 169-172], as well as can be obtained from using hybrid approach that distinguishes high and low-frequency oscillators [173, 174], or utilizing projected modes [175]. Besides, more realistically, the spectral density can be obtained directly from all-atom molecular dynamics simulations by taking Fourier transform of the energy-gap time correlation function [164, 176], $C_{UU}(t)$ = $\langle U(t) U(0) \rangle$$-$$\langle U \rangle^2$, where $U$ = $H_{00}$$-$$H_{11}$ is the instantaneous energy gap between the two electronic states:

The discretization of the spectral density into $N$ ($N$$\gg$1) nuclear normal modes can be realized by solving the following equation for the frequency $\omega_j$ of the $j$-th mode: $(j = 1,2,\ldots,N)$ [164]:

Once the frequency of the $j$-th mode is determined, the corresponding system-bath coupling coefficients, $c_j$, is obtained via the relation

This procedure generates a bath frequency set that corresponds to equal partition of the reorganization energy given by

Using the above-mentioned approach, the condensed phase system (bath) information is incorporated to the model through the spectral density decomposition resulting in discrete bath frequencies and the electron-vibration couplings. It is recently demonstrated that one can go beyond two-state models and construct multi-state harmonic models from all-atom simulation [166, 177]. However, casting into the set of normal modes also means that the explicit anharmonicity is effectively neglected, despite that the harmonic modes are the best fit to the energy-gap time correlation function from all-atom MD simulation. Although the harmonic model construction is essential for mapping the all-atom interaction onto model Hamiltonian that is easy to evaluate, it is desirable to test more realistic Hamiltonians that are generally anharmonic force fields or are obtained with on-the-fly electronic structure calculations.

Ⅵ.   EMERGING DIRECTIONS OF QUANTUM DYNAMICS RESEARCH USING GQME

Here, we curate a topical discussion on the emerging directions of future quantum dynamics research on different areas and directions in the development of GQME.

A   Machine learning for dynamics of open quantum systems

The first one is concerned with the integration of machine learning techniques for accelerating or improving GQME-based methods for simulating quantum dynamics, especially for that of a condensed phase system. Neural network (NN) is one of the foremost effective schemes for performing machine learning. Recently, attempts have been made by several investigators to utilize NN to directly learn the time-local generators [178] and to solve the quantum many-body problems by using Boltzmann machine with NN embedded quantum states [179, 180]. Moreover, researchers started to use variational quantum Monte Carlo techniques for integrating quantum master equation to solve for the reduced density matrix of large size with a preferable computational scaling [181]. Brief accounts of many other relevant works that integrate machine learning and quantum dynamics could also be found in Ref.[182].

The performance of NN machine learning is often affected significantly by the depth of the hidden layers, which is a set of mathematical activation functions. Deep NNs with many hidden layers have found several applications for solving problems of growing complexity. In relation to quantum dynamics simulation, method for solving quantum master equation with the help of a deep quantum neural network (i.e., deep neural network deployed on quantum computing platform) was recently proposed and demonstrated to accurately reproduce the known benchmark results for 1-dimensional and 2-dimensional Ising models [183]. Convolutional neural networks were also recently demonstrated to be capable of predicting long-time dynamics based on short-time dynamics training pieces sampled from long-time trajectories [184].

Machine learning techniques based on NN implementations for various aspects of chemical dynamics were recently summarized by Kulichenko et al. [185]. On this line of development, methods for solving equilibrium continuity dynamics that can be employed as an approximate solution toward solving quantum Liouville dynamics [186], and for transforming the non-Markovian system dynamics into a Markovian system dynamics with the effective bath [187] were recently demonstrated. Machine learning for optimizing the model construction has been applied to represent molecular liquids as system-bath model [188], exciton and electron transfer phenomenon [189], and the effective dimension-reduced bath [190].

B   Beyond the common choice of projection operator

Up to now, we have introduced two forms of projection operators, one of which is the thermal projection operator or the so-called Feynman-Vernon type [191] (Eq.(17)), and the other of which is the electronic-nuclear type (Eq.(80)). In the past decade or so, several other ways of how the projection operators can be defined and utilized have also started to emerge. Notably it is the projection operator to account for strong system-bath coupling, which is also called the correlated projection operator that does not factorize the system and the bath completely. Several correlated projection operator forms have been derived and proven to give better performance when dealing with a system with large coupling [15, 192-194]. Another projection operator is the self-consistent projection operator that is called the consistent Mori projection operator (c-MoP) [195, 196]. The c-MoP was developed with a focus on improving the efficiency and accuracy of the numerical implementation [195, 196].

C   Quantum dynamics for challenging situations

Although GQME can give the numerically exact reduced system dynamics if fed with exact projection-free inputs, it still has several challenges in realistic applications. The challenging situation for GQME includes no available exact SDBCFs, long-lived memory, and no straightforward partitioning of the Hamiltonian. For obtaining accurate SDBCFs, there are difficult situations like in the nonergodic systems [197], conical intersections [71, 72] or correlated vibrations in the non-Condon case [198], very large number of system DOFs, and no effective quantum dynamics to obtain converged SDBCFs, etc. For the long-lived memory, it could be not advantageous to use GQME if the memory is not short-lived, in which case, directly evaluating the system dynamics would be equivalent to using GQME if not simpler [148]. Similarly when the system has complicated memory effects [199], very long time dynamics [200], or long-lived bath correlation due to structured environment [201], GQME may not be easy to simulate. The existing of the inhomogeneous term tends to be the case when the system has a long memory thus not easy to get decoupled from the bath for the initial preparation. For the partitioning of Hamiltonian, it is not always straightforward to express the overall Hamiltonian into a certain form even if using the perturbation type, such as the nuclear Hamiltonians and the electronic couplings in metallic systems. In addition, the definition of the projection operator depends on the partitioning of the overall Hamiltonian. Thus, both of Hamiltonian partitioning and projection operator definition will also influence the inhomogeneous term. Work on such directions is promising and new forms of GQME or memory kernel may be needed to circumvent some of the above-mentioned issues.

D   Quantum dynamics on quantum computer

In order to accelerate the methodology development and broaden the application of numerically exact quantum dynamics for open quantum systems including GQME, researchers start to develop open-source software and build quantum computing platforms. Versatile open-source simulation and analysis toolkits and platforms such as the Quantum Toolbox in Python (QuTiP) [202] and Quantum Dynamics Hub [203] become available to the community for free. Another key avenue that might not only aid but may also disrupt the quantum dynamics community is the use of quantum computing technology for solving GQME or other exact quantum dynamics. Along with this trend, application-friendly quantum computing toolkits such as the open-source quantum software project [204] and robust quantum programming languages (e.g., Qiskit [205], Cirq [206], and qsim [207]) emerged.

On the hardware side, the availability of cloud-based quantum computing hardware [5] for simulating quantum master equations and other quantum dynamics or thermodynamics also start to democratize the access to quantum computing. The algorithm that might be able to simulate quantum dynamics on quantum computers has been around for about a decade [208, 209]. The general process for simulating quantum dynamics on a quantum computer can be summarized in three steps: system preparation, imposing propagation and operation gates on the system, and system measurement. However, the real implementation on actual quantum computing platforms is often limited by the physical constraint inherent to the quantum computing hardware (e.g. length of system decoherence, circuit depth, quantum error). A general proposal for many-body simulation on a quantum computer was put forward a few years ago [210]. A proof-of-concept demonstration for simulating the FMO complex using a quantum computer was also recently reported [211]. Currently, quantum dynamics simulations can be implemented on the noisy intermediate scale quantum (NISQ) [212] hardware by implementing NISQ algorithms [213].

Excitingly, now the combination of quantum computing software and hardware for simulating quantum dynamics has been shown to have preferable scaling [211, 214, 215] allowing investigators to probe problems that are previously intractable or inconceivable to be tackled by a classical supercomputer.

To accelerate the applications, several authors provided tutorial guidance on how to formulate and implement quantum dynamics on quantum computers [215]. With more people gaining interest in this new field, more quantum dynamics simulation on quantum computers will be possible in the near future.

Ⅶ.   CONCLUDING REMARKS AND FUTURE OUTLOOK

In this review, we present a pedagogical introduction to the topic of the generalized quantum master equation (GQME) and summarize recent advances. We believe this review will be a useful material toward lowing the barrier for learning about the theory and implementation of GQME. GQME gives rise to numerically exact reduced system dynamics if the projection-free inputs are evaluated accurately. GQME in combination with the projection-free inputs evaluated with approximate quantum dynamical methods can generate more accurate results efficiently compared with using the same approximate method alone. Usually the memory kernel is short-lived, which will enable us to describe the arbitrary long time system dynamics without the simulation of the actual long-time quantum dynamics. In the tutorial sections, we provided a step-by-step derivation for the Nakajima-Zwanzig GQME and the time-convolutionless GQME using the projection operator technique. We also presented various schemes that can be used to express the memory kernel or TCL generator into the projection-free inputs given the choice of specific projection operators. For calculating the memory kernel in NZ-GQME, we showed the projection-free inputs that can be expressed in terms of the system-dependent bath correlation functions. A practical take for the numerical evaluation of the projection-free inputs was also provided. So far, GQME techniques have been applied to many model Hamiltonians, such as nonadiabatic dynamics in condensed phases and the excitonic energy transfer process in the light-harvesting photosynthetic complexes.

Furthermore, we reviewed the emerging trends of the development of GQME-related techniques. Specifically, we highlighted the attempts of using machine learning methods to overcome the bottleneck of the current GQME implementations, and the emergence of quantum algorithms and quantum computing platforms for simulating exact quantum dynamics make it hopeful for understanding increasingly large and complex systems.

Ⅷ.   ACKNOWLEDGMENTS

Xiang Sun acknowledges support from NYU Shanghai, the National Natural Science Foundation of China (No.21903054), the Hefei National Laboratory for Physical Sciences at the Microscale (No.KF2020008), the Shanghai Sailing Program (No.19YF1435600), and the Program for Eastern Young Scholar at Shanghai Institutions of Higher Learning.

APPENDIX A.   Mathematical Transforms and Expansions

1   Laplace transform

Laplace transform is useful for expressing differential equations in terms of algebraic equations that we encounter in deriving the NZ-GQME. Laplace transform for a given function $f(t)$ of a real argument $t$ is defined as

where the variable $s$ is complex. For the first order derivative $f'(t)$ = $\mathrm d f(t)/\mathrm dt$, it is easy to show using integration by parts that its Laplace transform is given by

Since exponential function has Laplace transform

we have

Also, we use the convolution theorem for Laplace transform:

2   Dyson decomposition

The Dyson decomposition was initially derived for unification of quantum electrodynamic theory [216]. The utility of Dyson decomposition of propagators is now ubiquitous both in classical or quantum statistical mechanics[91, 217, 218].

Applying Dyson decomposition to quantum propagators, we have

For example, substituting $\mathcal{A}$ = $\mathcal{L}$ and $\mathcal{B}$ = $\mathcal{Q}\mathcal{L}$ = $\mathcal{L}$$-$$\mathcal{P}\mathcal{L}$ into Eq.(A9), we obtain

Also note that the Dyson identity is commonly used in its Laplace transform (or resolvent) representation such that

The above identity is general, in the sense that it is applied to any operators.

3   Volterra equation

Most of GQMEs require to solve Volterra equation of the second kind [1, 38, 217, 219]. For a function $f(t), g(t)$, and $h(t,\tau)$, where the $f(t)$ is unknown and the other two functions is known, we have

The equation is to be solved iteratively until convergence, typically starts with an initial guess $f(t)$ = $g(t)$. Typically, the numerical evaluation used is method based on the trapezoidal rule, or by constructing a series of Legendre polynomial and performing integration by parts.

For instance, using the extended trapezoidal rule, and denoting the elements of $f(t)$ at iteration $i$ as $f_{jk}^i(t)$ as a function of time steps $t_n$ = $n\Delta t ; n = 0,1,2,\ldots,N_{\mathrm{max}}$, we have

Here, $h_{jl}$ and $g_{jk}$ denote the matrix elements of the function $h(t,\tau)$ and $g(t)$, respectively. Taking the Volterra equation shown in Eq.(90) as an example, $f(t)$$\rightarrow$$\mathcal{K}^{\mathrm{M}}(t)$, $g(t)$$\rightarrow$$i\dot{\mathcal{F}}(t)$$-$$\frac{1}{\hbar}\mathcal{F}(t) \langle\mathcal{L}\rangle_{\mathrm{N}}^0$, and $h(t,\tau)$ = $i\mathcal{F}(t-\tau)$.

APPENDIX B.   Proof of Eq.(56)

First note that the LHS of Eq.(56) becomes

Then, we consider the three contributions of splitting the leftmost Liouville operator $\mathcal{L}$ = $\mathcal{L}_{ {\rm{S}}}$+$\mathcal{L}_{ {\rm{B}}}$+$\mathcal{L}_{ {\rm{BS}}}$. First, for the bath $\mathcal{L}_{ {\rm{B}}}$ we have

for any operator $\hat A$ = $\sum\limits_{j} \hat S_j \otimes \hat B_j$, where $\hat S_j$ and $\hat B_j$ are system and bath operators, respectively, since it can be shown that $\mathrm{Tr}_ {\rm{B}}\left[\mathcal L_ {\rm{B}} \hat A\right]$ = $\mathrm{Tr}_ {\rm{B}}\left\{\left[ \hat H_ {\rm{B}},\hat A\right]\right\}$ = $\sum\limits_j \hat S_j \otimes \mathrm{Tr}_ {\rm{B}}\left\{\left[ \hat H_ {\rm{B}},\hat B_j\right]\right\}$ = 0.

Second, for the system $\mathcal L_ {\rm{S}}$, we have for any operator $\hat A$

since

and thus using Eq.(B3) we prove Eq.(B2) as below

Thus, only the system-bath coupling $\mathcal L_{ {\rm{BS}}}$ on the leftmost survives, i.e.

Next, consider rightmost Liouville operator components in Eq.({1.5}), namely $\mathcal {QL} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma$ = $\mathcal {QL}_ {\rm{S}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma$+ $\mathcal {QL}_ {\rm{B}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma$+$\mathcal {QL}_{ {\rm{BS}}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma$. One can show that

since

1. $\mathcal {QL}_ {\rm{S}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma = \mathcal {L}_ {\rm{S}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma -\mathcal {PL}_ {\rm{S}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma \nonumber\\ = \mathcal {L}_ {\rm{S}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma - \hat \rho_ {\rm{B}}^{\mathrm{eq}} \mathrm{Tr}_ {\rm{B}}\left[ \hat \rho_ {\rm{B}}^{\mathrm{eq}}\right] \mathcal {L}_ {\rm{S}} \hat \sigma\\ = \mathcal {L}_ {\rm{S}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma - \hat \rho_B^{\mathrm{eq}} \mathcal {L}_ {\rm{S}} \hat \sigma \\ = \mathcal {L}_ {\rm{S}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma - \mathcal {L}_ {\rm{S}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma = 0$

2. $\mathcal {L}_ {\rm{B }} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma = \left(\mathcal {L}_ {\rm{B}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\right) \hat \sigma = \left[ \hat H_ {\rm{B}}, \hat \rho_ {\rm{B}}^{\mathrm{eq}}\right] \hat \sigma = 0$

3. $\mathcal {QL}_{ {\rm{BS}}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma = \mathcal {L}_{ {\rm{BS}}} \hat \rho_B^{\mathrm{eq}}\hat \sigma - \mathcal {PL}_{ {\rm{BS}}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma \\ = \mathcal {L}_{ {\rm{BS}}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma - \hat \rho_ {\rm{B}}^{\mathrm{eq}} \mathrm{Tr}_ {\rm{B}} \left[ \mathcal {L}_{BS} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\right] \hat \sigma \\ = \mathcal {L}_{ {\rm{BS}}} \hat \rho_ {\rm{B}}^{\mathrm{eq}}\hat \sigma$

where $\operatorname{Tr}_{\mathrm{B}}\left[\mathcal{L}_{\mathrm{BS}} \hat{\rho}_{\mathrm{B}}^{\mathrm{eq}}\right]=0$ using Eq.(53).

Therefore, we arrive at Eq.(56):