Chinese Journal of Polar Since  2017, Vol. 30 Issue (4): 395-403

#### The article information

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

Theory of Quantum Dissipation in a Class of Non-Gaussian Environments

Chinese Journal of Polar Since, 2017, 30(4): 395-403

http://dx.doi.org/10.1063/1674-0068/30/cjcp1706123

### Article history

Accepted on: June 27, 2017
Theory of Quantum Dissipation in a Class of Non-Gaussian Environments
Rui-xue Xu, Yang Liu, Hou-dao Zhang, Yan YiJing
Dated: Received on June 17, 2017; Accepted on June 27, 2017
Hefei National Laboratory for Physical Sciences at the Microscale and Department of Chemical Physics and Synergetic Innovation Center of Quantum Information and Quantum Physics and iChEM, University of Science and Technology of China, Hefei 230026, China
Author: YiJing Yan, E-mail: yanyj@ustc.edu.cn
Abstract: In this work we construct a novel dissipaton-equation-of-motion (DEOM) theory in quadratic bath coupling environment, based on an extended algebraic statistical quasi-particle approach. To validate the new ingredient of the underlying dissipaton algebra, we derive an extended Zusman equation via a totally different approach. We prove that the new theory, if it starts with the identical setup, constitutes the dynamical resolutions to the extended Zusman equation. Thus, we verify the generalized (non-Gaussian) Wick's theorem with dissipatons-pair added. This new algebraic ingredient enables the dissipaton approach being naturally extended to nonlinear coupling environments. Moreover, it is noticed that, unlike the linear bath coupling case, the influence of a non-Gaussian environment cannot be completely characterized with the linear response theory. The new theory has to take this fact into account. The developed DEOM theory manifests the dynamical interplay between dissipatons and nonlinear bath coupling descriptors that will be specified. Numerical demonstrations will be given with the optical line shapes in quadratic coupling environment.
Key words: Non-Gaussian environment     Quantum dissipation     Dissipaton theory
Ⅰ. INTRODUCTION

Quantum dissipation plays crucial roles in many fields of modern science. Exact theories include the Feynman-Vernon influence functional path integral approach [1], and its differential equivalence, the hierarchical-equations-of-motion (HEOM) formalism [2-6]. How-ever, almost all existing quantum dissipation theories exploit the Gaussian-Wick's thermodynamical statis-tics [7-9], which is strictly valid only for linear bath couplings. Intrinsically, a linear bath coupling implies a weak backaction of system on environment. The lowest non-Gaussian environment influence requires a quadratic bath coupling.

In this work, we extend the dissipaton equation of motion (DEOM) theory [10, 11] to treat the linear-plus-quadratic bath coupling environment. This theory goes with a statistical quasi-particle ("dissipaton") description for the hybrid environment that can be either bosonic or fermionic or excitonic. Dynamical variables in DEOM are the dissipaton density operators (DDOs), for both the reduced system and the hybrid bath dynamics [10, 11]. The latter could also be measured experimentally, via such as the Fano interference [12-16], vibronic spectroscopy with non-Condon polarized environment [17], and transport current noise spectrum [18].

Dissipaton algebra plays essential roles [10, 11].It consists of the generalized (non-Gaussian) Wick's theorem and the generalized diffusion equation. This noval algebra leads to the rules on how the DDOs evolves in time, and further on their relations to experimental measurable quantities that involve explicitly the hybrid bath dynamics [14-18]. From the algebraic construction point of view, the new DEOM theory in quest amounts to the establishment of the generalized Wick's theorem with dissipatons-pairs added. This will be the new ingredient of the dissipaton algebra for treating the quadratic bath coupling in study.

Another important issue is concerned with the characterization of nonlinear coupling bath. On top of the interacting bath correlation function description [7-9], additional information would be needed. This is a general concern in any non-Gaussian environment theories. To address this issue, we propose a polarization model to determine both the linear and nonlinear bath coupling strengths. This model resolves this issue, with a single additional parameter, on top of the conventional linear response theory.

In this work, we construct the DEOM formalism, via the dissipaton algebra, including the aforementioned new ingredient for treating quadratic bath coupling. We validate this new ingredient, the generalized Wick's theorem with dissipatons-pairs added. To do that we derive an extended Zusman equation via a totally different approach. We prove that the DEOM formalism, if it started with the same setup, constitutes the dynamical resolutions to the extended Zusman equation. Therefore, as the algebraic construction is concerned, we would have also confirmed the DEOM formulations. The linear and quadratic bath coupling strength parameters, will be discussed on the basis of the nonlinear polarization model. Numerical DEOM demonstrations are then carried out on the optical line shapes in the nonlinear coupling environment.

Ⅱ. DISSIPATON DYNAMICS THEORY A. Statistical quasi-particle description

 $$$\label{HT} H_{\rm{T}} =H_{\rm{S}} + h_{\rm{B}} + \hat Q_{\rm{S}}(\alpha_1\hat x_{\rm{B}}+\alpha_2\hat x_{\rm{B}}^2)$$$ (1)

The system Hamiltonian $H_{\rm{S}}$ and dissipative operator $\hat Q_{\rm{S}}$ are arbitrary. The latter is set to be dimensionless. The bath Hamiltonian and the hybridization bath operator (solvation mode) are given, respectively,

 ${h_{\rm{B}}} = \frac{1}{2}\sum\limits_j {{\omega _j}} (\hat p_j^2 + \hat q_j^2)$ (2)
 ${\hat x_{\rm{B}}} = \sum\limits_j {{c_j}} {\hat q_j}$ (3)

Throughout this work we set $\hbar$=1 and $\beta$=$1/(k_{\textrm{B}}T)$, with $k_\textrm{B}$ and $T$ being the Boltzmann constant and temperature. Let $\langle\hat O \rangle_{\rm{B}}$$\equiv$${\rm tr}_{\rm{B}}\big(\hat O \textrm{e}^{-\beta h_{\rm{B}}}\big)/{\rm tr}_{\rm{B}}\textrm{e}^{-\beta h_{\rm{B}}}$ be the bare bath ensemble average. Define $\hat x_{\rm{B}}(t)$$\equiv$$\textrm{e}^{ih_{\rm{B}}t}\hat x_{\rm{B}} \textrm{e}^{-ih_{\rm{B}}t}$. Set hereafter $t\geq 0$ for the time variable. We have [7, 9]

 $\begin{array}{l} {\langle {{\hat x}_{\rm{B}}}(t){{\hat x}_{\rm{B}}}(0)\rangle _{\rm{B}}} = \frac{1}{\pi }\int_{ - \infty }^\infty {\rm{d}} \omega \frac{{{{\rm{e}}^{ - i\omega t}}\chi _{\rm{B}}^{(i)}(\omega )}}{{1 - {{\rm{e}}^{ - \beta \omega }}}}\\ = \sum\limits_{k = 1}^K {{\eta _k}} {{\rm{e}}^{ - {\gamma _k}t}} \end{array}$ (4)

The first expression in Eq.(3) is the fluctuation-dissipation theorem [7, 9], with $\chi^{(i)}_{\rm{B}}(\omega)$ being the imaginary part of

 $$$\label{chiw_def} \chi_{\rm{B}}(\omega)\equiv i\int^{\infty}_{0}\!{\rm d}t\, \textrm{e}^{i\omega t} \langle[\hat x_{\rm{B}}(t), \hat x_{\rm{B}}(0)]\rangle_{\rm{B}}$$$ (5)

