Chinese Journal of Chemical Physics  2018, Vol. 31 Issue (3): 245-256

The article information

Yang Liu, Rui-xue Xu, Hou-dao Zhang, YiJing Yan

Dissipaton Equation of Motion Theory versus Fokker-Planck Quantum Master Equation

Chinese Journal of Chemical Physics, 2018, 31(3): 245-256

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

Article history

Accepted on: May 28, 2018
Dissipaton Equation of Motion Theory versus Fokker-Planck Quantum Master Equation
Yang Liu, Rui-xue Xu, Hou-dao Zhang, YiJing Yan
Dated: Received on April 27, 2018; Accepted on May 28, 2018
Hefei National Laboratory for Physical Sciences at the Microscale and Department of Chemical Physics and Collaborative Innovation Center of Chemistry for Energy Materials(iChEM), University of Science and Technology of China, Hefei 230026, China
*Author to whom correspondence should be addressed. YiJing Yan, yanyj@ustc.edu.cn
Abstract: The quest of exact and nonperturbative methods on quantum dissipation with nonlinear coupling environments remains in general a great challenge. In this review we present a comprehensive account on two approaches to the entangled system-and-environment dynamics, in the presence of linear-plus-quadratic coupling bath. One is the dissipaton-equation-ofmotion (DEOM) theory that has been extended recently to treat the nonlinear coupling environment. Another is the extended Fokker-Planck quantum master equation (FP-QME) approach that will be constructed in this work, based on its DEOM correspondence. We closely compare these two approaches, with the focus on the underlying quasi-particle picture, physical implications, and implementations.
Key words: Quadratic couping bath    Dissipaton equation of motion    Fokker-Planck quantum master equation    Quantum dissipation theory
Ⅰ. INTRODUCTION

Quantum dissipation plays crucial roles in many fields of science. Various theoretical methods have been proposed since 1950s, with the focus on reduced system dynamics of primary interest. The revival research witnessed recently in quantum dissipation theory owes partially to its close relation to quantum devices and nanomaterials functions. Accurate and efficient methods are in need nowadays. Exact approaches include the path-integral influence functional formalism [1], its differential equivalences via the hierarchical-equations-of-motion (HEOM) [2-8] and the related stochastic methods [9-12]. However, these Feynman-Vernon's influence functional based theories exploit the Gaussian-Wick's theorem that is strictly valid only for linear coupling bath [13-16]. This would have intrinsically assumed a weak backaction of system on environment.

Exact and nonperturbative treatments would be needed, at least up to the quadratic coupling bath level. The quest here is very challenging, as it goes beyond the Gaussian-Wick's theorem. Along the development of the dissipaton-equation-of-motion (DEOM) approach [17-19], we had recently extended this theory with nonlinear coupling bath [20, 21]. The DEOM theory is a second-quantization advancement of the HEOM formalism [2-7]. It offers a unified quasi-particle (dissipaton) description for coupling environments that can be either bosonic, fermionic, or hard-core bosonic [17]. Dynamical variables are now no longer mathematical auxiliaries, but physically well-defined quantities [17-19]. Beside the hierarchical dynamics equations governing evolutions, DEOM comprises also the novel dissipaton algebra, especially the generalized Wick's theorem. This enables DEOM a theory for not only the reduced system, but also the hybrid bath dynamics [17-19].

This review presents a comprehensive account on the extended DEOM theory, with the linear-plus-quadratic coupling environments [20, 21]. Moreover, we will also construct the Fokker-Planck quantum master equation (FP-QME) via the corresponding DEOM. The FP-QME approach treats explicitly the reduced "core-system" dynamics [22-28]. It starts with dividing the overall bath into the "first-shell" (solvation hereafter) and "secondary" parts. The core system comprises both the primary system and the directly coupled solvation modes that experience Brownian oscillators motion in the secondary bath environment. The most representing FP-QME is the Caldeira-Leggett's QME [27, 28]. However, it involves the high-temperature treatment on the secondary bath. This significantly compromises the range of applicability, as the underlying solvation dynamics is basically classical [22-28]. Caustic corrections are needed to bridge between the classical secondary bath and the primary quantum systems. An improved FP-QME, based on its DEOM equivalence, has been proposed recently [29-31], with the caustics being treated in a semiclassical manner.

Ⅱ. PRELUDE

We will be interested in both DEOM and FP-QME theories of an open quantum system in the presence of linear-plus-quadratic coupling bath environment. Without loss of generality, we focus on the single-dissipative mode case. The total system-and-bath composite Hamiltonian assumes the form of [21]

 $\begin{eqnarray}\label{HT} H_{\rm{T}}=H_{0} + h_{\rm{B}} + \hat Q_{\rm{S}}(\alpha_0+\alpha_1\hat x_{\rm{B}}+\alpha_2\hat x_{\rm{B}}^2) \end{eqnarray}$ (2.1)

While the Hamiltonian $H_{0}$ and dissipative mode $\hat Q_{\rm{S}}$ of system are arbitrary, the bath Hamiltonian $h_{\rm{B}}$ and solvation coordinate $\hat x_{\rm{B}}$ assume

 $\begin{eqnarray}\label{hB0} \begin{split} h_{\rm{B}}&= \frac{1}{2}\sum\limits_j \omega_j(\hat p_j^2 + \hat q_j^2)\\ \hat x_{\rm{B}}&=\sum\limits_j c_j \hat q_j \end{split} \end{eqnarray}$ (2.2)

Throughout this paper we set $\hbar$=1 and $\beta$=$1/(k_{\rm{B}}T)$, with $k_{\rm{B}}$ being the Boltzmann constant and $T$ the temperature. Adopt the dimensionless solvation coordinate, $\hat x_{\rm{B}}$ in Eq.(2.2). Let $\langle \hat O \rangle_{\rm{B}}\equiv$ ${\rm tr}_{\rm{B}} (\hat O \text{e}^{-\beta h_{\rm{B}}})/{\rm tr}_{\rm{B}}\text{e}^{-\beta h_{\rm{B}}}$, an average over the bare bath thermal equilibrium ensembles. The effective system Hamiltonian, as inferred from Eq.(2.1), would assume

 $\begin{eqnarray}\label{H_eff} H^{\rm eff}_{\rm{S}} &=& H_{0}+ (\alpha_0 + \alpha_2\langle\hat x^2_{\rm{B}}\rangle_{\rm{B}})\hat Q_{\rm{S}}\nonumber\\ &\equiv & H_{\rm{S}}+ \alpha_2 \langle\hat x^2_{\rm{B}}\rangle_{\rm{B}}\hat Q_{\rm{S}} \end{eqnarray}$ (2.3)

Denote $\hat x_{\rm{B}}(t)\equiv \text{e}^{ih_{\rm{B}}t}\hat x_{\rm{B}}\text{e}^{-ih_{\rm{B}}t}$. Set hereafter $t\geq$0 for the time variable, unless further specified.

It is worth noting that Eq.(2.1) is just one prescription of the total composite Hamiltonian, with the specified reference bath of $h_{\rm{B}}$. In the presence of nonlinear coupling bath, different reference bath Hamiltonians are no longer mirror image to each other. It requires that the total composition $H_{\rm{T}}$ be of the prescription invariance, with difference reference bath Hamiltonians. Moreover, the characterization on the coupling bath $\alpha$-parameters in Eq.(2.1), including $\alpha_0$ for the energy reorganization, goes beyond the linear response theory. Physically supported models would be needed to complete the characterization of the $\alpha$-parameters [20, 21]; see Section Ⅳ for further elaborations.

It is also noticed that the solvation coordinate $\hat x_{\rm{B}}$, as described in Eq.(2.2), is a macroscopic bath degree of freedom. It couples directly to the system dissipative mode, $\hat Q_{\rm{S}}$, and meanwhile it is subject to a Brownian motion in the secondary bath environment. In contract to the DEOM theory in Section Ⅲ, the extended FP-QME to be developed in Section Ⅳ describes the Brownian motion explicitly. In principle, one can have the "core-DEOM" theory, with the linear coupling secondary environment. However, this approach is much more expensive in general, as it accesses also the dynamics of the secondary dissipatons. We will see later that the effects of nonlinear coupling bath could be taken into account via the generalized Wick's theorem (cf.Ⅲ.B).

Turn to the characterization of dissipatons. These are statistical quasi-particles, defined via the linear solvation correlation function in the bare bath ensembles at thermal equilibrium [17-19]. Let $\chi_{\rm{B}}(\omega)$ be the Laplace frequency resolution of the solvation response function. Its real part is symmetric and imaginary part is antisymmetric. The latter reads [14]

 $\begin{eqnarray}\label{chiJw} {\rm Im}\chi_{\rm{B}}(\omega) = \frac{1}{2}\!\int^{\infty}_{-\infty}{\rm{d}} t\, \text{e}^{i\omega t}\langle [\hat x_{\rm{B}}(t), \hat x_{\rm{B}}(0)]\rangle_{\rm{B}} \end{eqnarray}$ (2.4)

This defines the solvation bath spectral density, in agree with ${\rm Im}\chi_{\rm{B}}(\omega>0)$=$\displaystyle\frac{\pi}{2}\sum_j c^2_j\delta(\omega-\omega_j)$, the microscopic formula following directly from Eq.(2.2). Together with the detailed-balance relation, Eq.(2.4) leads to [14]

 $\begin{eqnarray}\label{FDT} \langle\hat x_{\rm{B}}(t)\hat x_{\rm{B}}(0)\rangle_{\rm{B}} = \frac{1}{\pi} \int^{\infty}_{-\infty}\!\! {\rm d}\omega\, \frac{\text{e}^{-i\omega t}{\rm Im}\chi_{\rm{B}}(\omega)}{1-\text{e}^{-\beta\omega}} \end{eqnarray}$ (2.5)

This is the fluctuation-dissipation theorem [13-16]. The time-reversal relation for correlation functions reads [14]

 $\begin{eqnarray}\label{time-reversal} \langle\hat x_{\rm{B}}(0)\hat x_{\rm{B}}(t)\rangle_{\rm{B}} = \langle\hat x_{\rm{B}}(t)\hat x_{\rm{B}}(0)\rangle^{\ast}_{\rm{B}} \end{eqnarray}$ (2.6)

To proceed, we exploit a certain efficient and accurate sum-over-poles scheme [32, 33] to decompose the Fourier integrand in Eq.(2.5), followed by the Cauthy's contour integration in the lower-half plane. We obtain

 $\begin{eqnarray}\label{FF_exp_all} \langle\hat x_{\rm{B}}(t)\hat x_{\rm{B}}(0)\rangle_{\rm{B}} = \sum\limits_{k=1}^K \eta_k \text{e}^{-\gamma_k t} \end{eqnarray}$ (2.7)

It is easy to show that (ⅰ) the individual exponential term arising from the Bose function is real, and (ⅱ) those $\{\gamma_k\}$ arising from the spectral density function are either real or complex conjugate-paired [8, 18]. In other words, the associate index, $\bar k$, defined via $\gamma_{\bar k}\equiv \gamma^{\ast}_{k}$, must also appear in Eq.(2.7). We can therefore express the time-reversal counterpart to Eq.(2.7) via Eq.(2.6) as

 $\begin{eqnarray}\label{time-reversal_2} \langle\hat x_{\rm{B}}(0)\hat x_{\rm{B}}(t)\rangle_{\rm{B}} = \sum\limits_{k=1}^K\eta^{\ast}_{\bar k}\text{e}^{-\gamma_k t} \end{eqnarray}$ (2.8)

The dissipaton decomposition on the solvation coordinate goes as follows,

 $\begin{eqnarray}\label{xB_in_f} \hat x_{\rm{B}} = \sum\limits_{k=1}^K \hat f_{k} \end{eqnarray}$ (2.9)

with

 $\begin{eqnarray}\label{ff_t} \begin{split} \langle\hat f_{k}(t)\hat f_{j}(0)\rangle_{\rm{B}}&=\delta_{kj}\eta_{k}\text{e}^{-\gamma_k t} \\ \langle\hat f_{j}(0)\hat f_{k}(t)\rangle_{\rm{B}}&=\delta_{kj}\eta^{\ast}_{\bar k}\text{e}^{-\gamma_k t} \end{split} \end{eqnarray}$ (2.10)

Apparently, both Eqs. (2.7) and (2.8) are satisfied.

