The article information
 Ruixue Xu, Xuecheng Tao, Yao Wang, Yang Liu, Houdao Zhang, YiJing Yan
 徐瑞雪, 陶雪成, 王尧, 刘阳, 张厚道, 严以京
 A HierarchicalEquationofMotion Based Semiclassical Approach to Quantum Dissipation^{†}
 基于级联方程的量子耗散半经典方法
 Chinese Journal of Chemical Physics, 2018, 31(4): 608612
 化学物理学报, 2018, 31(4): 608612
 http://dx.doi.org/10.1063/16740068/31/cjcp1807172

Article history
 Received on: July 19, 2018
 Accepted on: July 22, 2018
Quantum dissipation plays crucial roles in many fields of modern science. It describes the relaxation dynamics of a centered system to thermal equilibrium embedded in a thermal environment. Assuming Gaussian property for the environment statistics, exact approaches have been developed including the pathintegral influence functional formalism [1], its differential equivalences via the hierarchicalequationofmotion (HEOM) [27] or the related stochastic methods [812], and a dissipatonequationofmotion approach as a secondquantization advancement of the HEOM formalism [1319]. Theoretical developments have also been extended to solventinduced polarization dynamics [1416] and some classes of nonGaussian environments [12, 14, 1719].
The HEOM approach has become nowadays a standard method to simulate nonperturbative and nonMarkovian quantum dissipative dynamics. However, since its dynamical components are a set of density matrices, the HEOM method encounters low numerical efficiency in simulating multilevel systems, even with the stateoftheart filtering propagator [20] and Padé spectrum decomposition [7, 21, 22] which have largely accelerated its time evolution. For example, the HEOM will be difficult to apply if the primarily interested systems are nuclear motions along certain potentials meanwhile strongly influenced by a thermal environment at a relatively high temperature. This is because it is hard to find a proper representation in advance to describe such a system with a small size of basis set.
On the other hand, occasionally our objective is just the physical observation of dynamic variables such as coordinates or momenta, or other phasespace properties as the nuclear degrees of freedom are concerned. Equations of motion of those interested expectation values can be derived to avoid the timeconsuming simulation of density matrices evolutions. In quantum mechanics for general systems, the resulting equations of motion will be coupled equations between different order phasespace moments. This gives at common sense a semiclassical (SC) approach if the order of phasespace moments are truncated at a certain level.
In this work, we apply this strategy to HEOM and construct a SCHEOM formalism. The dynamical components in SCHEOM are then different order phasespace moments of not only the primary reduced system density operator but also the auxiliary density operators (ADOs) of HEOM. We take the hierarchical quantum master equation in Ref.[23] as an example to demonstrate the methodology. Numerical illustrations are carried out on linear spectra of anharmonic oscillators with all involved calculations preformed within a minute on a single node normal computer.
Ⅱ. PRELUDE OF HEOMIn this section, we briefly introduce the background of bath statistics and the hierarchical quantum master equation in Ref.[23]. The total Hamiltonian assumes
$ \begin{eqnarray} H_{\rm{T}}(t)=H(t)+QF_{\rm{B}}(t)\ {\rm and}\ H(t)=H_{\rm{S}}\hat\mu\epsilon(t) \end{eqnarray} $  (2.1) 
Here
$ \begin{eqnarray} C(t)\equiv\langle F_{\rm{B}}(t)F_{\rm{B}}(0)\rangle_{\rm{B}} \end{eqnarray} $  (2.2) 
which is related to the bath spectral density
$ \begin{eqnarray}\label{FDT} C(t)=\frac{1}{\pi}\int_{\infty}^{\infty}\!\!\textrm{d}\omega \textrm{e}^{i\omega t} \frac{J(\omega)}{1\textrm{e}^{\beta\omega}} \end{eqnarray} $  (2.3) 
Throughout this work we set
In this work, the SCHEOM approach will be exemplified on the hierarchical quantum master equation in Ref.[23], where the Drude bath model is assumed, i.e.
$ \begin{eqnarray} J(\omega)=\frac{2\lambda\gamma\omega}{\omega^2+\gamma^2} \end{eqnarray} $  (2.4) 
and the bosonic function is expanded as
$ \begin{eqnarray}\label{exp_3} \frac{1}{1\textrm{e}^{\beta\omega}}=\frac{1}{\beta\omega}+\frac{1}{2}+\frac{\beta\omega}{12} +{\cal O}(\beta\omega)^3 \end{eqnarray} $  (2.5) 
Using the three lowest expansion terms in Eq.(2.5), the approximated correlation function
$ \begin{eqnarray}\label{Capp} C_{\rm{app}}(t)=(c_ric_i)\textrm{e}^{\gamma t}+2\Delta\delta(t) \end{eqnarray} $  (2.6) 
with
$ \begin{eqnarray} c_r=2\lambda k_\textrm{B}T\gamma\Delta, \quad c_i=\lambda\gamma, \quad \Delta=\frac{\lambda\gamma}{6k_\textrm{B}T} \end{eqnarray} $  (2.7) 
Denote [23]
$ \begin{eqnarray}\label{cripara} \Gamma=\big[6+\sqrt{12+(\beta\gamma)^2}\, \big]/\beta, \quad \kappa=\sqrt{\Gamma/\Delta} \end{eqnarray} $  (2.8) 
and
$ \begin{eqnarray}\label{criterion} \min\{\Gamma/\Omega_{\rm{S}}, \kappa\}>5 \end{eqnarray} $  (2.9) 
With Eq.(2.6), the final form of the hierarchical quantum master equation is obtained as [23]
$ \begin{eqnarray}\label{hqme} \dot\rho_n(t)\hspace{0.12cm}&=&\hspace{0.12cm}i[H(t), \rho_n(t)]n\gamma\rho_n(t)\Delta[Q, [Q, \rho_n(t)]]\nonumber\\ &&{}\hspace{0.12cm}i\sqrt{n/c_r}\, \big(c_r[Q, \rho_{n1}(t)]ic_i\{Q, \rho_{n1}(t)\}\big)\nonumber\\ &&{}\hspace{0.12cm}i\sqrt{(n+1)c_r}\, [Q, \rho_{n+1}(t)] \end{eqnarray} $  (2.10) 
Here,
The hierarchical construction of the equations of motion resolves not only the systembath coupling strength but also the memory time scale of bath correlation. Numerically the hierarchy has to be truncated at a certain tier
We now consider a general nuclear degree of freedom along a general anharmonic potential, under an external field action and embedded in a thermal bath. The system Hamiltonian is
$ \begin{eqnarray} && H_{\rm{S}}=\frac{p^2}{2M}+V(x) \\ \ \ &&\ \ {\rm with}\ \ V(x)=\sum\limits_{j=1} v_j x^j\nonumber \end{eqnarray} $  (3.1) 
The dipolar operator
$ \begin{eqnarray} \hat\mu(x) = \sum\limits_{j=1} u_j x^j \end{eqnarray} $  (3.2) 
$ \begin{eqnarray} Q(x) = \sum\limits_{j=1} w_j x^j \end{eqnarray} $  (3.3) 
The dynamic components in the SCHEOM to be constructed are introduced as
$ \begin{eqnarray}\label{pxn} \langle x^kp^m\rangle_n\equiv{\rm tr}_{\rm{S}}(x^kp^m\rho_n) \end{eqnarray} $  (3.4) 
For clarity of the final equation, their complex conjugate counterparts related formulations would not be written explicitly in the following text.
For a general dynamic operator
$ \begin{eqnarray}\label{Ahqme} \frac{\partial}{\partial t} \langle \hat A\rangle_n \hspace{0.12cm}&=&\hspace{0.12cm}i\langle[\hat A, H_{\rm{S}}]\rangle_n+i\langle[\hat A, \hat\mu]\rangle_n\epsilon(t)\nonumber\\ & &{}\hspace{0.12cm}n\gamma \langle \hat A\rangle_n\Delta\langle[[\hat A, Q], Q]\rangle_n\nonumber\\ & &{}\hspace{0.12cm}i\sqrt{n/c_r}\, \big(c_r\langle[\hat A, Q]\rangle_{n1}ic_i\langle\{\hat A, Q\}\rangle_{n1}\big)\nonumber\\ & &{}\hspace{0.12cm}i\sqrt{(n+1)c_r}\, \langle[\hat A, Q]\rangle_{n+1} \end{eqnarray} $  (3.5) 
The main task of this work is thus to derive the equations of motion of
To derive the SCHEOM, the following relation would be needed frequently:
$ \begin{eqnarray}\label{pmxk} [x^k, p^m]=\sum\limits_{l=1}^{\min(m, k)}(i)^lC^l_mP^l_k x^{kl}p^{ml} \end{eqnarray} $  (3.6) 
Here
$ \begin{eqnarray} \frac{\partial^m}{\partial x^m}[x^kf(x)] =\sum\limits_{l=0}^{\min(m, k)}C^l_mP^l_k x^{kl}\frac{\partial^{ml}}{\partial x^{ml}}f(x) \end{eqnarray} $  (3.7) 
Together with
$ \begin{eqnarray} \{x^kp^m, x^j\}=2x^{k+j}p^mx^k[x^j, p^m] \end{eqnarray} $  (3.8) 
the final form of SCHEOM on
$ \begin{eqnarray}\label{scheom} \frac{\partial}{\partial t}\langle x^kp^m\rangle_n\hspace{0.12cm} &=&\hspace{0.12cm} \frac{i}{2M}\sum\limits_{l=1}^{\min(2, k)}(i)^lC^l_2P^l_k\langle x^{kl}p^{2+ml}\rangle_n\nonumber\\ & &{}\hspace{0.12cm}n\gamma \langle x^kp^m\rangle_ni\sum\limits_j[v_ju_j\epsilon(t)]\times\nonumber\\ & &{}\hspace{0.12cm} \sum\limits_{l=1}^{\min(m, j)}(i)^lC^l_mP^l_j\langle x^{k+jl}p^{ml}\rangle_n\nonumber\\ & &{}\hspace{0.12cm} \Delta\sum\limits_{jj'}w_jw_{j'}\sum\limits_{l=1}^{\min(m, j)}\sum\limits_{l'=1}^{\min(ml, j')}(i)^{l+l'}\times\nonumber\\ & &{}\hspace{0.12cm} C^l_mP^l_jC^{l'}_{ml}P^{l'}_{j'}\langle x^{k+j+j'll'}p^{mll'}\rangle_n\nonumber\\ & &{}\hspace{0.12cm}i\sqrt{(n+1)c_r}\, \sum\limits_jw_j\times\nonumber\\ & &{}\hspace{0.12cm} \sum\limits_{l=1}^{\min(m, j)}(i)^lC^l_mP^l_j\langle x^{k+jl}p^{ml}\rangle_{n+1}\nonumber\\ & &{}\hspace{0.12cm} i\sqrt{n/c_r}\, (c_ric_i)\sum\limits_jw_j\times\nonumber\\ & &{}\hspace{0.12cm} \sum\limits_{l=1}^{\min(m, j)}(i)^lC^l_mP^l_j\langle x^{k+jl}p^{ml}\rangle_{n1}\nonumber\\ & &{}\hspace{0.12cm}2c_i\sqrt{n/c_r}\, \sum\limits_jw_j\langle x^{k+j}p^{m}\rangle_{n1} \end{eqnarray} $  (3.9) 
The total dimension of this equation is
For application, we consider the linear spectrum of an anharmonic oscillator, evaluated via the linear response approach on Eq.(3.9). To the first order of the external field, the linear response relation is written as
$ \begin{eqnarray} \delta P^{(1)}(t)\hspace{0.12cm} &=&\hspace{0.12cm} {\rm tr}_{\rm{S}}{\rm tr}_{\rm{B}}\{\hat\mu[\rho^{(1)}_{\rm{tot}}(t)\rho_{\rm{tot}}^{\rm eq}]\}\nonumber\\ \hspace{0.12cm} &=&\hspace{0.12cm} \langle\hat\mu\rangle^{(1)}_0\langle\hat\mu\rangle_{\rm eq}\nonumber\\ \hspace{0.12cm}&\equiv&\hspace{0.12cm} \int_{\infty}^t\!\!\textrm{d}\tau\, \chi(t\tau)\epsilon(\tau)\quad\quad \end{eqnarray} $  (4.1) 
Here,
To derive the perturbative solution of Eq.(3.9) with respect to the external field, we recast Eq.(3.9) in the matrixvector form as
$ \begin{eqnarray}\label{mateom} \dot{\mathbf X}(t)=[{\mathbf\Lambda}+{\mathbf D}\epsilon(t)]{\mathbf X}(t) \end{eqnarray} $  (4.2) 
with
$ \begin{eqnarray} {\mathbf X}(t)=\textrm{e}^{{\mathbf\Lambda}t}{\mathbf X}_I(t) \end{eqnarray} $  (4.3) 
we have
$ \begin{eqnarray} \dot{\mathbf X}_I(t)=\textrm{e}^{{\mathbf\Lambda}t}{\mathbf D}\textrm{e}^{{\mathbf\Lambda}t}{\mathbf X}_I(t)\epsilon(t) \end{eqnarray} $  (4.4) 
which leads to
$\begin{eqnarray}\label{linearresp} \delta{\mathbf X}^{(1)}(t) = \int_{\infty}^t\textrm{d}\tau \textrm{e}^{{\mathbf\Lambda}(t\tau)}{\mathbf D}{\mathbf X}_{\rm eq}\epsilon(\tau)\nonumber\\ \equiv \int_{\infty}^t\textrm{d}\tau{\chi}(t\tau)\epsilon(\tau) \end{eqnarray} $  (4.5) 
For numerical demonstration, the potential of the anharmonic oscillator is assumed as
$ \begin{eqnarray} V(x)=V_\textrm{e}\Big[\Big(\frac{x}{a}\Big)^2\theta\Big(\frac{x}{a}\Big)^3\Big] \end{eqnarray} $  (4.6) 
The parameters are chosen as
In this work, we propose a SCHEOM approach, which is constructed on the basis of the wellestablished HEOM formalism, with its dynamic components being phasespace moments instead of the density matrices. The objective is to achieve the numerical executability to simulate nonperturbative and nonMarkovian quantum dissipative dynamics of anharmonic nuclear systems strongly influenced by a thermal bath at a relatively high temperature. Such systems would require practically a rather large basis set in representation leading to the densitymatrix based methods pretty difficult to perform.
For clarity of the SCHEOM formalism and its construction, the present work is exemplified with the hierarchical quantum master equation of Eq.(2.10) for the case of single dissipative mode with single exponential term bath correlation. The numerical demonstration is made on the linear response of a nuclear oscillator with anharmonicity. More generalized SCHEOM construction and more complex applications will be carried out in future work. Also the truncation on higherorder phasespace moments need further elaboration considering related physical inherent properties.
Ⅵ. ACKNOWLEDGEMENTSThis work was supported by the National Natural Science Foundation of China (No.21373191, No.21573202, No.21633006, and No.21703225), and the Fundamental Research Funds for the Central Universities (No.2030020028, No.2060030025, and No.2340000074).
[1]  R. P. Feynman, and F. L. Vernon Jr, Ann. Phys. 24 , 118 (1963). DOI:10.1016/00034916(63)90068X 
[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]  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 
[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. S. Jin, X. Zheng, and Y. J. Yan, J. Chem. Phys. 128 , 234703 (2008). DOI:10.1063/1.2938087 
[7]  J. J. Ding, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 136 , 224103 (2012). DOI:10.1063/1.4724193 
[8]  J. S. Shao, J. Chem. Phys. 120 , 5053 (2004). DOI:10.1063/1.1647528 
[9]  J. M. Moix, and J. S. Cao, J. Chem. Phys. 139 , 134106 (2013). DOI:10.1063/1.4822043 
[10]  Y. L. Ke, and Y. Zhao, J. Chem. Phys. 145 , 024101 (2016). DOI:10.1063/1.4955107 
[11]  C. Y. Hsieh, and J. Cao, J. Chem. Phys. 148 , 014103 (2018). DOI:10.1063/1.5018725 
[12]  C. Y. Hsieh, and J. Cao, J. Chem. Phys. 148 , 014104 (2018). DOI:10.1063/1.5018726 
[13]  Y. J. Yan, J. Chem. Phys. 140 , 054105 (2014). DOI:10.1063/1.4863379 
[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/s1142601554992 
[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]  R. X. Xu, Y. Liu, H. D. Zhang, and Y. J. Yan, Chin. J. Chem. Phys. 30 , 395 (2017). DOI:10.1063/16740068/30/cjcp1706123 
[18]  R. X. Xu, Y. Liu, H. D. Zhang, and Y. J. Yan, J. Chem. Phys. 148 , 114103 (2018). DOI:10.1063/1.4991779 
[19]  Y. Liu, R. X. Xu, H. D. Zhang, and Y. J. Yan, Chin. J. Chem. Phys. 31 , 245 (2018). DOI:10.1063/16740068/31/cjcp1804083 
[20]  Q. Shi, L. P. Chen, G. J. Nan, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 130 , 084105 (2009). DOI:10.1063/1.3077918 
[21]  J. Hu, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 133 , 101106 (2010). DOI:10.1063/1.3484491 
[22]  J. Hu, M. Luo, F. Jiang, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 134 , 244106 (2011). DOI:10.1063/1.3602466 
[23]  R. X. Xu, B. L. Tian, J. Xu, Q. Shi, and Y. J. Yan, J. Chem. Phys. 131 , 214111 (2009). DOI:10.1063/1.3268922 