The second expression of Eq.(4) presents an exponential series expansion of the linear bath correlation function. Set hereafter $\hat x_{\rm{B}}$ to be dimensionless, so that the bath coupling parameters, $\alpha_1$ and $\alpha_2$ in Eq.(1), are of frequency unit. It is well known that, for a complete characterization of the nonlinear environment influence (\alpha_2$$\neq0), additional information, on top of Eq.(4), is needed. We will address this issue latter. Dissipatons, \{\hat f_k\}, arise strictly via the linear bath coupling part, as follows:  {\hat x_{\rm{B}}} = \sum\limits_{k = 1}^K {{{\hat f}_k}} (6) with  $$\label{ff_t} \begin{split} \langle \hat f_{k}(t)\hat f_{j}(0)\rangle_{\rm{B}}&=\delta_{kj}\eta_{k}\textrm{e}^{-\gamma_k t} \ \\ \langle \hat f_{j}(0)\hat f_{k}(t)\rangle_{\rm{B}}&=\delta_{kj}\eta^{\ast}_{\bar k}\textrm{e}^{-\gamma_k t} \end{split}$$ (7) The associated index \bar k in the second expression above is defined via \gamma_{\bar k}$$\equiv\gamma^{\ast}_{k}. It is easy to verify that both Eq.(4) and its time reversal are reproduced. Denote for later use  $$\label{ff_t0} \begin{split} \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} \\ \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}$$ (8) They differ from \langle \hat f_{k}\hat f_{j}\rangle_{\rm{B}}, which will also appear explicitly in the dissipaton algebra (cf. Eq.(13) and Eq.(14)). From Eq.(6), we have \langle\hat x^2_{\rm{B}}\rangle_{\rm{B}}=\displaystyle\sum_{kj} \langle \hat f_{k}\hat f_{j}\rangle_{\rm{B}}. As defined in Eq.(7), dissipatons \{\hat f_k\} are statistically independent quasi-particles. Each of them is a macroscopic linear combination of bath degrees of freedoms. Individual \hat f_k is characterized by a single damping parameter, \gamma_{k}, that can be complex, and a joint-pair of interacting strength parameters, \eta_k and \eta^{\ast}_{\bar k}. B. Dissipaton algebra and DEOM formalism Dynamical variables in DEOM are dissipaton density operators (DDOs) [10, 11]:  $$\label{DDOs} \rho^{(n)}_{\bf n}(t)\equiv \rho^{(n)}_{n_1\cdots n_K}(t)\equiv {\rm tr}_{\rm{B}}\big[ (\hat f_{K}^{n_K}\cdots\hat f_{1}^{n_1})^{\circ}\rho_{\rm{T}}(t)\big]$$ (9) The dissipatons product inside (\cdots)^{\circ} is irreducible, such that(\text{c-number})^{\circ}=0. Bosonic dissipatons satisfy the symmetric permutation, (\hat f_k\hat f_j)^{\circ}=(\hat f_j\hat f_k)^{\circ}. Physically, each DDO of Eq.(9) stands for a given configuration of the total n=n_1+\cdots +n_k dissipatons. Denote for the use below the associated DDO's index, {\bf n}^{\pm}_k, which differs from {\bf n}\equiv n_1\cdots n_K at the specified n_k by \pm 1. Similarly, {\bf n}^{\pm\pm}_{kj} differs from {\bf n} at the specified n_k and n_j that are replaced by n_k\pm 1 and n_j\pm 1, respectively. The DEOM formalism can be easily constructed via the algebraic dissipaton approach [10, 11]. The construction starts with applying the Liouville-von Neuman equation, \dot\rho_{\rm{T}}(t)=-i[H_{\rm{S}}+h_{\rm{B}}+H_{\rm{SB}}, \rho_{\rm{T}}(t)], for the total density operator in Eq.(9). The bath h_{\rm{B}}-action and the system-bath coupling H_{\rm{SB}}-action are then readily evaluated with the generalized diffusion equation and the generalized Wick's theorem, respectively. The generalized diffusion equation arises from the single-damping parameter characteristic, as shown in Eq.(7), for both its forward and backward correlation functions. This feature leads to [10, 11]  $$\label{diff0} {\rm tr}_{\rm{B}}\Big[\Big(\frac{\partial \hat f_k}{\partial t}\Big)_{\rm{B}}\rho_{\rm{T}}(t)\Big] = -\gamma_{k}\, {\rm tr}_{\rm{B}}\big[\hat f_k\rho_{\rm{T}}(t)\big]$$ (10) Together with \big(\partial \hat f_k/\partial t\big)_{\rm{B}}=-i[\hat f_k, h_{\rm{B}}], the generalized diffusion equation leads immediately to  \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}}(t)]\big\}\nonumber\\ &&\quad\quad = i\, {\rm tr}_{\rm{B}}\!\left\{\! \big[(\hat f_{K}^{n_K}\cdots\hat f_{1}^{n_1})^{\circ}, h_{\rm{B}}\big]\rho_{\rm{T}}(t)\right\}\nonumber\\ &&\quad\quad = \left(\sum\limits_{k=1}^{K} n_k\gamma_k\right)\rho^{(n)}_{\bf n}(t) \end{eqnarray} (11) This is the bath h_{\rm{B}}-action contribution to the DDOs dynamics. The generalized diffusion Eq.(10) or Eq.(11), with an arbitrary total composite \rho_{\rm{T}}(t), has been validated previously [10, 11]. The generalized Wick's theorem (GWT) deals with the system-bath coupling H_{\rm{SB}}-action [10, 11]. In this work this theorem has two ingredients, GWT-1 and GWT-2, in relation to the linear and quadratic bath couplings, respectively. The GWT-1 reads [10, 11]  \begin{eqnarray}\label{Wick1} {\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] = \rho^{(n+1)}_{{\bf n}^{+}_{j}}(t) +\nonumber\\ \sum\limits_{j} n_k \langle \hat f_{k}\hat f_{j}\rangle^{>}_{\rm{B}}\rho^{(n-1)}_{{\bf n}^{-}_{k}}(t) \end{eqnarray} (12) The expression for {\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 goes with \langle \hat f_{j}\hat f_{k}\rangle^{<}_{\rm{B}}. It together with \langle \hat f_{k}\hat f_{j}\rangle^{>}_{\rm{B}} is defined in Eq.(8). The GWT-1 has been well established and used in evaluating the commutator action of linear bath coupling terms [10, 11]. The GWT-2 is concerned with the quadratic bath couplings, where a pair of dissipatons (\hat f_j\hat f_{j'}) participate insimultaneously without time-ordering.This new ingredient of dissipaton algebra would read  \begin{align} & \rm{t}{{\rm{r}}_{\rm{B}}}[{{(\mathit{\hat{f}}_{\mathit{K}}^{{{\mathit{n}}_{\mathit{K}}}}\cdots \mathit{\hat{f}}_{1}^{{{\mathit{n}}_{1}}})}^{{}^\circ }}({{{\mathit{\hat{f}}}}_{\mathit{j}}}{{{\mathit{\hat{f}}}}_{{\mathit{{j}'}}}}){{\rho }_{\rm{T}}}(\mathit{t})]=\sum\limits_{\mathit{k}}{{{\mathit{n}}_{\mathit{k}}}}\langle {{{\mathit{\hat{f}}}}_{\mathit{k}}}{{{\mathit{\hat{f}}}}_{\mathit{j}}}\rangle _{\rm{B}}^{>}\cdot \\ & \rm{t}{{\rm{r}}_{\rm{B}}}[{{(\mathit{\hat{f}}_{\mathit{K}}^{{{\mathit{n}}_{\mathit{K}}}}\cdots \mathit{\hat{f}}_{k}^{{{\mathit{n}}_{\mathit{k}}}-1}\cdots \mathit{\hat{f}}_{1}^{{{\mathit{n}}_{1}}})}^{{}^\circ }}{{{\mathit{\hat{f}}}}_{{\mathit{{j}'}}}}{{\rho }_{\rm{T}}}(\mathit{t})]+ \\ & \rm{t}{{\rm{r}}_{\rm{B}}}[{{(\mathit{\hat{f}}_{\mathit{K}}^{{{\mathit{n}}_{\mathit{K}}}}\cdots \mathit{\hat{f}}_{1}^{{{\mathit{n}}_{1}}}{{{\mathit{\hat{f}}}}_{\mathit{j}}})}^{{}^\circ }}{{{\mathit{\hat{f}}}}_{{\mathit{{j}'}}}}{{\rho }_{\rm{T}}}(\mathit{t})] \\ \end{align} (13) with the last term being evaluated as  \begin{eqnarray}\label{Wick21} \ {\rm tr}_{\rm{B}}\big[(\hat f_{K}^{n_K}\cdots\hat f_{1}^{n_1}\hat f_j)^{\circ} \hat f_{j'}\rho_{\rm{T}}(t)\big]= \rho^{(n+2)}_{{\bf n}^{++}_{jj'}}(t) +\nonumber\\ \langle\hat f_j \hat f_{j'}\rangle_{\rm{B}}\rho^{(n)}_{{\bf n}}(t) + \sum\limits_{k} n_k \langle\hat f_k\hat f_{j'}\rangle^{>}_{\rm{B}}\rho^{(n)}_{{\bf n}^{-+}_{kj}}(t) \end{eqnarray} (14) Together with the first term in Eq.(13) being evaluated via GWT-1, we obtain  \begin{align} & \rm{t}{{\rm{r}}_{\rm{B}}}[{{(\hat{f}_{K}^{{{n}_{K}}}\cdots \hat{f}_{1}^{{{n}_{1}}})}^{{}^\circ }}({{{\hat{f}}}_{j}}{{{\hat{f}}}_{{{j}'}}}){{\rho }_{\rm{T}}}(t)]=\rho _{\bf{n}_{j{j}'}^{++}}^{(n+2)}(t)+ \\ & {{\langle {{{\hat{f}}}_{j}}{{{\hat{f}}}_{{{j}'}}}\rangle }_{\rm{B}}}\rho _{\bf{n}}^{(n)}(t)+\sum\limits_{k}{{{n}_{k}}}[\langle {{{\hat{f}}}_{k}}{{{\hat{f}}}_{j}}\rangle _{\rm{B}}^{>}\rho _{\bf{n}_{k{j}'}^{-+}}^{(n)}(t)+ \\ & \langle {{{\hat{f}}}_{k}}{{{\hat{f}}}_{{{j}'}}}\rangle _{\rm{B}}^{>}\rho _{\bf{n}_{kj}^{-+}}^{(n)}(t)]+\sum\limits_{k, {k}'}{{{n}_{k}}}({{n}_{{{k}'}}}-{{\delta }_{k{k}'}}) \\ & \langle {{{\hat{f}}}_{k}}{{{\hat{f}}}_{j}}\rangle _{\rm{B}}^{>}\langle {{{\hat{f}}}_{{{k}'}}}{{{\hat{f}}}_{{{j}'}}}\rangle _{\rm{B}}^{>}\rho _{{{\bf{n}}_{\overline{k}{{\ }_{\overline{{{k}'}}}}}}}^{(n-2)}(t) \\ \end{align} (15) The expression for {\rm tr}_{\rm{B}}\big[(\hat f_{K}^{n_K}\cdots\hat f_{1}^{n_1})^{\circ}\rho_{\rm{T}}(t)(\hat f_j\hat f_{j'})\big] is similar, but with each individual \langle \hat f_{k}\hat f_{j}\rangle^{>}_{\rm{B}} being replaced by \langle \hat f_{j}\hat f_{k}\rangle^{<}_{\rm{B}}. The GWT-2 is to be used in the evaluation of the quadratic bath coupling contribution. The DEOM formalism in the presence of both linear and quadratic bath couplings can now be readily constructed via the above dissipaton algebra. The final results read  \begin{array}{l} \dot \rho _{\bf{n}}^{(n)} = - (i{\mathcal{L}_{{\rm{eff}}}} + \sum\limits_k {{n_k}} {\gamma _k})\rho _{\bf{n}}^{(n)} - i2{\alpha _2}\sum\limits_{kj} {{n_k}} {C_k}\rho _{{\bf{n}}_{\mathit{kj}}^{ - + }}^{(n)} - \\ i{\alpha _2}\sum\limits_{kj} [A\rho _{{\bf{n}}_{\mathit{kj}}^{ + + }}^{(n + 2)} + {n_k}({n_j}-{\delta _{jk}}){B_{kj}}\rho _{{\bf{n}}_{\mathit{kj}}^{--}}^{(n -2)}] - \\ i{\alpha _1}\sum\limits_k ( A\rho _{{\bf{n}}_\mathit{k}^ + }^{(n + 1)} + {n_k}{C_k}\rho _{{\bf{n}}_\mathit{k}^ - }^{(n - 1)}) \end{array} (16) Here,  \begin{array}{*{20}{l}} {{\mathcal{L}_{{\rm{eff}}}}\hat O}&{ \equiv [{H_{{\rm{eff}}}}, \hat O]\;\;{\mkern 1mu} A\hat O \equiv \left[{{{\hat Q}_{\rm{S}}}, \hat O} \right]}\\ {{B_{kj}}\hat O}&{ \equiv {\eta _k}{\eta _j}{{\hat Q}_{\rm{S}}}\hat O - \eta _{\bar k}^ * \eta _{\bar j}^ * \hat O{{\hat Q}_{\rm{S}}}}\\ {{C_k}\hat O}&{ \equiv {\eta _k}{{\hat Q}_{\rm{S}}}\hat O - \eta _{\bar k}^ * \hat O{{\hat Q}_{\rm{S}}}} \end{array} (17) with  $$\label{Heff} H_{\rm eff} \equiv H_{\rm{S}}+\alpha_2\langle\hat x^2_{\rm{B}}\rangle_{\rm{B}}\hat Q_{\rm{S}}$$ (18) Evidently, the dissipatons defined in Eqs.(4)-(7) are strictly based on the linear bath coupling part that satisfies Gaussian-Wick's statistics [7-9]. The quadratic non-Gaussian bath influences are treated via the GWT-2, Eq.(15). This is the virtue of the DEOM theory that includes the powerful dissipaton algebra [10, 11]. In general, Eqs. (10)-(15) are all non-Gaussian operators in the system subspace. Ⅲ. VALIDATION ON DISSIPATON ALGEBRA WITH EXTENDED ZUSMAN EQUATION This section is devoted to validating the GWT-2, Eqs. (13) and (14), the new ingredient of the dissipaton algebra presented above. It together with the well-established Eqs. (10) and (12) lead immediately and unambiguously to the extended DEOM (Eq.(16)). Therefore, from the algebraic construction point of view, the required validation can be made with the dissipaton basis set of size K=1. This amounts to the formal setting Eq.(4) with  $$\label{xxt_single} \langle \hat x_{\rm{B}}(t)\hat x_{\rm{B}}(0)\rangle_{\rm{B}} \simeq \eta \textrm{e}^{-\gamma t}$$ (19) The dissipaton index {\bf n} can be omitted; i.e., \rho^{(n)}_{\bf n}(t)= \rho^{(n)}(t)\equiv \hat\rho_{n}(t), for the basis set size K=1 case, in which the DEOM (Eq.(16)) reads  \begin{eqnarray}\label{DEOM_single} \dot{\hat\rho}_{n}(t)= - ( i{\cal L}_{\rm eff} + n\gamma +i2\alpha_2n{\cal C}) \hat\rho_{n}(t) -\nonumber\\ i\alpha_2\left[{\cal A}\hat \rho_{n+2}(t)+n(n-1){\cal B}\hat\rho_{n-2}(t) \right] -\nonumber\\i\alpha_1\left[{\cal A}\hat \rho_{n+1}(t) + n{\cal C}\hat \rho_{n-1}(t) \right] \end{eqnarray} (20) The involving superoperators {\cal A}, {\cal B} and {\cal C}, are the same as those in Eq.(17), but without dissipaton indexes. In the following, on the basis of Eq.(19) which will be called the Zusman setup, we construct the extended Zusman equation via a totally different approach. By showing that the extended Zusman equation is identical to Eq.(20), we validate the dissipaton algebra and thus the extended DEOM theory, Eq.(16). The proof here is rigourous, due to the nature of the algebraic construction, despite of the fact that the Zusman setup, Eq.(19), itself could even be a bad approximation. It is well known that the Zusman setup, Eq.(19), is equivalent to the combination of the high-temperature (HT) and the Smoluchowski limits. The HT limit is characterized with  $$\label{higtT} \frac{1}{1-\textrm{e}^{-\beta\omega}} \simeq \frac{1}{\beta\omega}+\frac{1}{2}$$ (21) In this case, the solvation mode is a classical Brownian motion in the secondary bath environment. The latter exerts a stochastic force \widetilde F(t) and friction constant \zeta on the solvation mode. The corresponding Langevin equation reads  $$\label{Langevin} \dot{p}_{\rm{B}}(t)=-\omega_{\rm{B}}x_{\rm{B}}(t) -\zeta p_{\rm{B}}(t) +\widetilde F (t)$$ (22) with  $$\label{corr_HT} \langle \widetilde F(t)\widetilde F(0)\rangle^{\text{ HT}} = \frac{2\zeta}{\beta\omega_{\rm{B}}}\delta(t) + i \frac{\zeta}{\omega_{\rm{B}}}\, {\dot\delta}(t)$$ (23) This is the high-temperature fluctuation-dissipation theorem. The resultant Caldeira-Leggett's equation reads [19-21]  \begin{eqnarray}\label{rhoCL} \frac{\partial}{\partial t}{\hat\rho}^{\text{ HT}}_{\text{ W}} =-(i{\cal L}_{\rm{S}}+\hat L_{\text{ FP}}){\hat\rho}^{\text{ HT}}_{\text{ W}}+ \frac{\partial}{\partial p_{\rm{B}}}\Big(\frac{\alpha_1}{2}+\alpha_2x_{\rm{B}}\Big)\cdot\nonumber\\ \{\hat Q_{\rm{S}}, {\hat\rho}^{\text{ HT}}_{\text{ W}}\} -i\left(\alpha_1x_{\rm{B}} +\alpha_2x_{\rm{B}}^2-\frac{\alpha_2}{4} \frac{\partial^2}{\partial p_{\rm{B}}^2}\right) [\hat Q_{\rm{S}}, {\hat\rho}^{\text{ HT}}_{\text{ W}}] \end{eqnarray} (24) Here, {\hat\rho}^{\text{ HT}}_{\text{ W}}\equiv {\hat\rho}^{\text{ HT}}_{\text{ W}}(x_{\rm{B}}, p_{\rm{B}};t) denotes the reduced system-plus-solvation density operator in the HT limit, with the solvation mode in the Wigner representation. While {\cal L}_{\rm{S}}\, \cdot\,\equiv $$[H_{\rm{S}}\, \cdot\, ] is the system Liouvillian, \hat L_{\text{ FP}} in Eq.(24) is the Fokker-Planck operator [22],  $$\label{calFP} \hat L_{\text{ FP}}= \omega_{\rm{B}}\Big(\frac{\partial}{\partial x_{\rm{B}}}p_{\rm{B}} -\frac{\partial}{\partial p_{\rm{B}}}x_{\rm{B}}\Big) - \frac{\zeta}{\beta\omega_{\rm{B}}} \frac{\partial^2}{\partial p_{\rm{B}}^2} -\zeta\frac{\partial}{\partial p_{\rm{B}}}p_{\rm{B}}$$ (25) To complete the Zusman setup, Eq.(19), consider further the Smoluchowski (or strongly-overdamped) limit; i.e., \zeta$$\gg$$\omega_{\rm{B}}, whereas \omega_{\rm{B}}^2/\zeta=\gamma remains finite to be the exponent in Eq.(19). Moreover, it is easy to show that \langle x^2_{\rm{B}}\rangle_{\rm{B}}=\langle p^2_{\rm{B}}\rangle_{\rm{B}}\rightarrow(\beta\omega_{\rm{B}})^{-1}, in the high-temperature limit. The identities here will be used in eliminating the appearance of \beta\omega_{\rm{B}} in the formulations below. The pre-exponential coefficient in Eq.(19) reads then  $$\label{eta_Zusman} \eta \equiv \eta_r+i\eta_i = \langle x^2_{\rm{B}}\rangle_{\rm{B}}\left(1-\frac{i\beta\gamma}{2}\right)$$ (26) In the strongly-overdamped limit, the momentum p_{\rm{B}} would no longer be a correlated dynamical variable. The equation of motion for \hat\rho(x_{\rm{B}};t)=\int{\rm d}p_{\rm{B}}\hat\rho_\text{ W}(x_{\rm{B}}, p_{\rm{B}};t), which is closed now, can be obtained via the standard Fokker-Planck-Smoluchowski algorithm [22]. Another equivalent but much simpler approach is the so-called diffusion mapping method [23]; i.e., mapping each individual p_{\rm{B}}-space variable to its limiting diffusive x_{\rm{B}}-space correspondence. This method makes a simple use of the Langevin equation Eq.(22), which in the strongly-overdamped limit reduces to  $$\label{Smolimit0} 0 = -\dot p_{\rm{B}} = \omega_{\rm{B}} x_{\rm{B}} + \zeta p_{\rm{B}} - \widetilde F(t)$$ (27) Consider further the following two thermodynamic relations,  $$\label{partial_xp} \begin{split} \frac{\partial}{\partial p_{\rm{B}}}\hat\rho^{\text{ HT}}_{\text{ W}} &= -\frac{1}{\langle p_{\rm{B}}^2\rangle}_{\rm{B}} p_{\rm{B}} \hat\rho^{\text{ HT}}_{\text{ W}} \\ \frac{\partial}{\partial x_{\rm{B}}}\hat\rho^{\text{ HT}}_{\text{ W}} &= -\frac{1}{\langle x_{\rm{B}}^2\rangle}_{\rm{B}} x_{\rm{B}} \hat\rho^{\text{ HT}}_{\text{ W}} \end{split}$$ (28) Together with p_{\rm{B}}/x_{\rm{B}}$$\simeq$$-\omega_{\rm{B}}/\zeta, as implied in Eq.(27), and {\langle x_{\rm{B}}^2\rangle}_{\rm{B}}={\langle p_{\rm{B}}^2\rangle}_{\rm{B}} in study here, we would have then  $$\label{map0pm} \frac{\partial}{\partial p_{\rm{B}}} \hat\rho^{\text{ HT}}_{\text{ W}} \simeq \frac{\omega_{\rm{B}}}{\zeta} \frac{\partial}{\partial x_{\rm{B}}}\hat\rho^{\text{ HT}}_{\text{ W}} \,$$ (29) The above results of the Zusman setup lead to the following rules of diffusion mapping [23],  \frac{\partial }{{\partial {p_{\rm{B}}}}} \Rightarrow \frac{{{\omega _{\rm{B}}}}}{\zeta }\frac{\partial }{{\partial {x_{\rm{B}}}}} (30)  {p_{\rm{B}}} \Rightarrow - {\langle x_{\rm{B}}^2\rangle _{\rm{B}}}\frac{{{\omega _{\rm{B}}}}}{\zeta }\frac{\partial }{{\partial {x_{\rm{B}}}}} (31) The Smoluchowski limit to the Calderia-Leggett's master equation is now readily obtained by replacing all those p_{\rm{B}}-dependent operators in Eq.(24) and Eq.(25). In particular, the Fokker-Planck operator becomes the Smoluchowski or diffusion operator,  \begin{eqnarray}\label{LD} {\hat L}_{\text{ FP}} \Longrightarrow{\hat L}_{\text{ D}} \equiv-\gamma \Big(\eta_r\frac{\partial^2}{\partial x^2_{\rm{B}}} + \frac{\partial}{\partial x_{\rm{B}}}x_{\rm{B}}\Big) \end{eqnarray} (32) While \eta_r={\langle x_{\rm{B}}^2\rangle}_{\rm{B}} was defined in Eq.(26), \gamma=\omega^2_{\rm{B}}/\zeta assumes the diffusion constant here. The Smoluchowski limit of Eq.(24) has an extended Zusman equation form,  \begin{eqnarray}\label{exZE1} &&\frac{\partial}{\partial t}\hat\rho(x_{\rm{B}};t) + (i{\cal \mathcal{L}}_{\rm{S}}+{\hat L}_{\rm{D}})\hat\rho(x_{\rm{B}};t)= -i\Big[\alpha_1 x_{\rm{B}}+\nonumber\\ &&{}\quad\quad\alpha_2\Big(x_{\rm{B}}^2 -\eta^{2}_i \frac{\partial^2}{\partial x_{\rm{B}}^2}\Big)\Big] \big[\hat Q_{\rm{S}}, \hat\rho(x_{\rm{B}};t)\big]- 2\eta_i \frac{\partial}{\partial x_{\rm{B}}}\nonumber\\ &&{}\quad\quad \Big[\big(\frac{\alpha_1}{2}+\alpha_2 x_{\rm{B}}\big) \big\{\hat Q_{\rm{S}}, \hat\rho(x_{\rm{B}};t)\big\} \Big] \end{eqnarray} (33) This recovers the conventional Zusman equation [24-26], in the absence of the quadratic bath coupling (\alpha_2=0). We have also derived Eq.(33) via the standard Fokker-Planck-Smoluchowski approach [22]; however, the derivations are too mathematical and tedious, see Appendix for details. The present universal diffusion mapping approach with Eq.(30) and Eq.(31) is much simpler and physically more appealing. It is easy to verify that the DEOM formalism, Eq.(20), is just the dynamical resolution to the extended Zusman Eq.(33). More precisely,  $$\label{rhon_rho_x} \hat\rho_n(t)= \Big(\frac{\eta_r}{2}\Big)^{n/2} \!\!\int_{-\infty}^{\infty}\!\!{\rm d}x_{\rm{B}} H_n\Big(\frac{x_{\rm{B}}}{\sqrt{2\eta_r}}\Big)\hat\rho(x_{\rm{B}};t)$$ (34) or  \hat \rho ({x_{\rm{B}}};t) = \sum\limits_{n = 0}^\infty {{\phi _n}} ({x_{\rm{B}}}){\hat \rho _n}(t) (35) where H_n(x) is the nth-order Hermite polynomial, and  $$\label{phin} \phi_n(x_{\rm{B}})= \frac{(2\eta_r)^{-n/2} }{n!\sqrt{2\pi\eta_r}} \exp\Big(\!-\!\dfrac{x^2_{\rm{B}}}{2\eta_r}\Big)H_n\Big(\frac{x_{\rm{B}}}{\sqrt{2\eta_r}}\Big)$$ (36) We have thus validated the dissipaton algebra, Eqs.(10)-(15). This is the purpose of the above comparisons between the dissipaton approach and the present system-and-solvation composite description. The dissipaton algebra, including the new ingredient, Eq.(13) and Eq.(14), the generalized Wick's theorem with dissipatons-pairs added, is also by de facto established. Therefore, the extended DEOM Eq.(16) for general cases is also validated, due to its algebraic construction nature. Ⅳ. INTERPLAY BETWEEN DISSIPATONS AND ENVIRONMENT PARAMETERS Turn to the issue on the bath coupling parameters, \alpha_1 and \alpha_2. It is crucial to have a physical support on the nonlinear coupling bath descriptors. This issue is directly related to the extended DEOM theory, which should describe the dynamical interplay between dissipatons and nonlinear bath couplings. Erroneous descriptors of \alpha_1 and \alpha_2 would result in unphysical DEOM dynamics. In the following, we propose a polarization model to determine both \alpha_1 and \alpha_2. For clarity, we consider a chromophore system, with its ground |g\rangle and an excited |e\rangle states being engaged in optical excitations. The total system-and-bath composite Hamiltonian in the presence of external classical laser field E(t) assumes  $$\label{HT_spectrum} H_{\rm{T}}(t)=h_\textrm{g}|g\rangle\langle g|+(h_\textrm{e}+\omega_{\textrm{eg}})|e\rangle\langle e| - \hat \mu_{\rm{S}} E(t)$$ (37) with \hat \mu_{\rm{S}} = \mu(|g\rangle\langle e|+|g\rangle\langle e|) and  $$\label{delta_hB} h_\textrm{e}-h_\textrm{g} = \alpha_0 + \alpha_1 \hat x_{\rm{B}} + \alpha_2\hat x^2_{\rm{B}}$$ (38) Here, h_\textrm{g} and h_\textrm{e} denote the bath Hamiltonians associating with the ground and excited system states, respectively. Eq.(37) assumes the form of Eq.(1); i.e.,  $$\label{HT1} H_{\rm{T}}(t) = H_{\rm{S}}(t) + h_{\rm{B}} + \hat Q_{\rm{S}}(\alpha_1\hat x_{\rm{B}}+\alpha_2\hat x^2_{\rm{B}})$$ (39) with the system Hamiltonian and dissipative mode,  {H_{\rm{S}}}(t) = ({\omega _{eg}} + {\alpha _0})|e\rangle \langle e| - {\hat \mu _{\rm{S}}}E(t) (40)  {\hat Q_{\rm{S}}} = |e\rangle \langle e| (41) In Eq.(39) the bath Hamiltonian goes with h_{\rm{B}}=h_\textrm{g}. The polarization model assumes  $$\label{hB_solvA} \begin{split} h_{\textrm{g}} &= \frac{1}{2}\omega_{\rm{B}}(\hat p^2_{\rm{B}}+\hat x^2_{\rm{B}}) + h^{\rm int}_{\rm{B}}(\hat x_{\rm{B}};\widetilde{\bf q}) \\ h_{\textrm{e}} &= \frac{1}{2}\omega'_{\rm{B}}(\hat p^{\prime 2}_{\rm{B}}+\hat x^{\prime 2}_{\rm{B}}) + h^{\rm int}_{\rm{B}}(\hat x'_{\rm{B}};\widetilde{\bf q}) \end{split}$$ (42) where \hat p'_{\rm{B}}=(\omega_{\rm{B}}/\omega'_{\rm{B}})^{1/2}\hat p_{\rm{B}},  $$\label{xB_e} \hat x'_{\rm{B}} = \Big(\frac{\omega'_{\rm{B}}}{\omega_{\rm{B}}}\Big)^{1/2}(\hat x_{\rm{B}}-d_{\rm{B}})$$ (43) and  $$\label{hB_int_def} h^{\rm int}_{\rm{B}}(\hat x_{\rm{B}};\widetilde{\bf q}) =\frac{1}{2}\sum\limits_k {{{\widetilde \omega }_k}} [ \widetilde p^2_k+(\widetilde q_k-\widetilde c_k\hat x_{\rm{B}})^2 \big]$$ (44) The physical picture of this model is as follows. The system is initially in the ground |g\rangle state, with \hat x_{\rm{B}} describing its first solvation shell of frequency \omega_{\rm{B}}. Upon excitation, the system in the excited |e\rangle state experiences different solvation environment. The reorganized first-shell solvation is described with \hat x'_{\rm{B}}. It has different frequency, \omega'_{\rm{B}}, and is also linearly shifted by d_{\rm{B}}, with respective to the ground-state solvation shell. The secondary environment (\widetilde{\bf q}) remains unchanged, as described by the same h^{\rm int}_{\rm{B}}(X_{\rm{B}};\widetilde{\bf q}), for its interacting with either X_{\rm{B}}=\hat x_{\rm{B}} or \hat x'_{\rm{B}} solvation mode. Apparently, the quadractic bath coupling vanishes when \omega'_{\rm{B}}=\omega_{\rm{B}}. The coupling strength of secondary bath with the solvation mode is given by  $$\label{wti_eta} \widetilde\eta \equiv \sum\limits_k {{{\widetilde \omega }_k}} \widetilde c^2_k$$ (45) The renormalized frequencies for \hat x_{\rm{B}} and \hat x'_{\rm{B}} would be  \begin{eqnarray}\label{wti_wB} \widetilde\omega_{\rm{B}}=\omega_{\rm{B}}+\widetilde\eta\end{eqnarray} (46)  \begin{eqnarray}\widetilde\omega'_{\rm{B}} = \omega'_{\rm{B}} + \widetilde\eta \end{eqnarray} (47) respectively. To proceed, we consider first the linear--displacement--mapping (LDM) ansatz. Let h^{(1)}_{\textrm{e}} be the special case of h_{\textrm{e}} at \omega'_{\rm{B}}=\omega_{\rm{B}}. The LDM ansatz assumes h^{(1)}_{\textrm{e}}-h_\textrm{g}, due to a linear displacement d_{\rm{B}}, be the result within linear response theory.That is [9]  \begin{array}{l} h_{\rm{e}}^{(1)} - {h_{\rm{g}}} = \frac{1}{2}{\widetilde \omega _{\rm{B}}}d_{\rm{B}}^2 - {d_{\rm{B}}}({\widetilde \omega _{\rm{B}}}{{\hat x}_{\rm{B}}} - \widetilde F)\\ = \lambda - \sqrt {2\lambda {\omega _{\rm{B}}}} {{\hat x}_{\rm{B}}} \end{array} (48) Here,  \begin{eqnarray}\label{wtiF_def} \widetilde F = \sum\limits_k \widetilde\omega_k\widetilde c_k\widetilde q_k \end{eqnarray} (49) Consequently,  \begin{eqnarray}\label{lambda} \lambda = \frac{1}{2}\widetilde\omega_{\rm{B}}d^2_{\rm{B}} \end{eqnarray} (50) and also \sqrt{2\lambda\omega_{\rm{B}}} \hat x_{\rm{B}}=d_{\rm{B}}(\widetilde\omega_{\rm{B}} \hat x_{\rm{B}} - \widetilde F).Together with considering the secondary--bath--free limit, where \widetilde F=0 and \widetilde\omega_{\rm{B}} = \omega_{\rm{B}}, we obtain  \begin{eqnarray}\label{LDM_wtiX} \widetilde F = \widetilde\omega_{\rm{B}}\Big(1-\sqrt{\omega_{\rm{B}}/\widetilde\omega_{\rm{B}}}\Big)\hat x_{\rm{B}} \end{eqnarray} (51) This relates the secondary bath, \widetilde F of Eq.(49), to the solvation coordinate \hat x_{\rm{B}}. Substituting Eqs.(42)-(51) for Eq.(38), followed by some elementary algebra, we obtain  \begin{eqnarray}\label{alpha_all} \begin{array}{l} \alpha_0 = \theta\widetilde\theta\lambda\\ \alpha_1 =\displaystyle\sqrt{2\theta\lambda\omega_{\rm{B}}}\Big(1+\frac{\sqrt{\theta}\widetilde\theta-1}{r}\Big) \\[2ex] \alpha_2 =\displaystyle\frac{\widetilde\omega_{\rm{B}}}{2} \left[2(\sqrt{\theta}-1)r +(\theta\widetilde\theta-2\sqrt{\theta}+1) \right] \end{array} \end{eqnarray} (52) where  $$\label{theta_phi} \theta \equiv \frac{\omega'_{\rm{B}}}{\omega_{\rm{B}}}, \ \ \, \widetilde\theta \equiv \frac{\widetilde\omega'_{\rm{B}}}{\widetilde\omega_{\rm{B}}}, \ \ \, r\equiv \sqrt{ \frac{\omega_{\rm{B}}}{\widetilde\omega_{\rm{B}}} }$$ (53) Apparently, \alpha_2=0 if \omega'_{\rm{B}}=\omega_{\rm{B}}. Moreover, the bath coupling \alpha-parameters of Eq.(52) should go along with the underlying interplay with the dissipatons defined in Eqs.(4)-(7). Given the temperature, dissipatons are determined by \chi^{(i)}_{\rm{B}}(\omega). On the other hand, the \alpha-parameters in Eq.(52) are functions of \widetilde\eta, \omega_{\rm{B}}, \omega'_{\rm{B}} and \lambda. While the latter two are free variables, \widetilde\eta and \omega_{\rm{B}} are dictated by the same \chi^{(i)}_{\rm{B}}(\omega) that determines the dissipatons. For the secondary bath coupling strength parameter (Eq.(45)), we have [9]  $$\label{wti_eta_wtiJ} \widetilde\eta = \frac{1}{\pi}\int^{\infty}_{-\infty}\!\! {\rm d}\omega \frac{\widetilde J(\omega)}{\omega}$$ (54) Here, \widetilde J(\omega) denotes the interacting secondary bath spectral density that can be expressed in terms of [9]  $$\label{tiJw_in_chiBw} \widetilde J (\omega) = \frac{\chi^{(i)}_{\rm{B}}(\omega)}{|\chi_{\rm{B}}(\omega)|^2}$$ (55) Note that \chi_{\rm{B}}(\omega)$$\equiv$$\chi^{(r)}_{\rm{B}}(\omega)+i\chi^{(i)}_{\rm{B}}(\omega) is completely determined by either the real or imaginary part via the Kramers-Kronig relation [7-9]. For the solvation mode frequency, we have [9]  $$\label{wB_via_chi} \omega_{\rm{B}}= \frac{1}{\chi_{\rm{B}}(0)} = \frac{1}{\pi}\int^{\infty}_{-\infty}\!\! {\rm d}\omega\, \omega \chi^{(i)}_{\rm{B}}(\omega)$$ (56) Here, \chi_{\rm{B}}(0)=\chi^{(r)}_{\rm{B}}(0), as \chi^{(i)}_{\rm{B}}(\omega=0)=0. The Kramers-Kronig relation reads here [7-9]  $$\label{KK0} \chi_{\rm{B}}(0) = \frac{1}{\pi}\int^{\infty}_{-\infty}\!\! {\rm d}\omega\, \frac{\chi^{(i)}_{\rm{B}}(\omega)}{\omega}$$ (57) The above identities describe the determination of \omega_{\rm{B}} via any given \chi^{(i)}_{\rm{B}}(\omega). Actually, Eq.(56) follows in line with the definition of \chi_{\rm{B}}(\omega) in Eq.(5). For numerical illustrations, we adopt the Drude model for the secondary bath. The resultant \chi_{\rm{B}}(\omega), which supports Eqs.(54)-(57), reads  $$\label{chiBw_model} \chi_{\rm{B}}(\omega) = \frac{\omega_{\rm{B}}}{ \omega^2_{\rm{B}}-\omega^2+\dfrac{\widetilde\eta\omega_{\rm{B}}\omega}{\omega+i\widetilde\gamma} }$$ (58) We set the parameters (in unit of \omega_{\rm{B}}) \widetilde\gamma=10 and \widetilde\eta=20, at k_{\rm{B}}T=1; and also \lambda=1 for nonzero linear bath coupling strength. FIG. fig 1 depicts the evaluated linear absorption spectra, at four representing values of \theta=\omega'_{\rm{B}}/\omega_{\rm{B}}. When \theta=1 there is only the linear bath coupling. The other three (solid) curves with \theta$$\neq1 are of both the linear and the quadratic bath couplings, which are presented in parallel with their linear-free (thin) counterparts. In contrast with the pure linear bath coupling ($\theta$=1) case, the spectrum lineshape is generally asymmetric. The observed skews in individual lineshape profiles, which show non-monotonic dependence on $\theta$, are all in qualitative agreements with the secondary-bath-free but analytical results [27].

 FIG. 1 Evaluated absorption lineshapes. See text for the details.
Ⅴ. CONCLUSION

Evidently, the dissipaton algebra leads readily to the DEOM formalism. The key contribution of this work is the establishment of the generalized Wick's theorem with a pair of dissipatons added; i.e., the GWT-2, Eq.(13) and Eq.(14). The other ingredients of the dissipaton algebra had all been well established in our previous work [10, 11].

The new ingredient, Eq.(13) and Eq.(14), which is now verified unambiguously, can be used consecutively to treat further higher-order nonlinear bath couplings. For example, the GWT-3, illustrated with the dissipaton-basis-set size of one, would go with [cf.Eq.(13)]

 \begin{align} &\text{t}{{\text{r}}_{\text{B}}}[{{({{{\hat{f}}}^{n}})}^{{}^\circ }}(\hat{f}\hat{f}\hat{f}){{\rho }_{\text{T}}}]=n\langle \hat{f}\hat{f}\rangle _{\text{B}}^{>}\text{t}{{\text{r}}_{\text{B}}}[{{({{{\hat{f}}}^{n-1}})}^{{}^\circ }}(\hat{f}\hat{f}){{\rho }_{\text{T}}}]+ \\ &\text{t}{{\text{r}}_{\text{B}}}[{{({{{\hat{f}}}^{n}}\hat{f})}^{{}^\circ }}(\hat{f}\hat{f}){{\rho }_{\text{T}}}] \\ \end{align} (59)

While the first quantity is evaluated by using Eq.(13) and Eq.(14), the second quantity would be

 \begin{align} &\text{t}{{\text{r}}_{\text{B}}}[{{({{{\hat{f}}}^{n}}\hat{f})}^{{}^\circ }}(\hat{f}\hat{f}){{\rho }_{\text{T}}}]={{\langle {{{\hat{f}}}^{2}}\rangle }_{\text{B}}}\text{t}{{\text{r}}_{\text{B}}}[{{({{{\hat{f}}}^{n}})}^{{}^\circ }}\hat{f}{{\rho }_{\text{T}}}]+ \\ &n\langle \hat{f}\hat{f}\rangle _{\text{B}}^{>}\text{t}{{\text{r}}_{\text{B}}}[{{({{{\hat{f}}}^{n-1}}\hat{f})}^{{}^\circ }}\hat{f}{{\rho }_{\text{T}}}]+ \\ &\text{t}{{\text{r}}_{\text{B}}}[{{({{{\hat{f}}}^{n}}\hat{f}\hat{f})}^{{}^\circ }}\hat{f}{{\rho }_{\text{T}}}] \\ \end{align} (60)

The first two quantities are evaluated via the GWT-1 and the GWT-2 of Eq.(14), respectively. The last quantity above goes with

 \begin{align} &\text{t}{{\text{r}}_{\text{B}}}[{{({{{\hat{f}}}^{n}}\hat{f}\hat{f})}^{{}^\circ }}\hat{f}{{\rho }_{\text{T}}}]={{\rho }^{(n+3)}}+(2{{\langle {{{\hat{f}}}^{2}}\rangle }_{\text{B}}}+ \\ &n\langle \hat{f}\hat{f}\rangle _{\text{B}}^{>}){{\rho }^{(n+1)}} \\ \end{align} (61)

The GWT-3 is then completed. The GWT-$n$ follows the same recursive procedure. Thus, the present work represents a major advancement in the DEOM theory, with the specified class of non-Gaussian coupling environments that could be physically characterized, as illustrated in this work.

Ⅵ. ACKNOWLEDGEMENTS

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

APPENDIX A: The Smoluchowski limit: Conventional approach

This appendix utilizes the standard textbook approach of Ref.[22]to derive the extended Zusman Eq.(33).This is the application of the Smoluchowski limit to the Caldeira-Leggett's equation (Eq.(24)), and derive the closed equation for

 $\begin{eqnarray}\label{rhoZusman} \hat\rho(x_{\rm{B}};t)=\int_{-\infty}^\infty\!\! {\rm d} p_{\rm{B}}\, {\hat\rho}^{ \textrm{HT}}_{\textrm{W}}(x_{\rm{B}}, p_{\rm{B}};t) \end{eqnarray}$ (A1)

We start with the FP operator, Eq.(25), which has the coherent and incoherent contributions:

 $\begin{eqnarray}\label{hatL_FP_app} {\hat L}_{\rm{FP}}={\hat L}^{\textrm{coh}}_{\rm{FP}}+{\hat L}^{\textrm{incoh}}_{\rm{FP}} \end{eqnarray}$ (A2)

where

 $\begin{eqnarray}\label{hatL_coh_inc} \begin{array}{l} \displaystyle{\hat L}^{\textrm{coh}}_{\rm{FP}}= \omega_{\rm{B}}\left(\frac{\partial}{\partial x_{\rm{B}}}p_{\rm{B}} -\frac{\partial}{\partial p_{\rm{B}}}x_{\rm{B}}\right) \\[2ex] \displaystyle {\hat L}^{\textrm{incoh}}_{\rm{FP}}=-\zeta\frac{\partial}{\partial p_{\rm{B}}} \left(\frac{1}{\beta\omega_{\rm{B}}}\frac{\partial}{\partial p_{\rm{B}}}+p_{\rm{B}} \right) \end{array} \end{eqnarray}$ (A3)

It is easy to obtain [22]

 $\begin{eqnarray}\label{wtiL_FP} {\widetilde L}^{\textrm{incoh}}_{\rm{FP}} \equiv \psi^{-1}_0(p_{\rm{B}}){\hat L}^{\textrm{incoh}}_{\rm{FP}}\psi_0(p_{\rm{B}}) =\zeta \hat a^{†}\hat a \end{eqnarray}$ (A4)

where

 $\begin{eqnarray}\label{psi_0} \psi_0(p_{\rm{B}})=\left(\frac{\beta\omega_{\rm{B}}}{2\pi}\right)^{1/4} \exp\left(-\frac{\beta\omega_{\rm{B}}}{4}p_{\rm{B}}^2\right) \end{eqnarray}$ (A5)

and

 $\begin{array}{l} \hat a = \frac{{\sqrt {\beta {\omega _{\rm{B}}}} }}{2}{p_{\rm{B}}} + \frac{1}{{\sqrt {\beta {\omega _{\rm{B}}}} }}\frac{\partial }{{\partial {p_{\rm{B}}}}}\\ {{\hat a}^† } = \frac{{\sqrt {\beta {\omega _{\rm{B}}}} }}{2}{p_{\rm{B}}} - \frac{1}{{\sqrt {\beta {\omega _{\rm{B}}}} }}\frac{\partial }{{\partial {p_{\rm{B}}}}} \end{array}$ (A6)

These are the bosonic annihilation and creation operators, satisfying $[\hat a, \hat a^{†}]$=1. The normalized eigen solutions to Eq.(A4) are therefore

 $\begin{eqnarray} {\widetilde L}^{\textrm{incoh}}_{\rm{FP}}\psi_{n}(p_{\rm{B}})=n\psi_{n}(p_{\rm{B}}), \ \ \, n=0, 1, \cdots, \end{eqnarray}$ (A7)

with the ground state $\psi_{0}(p_{\rm{B}})$ of Eq.(A5) and

 $\begin{eqnarray}\label{psi_n} \hat a^{†}\psi_{n}(p_{\rm{B}})= \sqrt{n+1}\, \psi_{n+1}(p_{\rm{B}}) \end{eqnarray}$ (A8)

Apparently, $\{\psi_n(p_{\rm{B}})\}$ are all real.

On the other hand, from Eq.(A6), we evaluate

 $\begin{array}{l} \frac{{\widetilde \partial }}{{\partial {p_{\rm{B}}}}} \equiv \psi _0^{ - 1}({p_{\rm{B}}})\frac{\partial }{{\partial {p_{\rm{B}}}}}{\psi _0}({p_{\rm{B}}}) = - \sqrt {\beta {\omega _{\rm{B}}}} {{\hat a}^† }\\ {\widetilde p_{\rm{B}}} \equiv {p_{\rm{B}}} = \frac{1}{{\sqrt {\beta {\omega _{\rm{B}}}} }}(\hat a + {{\hat a}^† }) \end{array}$ (A9)

The transformed ${\widetilde L}^{\textrm{coh}}_{\rm{FP}}$$\equiv$$\psi^{-1}_0(p_{\rm{B}}) {\hat L}^{\textrm{coh}}_{\rm{FP}} \psi_0(p_{\rm{B}})$ in Eq.(A3) reads then

 $\widetilde L_{{\rm{FP}}}^{{\rm{coh}}} = {\omega _{\rm{B}}}\left[{\frac{1}{{\sqrt {\beta {\omega _{\rm{B}}}} }}\frac{\partial }{{\partial {x_{\rm{B}}}}}(\hat a + {{\hat a}^† }) + \sqrt {\beta {\omega _{\rm{B}}}} {x_{\rm{B}}}{{\hat a}^† }} \right]$ (A10)

We have also

 $\begin{eqnarray}\label{wti_partial_p2} \widetilde{\frac{\partial^2}{\partial p_{\rm{B}}^2}}=\Big(\widetilde{\frac{\partial}{\partial p_{\rm{B}}}}\Big)^2 =\beta\omega_{\rm{B}}\hat a^{†\, 2} \end{eqnarray}$ (A11)

Turn now to the transformed Caldeira-Leggett's equation (Eq.(24)),

 $\begin{array}{l} \frac{\partial }{{\partial t}}\widetilde \rho _{\rm{W}}^{{\rm{HT}}} = - [i{\mathcal{L}_{\rm{S}}} + (\widetilde L_{{\rm{FP}}}^{{\rm{coh}}} + \widetilde L_{{\rm{FP}}}^{{\rm{incoh}}})]\widetilde \rho _{\rm{W}}^{{\rm{HT}}} + (\frac{{{\alpha _1}}}{2} + \\ {\alpha _2}{x_{\rm{B}}})\{ {{\hat Q}_{\rm{S}}}, \frac{{\widetilde \partial }}{{\partial {p_{\rm{B}}}}}\widetilde \rho _{\rm{W}}^{{\rm{HT}}}\} - i({\alpha _1}{x_{\rm{B}}} + \\ {\alpha _2}x_{\rm{B}}^2 - \frac{{{\alpha _2}}}{4}\frac{{\widetilde {{\partial ^2}}}}{{\partial p_{\rm{B}}^2}})[{{\hat Q}_{\rm{S}}}, \widetilde \rho _{\rm{W}}^{{\rm{HT}}}] \end{array}$ (A12)

for

 $\begin{array}{l} \widetilde \rho _{\rm{W}}^{{\rm{HT}}}({x_{\rm{B}}}, {p_{\rm{B}}};t) \equiv \psi _0^{ - 1}({p_{\rm{B}}})\hat \rho _{\rm{W}}^{{\rm{HT}}}({x_{\rm{B}}}, {p_{\rm{B}}};t)\\ = \sum\limits_{n = 0}^\infty {{{\widetilde \rho }_n}} ({x_{\rm{B}}};t){\psi _n}({p_{\rm{B}}}) \end{array}$ (A13)

The second expression goes with the complete and orthonormal basis set of $\{\psi_n(p_{\rm{B}}); n=0, 1, \cdots\}$. Therefore,

 $\begin{eqnarray}\label{wtirhon} \widetilde\rho_n(x_{\rm{B}};t)= \int^{\infty}_{-\infty}\!\!{\rm d}p_{\rm{B}}\, \psi_n(p_{\rm{B}}){\widetilde\rho}^{\textrm{ HT}}_{\textrm{W}}(x_{\rm{B}}, p_{\rm{B}};t) \end{eqnarray}$ (A14)

Together with Eqs.(A4)-(A11), we obtain

 $\begin{array}{l} \frac{\partial }{{\partial t}}{\widetilde \rho _n} = - (i{\mathcal{L}_{\rm{S}}} + n\zeta ){\widetilde \rho _n} - i({\alpha _1}{x_{\rm{B}}} + {\alpha _2}x_{\rm{B}}^2)[{{\hat Q}_{\rm{S}}}, {\widetilde \rho _n}] - \\ {\omega _{\rm{B}}}\sqrt {n + 1} \sqrt {\beta {\omega _{\rm{B}}}} (\hat D - {x_{\rm{B}}}){\widetilde \rho _{n + 1}} - \\ {\omega _{\rm{B}}}\sqrt n \sqrt {\beta {\omega _{\rm{B}}}} \hat D{\widetilde \rho _{n - 1}} - \\ \sqrt n \sqrt {\beta {\omega _{\rm{B}}}} (\frac{{{\alpha _1}}}{2} + {\alpha _2}{x_{\rm{B}}})\{ {{\hat Q}_{\rm{S}}}, {\widetilde \rho _{n - 1}}\} + \\ i\sqrt {n(n - 1)} {\mkern 1mu} \frac{{{\alpha _2}}}{4}\beta {\omega _{\rm{B}}}[{{\hat Q}_{\rm{S}}}, {\rm{ }}{\widetilde \rho _{n-2}}] \end{array}$ (A15)

where

 $\begin{eqnarray}\label{app_hatD} \hat D\equiv (\beta\omega_{\rm{B}})^{-1} \frac{\partial}{\partial x_{\rm{B}}}+x_{\rm{B}} \end{eqnarray}$ (A16)

Consider now the Smoluchowski limit, where $\zeta$$\gg$$\omega_{\rm{B}}$.To derive a single closed equation in this limit, it requires further [22] (ⅰ) the core Hamiltonian, $H_{\rm{S}}$+ $\hat Q_{\rm{S}}(\alpha_1x_{\rm{B}}$+$\alpha_2x_{\rm{B}}^2)$, could be neglected compared to the friction, $\zeta$; (ⅱ) $\partial_t {\widetilde\rho}_{n>0}$=0; and (ⅲ) ${\widetilde\rho}_{n>2}$=0, as the bath coupling is considered up to the quadratic level, cf.Eq.(A15). Consequently, while

 $\begin{array}{l} \frac{\partial }{{\partial t}}{\widetilde \rho _0} = - i{\mathcal{L}_{\rm{S}}}{\widetilde \rho _0} - i({\alpha _1}{x_{\rm{B}}} + {\alpha _2}x_{\rm{B}}^2)[{{\hat Q}_{\rm{S}}}, {\widetilde \rho _0}] - \\ {\omega _{\rm{B}}}\sqrt {\beta {\omega _{\rm{B}}}} (\hat D - {x_{\rm{B}}}){\widetilde \rho _1} \end{array}$ (A17)

Eq.(A15) with $n=1$ and 2 becomes, respectively,

 $\begin{array}{l} 0 = - \zeta {\widetilde \rho _1} - {\omega _{\rm{B}}}\sqrt 2 \sqrt {\beta {\omega _{\rm{B}}}} (\hat D - {x_{\rm{B}}}){\widetilde \rho _2} - \\ {\omega _{\rm{B}}}\sqrt {\beta {\omega _{\rm{B}}}} \hat D{\widetilde \rho _0} - \sqrt {\beta {\omega _{\rm{B}}}} (\frac{{{\alpha _1}}}{2} + {\alpha _2}{x_{\rm{B}}})\{ {{\hat Q}_{\rm{S}}}, {\widetilde \rho _0}\} \end{array}$ (A18)
 $\begin{array}{l} 0 = - 2\zeta {\widetilde \rho _2} - {\omega _{\rm{B}}}\sqrt 2 \sqrt {\beta {\omega _{\rm{B}}}} \hat D{\widetilde \rho _1} - \sqrt 2 \sqrt {\beta {\omega _{\rm{B}}}} \cdot \\ (\frac{{{\alpha _1}}}{2} + {\alpha _2}{x_{\rm{B}}})\{ {{\hat Q}_{\rm{S}}}, {\widetilde \rho _1}\} + i\sqrt 2 {\mkern 1mu} \frac{{{\alpha _2}}}{4}\beta {\omega _{\rm{B}}}[{{\hat Q}_{\rm{S}}}, {\rho _0}] \end{array}$ (19)

The latter two, together with Eq.(A16) and $\gamma$$\equiv$$\omega_{\rm{B}}^2/\zeta$, result in

 $\begin{array}{l} \zeta {\widetilde \rho _1} = \frac{{{\omega _{\rm{B}}}}}{\zeta }\frac{\partial }{{\partial {x_{\rm{B}}}}}[{\omega _{\rm{B}}}\hat D{\widetilde \rho _1} + (\frac{{{\alpha _1}}}{2} + {\alpha _2}{x_{\rm{B}}})\{ {{\hat Q}_{\rm{S}}}, {\widetilde \rho _1}\}] - \\ i\frac{{{\alpha _2}}}{4}\frac{{\beta \gamma }}{{\sqrt {\beta {\omega _{\rm{B}}}} }}\frac{\partial }{{\partial {x_{\rm{B}}}}}[{{\hat Q}_{\rm{S}}}, {\widetilde \rho _0}] - {\omega _{\rm{B}}}\sqrt {\beta {\omega _{\rm{B}}}} \hat D{\widetilde \rho _0} - \\ \sqrt {\beta {\omega _{\rm{B}}}} (\frac{{{\alpha _1}}}{2} + {\alpha _2}{x_{\rm{B}}})\{ {{\hat Q}_{\rm{S}}}, {\widetilde \rho _0}\} \\ \approx -i\frac{{{\alpha _2}}}{4}\frac{{\beta \gamma }}{{\sqrt {\beta {\omega _{\rm{B}}}} }}\frac{\partial }{{\partial {x_{\rm{B}}}}}[{{\hat Q}_{\rm{S}}}, {\widetilde \rho _0}] - {\omega _{\rm{B}}}\sqrt {\beta {\omega _{\rm{B}}}} \hat D{\widetilde \rho _0}\\ - \sqrt {\beta {\omega _{\rm{B}}}} (\frac{{{\alpha _1}}}{2} + {\alpha _2}{x_{\rm{B}}})\{ {{\hat Q}_{\rm{S}}}, {\widetilde \rho _0}\} \end{array}$ (A20)

The second expression is obtained by considering the Smoluchowski limit where $\zeta$$\gg$$\omega_{\rm{B}}$. Note also (32) and (A16))

 $\begin{eqnarray}\label{app_LD} \hat L_{\rm{D}}=-\gamma\frac{\partial}{\partial x_{\rm{B}}}\hat D \end{eqnarray}$ (A21)

and (cf. Eq.(26))

 $\begin{eqnarray}\label{app_eta} \begin{split} \eta_r&=\langle x_{\rm{B}}^2\rangle_{\rm{B}}=(\beta \omega_{\rm{B}})^{-1} \\ \eta_i&=-\frac{\beta\gamma\langle x_{\rm{B}}^2\rangle_{\rm{B}}}{2}=-\frac{\gamma}{2\omega_{\rm{B}}} \end{split} \end{eqnarray}$ (A22)

The high-temperature relation is used here; see comments after Eq.(25).Substituting Eqs.(A20)-(A22) into Eq.(A17), followed by some simple algebra, we obtain

 $\begin{array}{l} \frac{\partial }{{\partial t}}{\widetilde \rho _0} = - (i{\mathcal{L}_{\rm{S}}} + {{\hat L}_{\rm{D}}}){\widetilde \rho _0} - i[{\alpha _1}{x_{\rm{B}}} + {\alpha _2}(x_{\rm{B}}^2-\eta _i^2\frac{{{\partial ^2}}}{{\partial x_{\rm{B}}^2}})] \cdot \\ [{{\hat Q}_{\rm{S}}}, {\widetilde \rho _0}] - 2{\eta _i}\frac{\partial }{{\partial {x_{\rm{B}}}}}\left[{(\frac{{{\alpha _1}}}{2} + {\alpha _2}{x_{\rm{B}}})\{ {{\hat Q}_{\rm{S}}}, {{\widetilde \rho }_0}\} } \right] \end{array}$ (A23)

On the other hand, from Eq.(A1), Eq.(A13), and Eq.(A14), we have

 $\begin{eqnarray} \hat\rho(x_{\rm{B}};t)=\int_{-\infty}^\infty\!\! {\rm d} p_{\rm{B}}\, \psi_0(p_{\rm{B}}){\widetilde\rho}^{\textrm{ HT}}_{\textrm{W}}(x_{\rm{B}}, p_{\rm{B}};t) =\widetilde\rho_0(x_{\rm{B}};t)\nonumber\\ \end{eqnarray}$ (A24)

Therefore, Eq.(A23) is equivalent to Eq.(33). We have thus completed the standard Fokker-Planck-Smoluchowski approach [22] to the construction of the extended Zusman equation. Apparently, the novel method of construction, on the basis of Eqs.(27)-(29) that result in the rules of diffusion mapping, Eq.(30) and Eq.(31), is much simpler and physically more ap-pealing.

Reference
 [1] R. P. Feynman, and F. L. Vernon Jr, Ann. Phys. 24 , 118 (1963). DOI:10.1016/0003-4916(63)90068-X [2] Y. Tanimura, Phys. Rev. A 41 , 6676 (1990). DOI:10.1103/PhysRevA.41.6676 [3] Y. Tanimura, J. Phys. Soc. Jpn. 75 , 082001 (2006). DOI:10.1143/JPSJ.75.082001 [4] 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 [5] 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 [6] J. S. Jin, X. Zheng, and Y. J. Yan, J. Chem. Phys. 128 , 234703 (2008). DOI:10.1063/1.2938087 [7] U. Weiss, Quantum Dissipative Systems, 3rd Edn. , Series in Modern Condensed Matter Physics, Vol. 13, Singapore: World Scientific, (2008). [8] H. Kleinert, Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets, 5th Edn. , Singapore: World Scientific, (2009). [9] Y. J. Yan, and R. X. Xu, Annu. Rev. Phys. Chem. 56 , 187 (2005). DOI:10.1146/annurev.physchem.55.091602.094425 [10] Y. J. Yan, J. Chem. Phys. 140 , 054105 (2014). DOI:10.1063/1.4863379 [11] Y. J. Yan, J. S. Jin, R. X. Xu, and X. Zheng, Frontiers Phys. 11 , 110306 (2016). DOI:10.1007/s11467-016-0513-5 [12] U. Fano, Phys. Rev. 124 , 1866 (1961). DOI:10.1103/PhysRev.124.1866 [13] A. E. Miroshnichenko, S. Flach, and Y. S. Kivshar, Rev. Mod. Phys. 82 , 2257 (2010). DOI:10.1103/RevModPhys.82.2257 [14] H. D. Zhang, R. X. Xu, X. Zheng, and Y. J. Yan, J. Chem. Phys. 142 , 024112 (2015). DOI:10.1063/1.4905494 [15] R. X. Xu, H. D. Zhang, X. Zheng, and Y. J. Yan, Sci. China Chem. 58 , 1816 (2015). DOI:10.1007/s11426-015-5499-2 [16] 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 [17] H. D. Zhang, Q. Qiao, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 145 , 204109 (2016). DOI:10.1063/1.4968031 [18] J. S. Jin, S. K. Wang, X. Zheng, and Y. J. Yan, J. Chem. Phys. 142 , 234108 (2015). DOI:10.1063/1.4922712 [19] A. O. Caldeira, and A. J. Leggett, Physica A 121 , 587 (1983). DOI:10.1016/0378-4371(83)90013-4 [20] A. Garg, J. N. Onuchic, and V. Ambegaokar, J. Chem. Phys. 83 , 4491 (1985). DOI:10.1063/1.449017 [21] M. Thoss, H. B. Wang, and W. H. Miller, J. Chem. Phys. 115 , 2991 (2001). DOI:10.1063/1.1385562 [22] H. Risken, The Fokker-Planck Equation, Methods of Solution and Applications, 2nd Edn. , Berlin: SpringerVerlag, (1989). [23] H. D. Zhang, J. Xu, R. X. Xu, and Y. J. Yan, “Modified Zusman equation for quantum solvation dynamics and rate processes, ” in Reaction Rate Constant Computations: Theories and Applications, edited by K. L. Han and T. S. Chu, Ch. 13, RSC Theoretical and Computational Chemistry Series No. 6, London, (2014). http://dx.doi.org/10.1039/9781849737753-00319. [24] L. D. Zusman, Chem. Phys. 49 , 295 (1980). DOI:10.1016/0301-0104(80)85267-0 [25] L. D. Zusman, Chem. Phys. 80 , 29 (1983). DOI:10.1016/0301-0104(83)85166-0 [26] D. Y. Yang, and R. I. Cukier, J. Chem. Phys. 91 , 281 (1989). DOI:10.1063/1.457514 [27] Y. J. Yan, and S. Mukamel, J. Chem. Phys. 85 , 5908 (1986). DOI:10.1063/1.451502