As highlighted above, dissipatons consist of linear and statistically independent macroscopic degree of freedom, arising from the solvation coordinate. Each dissipaton goes by a single exponent, for its forward and backward correlation functions. This feature leads to [17, 18]

 $\begin{eqnarray}\label{diff_0} {\rm tr}_{\rm{B}}\left[\Big(\frac{\partial \hat f_k}{\partial t}\Big)_{\rm{B}}\rho_{\rm{T}}(t)\right] = -\gamma_{k}\, {\rm tr}_{\rm{B}}\big[\hat f_k\rho_{\rm{T}}(t)\big] \end{eqnarray}$ (2.11)

It will lead to the generalized diffusion equation in the DEOM theory, see Eq.(3.2).

Ⅲ. DISSIPATON-EQUATION-OF-MOTION THEORY A. Second quantization description

Dynamical variables in the DEOM theory are the so-called dissipaton density operators (DDOs) [17, 18]:

 $\begin{eqnarray}\label{DDO_def} \rho^{(n)}_{\bf n}(t)&\equiv & \rho^{(n)}_{n_1\cdots n_K}(t)\nonumber\\ &\equiv & {\rm tr}_{\rm{B}}\big[ (\hat f_{K}^{n_K}\cdots\hat f_{1}^{n_1})^{\circ}\rho_{\rm{T}}(t)\big] \end{eqnarray}$ (3.1)

Here, $n$=$n_1$+$\cdots$+$n_K$, with $n_k$=0, 1, 2, $\cdots$. The reduced system density operator is just $\rho_{\rm{S}}(t)\equiv \rho^{(0)}(t)$. In general, $\rho^{(n)}_{\bf n}(t)$ is an $n$-dissipatons DDO, with the configuration of ${\bf n}\equiv \{n_1\cdots n_K\}$ that is an ordered set of occupation numbers on individual dissipatons. More precisely, the product of $n$-dissipatons inside the circled parentheses, $(\cdots)^{\circ}$, is irreducible. It goes by $(\hat f_k\hat f_j)^{\circ}$=$(\hat f_j\hat f_k)^{\circ}$, for Bosonic dissipatons that satisfy symmetric permutation. In other words, Eq.(3.1) resembles the second-quantization representation of a bosonic permanent. The DDOs with fermionic environment are similar, but resemble the Slater determinants, having the occupation numbers of $0$ or $1$ only, due to the antisymmetric permutation relation [17, 18]. The irreducibility notion, $(\cdots)^{\circ}$, is also closely related to the generalized Wick's theorem (cf. Section Ⅲ.C). This is the most useful ingredient of the dissipaton algebra, which enables DEOM truly a type of "core-system" theory. Apparently, the DEOM space resembles a Fock space in quantum mechanics, with individual elements being the form of ${\mathit{\boldsymbol{\rho}}}(t)$=$\big\{\rho^{(n)}_{\textbf{n}}(t); n$=0, 1, $\cdots, L\big\}$. Here, $L$ denotes the converged many-dissipaton level.

The generalized diffusion equation, which dictates the influence of $h_{\rm{B}}$ on DDOs without approximation, reads [17, 18]

 $\begin{eqnarray}\label{diff} i\, {\rm tr}_{\rm{B}}\big\{ (\hat f_{K}^{n_K}\cdots\hat f_{1}^{n_1})^{\circ} [h_{\rm{B}}, \rho_{\rm{T}}]\big\} = \Big(\sum\limits_{k=1}^{K} n_k\gamma_k\Big)\rho^{(n)}_{\bf n} \end{eqnarray}$ (3.2)

This is obtained via Eq.(2.11), together with the Heisenberg equation, $(\partial \hat f_k/\partial t)_{\rm{B}}$=$-i[\hat f_k, h_{\rm{B}}]$ and the trace cyclic invariance. It is noticed that in FP-QMEs, the influence of $h_{\rm{B}}$ is described by the FP operator that involves a certain bi-dissipaton ansatz (cf. Section Ⅴ.B). Equation (3.2) does generalize the eigenequation of the FP operator (cf. Eq.(6.18)). It resembles the second-quantization result on the influence of bare-bath $h_{\rm{B}}$.

B. The generalized Wick's theorem

The generalized Wick's theorem is the most important ingredient of the dissipaton algebra [17-21]. It enables DEOM a theory for entangled system-and-bath dynamics [34-37], on top of the hierarchical dynamics an algebraic construction (cf. Section Ⅲ.C). For bookkeeping, we denote the values of Eq.(2.10) at $t\rightarrow$0+ by

 $\begin{eqnarray}\label{ff0} \begin{split} %\label{ff0_a} \langle\hat f_{k}\hat f_{j}\rangle^{>}_{\rm{B}} &\equiv \langle\hat f_{k}(0+)\hat f_{j}(0)\rangle_{\rm{B}} =\delta_{kj}\eta_{k} \\ %\label{ff0_b} \langle\hat f_{j}\hat f_{k}\rangle^{<}_{\rm{B}} &\equiv \langle\hat f_{j}(0)\hat f_{k}(0+)\rangle_{\rm{B}} =\delta_{kj}\eta^{\ast}_{\bar k} \end{split} \end{eqnarray}$ (3.3)

Moreover, the irreducibility notion, $(\cdots)^{\circ}$ in Eq.(3.1), facilitates the presentations below.

Generalized Wick's Theorem-1 (GWT-1) deals with the case of adding one dissipaton each time to DDOs of Eq.(3.1). For the forward action of a dissipaton, it reads [17, 18]:

 $\begin{eqnarray}\label{Wick10} &&{\rm tr}_{\rm{B}}\big[(\hat f_{K}^{n_K}\cdots\hat f_{1}^{n_1})^{\circ} \hat f_j\rho_{\rm{T}}(t)\big]\nonumber\\ &&= \rho^{(n+1)}_{{\bf n}^{+}_{j}}(t) + \sum\limits_{k=1}^{K} n_k \langle\hat f_{k}\hat f_{j}\rangle^{>}_{\rm{B}}\rho^{(n-1)}_{{\bf n}^{-}_{k}}(t) \end{eqnarray}$ (3.4)

The associated($n\pm$1)-particle DDO, $\rho^{(n \pm 1)}_{{\bf n}^{\pm}_{k}}(t)$, has the configuration ${\bf n}^{\pm}_{k}$ that differs from ${\bf n}\equiv \{n_1\cdots n_K\}$ by replacing the specified $n_k$ with $n_k\pm$1. The backward-action counterpart, ${\rm tr}_{\rm{B}}\big[(\hat f_{K}^{n_K}\cdots\hat f_{1}^{n_1})^{\circ} \rho_{\rm{T}}(t)\hat f_j\big]$, is similar, but with $\langle\hat f_{k}\hat f_{j}\rangle^{>}_{\rm{B}}$ being replaced by $\langle\hat f_{j}\hat f_{k}\rangle^{<}_{\rm{B}}$; see Eq.(3.3) for their values. By using Eq.(2.9), we report GWT-1 as

 $\begin{eqnarray}\label{Wick_xB} &&{\rm tr}_{\rm{B}}\big[(\hat f_{K}^{n_K}\cdots\hat f_{1}^{n_1})^{\circ} \hat x_{\rm{B}}\rho_{\rm{T}}(t)\big] =\nonumber\\ &&\sum\limits_{k=1}^K \rho^{(n+1)}_{{\bf n}^{+}_{k}}(t) + \sum\limits_{k=1}^{K} n_k \eta_{k} \rho^{(n-1)}_{{\bf n}^{-}_{k}}(t) \end{eqnarray}$ (3.5)

The backward-action counterpart is similar, but with $\eta_{k}$ being replaced by $\eta^{\ast}_{\bar k}$. Apparently, GWT-1 deals with the linear coupling bath.

GWT-2 is concerned with the quadratic coupling bath, where a pair of dissipatons are added simultaneously each time. This would suggest that GWT-2 read [20, 21]:

 $\begin{eqnarray} \label{Wick2} &&{\rm tr}_{\rm{B}}\big[(\hat f_{K}^{n_K}\cdots\hat f_{1}^{n_1})^{\circ} (\hat f_j\hat f_{j'})\rho_{\rm{T}}(t)\big] \nonumber\\ &&=\langle\hat f_j \hat f_{j'}\rangle_{\rm{B}}\rho^{(n)}_{{\bf n}}(t) +\rho^{(n+2)}_{{\bf n}^{++}_{jj'}}(t)+\nonumber\\ && \sum\limits_{k} n_k \langle\hat f_k\hat f_{j'}\rangle^{>}_{\rm{B}}\rho^{(n)}_{{\bf n}^{-+}_{kj}}(t) +\sum\limits_{k} n_k \langle\hat f_k\hat f_j\rangle^{>}_{\rm{B}}\rho^{(n)}_{{\bf n}^{-+}_{kj'}}(t) +\nonumber\\ && \sum\limits_{k, k'} n_k(n_{k'}-\delta_{kk'})\langle\hat f_k\hat f_j\rangle^{>}_{\rm{B}}\langle\hat f_{k'}\hat f_{j'}\rangle^{>}_{\rm{B}} \rho^{(n-2)}_{{\bf n}^{-}_{k k'}}(t) \end{eqnarray}$ (3.6)

The $\langle\hat f_j \hat f_{j'}\rangle_{\rm{B}}$ in the first term takes at the same time for the specified dissipatons pair that are added simultaneously. The associated DDO index, ${\bf n}^{\pm\pm}_{kj}$, differs from ${\bf n}\equiv \{n_1\cdots n_K\}$ on the specified subindexes, $n_k$ and $n_j$, that are replaced by $n_k\pm$1 and $n_j\pm$1, respectively. By using Eqs. (2.9) and (3.3), we can recast GWT-2 as

 $\begin{eqnarray}\label{Wick_xB2} &&{\rm tr}_{\rm{B}}\big[(\hat f_{K}^{n_K}\cdots\hat f_{1}^{n_1})^{\circ} \hat x^2_{\rm{B}}\rho_{\rm{T}}(t)\big]\nonumber\\ &&=\langle\hat x^2_{\rm{B}}\rangle_{\rm{B}}\rho^{(n)}_{\bf n}(t) +\sum\limits_{kj}\rho^{(n+2)}_{{\bf n}^{++}_{kj}}(t) +2\sum\limits_{kj} n_k\eta_k\rho^{(n)}_{{\bf n}^{-+}_{kj}}(t)+ \nonumber\\ &&\quad\sum\limits_{kj} n_k(n_j-\delta_{kj})\eta_k\eta_j\rho^{(n-2)}_{{\bf n}^{-}_{kj}}(t) \end{eqnarray}$ (3.7)

The backward-action counterpart is similar, but with $\eta_k$ and $\eta_j$ being replaced by $\eta^{\ast}_{\bar k}$ and $\eta^{\ast}_{\bar j}$, respectively.

C. Hierarchical dynamics equations

The complete DEOM theory consists also the hierarchical dynamics formalism for DDOs. Let us start with the time derivative on Eq.(3.1),

 $\begin{eqnarray}\label{dot_DDO} \dot\rho^{(n)}_{\bf n}(t)&\equiv & \rho^{(n)}_{n_1\cdots n_K}(t)\nonumber\\ &\equiv & {\rm tr}_{\rm{B}}\big[ (\hat f_{K}^{n_K}\cdots\hat f_{1}^{n_1})^{\circ} {\dot\rho}_{\rm{T}}(t)\big] \end{eqnarray}$ (3.8)

where $\dot\rho_{\rm{T}}(t)=-i[H_{\rm{T}}, \rho_{\rm{T}}(t)]$. We evaluate Eq.(3.8) with the each individual term in the total system-and-bath composite $H_{\rm{T}}$ of Eq.(2.1). While the generalized diffusion Eq.(3.2) has taken care of the bath $h_{\rm{B}}$ action, the generalized Wick's theorem is the vehicle to the algebraic evaluation on the influence of system-coupling bath. In particular, GWT-1 and GWT-2, Eqs. (3.5) and (3.7) and their backward-action counterparts, evaluate the influences of linear and quadractic coupling bath, respectively. Apparently, $\langle\hat x^2_{\rm{B}}\rangle_{\rm{B}}$ in the first term of Eq.(3.7) enters into the effective system Hamiltonian, $H^{\rm eff}_{\rm{S}}$ of Eq.(2.3).

The final hierarchical dynamics formalism reads

 $\begin{eqnarray}\label{DEOM} \dot\rho^{(n)}_{\bf n} \hspace{-0.12cm}& =&\hspace{-0.12cm} - \Big(i{\cal L}^{\rm eff}_{\rm{S}} + \sum\limits_{k}n_k\gamma_k\Big)\rho^{(n)}_{\bf n} -i2\alpha_2 \sum\limits_{kj} n_k {\cal C}_k \rho^{(n)}_{{\bf n}^{-+}_{kj}} -\nonumber\\ &&\hspace{-0.12cm}i\alpha_2\sum\limits_{kj} \Big[{\cal A}\rho^{(n+2)}_{{\bf n}^{++}_{kj}} +n_k(n_{j}-\delta_{jk}){\cal B}_{kj}\rho^{(n-2)}_{{\bf n}^{-}_{kj}} \Big]-\nonumber\\ &&\hspace{-0.12cm}i\alpha_1\sum\limits_{k}\Big({\cal A}\rho^{(n+1)}_{{\bf n}^{+}_{k}} +n_k{\cal C}_k\rho^{(n-1)}_{{\bf n}^{-}_{k}}\Big) \end{eqnarray}$ (3.9)

Here ${\cal L}^{\rm eff}_{\rm{S}}\hat O$=$[H^{\rm eff}_{\rm{S}}, \hat O]$, and

 $\begin{eqnarray}\label{calABC} \begin{split} &{\cal A}\hat O \equiv \big[\hat Q_{\rm{S}}, \hat O]\\ &{\cal C}_k\hat O \equiv \eta_k \hat Q_{\rm{S}}\hat O -\eta^{\ast}_{\bar k}\hat O\hat Q_{\rm{S}}\\ &{\cal B}_{kj} \hat O \equiv \eta_k\eta_j \hat Q_{\rm{S}} \hat O - \eta^{\ast}_{\bar k}\eta^{\ast}_{\bar j}\hat O \hat Q_{\rm{S}} \end{split} \end{eqnarray}$ (3.10)

In the absence of quadratic coupling ($\alpha_2$=0), Eq.(3.9) reduces to the HEOM formalism [3-6]. The latter is a path integral influence functional based theory [1], with $\{\rho^{(n>0)}_{\bf n}\}$ being considered as mathematical auxiliaries. The above observations validate the generalized diffusion Eq.(2.11) that leads to Eq.(3.2), and GWT-1, Eq.(3.4) or Eq.(3.5), and the time-reversal counterpart. While the conventional diffusion equation [38] and Wick's theorem [14, 15] are concerned only with $\cal C$-number properties of Gaussian statistics, those two generalizations, Eqs. (2.11) and (3.4), go with reduced operators in the system-subspace, with arbitrary Hamiltonian $H_{\rm{S}}$ and arbitrary dissipative mode $\hat Q_{\rm{S}}$.

Nevertheless, GWT-2, Eq.(3.6) or Eq.(3.7), is subject to further scrutiny. Analytically, we had validated it in the high-temperature regime [20]; however, numerically, we found it would be only applicable up to a moderate strength of quadratic coupling bath [21]. The above observations would suggest that Eq.(3.7) would be a sort of Ehrenfest mean-field treatment on quadratic coupling bath [21]. We will defer this issue to future study.

Ⅳ. THEORY OF NONLINEAR COUPLING BATH DESCRIPTORS

To complete the quantum dissipation theory with nonlinear coupling environments, we shall further have physically supported $\alpha$-parameters in Eq.(2.1). It is well known that in the presence of quadratic coupling bath, the $\alpha$-parameters cannot be characterized solely via the linear response theory. Additional information is needed.

In the follow we present an extended solvation model approach to the required linear-plus-quadratic coupling bath descriptors [21]. Let us examine the system-bath interaction, the last term of Eq.(2.1). Physically it arises from the surrounding environment rearrangements in response to the dissipative system operator $\hat Q_{\rm{S}}$. The solvation model approach is to rewrite Eq.(2.1) as

 $\begin{eqnarray}\label{HT0} H_{\rm{T}}=H_{0} + h_{\rm{B}}+ \hat Q_{\rm{S}}(h'_{\rm{B}}-h_{\rm{B}}) \end{eqnarray}$ (4.1)

with

 $\begin{eqnarray}\label{del_hB_def} \delta h_{\rm{B}} &\equiv & h'_{\rm{B}}-h_{\rm{B}}\nonumber\\ &=& \alpha_0 + \alpha_1\hat x_{\rm{B}}+\alpha_2\hat x_{\rm{B}}^2 \end{eqnarray}$ (4.2)

Meanwhile, we examine the reference-bath invariance requirement by recasting the last two terms of Eq.(4.1) as

 $\begin{eqnarray}\label{hBhB} h_{\rm{B}} + \hat Q_{\rm{S}}(h'_{\rm{B}}-h_{\rm{B}}) =h'_{\rm{B}} + (1-\hat Q_{\rm{S}})(h_{\rm{B}}-h'_{\rm{B}}) \end{eqnarray}$ (4.3)

This identity relates between the reference-$h_{\rm{B}}$ and the reference-$h'_{\rm{B}}$ based descriptions, with $\delta h_{\rm{B}}\equiv h'_{\rm{B}}-h_{\rm{B}}$ being given by Eq.(4.2), whereas

 $\begin{eqnarray}\label{del_hBe_def} \delta h'_{\rm{B}}&\equiv &h_{\rm{B}}-h'_{\rm{B}}\nonumber\\ &=& \alpha'_0+\alpha'_1\hat x'_{\rm{B}}+\alpha'_2\hat x^{\prime 2}_{\rm{B}} \end{eqnarray}$ (4.4)

We can then recast Eq.(4.3) in terms of

 $\begin{eqnarray}\label{alpha_alpha} \alpha_0+\alpha_1\hat x_{\rm{B}}+\alpha_2\hat x^{2}_{\rm{B}} =-(\alpha'_0+\alpha'_1\hat x'_{\rm{B}}+\alpha'_2\hat x^{\prime 2}_{\rm{B}}) \end{eqnarray}$ (4.5)

Here the primed quantities and unprimed ones are associated with the $h'_{\rm{B}}$- and the $h_{\rm{B}}$-based prescriptions, respectively. This follows the spectroscopic notation, as if $h'_{\rm{B}}$ and $h_{\rm{B}}$ were the molecular nuclear Hamiltonians associated with electronic excited and ground states, respectively [39, 40]. Apparently, the linear and quadratic terms arise from the linear displacement and frequency shift, respectively; i.e.,

 $\begin{eqnarray}\label{hatx'} \hat x' = (\Omega'/\Omega)^{1/2}(\hat x-d) \end{eqnarray}$ (4.6)

To proceed, let us divide the overall bath, $h_{\rm{B}}$, into the solvation mode and secondary bath parts. In line with the complete square-form of $h_{\rm{B}}$ in Eq.(2.2), its partition would read [27, 28]

 $\begin{eqnarray}\label{hB} h_{\rm{B}} = \frac{1}{2}\Omega\big(\hat p^2_{\rm{B}}+\hat x^2_{\rm{B}}\big) + \frac{1}{2}\sum\limits_{k} \tilde\omega_k \Big[ \tilde p_k^{2} + \Big(\tilde x_k-\frac{\tilde c_k} {\tilde\omega_k}\hat x_{\rm{B}}\Big)^2\Big]\nonumber\\ \end{eqnarray}$ (4.7)

As implied here the solvation mode is subject to a Brownian oscillator motion in the secondary bath, with $\tilde h_{\rm{B}}$=$\displaystyle\frac{1}{2}\sum_j \tilde\omega_k(\tilde p^2_k+\tilde x^2_k)$. The latter gives rise to the stochastic random force:

 $\begin{eqnarray} \tilde F(t)=\text{e}^{i\tilde h_{\rm{B}}t}\Big(\sum\limits_{k}\tilde c_k\tilde x_{k}\Big)\text{e}^{-i\tilde h_{\rm{B}}t} \end{eqnarray}$ (4.8)

and also the associated friction kernel:

 $\begin{eqnarray}\label{wtizeta_t_def} \tilde\zeta(t)=\Omega\sum\limits_{j}(\tilde c^2_k/\tilde\omega_k)\cos(\tilde\omega_k t) \end{eqnarray}$ (4.9)

on the solvation mode. The Langevin equation reads

 $\begin{eqnarray}\label{Langevin} \dot{\hat p}_{\rm{B}}(t) = -\Omega \hat x_{\rm{B}}(t) - \int^t_0{\rm{d}}\tau\tilde\zeta(t-\tau)\hat p_{\rm{B}}(\tau) +\tilde F(t)\quad \end{eqnarray}$ (4.10)

Together with $\dot{\hat x}_{\rm{B}}$=$\Omega\hat p_{\rm{B}}$, we obtain $\chi_{\rm{B}}(\omega)$ in Eq.(2.4) the expression [14, 15],

 $\begin{eqnarray}\label{chiBw} \chi_{\rm{B}}(\omega) = \frac{\Omega}{\Omega^2-\omega^2 -i\omega\zeta(\omega)} \end{eqnarray}$ (4.11)

with the complex friction function,

 $\begin{eqnarray}\label{zeta_w} \zeta(\omega) \equiv \int^{\infty}_{0}\!{\rm d}t\, \text{e}^{i\omega t}\tilde\zeta(t) \end{eqnarray}$ (4.12)

It is also easy to obtain [14]

 $\begin{eqnarray}\label{wtiJw_zetaw} \tilde J(\omega) \equiv \frac{1}{2}\int^{\infty}_{-\infty}{\rm{d}} t\, \text{e}^{i\omega t}\langle[\tilde F(t), \tilde F(0)]\rangle_{\tilde{\rm{B}}} = \frac{\omega}{\Omega}{\rm Re}\, \zeta(\omega)\nonumber\\ \end{eqnarray}$ (4.13)

The first identity defines the secondary interacting bath spectral density (cf. Eq.(2.4)). Applying Eq.(4.11) for Eq.(4.13) leads to $\tilde\omega J(\omega)={\rm Im}\chi_{\rm{B}}(\omega)/|\chi_{\rm{B}}(\omega)|^2$. This is the relation between the secondary interacting bath spectral density and that of the overall one [14, 23]. The interacting secondary bath correlation function reads [cf. Eq.(2.5)]

 $\begin{eqnarray}\label{FDT2} \langle\tilde F(t)\tilde F(0)\rangle_{\tilde\beta} = \frac{1}{\pi}\! \int^{\infty}_{-\infty}\!\! {\rm d}\omega\, \frac{\text{e}^{-i\omega t}\tilde J(\omega)} {1-\text{e}^{-\beta\omega}} \end{eqnarray}$ (4.14)

Apparently, the complete characterization requires not only the $\alpha$-parameters, but also $\chi_{\rm{B}}(\omega)$ of Eq.(4.11), as inferred from Eq.(2.5).

The expression of $h'_{\rm{B}}$ is similar to Eq.(4.7), but with primed quantities. The secondary bath coordinates are given by $\tilde x'_k$=$(\tilde\omega'_k/\tilde\omega_k)^{1/2}(\tilde x_k-\tilde d_k)$, with respect to their reference $h_{\rm{B}}$-based counterparts. However, as elaborated in Ref.[21], the physical model for Eqs.(4.2) or (4.4) implies that the secondary bath environments do not contribute the quadratic coupling, i.e., $\tilde\omega'_k$=$\tilde\omega_k$, but only the linear displacements $\{\tilde d_k\}$. It also requires that the effects of secondary bath be representable with the solvation coordinate $\hat x$ only, on the basis of a set of well-defined linear displacement mapping rules [21]. This is achievable when the radio, $\tilde c'_k/\tilde c_k$, between the coupling bath strength parameter in $h'_{\rm{B}}$ and that in $h_{\rm{B}}$ of Eq.(4.7), can be effectively replaced by a single $k$-independent polarization constant [21]. This undetermined constant actually dictates the reference $h'_{\rm{B}}$-based $\zeta'(\omega)$, with respect to the $h_{\rm{B}}$-based $\zeta(\omega)$ of Eq.(4.12) with Eq.(4.9). The characterization is then carried out by using explicitly the prescription requirement, as described by Eqs. (4.2)-(4.6).

The final results in the $h_{\rm{B}}$-based prescription read [21]

 $\begin{eqnarray}\label{alpha_all} \begin{split} \alpha_0 &= \lambda \left(\frac{\Omega'}{\Omega}\right)^2 \\ \alpha_1 &= -(2\lambda\Omega)^{1/2}\left(\frac{\Omega'}{\Omega}\right)^2\\ \alpha_2 &= \frac{\Omega}{2}\left[\left(\frac{\Omega'}{\Omega}\right)^2-1\right] \end{split} \end{eqnarray}$ (4.15)

Here, $\lambda\equiv \Omega d^2/2$ denotes the reorganization energy, due to the linear displacement in Eq.(4.6), whereas $d^2/2$ amounts to the Huang-Rhys factor.

The $h'_{\rm{B}}$-based counterparts read [21]

 $\begin{eqnarray}\label{alphap_all} \begin{split} \alpha'_0 &= \lambda \\ \alpha'_1&= \left(2\lambda \frac{\Omega^2}{\Omega'}\right)^{1/2} \\ \alpha'_2 &= \frac{\Omega^2}{2\Omega'}\left[1-\left(\frac{\Omega'}{\Omega}\right)^2\right] \end{split} \end{eqnarray}$ (4.16)

Meanwhile we have also identified the aforementioned polarization constant to be of $(\tilde c'_{k}/\tilde c_{k})^2$=$\Omega/\Omega'$. It gives rise to have $\zeta'(\omega)$=$\zeta(\omega)$, the same as that of the reference $h_{\rm{B}}$-based prescription. The $h'_{\rm{B}}$-based counterpart to Eq.(4.11) is therefore [21]

 $\begin{eqnarray}\label{chiBw_e_final} \hat\chi'(\omega) = \frac{\Omega'}{\Omega^{\prime 2}-\omega^2 -i\omega\zeta(\omega)} \end{eqnarray}$ (4.17)

The required nonlinear coupling bath descriptors are now completed.

Ⅴ. EXTENDED FOKKER-PLANCK QUANTUM MASTER EQUATION A. Remarks on quantum versus semiclassical bath

It is worth reemphasizing that both the DEOM and FP-QME are "core-system" dynamics theories. While the former reveals the solvation dynamics via generalized Wick's theorem, the latter deals explicitly with the reduced core-system density operator. Here, the overall bath $h_{\rm{B}}$ of Eq.(2.2) is divided into two parts as Eq.(4.7): the solvation mode, $\Omega(\hat p^2_{\rm{B}}+\hat x^2_{\rm{B}})/2$, and the interacting secondary bath part, denoted as $h_{\rm{B}\tilde{\rm{B}}}$ below. The total composite Hamiltonian, Eq.(2.1), can be recast as

 $\begin{eqnarray}\label{HT1} H_{\rm{T}} = H_\text{core} + h_{\rm{B}\tilde{\rm{B}}} = H'_\text{core} + h_{\rm{B}} \end{eqnarray}$ (5.1)

with (noting that $H_{\rm{S}}\equiv H_0$+$\alpha_0\hat Q_{\rm{S}}$ as Eq.(2.3))

 \begin{align} H'_\text{core} &\equiv H_{\rm{S}} + \hat Q_{\rm{S}}(\alpha_1\hat x_{\rm{B}}+\alpha_2\hat x_{\rm{B}}^2)\nonumber\\ & \equiv H_{\rm{S}} + V_{\rm{SB}} \end{align} (5.2)
 \begin{align} H_\text{core} &\equiv H'_\text{core} +\frac{1}{2}\Omega(\hat p^2_{\rm{B}}+\hat x^2_{\rm{B}}) \label{Hcore_full} \end{align} (5.3)

Denote also ${\cal L}_\text{core}\hat O\equiv [H_\text{core}, \hat O]$ and ${\cal L}'_\text{core}\hat O\equiv [H'_\text{core}, \hat O]$.

An exact FP-QME for the reduced core-system density operator, denoted as $\hat\varrho(t)$ hereafter, would read

 $\begin{eqnarray}\label{FP-QME_exa} \dot{\hat\varrho}(t) =-(i{\cal L}'_\text{core}\! + {\cal L}^{\rm exa})\hat\varrho(t) \end{eqnarray}$ (5.4)

with the quantum Fokker-Planck operator of

 $\begin{eqnarray}\label{calFP_exa} {\cal L}^{\rm exa}_{\rm{FP}}\hat O \equiv \frac{i}{2}\Omega[\hat p^2_{\rm{B}}+\hat x^2_{\rm{B}}, \hat O]+{\cal R}^{\rm exa}_{\rm{B}}\hat O \end{eqnarray}$ (5.5)

One can investigate ${\cal L}^{\rm exa}$ in the absence of the primary system. This would be equivalent to the Langevin equation (Eq.(4.10)) for the Brownian dynamics of the solvation mode in the secondary bath. In general, the dissipative superoperator, ${\cal R}^{\rm exa}_{\rm{B}}$, depends on the initial condition on the total system-and-bath composite density operator [41, 42]. The resultant ${\cal R}^{\rm exa}_{\rm{B}}$ would be quite complicated and also time-dependent [41, 42]. Nevertheless, the general expression goes beyond the scope of this paper, as our main interest is the underlying quasi-particle picture. Simplification arises for the initial $\rho_{\rm{B}}(0)\simeq\rho^{\rm eq}_{\rm{B}}$ or the factorized $\rho_{\rm{T}}(0)$=$\rho_{\rm{S}}(0) \rho^{\rm eq}_{\rm{B}}$. In this case we have [43]

 $\begin{eqnarray}\label{calRexa} {\cal R}^{\rm exa}_{\rm{B}}\hat O &\simeq& \zeta \langle p^2_{\rm{B}}\rangle_{\rm{B}}\big[\hat x_{\rm{B}}, [\hat x_{\rm{B}}, \hat O]\big] +\frac{\zeta}{2}\big[\hat x_{\rm{B}}, \{\hat p_{\rm{B}}, \hat O\}\big] +\nonumber\\ &&\Omega\big(\langle p^2_{\rm{B}}\rangle_{\rm{B}}-\langle x^2_{\rm{B}}\rangle_{\rm{B}}\big) \big[\hat x_{\rm{B}}, [\hat p_{\rm{B}}, \hat O]\big] \end{eqnarray}$ (5.6)

Here, $\zeta\simeq\zeta(\omega$=0) denotes the effective friction constant; whereas $\langle x^2_{\rm{B}}\rangle_{\rm{B}}$ and $\langle p^2_{\rm{B}}\rangle_{\rm{B}}$, the thermal equilibrium phase-space variances of the solvation mode in bare bath $h_{\rm{B}}$, are evaluated via Eq.(2.5) as $\langle x^2_{\rm{B}}\rangle_{\rm{B}}$=$\langle\hat x(t$=0)$\hat x(0)\rangle_{\rm{B}}$ and $\Omega^{2}\langle p^2_{\rm{B}}\rangle_{\rm{B}}$=$\langle\dot{\hat x}(t$=0)$\dot{\hat x}(0)\rangle_{\rm{B}}$= $-\langle\ddot{\hat x}(t$=0)$\hat x(0)\rangle_{\rm{B}}$.

Some subtle issues arise in comparing the exact FP-QME, Eqs. (5.4)-(5.6), to the approximated FP-QMEs that assume a Markov secondary bath, with $\zeta(\omega)$ in Eq.(4.11) being replaced by the constant $\zeta$; i.e.,

 $\begin{eqnarray}\label{chiwBO} \chi_{\rm{B}}(\omega) = \frac{\Omega}{\Omega^2-\omega^2 -i\omega\zeta} \end{eqnarray}$ (5.7)

It is easy to show that in this case $\langle p^2_{\rm{B}}\rangle_{\rm{B}}$ diverges [14, 15]. As inferred from Eq.(4.13), the Markov secondary bath is of $\tilde J(\omega)$=$\omega\zeta/\Omega$. However, the resultant correlation funciton, $\langle\tilde F(t)\tilde F(0)\rangle_{\tilde\beta}$ of Eq.(4.14), is not a white noise, unless assuming further the classical high-temperature (HT) limit of

 $\begin{eqnarray}\label{HT_Bose} \frac{1}{1-\text{e}^{-\beta\omega}} \approx \frac{1}{2} + \frac{1}{\beta\omega} \qquad (\text{classical limit}) \end{eqnarray}$ (5.8)

In this limit,

 $\begin{eqnarray}\label{wtiFFt_HT} \langle\tilde F(t)\tilde F(0)\rangle_{\tilde\beta} \approx \frac{2\zeta}{\beta\Omega}\delta(t) + i\frac{\zeta}{\Omega}\dot\delta(t) \end{eqnarray}$ (5.9)

or

 $\begin{eqnarray}\label{variances_HT} \langle x^2_{\rm{B}}\rangle^{\rm{HT}}_{\rm{B}} = \langle p^2_{\rm{B}}\rangle^{\rm{HT}}_{\rm{B}} = \frac{1}{\beta\Omega} \end{eqnarray}$ (5.10)

In this case, Eq.(5.6) reduces to

 \begin{align}\label{calL_FP0} {\cal R}^{\rm{HT}}_{\rm{B}}\hat O = \frac{\zeta}{\beta\Omega} \big[\hat x_{\rm{B}}, [\hat x_{\rm{B}}, \hat O]\big] +i\frac{\zeta}{2}\, \big[\hat x_{\rm{B}}, \{\hat p_{\rm{B}}, \hat O\}\big] \end{align} (5.11)

It leads to Eq.(5.4) the Caldeira-Leggett's QME [27, 28],

 $\begin{eqnarray}\label{FP-QME_HT} \dot{\hat\varrho}(t) \approx -(i{\cal L}'_\text{core} + {\cal L}^{\rm{HT}}_{\rm{FP}})\hat\varrho(t) \end{eqnarray}$ (5.12)

with (cf. Eq.(5.5))

 $\begin{eqnarray}\label{calL_FP_HT} {\cal L}^{\rm{HT}}_{\rm{FP}}\hat\varrho(t) = i\frac{\Omega}{2}\big[\hat p^2_{\rm{B}}+\hat x^2_{\rm{B}}, \hat\varrho(t)\big] + {\cal R}^{\rm{HT}}_{\rm{B}} \end{eqnarray}$ (5.13)

This is a mixed quantum-classical theory, with $i{\cal L}_{\rm{core}}$ for the coherent core-system dynamics, and ${\cal R}^{\rm{HT}}_{\rm{B}}$ for the classical dissipative solvation dynamics in the secondary bath environment. However, as a classical bath theory, Eq.(5.12) is of a very limited applicability range in the parameters space.

It is anticipated that, as long as the coherent $i{\cal L}_\text{core}$ and the incoherent ${\cal R}_{\rm{B}}$ are not treated at the same level, a certain caustic compensation would be required. This would give rise to an additional dissipative superoperator, $\delta{\cal R}_{\rm{SB}}$, that depends explicitly on both the dissipative system and solvation bath modes. The resultant semiclassical FP-QME would assume the form of (see Section Ⅴ.B)

 $\begin{eqnarray}\label{FP-QME0} \dot{\hat\varrho}(t) = -(i{\cal L}'_\text{core} + {\cal L}_{\rm{FP}} + \delta{\cal R}_{\rm{SB}})\hat\varrho(t) \end{eqnarray}$ (5.14)

The high-temperature limit, Eq.(5.12), is associated with ${\cal L}^{\rm{HT}}_{\rm{FP}}$ of Eq.(5.13) and $\delta{\cal R}^{\rm{HT}}_{\rm{SB}}$=0.

It would be insightful to have a close inspection on the incoherent $\delta{\cal R}_{\rm{SB}}$ caustics versus the coherent interacting component of ${\cal L}'_{\rm{core}}$ that reads (cf. Eq.(5.2))

 $\begin{eqnarray}\label{calLcore} {\cal L}'_{\rm{core}} = {\cal L}_{\rm{S}} +\alpha_1 {\cal V}^{(1)}_{\rm{SB}} + \alpha_2 {\cal V}^{(2)}_{\rm{SB}} \end{eqnarray}$ (5.15)

Here,

 $\begin{eqnarray}\label{calV_SB} \begin{split} {\cal V}^{(1)}_{\rm{SB}}\hat O \equiv [\hat Q_{\rm{S}}\hat x_{\rm{B}}, \hat O]\\ {\cal V}^{(2)}_{\rm{SB}}\hat O \equiv [\hat Q_{\rm{S}}\hat x^2_{\rm{B}}, \hat O] \end{split} \end{eqnarray}$ (5.16)

To facilitate the later construction, we introduce

 $\begin{eqnarray}\label{oplusminus_def} \begin{split} x^\oplus_{\rm{B}}\hat O \equiv \frac{1}{2}\{\hat x_{\rm{B}}, \hat O\}\\ x^\ominus_{\rm{B}}\hat O \equiv -i[\hat x_{\rm{B}}, \hat O] \end{split} \end{eqnarray}$ (5.17)

and similarly $p^\oplus_{\rm{B}}$ and $p^\ominus_{\rm{B}}$, such that

 $\begin{eqnarray}\label{oplusminus} \frac{1}{2}\big[\hat x_{\rm{B}}, \{\hat p_{\rm{B}}, \hat O\}\big] \equiv ix^\ominus_{\rm{B}}p^\oplus_{\rm{B}}\hat O \end{eqnarray}$ (5.18)

We can then express Eq.(5.16) as

 $\begin{eqnarray} \begin{split} {\cal V}^{(1)}_{\rm{SB}} \hat O &=\big[\hat Q_{\rm{S}}, x^\oplus_{\rm{B}}\hat O\big] +\frac{i}{2}\big\{\hat Q_{\rm{S}}, x^\ominus_{\rm{B}}\hat O \big\} \\ {\cal V}^{(2)}_{\rm{SB}} \hat O &= \big[\hat Q_{\rm{S}}, (x^{\oplus 2}_{\rm{B}}-\frac{1}{4}x^{\ominus 2}_{\rm{B}})\hat O \big] +i\big\{\hat Q_{\rm{S}}, x^\oplus_{\rm{B}} x^\ominus_{\rm{B}}\hat O \big\} \end{split} \end{eqnarray}$ (5.19)

The caustic dissipation superoperator assumes also the form of $\delta{\cal R}_{\rm{SB}}$=$\alpha_1 \delta{\cal R}^{(1)}_{\rm{SB}} + \alpha_2 \delta{\cal R}^{(2)}_{\rm{SB}}$. It will be evident later that the final expressions (cf. Eqs. (5.30)-(5.32)) can be conveniently expressed in terms of

 $\begin{eqnarray} \begin{split} i{\cal V}^{(1)}_{\rm eff} \equiv i{\cal V}^{(1)}_{\rm{SB}} + \delta{\cal R}^{(1)}_{\rm{SB}}\\ i{\cal V}^{(2)}_{\rm eff} \equiv i{\cal V}^{(2)}_{\rm{SB}} + \delta{\cal R}^{(2)}_{\rm{SB}} \end{split} \end{eqnarray}$ (5.20)

Recast now Eq.(5.14), with Eq.(5.15), as

 $\begin{eqnarray}\label{FP-QME} \dot{\hat\varrho}(t) = -\big(i{\cal L}^{\rm eff}_{\rm{S}} + {\cal L}_{\rm{FP}} + i\alpha_1{\cal V}^{(1)}_{\rm eff}+i\alpha_2{\cal V}^{(2)}_{\rm eff}\big) \hat\varrho(t)\quad \end{eqnarray}$ (5.21)

It is worth to mention that we use ${\cal L}^{\rm eff}_{\rm{S}}$ here, instead of the presumed ${\cal L}_{\rm{S}}$, as it would be if the exact ${\cal L}^{\rm exa}_{\rm{FP}}$ or ${\cal R}^{\rm exa}_{\rm{B}}$, Eq.(5.6), were used. Actually, the adopted ${\cal L}_{\rm{FP}}$ below (cf. Eq.(5.22)) resembles formally as the high-temperature one, Eq.(5.13), with ${\cal R}^\text{HT}_{\rm{B}}$ of Eq.(5.11). It does not have the $\big[\hat x_{\rm{B}}, [\hat p_{\rm{B}}, \hat O]\big]$-term, as in Eq.(5.6). A proper ${\cal L}^{\rm eff}_{\rm{S}}$ or $H^{\rm eff}_{\rm{S}}$ may partially makeup this loss. We adopt $H^{\rm eff}_{\rm{S}}$ of Eq.(2.3), in line with the DEOM (Eq.(3.9)).

B. FP-QME with semiclassical bath

We will be interested in a semiclassical bath, having [30]

 $\begin{eqnarray}\label{calL_FP} {\cal L}_{\rm{FP}}\hat\varrho(t) &=& i\frac{\bar\Omega}{2}\big[\hat p^2_{\rm{B}}+\hat x^2_{\rm{B}}, \hat\varrho(t)\big] +\frac{\bar\zeta}{\beta\bar\Omega} \big[\hat x_{\rm{B}}, [\hat x_{\rm{B}}, \hat\varrho(t)]\big] +\nonumber\\&&i\frac{\bar\zeta}{2}\, \big[\hat x_{\rm{B}}, \{\hat p_{\rm{B}}, \hat\varrho(t)\}\big] \end{eqnarray}$ (5.22)

Formally it resembles Eq.(5.13) with Eq.(5.11); thus the well-established Fokker-Planck operator algebra [38] is directly applicable; see Section Ⅵ.A. Nevertheless, it replaces the solvation frequency and friction constant in Eq.(5.13) with temperature-dependent effective $\bar\Omega$ and $\bar\zeta$.

More specifically, ${\cal L}_{\rm{FP}}$ of Eq.(5.22) goes by a bi-dissipaton (exponential) approximant of the solvation correlation function, Eq.(2.7),

 $\begin{eqnarray}\label{xB_corr_FP_2} \langle\hat x_{\rm{B}}(t)\hat x_{\rm{B}}(0)\rangle_{\rm{B}} \approx \eta_{+} \text{e}^{-\gamma_{+}t}+\eta_{-} \text{e}^{-\gamma_{-}t} \end{eqnarray}$ (5.23)

where

 $\begin{eqnarray}\label{gam_pm} \gamma_{\pm}\equiv \frac{1}{2}\big[\bar\zeta\pm (\bar\zeta^2-4\bar\Omega^2)^{1/2}\big] \end{eqnarray}$ (5.24)

While $\bar\Omega$ and $\bar\zeta$ are positive, $\eta_+$ and $\eta_-$ in Eq.(5.23) are complex. Parametrization goes with the so-called minimum-dissipaton ansatz [29]:

 $\begin{split} &\eta_{+} + \eta_{-} = \langle\hat x^2_{\rm{B}}\rangle_{\rm{B}} \\ &\frac{\eta_{+}}{\gamma_{+}} + \frac{\eta_{-}}{\gamma_{-}} = \hat C_{\rm{B}}(\omega=0) \\ &\frac{\eta_{+}}{\gamma^2_{+}} + \frac{\eta_{-}}{\gamma^2_{-}} = -i\frac{{\rm{d}}\hat C_{\rm{B}}(\omega=0)}{{\rm{d}}\omega} \end{split}$ (5.25)

with the exact values via Eq.(2.5) on the left-hand-sides, where

 $\begin{eqnarray}\label{hatCwB_def} \hat C_{\rm{B}}(\omega) \equiv \int^{\infty}_0\!\!{\rm{d}} t\, \text{e}^{i\omega t} \langle\hat x_{\rm{B}}(t)\hat x_{\rm{B}}(0)\rangle_{\rm{B}} \end{eqnarray}$ (5.26)

Remarkably, for the Markovian Brownian osciallator of Eq.(5.7), one can analytically solve Eq.(5.25) and further identify its impressive accuracy range in the parameter space [29, 30]. One interesting result is that

 $\begin{eqnarray}\label{ratio_bar} \frac{\bar\zeta}{\zeta} = \frac{\bar\Omega^2}{\Omega^2} \leq 1 \end{eqnarray}$ (5.27)

which can be recast in terms of $\bar r_{\rm{BO}}\equiv \bar\zeta/(2\bar\Omega)$ versus $r_{\rm{BO}}\equiv \zeta/(2\Omega)$ as $\bar\Omega/\Omega$=$\bar r_{\rm{BO}}/r_{\rm{BO}}\leq$1. The lower the temperature the smaller the ratio is, and the equal sign holds in the HT limit of $T\rightarrow\infty$ only [29].

Introduce three real parameters, $\delta_0$, $\delta_1$ and $\delta_2$, via [29, 30]

 $\begin{eqnarray}\label{delta12_def} \begin{split} \eta_+ + \eta_- &\equiv \frac{1}{\beta\bar\Omega} + \delta_0 \, \\ \eta_+\gamma_+ + \eta_-\gamma_- &\equiv \bar\Omega\Big[\delta_1 +\frac{i}{2}(1+\delta_2)\Big]\, \end{split} \end{eqnarray}$ (5.28)

Together with Eqs.(5.25) and (5.27), we obtain

 $\begin{eqnarray}\label{MDA_all_delta0} \begin{split} \delta_0 &= \langle\hat x^2_{\rm{B}}\rangle_{\rm{B}} - \frac{1}{\beta{\bar\Omega}}\\ \delta_1 &= \frac{\bar\zeta}{\bar\Omega}\delta_0\\ \delta_2 &= \frac{\bar\Omega}{\Omega}-1 \end{split} \end{eqnarray}$ (5.29)

Physically, $\delta_0$ and $\delta_1$ are related to the short-time and long-time caustics, respectively, whereas $\delta_2$ measures the caustic renormalization. These caustic parameters all vanish at the high-temperature limit.

The above semiclassical parametrization scheme was originally proposed in Ref.[29]. It has also been used together with Eq.(5.22) in the construction of FP-QME (5.21) via its DEOM correspondence, in the absence of quadratic coupling bath [30].

In this work, we apply the same construction method on Eq.(5.21) with the linear-plus-quadratic coupling bath. The final results read

 $\begin{eqnarray}\label{FP_QME_termI} {\cal V}^{(1)}_{\rm eff}\hat O =[\hat Q_{\rm{S}}, {\cal W}_{\rm{B}} \hat O] +\frac{i}{2} \{\hat Q_{\rm{S}}, {\cal U}_{\rm{B}}\hat O \} \end{eqnarray}$ (5.30)
 $\begin{eqnarray}\label{FP_QME_termII} {\cal V}^{(2)}_{\rm eff}\hat O =\Big[\hat Q_{\rm{S}}, \Big({\cal W}^2_{\rm{B}}-\frac{1}{4}{\cal U}^2_{\rm{B}}\Big)\hat O\Big] +i \big\{\hat Q_{\rm{S}}, {\cal W}_{\rm{B}}{\cal U}_{\rm{B}}\hat O\big\}\nonumber\\ \end{eqnarray}$ (5.31)

where

 $\begin{eqnarray}\label{calWV_def} \begin{split} {\cal W}_{\rm{B}}\hat O &\equiv \big( x^\oplus_{\rm{B}}+\delta_0 p^\ominus_{\rm{B}} +\delta_1 x^\ominus_{\rm{B}}\big)\hat O \\ {\cal U}_{\rm{B}}\hat O &\equiv (1+\delta_2)x^\ominus_{\rm{B}}\hat O \end{split} \end{eqnarray}$ (5.32)
Ⅵ. FOKKER-PLANCK OPERATOR ALGEBRA VERSUS DISSIPATONS DYNAMICS A. Fokker-Planck operator algebra

In this section, we derive the above FP-QME formalism, via its DEOM correspondence. The latter is just the bi-dissipaton version of Eq.(3.9); see Section Ⅵ.B for the details. Exploited will also be the Fokker-Planck operator algebra [38]. The underlying quasi-particle picture will be elaborated in relation to that in Section Ⅲ.

The derivations will be carried out in the Wigner representation for the solvation mode. While reduced core-system $\hat\varrho_{\rm{W}}(t)\equiv \hat\varrho(x_{\rm{B}}, p_{\rm{B}};t)$ remains an operator in the system subspace, the superoperators ${\cal W}_{\rm{B}}$ and ${\cal U}_{\rm{B}}$ in Eq.(5.32) become solvation phase-space operators; i.e.,

 $\begin{eqnarray}\label{hatWU} \begin{split} \hat W_{\rm{B}} &\equiv [{\cal W}_{\rm{B}}]_{\rm{W}} = x_{\rm{B}} - \delta_0 \partial_x + \delta_1 \partial_p \, \\ \hat U_{\rm{B}} &\equiv [{\cal U}_{\rm{B}}]_{\rm{W}} = (1+\delta_2)\partial_p \, \end{split} \end{eqnarray}$ (6.1)

Here

 $\begin{eqnarray}\label{app_partial_notation} \partial_x \equiv \frac{\partial}{\partial x_{\rm{B}}}, \ \ \ \partial_p \equiv \frac{\partial}{\partial p_{\rm{B}}} \end{eqnarray}$ (6.2)

Meanwhile, ${\cal L}_{\rm{FP}}$ (Eq.(5.22)) assumes the FP operator,

 $\begin{eqnarray}\label{calL_FP_Wigner} \hat L_{\rm{FP}} &\equiv & [{\cal L}_{\rm{FP}}]_{\rm{W}}\nonumber\\ &=& \bar\Omega (\partial_{x} p_{\rm{B}} - \partial_{p} x_{\rm{B}}) - \frac{\bar\zeta}{\beta\bar\Omega}\partial^2_{pp} -\bar\zeta \partial_{p} p_{\rm{B}} \end{eqnarray}$ (6.3)

Its eigen solutions are analyzed via the similarity transformation, as follows [38].

 $\begin{eqnarray}\label{tiL_FP_ab0}\begin{split} \tilde L_{\rm{FP}} &\equiv \text{e}^{S}\hat L_{\rm{FP}}\text{e}^{-S} \ \ \\ S&=\frac{1}{4}\beta\bar\Omega(x^2_{\rm{B}}+p^2_{\rm{B}}) \end{split} \end{eqnarray}$ (6.4)

Under this transformation,

 $\begin{eqnarray}\label{ti_partial_x_p} \begin{split} \tilde\partial_{x} &= \text{e}^{S}\partial_x \text{e}^{-S} = -\frac{\beta\bar\Omega}{2} x_{\rm{B}}+\partial_{x} \equiv -\sqrt{\beta\bar\Omega}\, \hat a^{†}_{x} \\ \tilde\partial_{p} &= \text{e}^{S}\partial_x \text{e}^{-S} = -\frac{\beta\bar\Omega}{2} p_{\rm{B}}+\partial_{p} \equiv -\sqrt{\beta\bar\Omega}\, \hat a^{†}_{p} \end{split} \end{eqnarray}$ (6.5)

where

 $\begin{eqnarray} \qquad \begin{matrix} \hat a^{†}_{\mu} = \dfrac{\sqrt{\beta\bar\Omega}}{2} \mu_{\rm{B}} - \dfrac{1}{\sqrt{\beta\bar\Omega}} \partial_{\mu} \\ \hat a_{\mu} = \dfrac{\sqrt{\beta\bar\Omega}}{2} \mu_{\rm{B}} +\dfrac{1}{\sqrt{\beta\bar\Omega}} \partial_{\mu} \end{matrix} \ \ \, \text{with} \ \ \mu = x, p \end{eqnarray}$ (6.6)

We have

 $\begin{eqnarray}\label{ti_x_p} \begin{split} \tilde x_{\rm{B}} &= x_{\rm{B}} =\frac{1}{\sqrt{\beta\bar\Omega}}(\hat a_{x}+\hat a^{†}_{x}) \\ \tilde p_{\rm{B}} &=p_{\rm{B}} =\frac{1}{\sqrt{\beta\bar\Omega}}(\hat a_{p}+\hat a^{†}_{p}) \end{split} \end{eqnarray}$ (6.7)

Apparently, $[\hat a_{\mu}, \hat a^{†}_{\nu}]$=$\delta_{\mu\nu}$. It together with Eqs.(6.5) and (6.7) leads to Eq.(6.3) the $S$-transformed expression,

 $\begin{eqnarray}\label{app_tiL_FP} \tilde L_{\rm{FP}} =\bar \Omega(\hat a^†_{p}\hat a_{x}-\hat a^†_{x}\hat a_{p}) +\bar\zeta\hat a^†_{p}\hat a_{p} \end{eqnarray}$ (6.8)

Note that $\bar\Omega^2$=$\gamma_{+}\gamma_{-}$ and $\bar\zeta$=$\gamma_{+}\gamma_{-}$ (cf. Eq.(5.24)). Let

 $\begin{eqnarray}\label{app_rj_def} \begin{split} r_{1} \equiv \sqrt{\gamma_{+}/(\gamma_{+}-\gamma_{-})}\, \\ r_{2} \equiv \sqrt{\gamma_{-}/(\gamma_{+}-\gamma_{-})}\, \end{split} \end{eqnarray}$ (6.9)

Introduce now the operators, $\hat b^{\pm}_{j}; j$=1, 2, via

 $\hat a^{†}_{x} =r_{2}\hat b^{+}_1+r_{1}\hat b^{+}_2, \ \ \ \ \, \hat a^{†}_{p} = r_{1}\hat b^{+}_1 + r_{2}\hat b^{+}_2$ (6.10)
 $\hat a_{x} = -r_{2}\hat b^{-}_1+r_{1}\hat b^{-}_2 , \ \ \, \hat a_{p} = r_{1}\hat b^{-}_1 - r_{2}\hat b^{-}_2$ (6.11)

These are equivalent to

 $\begin{eqnarray}\label{c_pm_def} \begin{split} %\label{c_m_def} \hat b^{-}_1&\equiv r_{2}\hat a_{x}+r_{1}\hat a_{p}, \quad\ \ \hat b^{-}_2 \equiv r_{1}\hat a_{x}+r_{2}\hat a_{p} \\ \hat b^{+}_1&\equiv -r_{2}\hat a^{†}_{x}+r_{1}\hat a^{†}_{p} , \quad \hat b^{+}_2 \equiv r_{1}\hat a^{†}_{x}-r_{2}\hat a^{†}_{p} \end{split} \end{eqnarray}$ (6.12)

Despite of $\hat b^{\pm}_{j}\neq (\hat b^{\mp}_{j})^{†}$, they satisfy

 $\begin{eqnarray}\label{c_pm_commutators} [\hat b^{-}_{j}, \hat b^{+}_{k}]=\delta_{jk} ~\text{ and }~ [\hat b^{\pm}_{j}, \hat b^{\pm}_{k}]=0 \end{eqnarray}$ (6.13)

These lead to Eq.(6.8) the diagonalized form of

 $\begin{eqnarray}\label{tiL_FP_c} \tilde L_{\rm{FP}} =\gamma_{+}\hat b^{+}_{1}\hat b^{-}_{1}+\gamma_{-}\hat b^{+}_{2}\hat b^{-}_{2} \end{eqnarray}$ (6.14)

Due to the boson-like relations of Eq.(6.13), the standard harmonic oscillator algebra is applicable. We obtain

 $\hat b^{+}_j\psi_{n_j}\!(x_{\rm{B}}, p_{\rm{B}})=\sqrt{n_j+1}\, \psi_{n_j+1}(x_{\rm{B}}, p_{\rm{B}})$ (6.15)
 $\hat b^{-}_j\psi_{n_j}\!(x_{\rm{B}}, p_{\rm{B}})=\sqrt{n_j}\, \psi_{n_j-1}(x_{\rm{B}}, p_{\rm{B}})$ (6.16)

where $n_j$=0, 1, $\cdots$, with the normalized ground state,

 $\begin{eqnarray}\label{psi0_xpB} \psi_{0}(x_{\rm{B}}, p_{\rm{B}})=\Big(\frac{\beta\bar\Omega}{2\pi}\Big)^{1/2} \text{e}^{-\beta\bar\Omega/4(x^2_{\rm{B}}+p^2_{\rm{B}})} \end{eqnarray}$ (6.17)

Further, from $(\hat b^{+}_{j}\hat b^{-}_{j})\psi_{n_j}(x_{\rm{B}}, p_{\rm{B}})$=$n_j\psi_{n_j}(x_{\rm{B}}, p_{\rm{B}})$, we obtain Eq.(6.14) the eigen solutions,

 $\begin{eqnarray}\label{eig_tiL} \tilde L_{\rm{FP}}\Psi_{n_1 n_2} = (n_1\gamma_{+}+n_2\gamma_{-}) \Psi_{n_1 n_2} \end{eqnarray}$ (6.18)

where

 $\begin{eqnarray}\label{Psi_n1n2} \Psi_{n_1 n_2}(x_{\rm{B}}, p_{\rm{B}})\equiv \psi_{n_1}(x_{\rm{B}}, p_{\rm{B}})\psi_{n_2}(x_{\rm{B}}, p_{\rm{B}}) \end{eqnarray}$ (6.19)

The physical implications on the above diagonalization will be discussed at the end of this subsection.

Introduce for the later use the scaled (unnormalized) eigenfunctions [30],

 $\begin{eqnarray}\label{wtiPsi_def} \tilde\Psi_{n_1n_2}(x_{\rm{B}}, p_{\rm{B}}) = \frac{(-1)^{n_1}}{\sqrt{n_1!\, n_2!}} \frac{\big(\sqrt{\beta\bar\Omega}\big)^{n_1+n_2}}{r^{n_1}_{2}r^{n_2}_{1} } \Psi_{n_1n_2}\nonumber\\ \end{eqnarray}$ (6.20)

It is evident below that the phase-space actions on the scaled wavefunctions assume rather compact expressions. For $\tilde\partial_x$ and $\tilde\partial_p$ (Eq.(6.5)), applying Eq.(6.10) and then Eq.(6.15), we obtain

 $\tilde\partial_x\tilde\Psi_{n_1n_2} = r^2_2(n_1+1)\tilde\Psi_{n_1 +1, n_2} \! -r^2_1(n_2+1)\tilde\Psi_{n_1, n_2+1}$ (6.21)
 $\tilde\partial_p\tilde\Psi_{n_1n_2} =r_1r_2[(n_1 + 1)\tilde\Psi_{n_1+1, n_2} \! -(n_2+1)\tilde\Psi_{n_1, n_2+1}]$ (6.22)

Moreover, as the actions on the scaled eigenfunction are concerned,

 $\begin{eqnarray}\label{app_p+partialp} \tilde p_{\rm{B}}+\frac{1}{\beta\bar\Omega}\tilde\partial_p = -r_{1}r_2\Big(\tilde x_{\rm{B}}+\frac{1}{\beta\bar\Omega}\tilde\partial_x\Big) \end{eqnarray}$ (6.23)

where

 $\begin{eqnarray}\label{app_x+partialx} \Big(\tilde x_{\rm{B}}+\frac{1}{\beta\bar\Omega}\tilde\partial_x\Big)\tilde\Psi_{n_1n_2} =\tilde\Psi_{n_1-1, n_2}+\tilde\Psi_{n_1, n_2-1} \end{eqnarray}$ (6.24)

To conclude this subsection, we like to point out the implications on the above FP operator algebra. First of all, the eigen Eq.(6.18) resembles the generalized diffusion Eq.(3.2). It is noticed that the FP operator itself arises from the bare bath Hamiltonian $h_{\rm{B}}$; see Eq.(4.7) and also the comments following Eq.(5.5). Moreover, one can view Eq.(6.15) and Eq.(6.16) as the manifestation of the generalized Wick's theorem, Eq.(3.4) the GWT-1. The effect of quadratic coupling bath on the eigenfunctions will resemble the GWT-2 Eq.(3.6). The above observations suggest the quasi-particle or dissipaton picture in the FP-QME. Nevertheless, the dissipaton approach provides much simpler algebra also clearer physical picture.

B. The corresponding dissipatons dynamics

To facility the later construction of FP-QME Eq.(5.21) via its DEOM correspondence, we present explicitly the bi-dissipaton version of Eq.(3.9), as follows.

 $\begin{eqnarray}\label{DEOM2} \dot\rho_{n_1n_2} \! & =& - \big(i{\cal L}_{\rm eff} + n_1\gamma_{+} + n_2\gamma_{-}\big)\rho_{n_1n_2}-\nonumber\\ &&i\alpha_2 \big({\cal B}_{20}\rho^{\{n-2\}}_{20} +{\cal B}_{02}\rho^{\{n-2\}}_{02} +2{\cal B}_{11}\rho^{\{n-2\}}_{11} \big)-\nonumber\\ &&i\alpha_2\big({\cal A}\rho^{\{n+2\}}_{n_1n_2} +2{\cal C}_1\rho^{\{n\}}_{1}+2{\cal C}_2\rho^{\{n\}}_{2}\big)-\nonumber\\ &&i\alpha_1\big({\cal A}\rho^{\{n+1\}}_{n_1n_2} +{\cal C}_1\rho^{\{n-1\}}_{1}+{\cal C}_2\rho^{\{n-1\}}_{2}\big) \end{eqnarray}$ (6.25)

The linear coupling bath related quantities are

 $\begin{eqnarray}\label{varrho+1} \rho^{\{n+1\}}_{n_1n_2}(t) \equiv \rho_{n_1+1, n_2}(t) +\rho_{n_1, n_2+1}(t) \end{eqnarray}$ (6.26)

and

 $\begin{eqnarray}\label{varrho-1} \begin{split} &\rho^{\{n-1\}}_{1}(t) \equiv n_1\rho_{n_1-1, n_2}(t) \\ &\rho^{\{n-1\}}_{2}(t) \equiv n_2\rho_{n_1, n_2-1}(t) \end{split} \end{eqnarray}$ (6.27)

 $\begin{eqnarray}\label{varrho+2} \rho^{\{n+2\}}_{n_1n_2}(t) &\equiv & \rho_{n_1+2, n_2}(t) +\rho_{n_1, n_2+2}(t)+\nonumber\\&&2\rho_{n_1+1, n_2+1}(t) \end{eqnarray}$ (6.28)
 $\begin{eqnarray}\label{varrho+0} \begin{split} \rho^{\{n\}}_1(t) &\equiv n_1[\rho_{n_1-1, n_2+1}(t)+\rho_{n_1n_2}(t)]\\ \rho^{\{n\}}_2(t) &\equiv n_2[\rho_{n_1+1, n_2-1}(t)+\rho_{n_1n_2}(t)] \end{split} \end{eqnarray}$ (6.29)

and

 $\begin{eqnarray}\label{varrho-2} \begin{split} \rho^{\{n-2\}}_{20}(t) &\equiv n_1(n_1-1) \rho_{n_1-2, n_2}(t) \\ \rho^{\{n-2\}}_{02}(t) &\equiv n_2(n_2-1) \rho_{n_1, n_2-2}(t) \\ \rho^{\{n-2\}}_{11}(t) &\equiv n_1 n_2\rho_{n_1-1, n_2-1}(t) \end{split} \end{eqnarray}$ (6.30)

The superoperators involved in Eq.(6.25) follow those of Eq.(3.10); i.e., ${\cal A}\hat O\equiv [\hat Q_{\rm{S}}, \hat O]$

 $\begin{eqnarray}\label{calC_def} \begin{split} {\cal C}_1 \hat O &\equiv \eta_{+}\hat Q_{\rm{S}}\hat O - \bar\eta^{\ast}_{+} \hat O \hat Q_{\rm{S}} \\ {\cal C}_2 \hat O &\equiv \eta_{-}\hat Q_{\rm{S}}\hat O - \bar\eta^{\ast}_{-} \hat O \hat Q_{\rm{S}} \end{split} \end{eqnarray}$ (6.31)

and

 $\begin{eqnarray}\label{calB_def} \begin{split} {\cal B}_{20} &\equiv \eta^2_{+}\hat Q_{\rm{S}}\hat O - (\bar\eta^{\ast}_{+})^2 \hat O \hat Q_{\rm{S}} \\ {\cal B}_{02} &\equiv \eta^2_{-}\hat Q_{\rm{S}}\hat O - (\bar\eta^{\ast}_{-})^2 \hat O \hat Q_{\rm{S}} \\ {\cal B}_{11} &\equiv \eta_{+}\eta_{-}\hat Q_{\rm{S}}\hat O - \bar\eta^{\ast}_{+}\bar\eta^{\ast}_2 \hat O \hat Q_{\rm{S}} \end{split} \end{eqnarray}$ (6.32)

Note that $\bar\eta_{\pm}$ amounts to the pre-exponential coefficient $\eta_{\bar k}$ of Eq.(2.8). As specified there, it can now be recast as

 $\begin{eqnarray}\label{bar_eta_pm} \bar\eta_{\pm} \equiv \begin{cases} \eta_{\pm};& \text{if$\gamma_{\pm}$is real} \\ \eta_{\mp};&\text{if$\gamma_+ = \gamma^{\ast}_{-}$} \end{cases} \end{eqnarray}$ (6.33)

Moreover, as inferred from Eq.(5.23) and the time-reversal relation, Eq.(2.6), we have

 $\begin{eqnarray}\label{time_rev} \bar\eta^{\ast}_{+}\text{e}^{-\gamma_{+} t}+ \bar\eta^{\ast}_{-}\text{e}^{-\gamma_{-} t} =(\eta_{+}\text{e}^{-\gamma_{+} t}+\eta_{-}\text{e}^{-\gamma_{-} t})^{\ast} \end{eqnarray}$ (6.34)

Taking the time derivative at $t$=0, we obtain

 $\begin{eqnarray}\label{eqB1} \bar\eta^{\ast}_{+}\gamma_{+} + \bar\eta^{\ast}_{-}\gamma_{-} =(\eta_{+}\gamma_{+} + \eta_{-}\gamma_{-})^{\ast} \end{eqnarray}$ (6.35)

Together with $r_1/r_2$=$\gamma_+/\bar\Omega$ and $r_2/r_1=\gamma_-/\bar\Omega$ (cf. Eqs. (5.24) and (6.9); see also Eq.(5.28)), we have

 $\begin{eqnarray}\label{eqB2} \bar\eta^{\ast}_{+}\frac{r_1}{r_2} + \bar\eta^{\ast}_{-}\frac{r_2}{r_1} &=& \Big(\eta_{+}\dfrac{r_{1}}{r_{2}} + \eta_{-}\dfrac{r_{2}}{r_{1}}\Big)^{\ast}\nonumber\\ &=&\delta_1-\frac{i}{2}(1+\delta_2) \end{eqnarray}$ (6.36)

Together with Eq.(5.28), the above identities relate the complex $\eta$-parameters, in Eq.(6.31) and Eq.(6.32) of DEOM, to the real $\delta$-parameters in Eq.(5.32) of FP-QME.

C. The FP-QME space verse the DEOM space

We can now complete the construction of FP-QME (5.21), with Eq.(5.30)-Eq.(5.32) based on its DEOM counterpart, Eq.(6.25)-Eq.(6.32). Let $\tilde\varrho(x_{\rm{B}}, p_{\rm{B}}, t)\equiv e^S\hat\varrho(x_{\rm{B}}, p_{\rm{B}}, t)$, the same similarity transformation used in the LP operator diagonalization; cf. Eq.(6.4). The construction starts with

 \begin{align}\label{app_rhoW_in_DDO1} \tilde\varrho(x_{\rm{B}}, p_{\rm{B}}, t) &= \sum^{\infty}_{n_1, n_2} \rho_{n_1n_2}(t)\tilde\Psi_{n_1n_2}(x_{\rm{B}}, p_{\rm{B}}) \end{align} (6.37)

This variables-separation expression specifies the fact that the FP-QME dynamics is identical to its DEOM counterpart DEOM2. The expansion basis set, $\{\tilde\Psi_{n_1n_2}(x_{\rm{B}}, p_{\rm{B}})\}$, specified in Eq.(6.20), are those properly scaled eigenfunctions of the LP operators, as we did before for the linear coupling case [30].

In the following derivations, we adopt the notation of correspondence, "$\Leftrightarrow$", in which Eq.(6.37) reads

 \begin{align}\label{app_rhoW_in_DDO2} \tilde\varrho(x_{\rm{B}}, p_{\rm{B}}, t) \Leftrightarrow \rho_{n_1n_2}(t) \end{align} (6.38)

Apparently, it would also lead to the operator-level correspondence of $\hat\varrho(t)\Leftrightarrow\rho_{n_1n_2}(t)$.

The phase-space actions on $\tilde\varrho(x_{\rm{B}}, p_{\rm{B}}, t)$, as involved in Eq.(6.1), can now be carried out easily, by using Eq.(6.21)-Eq.(6.24). For example, Eq.(6.24) results in

 $\begin{eqnarray} \Big(\tilde x_{\rm{B}}+\frac{1}{\beta\bar\Omega}\tilde\partial_x\Big)\tilde\varrho &=& \sum\limits_{n_1, n_2}\rho_{n_1n_2}(\tilde\Psi_{n_1-1, n_2}+\tilde\Psi_{n_1, n_2-1})\nonumber\\ &=&\sum\limits_{n_1, n_2}(\rho_{n_1+1, n_2}+\rho_{n_1, n_2+1})\tilde\Psi_{n_1n_2}\nonumber\\ \end{eqnarray}$ (6.39)

Following the notation that defines Eq.(6.38) for Eq.(6.37), the last expression above reads

 $\begin{eqnarray}\label{app_x+partialx_DDO} \Big(\tilde x_{\rm{B}}+\frac{1}{\beta\bar\Omega}\tilde\partial_x\Big) \tilde\varrho \Leftrightarrow \rho_{n_1+1, n_2}+\rho_{n_1, n_2+1} \end{eqnarray}$ (6.40)

Similarly, we have

 $\begin{eqnarray}\label{app_partialxandp_DDO} \begin{split} &\tilde\partial_{x}\tilde\varrho \Leftrightarrow r^2_2 \rho^{\{n-1\}}_{1}- r^2_1 \rho^{\{n-1\}}_{2} \\ &\tilde\partial_{p}\tilde\varrho \Leftrightarrow r_1r_2 \big(\rho^{\{n-1\}}_{1}-\rho^{\{n-1\}}_{2}\big) \end{split} \end{eqnarray}$ (6.41)

The involved $\rho^{\{n-1\}}_{1}$ and $\rho^{\{n-1\}}_{2}$ were defined in Eq.(6.27). The solutions here are

 $\begin{eqnarray}\label{app_2} \begin{split} &\rho^{\{n-1\}}_1\Leftrightarrow \Big( \frac{r_1}{r_2}\tilde\partial_{p}-\tilde\partial_{x}\Big) \tilde\varrho \\ &\rho^{\{n-1\}}_2\Leftrightarrow \Big( \frac{r_2}{r_1}\tilde\partial_{p}-\tilde\partial_{x} \Big)\tilde\varrho \end{split} \end{eqnarray}$ (6.42)

Moreover, by using Eq.(6.26), we recast Eq.(6.40) as

 $\begin{eqnarray}\label{app_1} \rho^{\{n+1\}}_{n_1n_2} \Leftrightarrow \Big(\tilde x_{\rm{B}}+\frac{1}{\beta\bar\Omega}\tilde\partial_x\Big)\tilde\varrho \end{eqnarray}$ (6.43)

Now comparing the linear coupling bath contribution, the $\alpha_1$-term in Eq.(6.25), with that in Eq.(5.21), we have

 $\begin{eqnarray}\label{vs1} {\cal V}^{(1)}_{\rm eff}\hat\varrho \Leftrightarrow {\cal A}\rho^{\{n+1\}}_{n_1n_2} +{\cal C}_1\rho^{\{n-1\}}_{1}+{\cal C}_2\rho^{\{n-1\}}_{2} \end{eqnarray}$ (6.44)

Applying ${\cal A}\hat O\equiv [\hat Q_{\rm{S}}, \hat O]$ and Eq.(6.31) for ${\cal C}$-actions on the specified operators, with Eq.(6.42) and Eq.(6.43), we obtain

 $\begin{eqnarray}\label{vs1_deom} {\cal V}^{(1)}_{\rm eff}\hat\varrho \Leftrightarrow \Big(\tilde x_{\rm{B}} -\frac{1}{\beta\bar\Omega}\tilde\partial_{x} + \tilde\Gamma\Big) \hat Q_{\rm{S}} \tilde\varrho-\nonumber\\ \Big(\tilde x_{\rm{B}} -\frac{1}{\beta\bar\Omega}\tilde\partial_{x} +\widehat{\Gamma}\Big) \tilde\varrho \hat Q_{\rm{S}} \end{eqnarray}$ (6.45)

where

 $\begin{eqnarray} \tilde\Gamma &\equiv & \Big(\eta_{+}\frac{r_1}{r_2}+\eta_{-}\frac{r_2}{r_1}\Big)\tilde\partial_{p} +(\eta_{+}+\eta_{-})\tilde\partial_{x}\nonumber\\ &=& \Big[\delta_1+\frac{i}{2}(1+\delta_2)\Big]\tilde\partial_{p} +\Big(\delta_0+\frac{1}{\beta\bar\Omega}\Big)\tilde\partial_{x} \end{eqnarray}$ (6.46)

and

 $\begin{eqnarray} \widehat{\Gamma} &\equiv & \Big(\bar\eta^{\ast}_{+}\frac{r_1}{r_2} +\bar\eta^{\ast}_{-}\frac{r_2}{r_1}\Big)\tilde\partial_{p} +(\bar\eta^{\ast}_{+}+\bar\eta^{\ast}_{-})\tilde\partial_{x}\nonumber\\ &=& \Big[\delta_1-\frac{i}{2}(1+\delta_2)\Big]\tilde\partial_{p} +\Big(\delta_0+\frac{1}{\beta\bar\Omega}\Big)\tilde\partial_{x} \end{eqnarray}$ (6.47)

The second identity in each of the above two equations is obtained by using Eq.(6.36) and also $\bar\eta^{\ast}_{+}+\bar\eta^{\ast}_{-}=\eta_{+}+\eta_{-}$; see Eq.(5.28). We can then rearrange Eq.(6.45) as

 $\begin{eqnarray} {\cal V}^{(1)}_{\rm eff}\hat\varrho &\Leftrightarrow&[\hat Q_{\rm{S}}, \big(\tilde x_{\rm{B}}+\delta_0 \tilde\partial_{x} +\delta_1 \tilde\partial_{p}\big) \tilde\varrho\, ] +\nonumber\\ &&\frac{i}{2}(1+\delta_2)\{\hat Q_{\rm{S}}, \tilde\partial_{p} \tilde\varrho\, \} \end{eqnarray}$ (6.48)

By applying Eq.(6.1), we obtain

 $\begin{eqnarray}\label{FP_QME_termI_final} {\cal V}^{(1)}_{\rm eff}\hat\varrho =[\hat Q_{\rm{S}}, {\cal W}_{\rm{B}} \hat\varrho] +\frac{i}{2} \{\hat Q_{\rm{S}}, {\cal U}_{\rm{B}}\hat\varrho \} \end{eqnarray}$ (6.49)

This is just Eq.(5.30).

The influence of the quadratic coupling bath can be treated similarly, as follows. By comparing between the $\alpha_2$-term in Eq.(6.25), and that in Eq.(5.21), we have

 $\begin{eqnarray}\label{vs2} {\cal V}^{(2)}_{\rm eff}\hat\varrho(t) &\Leftrightarrow& {\cal B}_{20}\rho^{\{n-2\}}_{20} +{\cal B}_{02}\rho^{\{n-2\}}_{02} +2{\cal B}_{11}\rho^{\{n-2\}}_{11} +\nonumber\\&& {\cal A}\rho^{\{n+2\}}_{n_1n_2} +2{\cal C}_1\rho^{\{n\}}_{1}+2{\cal C}_2\rho^{\{n\}}_{2} \end{eqnarray}$ (6.50)

Apparently, the evaluations involve double-actions of linear operators in the Wigner phase-space. In specific, the double-action of Eq.(6.40) or Eq.(6.43) would lead to

 $\begin{eqnarray}\label{app_varrho+2_FP} \rho^{\{n+2\}}_{n_1n_2} \Leftrightarrow \Big(\tilde x_{\rm{B}}+\frac{1}{\beta\bar\Omega}\tilde\partial_x\Big)^2 \tilde\varrho \end{eqnarray}$ (6.51)

Moreover, by applying the phase-space action in Eq.(6.43), followed by each individual in Eq.(6.42), we obtain

 $\begin{eqnarray} \begin{split} &\rho^{\{n\}}_1\Leftrightarrow \Big[ \Big(\tilde x_{\rm{B}} + \frac{1}{\beta\bar\Omega}\tilde\partial_{x} \Big)\Big( \frac{r_1}{r_2}\tilde\partial_{p}-\tilde\partial_{x} \Big) \Big]\tilde\varrho \\ &\rho^{\{n\}}_2\Leftrightarrow \Big[ \Big(\tilde x_{\rm{B}} + \frac{1}{\beta\bar\Omega}\tilde\partial_{x} \Big)\Big( \frac{r_2}{r_1}\tilde\partial_{p}-\tilde\partial_{x} \Big) \Big]\tilde\varrho \end{split} \end{eqnarray}$ (6.52)

Furthermore, the pairs action in Eq.(6.42) result in

 $\begin{eqnarray} \begin{split} &\rho^{\{n-2\}}_{20}\Leftrightarrow \Big( \frac{r_1}{r_2}\tilde\partial_p-\tilde\partial_x \Big)^2\tilde\varrho \\ &\rho^{\{n-2\}}_{02}\Leftrightarrow \Big( \frac{r_2}{r_1}\tilde\partial_p-\tilde\partial_x \Big)^2\tilde\varrho \\ &\rho^{\{n-2\}}_{11}\Leftrightarrow \Big( \frac{r_1}{r_2}\tilde\partial_p-\tilde\partial_x \Big) \Big( \frac{r_2}{r_1}\tilde\partial_p-\tilde\partial_x \Big) \tilde\varrho \end{split} \end{eqnarray}$ (6.53)

By applying the same procedure, as from Eq.(6.44) to Eq.(6.49), we obtain Eq.(6.50) the expression,

 $\begin{eqnarray}\label{FP_QME_termII_final} {\cal V}^{(2)}_{\rm eff}\hat\varrho =\Big[\hat Q_{\rm{S}}, \Big({\cal W}^2_{\rm{B}}-\frac{1}{4}{\cal U}^2_{\rm{B}}\Big)\hat\varrho\, \Big] +i \big\{\hat Q_{\rm{S}}, {\cal W}_{\rm{B}}{\cal U}_{\rm{B}}\hat\varrho\big\}\nonumber\\ \end{eqnarray}$ (6.54)

This is just Eq.(5.31). We have thus completed the construction of FP-QME (5.21), with Eqs. (5.30)-(5.32).

Ⅶ. CONCLUDING REMARKS

In summary, we present a comprehensive account on the theories of quantum open system, with quadratic coupling bath environments. It is worth to reemphasize that even the nonlinear coupling bath characterization itself is nontrivial. The results presented in Section Ⅳ just represent one of physically supported schemes.

Both the DEOM and FP-QME approaches to the entangled system-and-environment dynamics are presented and scrutinized, with the focus on the underlying quasi-particle picture. Apparently, the DEOM approach enjoys a much clearer physical picture and much more friendly use in evaluating various systems. Nevertheless, the extended FP-QME, which is developed in this work, is carried out with a rather straightforward algorithm. The underlying equivalence between these two dynamics approaches would shed some light on the advancement of the DEOM theory with nonlinear coupling bath environments.

Ⅷ. ACKNOWLEDGEMENTS

This work was supported from the Ministry of Science and Technology (No.2016YFA0400900), the National Natural Science Foundation of China (No.21373191, No.21633006, and No.21303090), and the Fundamental Research Funds for the Central Universities (No.2030020028).

Reference
 [1] R. P. Feynman, F. L. Vernon Jr., Ann. Phys. 24 , 118 (1963). [2] J. S. Jin, X. Zheng, and Y. J. Yan, J. Chem. Phys. 128 , 234703 (2008). DOI:10.1063/1.2938087 [3] Y. Tanimura, Phys. Rev. A 41 , 6676 (1990). DOI:10.1103/PhysRevA.41.6676 [4] Y. Tanimura, J. Phys. Soc. Jpn. 75 , 082001 (2006). DOI:10.1143/JPSJ.75.082001 [5] R. X. Xu, P. Cui, X. Q. Li, Y. Mo, and Y. J. Yan, J. Chem. Phys. 122 , 041103 (2005). DOI:10.1063/1.1850899 [6] J. J. Ding, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 136 , 224103 (2012). DOI:10.1063/1.4724193 [7] 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 [8] R. X. Xu, and Y. J. Yan, Phys. Rev. E 75 , 031107 (2007). DOI:10.1103/PhysRevE.75.031107 [9] J. S. Shao, J. Chem. Phys. 120 , 5053 (2004). DOI:10.1063/1.1647528 [10] J. M. Moix, and J. S. Cao, J. Chem. Phys. 139 , 134106 (2013). DOI:10.1063/1.4822043 [11] Y. L. Ke, and Y. Zhao, J. Chem. Phys. 145 , 024101 (2016). DOI:10.1063/1.4955107 [12] Y. L. Ke, and Y. Zhao, J. Chem. Phys. 146 , 174105 (2017). DOI:10.1063/1.4982230 [13] R. Kubo, M. Toda, and N. Hashitsume, Statistical Physics Ⅱ: Nonequilibrium Statistical Mechanics, Berlin Heidelberg: Springer-Verlag, (1991). https://www.researchgate.net/publication/234291363_Statistical_Physics_II_Nonequilibrium_Statistical_Mechanics?ev=auth_pub [14] Y. J. Yan, and R. X. Xu, Annu. Rev. Phys. Chem. 56 , 187 (2005). DOI:10.1146/annurev.physchem.55.091602.094425 [15] U. Weiss, Quantum Dissipative Systems, 4th Edn. Singapore: World Scientific, (2012) [16] H. Kleinert, Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets, 5th Edn., Singapore: World Scientific, (2009). http://www.ams.org/mathscinet-getitem?mr=2518082 [17] Y. J. Yan, J. Chem. Phys. 140 , 054105 (2014). DOI:10.1063/1.4863379 [18] Y. J. Yan, J. S. Jin, R. X. Xu, and X. Zheng, Frontiers Phys. 11 , 110306 (2016). DOI:10.1007/s11467-016-0513-5 [19] H. D. Zhang, R. X. Xu, X. Zheng, and Y. J. Yan, Mol. Phys. 116 , 780 (2018). DOI:10.1080/00268976.2018.1431407 [20] R. X. Xu, Y. Liu, H. D. Zhang, and Y. J. Yan, Chin. J. Chem. Phys. 30 , 395 (2017). DOI:10.1063/1674-0068/30/cjcp1706123 [21] R. X. Xu, Y. Liu, H. D. Zhang, and Y. J. Yan, J. Chem. Phys. 148 , 114103 (2018). DOI:10.1063/1.4991779 [22] A. A. Golosov, R. A. Friesner, and P. Pechukas, J. Chem. Phys. 112 , 2095 (2000). DOI:10.1063/1.480888 [23] M. Thoss, H. B. Wang, and W. H. Miller, J. Chem. Phys. 115 , 2991 (2001). DOI:10.1063/1.1385562 [24] A. Garg, J. N. Onuchic, and V. Ambegaokar, J. Chem. Phys. 83 , 4491 (1985). DOI:10.1063/1.449017 [25] L. D. Zusman, Chem. Phys. 49 , 295 (1980). DOI:10.1016/0301-0104(80)85267-0 [26] L. D. Zusman, Chem. Phys. 80 , 29 (1983). DOI:10.1016/0301-0104(83)85166-0 [27] A. O. Caldeira, and A. J. Leggett, Ann. Phys. 149 , 374 (1983). DOI:10.1016/0003-4916(83)90202-6 [28] A. O. Caldeira, and A. J. Leggett, Physica A 121 , 587 (1983). DOI:10.1016/0378-4371(83)90013-4 [29] J. J. Ding, H. D. Zhang, Y. Wang, R. X. Xu, X. Zheng, and Y. J. Yan, J. Chem. Phys. 145 , 204110 (2016). DOI:10.1063/1.4967964 [30] J. J. Ding, Y. Wang, H. D. Zhang, R. X. Xu, X. Zheng, and Y. J. Yan, J. Chem. Phys. 146 , 024104 (2017). DOI:10.1063/1.4973610 [31] H. Liu, L. L. Zhu, S. M. Bai, and Q. Shi, J. Chem. Phys. 140 , 134106 (2014). DOI:10.1063/1.4870035 [32] J. Hu, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 133 , 101106 (2010). DOI:10.1063/1.3484491 [33] J. Hu, M. Luo, F. Jiang, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 134 , 244106 (2011). DOI:10.1063/1.3602466 [34] J. S. Jin, S. K. Wang, X. Zheng, and Y. J. Yan, J. Chem. Phys. 142 , 234108 (2015). DOI:10.1063/1.4922712 [35] H. D. Zhang, R. X. Xu, X. Zheng, and Y. J. Yan, J. Chem. Phys. 142 , 024112 (2015). DOI:10.1063/1.4905494 [36] H. D. Zhang, Q. Qiao, R. X. Xu, and Y. J. Yan, Chem. Phys. 481 , 237 (2016). DOI:10.1016/j.chemphys.2016.07.005 [37] H. D. Zhang, Q. Qiao, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 145 , 204109 (2016). DOI:10.1063/1.4968031 [38] H. Risken, The Fokker-Planck Equation, Methods of Solution and Applications, 2nd Edn., Berlin: SpringerVerlag, (1989). http://okmd.dyndns.co.za/1310318317313317317/The-Fokker-Planck-Equation-Methods-of-Solution-and-Applications-by-Hannes-Risken.pdf [39] S. Mukamel, The Principles of Nonlinear Optical Spectroscopy. New York: Oxford University Press (1995). [40] Y. J. Yan, and S. Mukamel, J. Chem. Phys. 85 , 5908 (1986). DOI:10.1063/1.451502 [41] B. L. Hu, J. P. Paz, and Y. Zhang, Phys. Rev. D 45 , 2843 (1992). DOI:10.1103/PhysRevD.45.2843 [42] R. X. Xu, B. L. Tian, J. Xu, and Y. J. Yan, J. Chem. Phys. 130 , 074107 (2009). DOI:10.1063/1.3078024 [43] R. X. Xu, Y. Mo, P. Cui, S. H. Lin, and Y. J. Yan, in Progress in Theoretical Chemistry and Physics, Vol. 12: Advanced Topics in Theoretical Chemical Physics, J. Maruani, R. Lefebvre, and E. Brändas, Eds. Dordrecht: Kluwer, 7-40(2003).