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

#### The article information

Haobin Wang, Xinzijian Liu, Jian Liu

Accurate Calculation of Equilibrium Reduced Density Matrix for the System-Bath Model: a Multilayer Multiconfiguration Time-Dependent Hartree Approach and its Comparison to a Multi-Electronic-State Path Integral Molecular Dynamics Approach

Chinese Journal of Chemical Physics, 2018, 31(4): 446-456

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

### Article history

Accepted on: June 15, 2018
Accurate Calculation of Equilibrium Reduced Density Matrix for the System-Bath Model: a Multilayer Multiconfiguration Time-Dependent Hartree Approach and its Comparison to a Multi-Electronic-State Path Integral Molecular Dynamics Approach
Haobin Wanga, Xinzijian Liub, Jian Liub
Dated: Received on May 29, 2018; Accepted on June 15, 2018
a. Department of Chemistry, University of Colorado Denver, Denver CO 80217-3364, USA;
b. Beijing National Laboratory for Molecular Sciences, Institute of Theoretical and Computational Chemistry, College of Chemistry and Molecular Engineering, Peking University, Beijing 100871, China
*Author to whom correspondence should be addressed. Haobin Wang, E-mail:haobin.wang@ucdenver.edu; Jian Liu, E-mail:jianliupku@pku.edu.cn
Abstract: An efficient and accurate method for computing the equilibrium reduced density matrix is presented for treating open quantum systems characterized by the system-bath model. The method employs the multilayer multiconfiguration time-dependent Hartree theory for imaginary time propagation and an importance sampling procedure for calculating the quantum mechanical trace. The method is applied to the spin-boson Hamiltonian, which leads to accurate results in agreement with those produced by the multi-electronic-state path integral molecular dynamics method.
Key words: Multilayer multiconfiguration time-dependent Hartree    Path integral    Equilibrium reduced density matrix    Imaginary time propagation
Ⅰ. INTRODUCTION

A challenging issue in theoretical chemistry is the accurate description of many-body quantum effects for complex systems. Some methods are "numerically exact" in principle, e.g., approaches that solve Schrödinger equations or quantum Liouville equations, or approaches employing Feynman path integral formalism. However, not all these methods give reliable results in practice. The primary difficulty is the rapid growth of the size of the Hilbert space (or Liouville space) as the number of degrees of freedom increases, which renders a brute-force approach (complete quadrature or full configuration-interaction) impractical. One thus needs to seek some smart algorithms to efficiently explore the relevant subspace in which the solution, though approximate, essentially converges to the true result. If a method fails to find (all or part of) such subspace, it will not give correct results even if this algorithm is based on an exact principle. Therefore, the development and testing of reliable numerically exact methods is essential to obtaining truly accurate results.

The class of problems we would like to investigate in this work is represented by the system-bath type of Hamiltonian, e.g. many models in condensed phase physics and chemistry for describing open quantum systems. Techniques developed to solve this class of problems include, but are not limited to, the numerical path integral approach [1-5], the numerical renormalization group method [6-8], the hierarchical equation of motion method [9-14], and the multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) theory [15-17]. Different from approximate approaches such as quantum perturbation theories or semiclassical methods, numerically accurate simulations on these model systems serve the same role as experiments, which is often the first step in developing a reliable, general theory for describing the underlying physical processes.

In a system-bath model, the latter represents the environment that influences the dynamics or thermodynamics of the former. Within linear response approximation the environment can often be modeled by a harmonic bath that is coupled linearly to the system [18, 19]. This makes it possible to integrate out the bath completely and concentrate only on the system dynamics. The price paid for this reduction is that the influence functionals or relevant kernels are nonlocal in time, which makes the study prohibitively expensive when the system becomes large. A complementary strategy is to represent the bath by discrete degrees of freedom, i.e., the same way as in a discrete numerical quadrature. One converges the dynamics/thermodynamics for this finite system-bath composite system by increasing the number of bath modes. The ML-MCTDH theory [15-17] employed in this work adopts such a strategy. As an added bonus, methods that belong to this group are often applicable to treating anharmonic baths and other types of potential energy functions. The practicality of the ML-MCTDH method lies on its efficient exploration of the Hilbert subspace, which is achieved by using a flexible tensor contraction scheme (vide infra).

The purpose of this paper is to present our implementation of the ML-MCTDH theory to evaluate equilibrium reduced density matrix (ERDM) for the system-bath model and compare the results to those obtained by the multi-electronic-state path integral molecular dynamics (MES-PIMD) approach [20]. The ML-MCTDH work is a generalization of the previous work on calculating energy eigenstates and other thermodynamic quantities [21-25]. In Section Ⅱ we will summarize the ML-MCTDH theory, followed by a discussion on practical calculations of ERDM in Section Ⅲ. Then we will provide in Section Ⅳ a brief introduction to the recently developed MES-PIMD approach [20]. In Section Ⅴ we will present examples on the spin-boson model. Finally, we conclude with an outlook to future extensions.

Ⅱ. VARIATIONAL APPROACH TO QUANTUM DYNAMICS: FROM CONVENTIONAL METHOD TO ML-MCTDH A. Conventional approach using static configurations

The numerical solution of the time-dependent Schrödinger equation using basis sets can be achieved by employing Dirac-Frenkel variational principle [26] (in this paper we use atomic units in which $\hbar$=1)

 $\begin{eqnarray} \label{dirac} \langle \delta\Psi(t) | i \frac{\partial} {\partial t} - \hat{H} | \Psi(t) \rangle = 0 \end{eqnarray}$ (2.1)

In conventional methods the wave function is expressed as a linear combination of time-independent configurations

 $\begin{eqnarray} \label{expan} |\Psi (t) \rangle &=& \sum\limits_{j_1} \sum\limits_{j_2} ... \sum\limits_{j_f} A_{j_1j_2...j_f}(t) \prod\limits_{\nu=1}^{f} |\phi_{j_\nu}^{(\nu)} \rangle \nonumber\\ &\equiv& \sum\limits_{J} A_J(t) |\Phi_J \rangle \end{eqnarray}$ (2.2)

where $J$ is a multi-dimensional index that runs through all the combinations of basis functions $\phi^{(\nu)}$ over $\nu$=1, 2, $..., f$-degrees of freedom. If for each degree of freedom $\nu$ the basis functions $\phi_n^{(\nu)}$'s are orthonormal, the equations of motion are obtained upon substitution of Eq.(2.1)

 $\begin{eqnarray} \label{fci} i \dot{A}_J(t) = \langle \Phi_J | \hat{H} | \Psi(t) \rangle = \sum\limits_L \langle \Phi_J | \hat{H} | \Phi_L \rangle A_L(t) \end{eqnarray}$ (2.3)

In this brute-force, full configuration interaction (CI) approach, the number of expansion terms grows exponentially with respect to the number of degrees of freedom.

The limitation of the simple expansion scheme in Eq.(2.2) is due to its unbiased exploration of the entire Hilbert space, which can be demonstrated by considering a special case where $f$ degrees of freedom are completely separable. In this case the solution of the time-dependent Schrödinger equation is given by a direct or Hartree product, $|\Psi(t) \rangle$=$\prod_{\nu=1}^{f} |\varphi^{(\nu)} (t) \rangle$. However, even for this uncorrelated case, the expansion in Eq.(2.2) still requires exponentially many terms, suggesting that most parameters are redundant. That is, the fact that tensor $A$ with $n^f$ terms in the full CI expansion can be contracted to $n$$\times$$f$ terms is completely ignored in Eq.(2.2).

Thus, it is natural to consider a more flexible approach to describing the dynamics in a reduced, but nonetheless numerically converged, parameter space. One may start from the single Hartree description for the completely uncorrelated dynamics above and gradually build in more time-dependent configurations as correlation becomes stronger. This is the multiconfiguration time-dependent Hartree (MCTDH) approach proposed by Meyer, Manthe, and Cederbaum [27-30].

B. The multiconfiguration time-dependent Hartree method

Within the MCTDH method the wave function has the form [27-30]

 $\begin{eqnarray} \label{newmc} |\Psi (t) \rangle = \sum\limits_{j_1} \sum\limits_{j_2} ... \sum\limits_{j_p} A_{j_1j_2...j_p}(t) \prod\limits_{\kappa=1}^{p} |\varphi_{j_\kappa}^{(\kappa)} (t) \rangle \end{eqnarray}$ (2.4)

where the summation is over all combinations of the single particle functions (SPFs) $|\varphi_{j_\kappa}^{(\kappa)} (t) \rangle$ for $\kappa$=1, 2, $..., p$. This tensor decomposition is usually called the Tucker form [31] or N-way singular value decomposition (SVD) in mathematics, which finds a wide range of applications [32] apart from MCTDH. One may impose the orthogonality condition

 $\begin{eqnarray} \label{ortho2} \langle \varphi_n^{(\kappa)} (t) | \varphi_m^{(\kappa)} (t) \rangle = \delta_{nm} \end{eqnarray}$ (2.5)

via the standard gauge condition

 $\langle \varphi_n^{(\kappa)} (0) | \varphi_m^{(\kappa)} (0) \rangle \hspace{-0.15cm} = \hspace{-0.15cm} \delta_{nm}$ (2.6a)
 $\begin{array}{l} \left\langle {\varphi _n^{(\kappa )}(t)|\frac{\partial }{{\partial t}}\varphi _m^{(\kappa )}(t)} \right\rangle \equiv \left\langle {\varphi _n^{(\kappa )}(t)| \dot \varphi _m^{(\kappa )}(t)} \right\rangle \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;= 0 \end{array}$ (2.6b)

Variations then lead to two sets of coupled MCTDH equations of motion [27-30]

 $i \dot{A}_J(t) \hspace{-0.15cm} = \hspace{-0.15cm} \left\langle \Phi_J(t) \left| \hat{H} \right| \Psi(t) \right\rangle\\ \ \ \ \ \ \ \ \ \ \ = \sum\limits_L \left\langle \Phi_J(t) \left| \hat{H} \right| \Phi_L(t) \right\rangle A_L(t)$ (2.7a)
 $i |\underline{\dot{\varphi}}^{(\kappa)}(t)\rangle \hspace{-0.15cm} = \hspace{-0.15cm} [1-\hat{P}^{(\kappa)}(t)]\nonumber\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [\hat{\rho}^{(\kappa)}(t)]^{-1} \langle \hat{H} \rangle^{(\kappa)} (t) |\underline{\varphi}^{(\kappa)}(t)\rangle$ (2.7b)

Eq.(2.7a) bears some similarity to Eq.(2.3), except that the Hamiltonian matrix $H_{JL}(t) $$\equiv$$\langle \Phi_J(t) | \hat{H} | \Phi_L(t) \rangle$ is now time-dependent. Eq.(2.7b) describes the time-dependence of the SPFs that does not appear in the conventional method. Other quantities are defined as follows: $|\underline{\varphi}^{(\kappa)}(t)\rangle$=$\{|\varphi_1^{(\kappa)}(t)\rangle, |\varphi_2^{(\kappa)}(t)\rangle, ...\}^T$ denotes the symbolic column vector of the SPFs for the $\kappa{\rm th}$ single particle group, $\langle \hat{H}\rangle^{(\kappa)} (t)$ denotes the mean-field operator, $\hat{\rho}^{(\kappa)}(t)$ is the reduced density matrix,

 $\langle \hat{H} \rangle^{(\kappa)}_{nm}(t) \hspace{-0.15cm} = \hspace{-0.15cm} \langle \Psi_n^{(\kappa)}(t) | \hat{H} | \Psi_m^{(\kappa)}(t) \rangle$ (2.8a)
 $\rho^{(\kappa)}_{nm} (t) \hspace{-0.15cm} = \hspace{-0.15cm}\langle \Psi_n^{(\kappa)}(t) | \Psi_m^{(\kappa)}(t) \rangle$ (2.8b)

and $P^{(\kappa)}(t)$ is the single particle space projector for the $\kappa{\rm th}$ single particle group,

 $\begin{eqnarray} P^{(\kappa)}(t) = \sum\limits_m |\varphi_m^{(\kappa)}(t)\rangle\langle \varphi_m^{(\kappa)}(t)| \end{eqnarray}$ (2.9)

In the expressions above, the single hole function $|\Psi_n^{(\kappa)}(t)\rangle$ for the $\kappa{\rm th}$ single particle group is defined as

 $\begin{eqnarray} |\Psi_n^{(\kappa)}(t)\rangle &=&\hspace{-0.15cm}\sum\limits_{j_1}...\sum\limits_{j_{\kappa-1}}\sum\limits_{j_{\kappa+1}}... \sum\limits_{j_p} \Big(A_{j_1...j_{k-1}nj_{k+1}...j_p}(t) \;\nonumber\\ && \prod\limits_{\lambda \neq \kappa }^{p} |\varphi_{j_\lambda}^{(\lambda)} (t) \rangle\Big) \end{eqnarray}$ (2.10a)

so that

 $\begin{eqnarray} |\Psi(t)\rangle = \sum\limits_n |\varphi_n^{(\kappa)}(t)\rangle |\Psi_n^{(\kappa)}(t)\rangle \end{eqnarray}$ (2.10b)

The expansion in the MCTDH approach, Eq.(2.4), resembles the full CI expression in Eq.(2.2), i.e., the total number of time-dependent configurations still scales exponentially versus the total number of single particle groups $p$. However the MCTDH method is applicable to more complex systems for two reasons: (ⅰ) the base of the exponential in the MCTDH approach, i.e. the number of physically important SPFs, is usually much smaller than the number of time-independent basis functions in the conventional approach; and (ⅱ) each single particle group may contain several physical degrees of freedom so that the number of the single particle groups ($p$) is usually much less than the number of physical degrees of freedom ($f$). The method recovers the single Hartree limit naturally for separable systems and systematically requires more configurations as correlation effects become more important.

The main limitation of the MCTDH approach outlined above lies in its way of constructing the SPFs, which is based on another full CI expansion employing the static Hartree products in the single particle subspace

 $\begin{eqnarray} \label{basism} |\varphi_n^{(\kappa)}(t)\rangle \hspace{-0.15cm} &=&\hspace{-0.15cm} \sum\limits_I B_{I}^{\kappa, n}(t) |u_I^{\kappa}\rangle \nonumber\\ &\equiv& \sum\limits_{i_1} \sum\limits_{i_2} ... \sum\limits_{i_{F(\kappa)}} B_{i_1i_2...i_{F(\kappa)}}^{\kappa, n}(t) \prod\limits_{q=1}^{F(\kappa)} |\phi_{i_q}^{\kappa, q} \rangle\nonumber\\ \end{eqnarray}$ (2.11)

The ML-MCTDH theory [15] successfully circumvents this bottleneck by using a hierarchical, dynamic contraction of the basis functions that constitute the original SPFs. This gives a flexible representation of the overall wave function and significantly increases the number of degrees that can be treated.

C. Multilayer multiconfiguration time-dependent Hartree theory

To begin with, the SPFs $|\varphi_{j_\kappa}^{(\kappa)} (t) \rangle$ in MCTDH can again be cast into a time-dependent multiconfiguration expansion

 $\begin{eqnarray} \label{L2SP} |\varphi_n^{(\kappa)}(t)\rangle \hspace{-0.15cm}&=&\hspace{-0.15cm} \sum\limits_I B_{I}^{\kappa, n}(t)\; |u_I^{\kappa}(t)\rangle\nonumber \\ &\equiv &\hspace{-0.15cm} \sum\limits_{i_1} \sum\limits_{i_2} ... \sum\limits_{i_{Q(\kappa)}} B_{i_1i_2...i_{Q(\kappa)}}^{\kappa, n}(t) \prod\limits_{q=1}^{Q(\kappa)} |v_{i_q}^{(\kappa, q)}(t) \rangle\nonumber \\ \end{eqnarray}$ (2.12)

For notation purpose we refer to the single particles introduced in the previous MCTDH section as the level one (L1), which in turn contains several level two (L2) single particles. An L1-SPF $|\varphi_n^{(\kappa)}(t)\rangle$ is expanded in the time-dependent L2-SPFs as in Eq.(2.12). The expansion of a two-layer ML-MCTDH wave function can thus be written in the form

 $\begin{array}{l} |\Psi (t)\rangle = \sum\limits_{{j_1}} {\sum\limits_{{j_2}} . } ..\sum\limits_{{j_p}} {{A_{{j_1}{j_2}...{j_p}}}} (t)\prod\limits_{\kappa = 1}^p {\left[ {\sum\limits_{{i_1}} {\sum\limits_{{i_2}} {...} } } \right.} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{{i_{Q(\kappa )}}} {B_{{i_1}{i_2}...{i_{Q(\kappa )}}}^{\kappa ,{j_\kappa }}} (t)\left. {\prod\limits_{q = 1}^{Q(\kappa )} | v_{{i_q}}^{(\kappa ,q)}(t)\rangle } \right] \end{array}$ (2.13)

for which Dirac-Frenkel variational principle can be split into

 $\langle \delta \Psi(t) | i \frac{\partial} {\partial t} - \hat{H} | \Psi(t) \rangle_{\rm top\; coefficients} = 0$ (2.14a)
 $\langle\delta \Psi(t) | i \frac{\partial}{\partial t} - \hat{H} | \Psi(t) \rangle_{\rm L1\; SPFs} = 0$ (2.14b)
 $\langle\delta \Psi(t) | i \frac{\partial}{\partial t} - \hat{H} | \Psi(t) \rangle_{\rm L2\; SPFs} = 0$ (2.14c)

The first two parts give equations that are in the same form as Eqs. (2.7a) and (2.7b),

 $i \dot{A}_J(t) = \sum\limits_L \left\langle \Phi_J(t) \left| \hat{H} \right| \Phi_L(t) \right\rangle A_L(t)$ (2.15a)
 $i |\underline{\dot{\varphi}}^{(\kappa)}(t) \rangle_{\rm L2\; coefficients} = [1-\hat{P}^{(\kappa)}(t) ][\hat{\rho}^{(\kappa)}(t)]^{-1} \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \langle \hat{H} \rangle^{(\kappa)} (t) |\underline{\varphi}^{(\kappa)}(t) \rangle$ (2.15b)

Here, the symbolic notation on the left hand side of Eq.(2.15b) means that the time derivative of the L1-SPFs is only taken with respect to the L2 expansion coefficients $B_{I}^{\kappa, n}$ and does not act on the L2 configuration $|u_I^{\kappa}(t)\rangle$=$\displaystyle{\prod_{q=1}^{Q(\kappa)} } |v_{i_q}^{(\kappa, q)}(t)\rangle$ in Eq.(2.12). In ML-MCTDH the configurations $|\Phi_J(t)\rangle$, the Hamiltonian matrix $H_{JL}(t)$=$\langle \Phi_J(t) | \hat{H} | \Phi_L(t) \rangle$, the L1 mean-field operator $\langle \hat{H} \rangle^k (t)$, and the L1-SPFs $|\underline{\varphi}^{(\kappa)}(t)\rangle$, all depend on the L2-SPFs $|v_{i_q}^{(\kappa, q)}(t) \rangle$. These quantities need to be built explicitly from the bottom layer SPFs and basis functions.

Variation with respect to the L2-SPFs leads to an expression similar to Eq.(2.15b)

 $\begin{eqnarray} \label{L2or} i |\underline{\dot{v}}^{(\kappa, q)}(t) \rangle &=&\hspace{-0.15cm} [1-\hat{P}^{(\kappa, q)}_{\rm L2}(t)] [\hat{\varrho}^{(\kappa, q)}(t)]^{-1} \nonumber\\&&\hspace{-0.15cm}\langle \hat{\mathcal{H}}\rangle^{(\kappa, q)} (t) |\underline{v}^{(\kappa, q)}(t) \rangle \end{eqnarray}$ (2.16)

Again $|\underline{v}^{(\kappa, q)}(t) \rangle$=$\{|v_1^{(\kappa, q)}(t)\rangle, \; |v_2^{(\kappa, q)}(t)\rangle, ...\}^T$ denotes the symbolic column vector of (the coefficients of) the L2-SPFs and \hat{P}_{\rm L2}^{(\kappa, q)}(t) $$\equiv \displaystyle{\sum_l} |v_l^{(\kappa, q)}(t)\rangle \langle v_l^{(\kappa, q)}(t) | the projection operator in the L2-SP space. The second layer reduced density matrix and mean-field operator are given with the aid of the second layer hole functions defined by  \begin{eqnarray} \label{L2hole} |\varphi_n^k(t)\rangle = \sum\limits_r |v_{r}^{(\kappa, q)}(t) \rangle \; |g_{n, r}^{(\kappa, q)}(t)\rangle \end{eqnarray} (2.17) and the first layer reduced density matrix and mean-field operator in Eq.(2.8). The expressions are  \varrho_{rs}^{(\kappa, q)}(t)\hspace{-0.15cm} = \hspace{-0.15cm} \sum\limits_n \sum\limits_m \rho^{(\kappa)}_{nm}(t) \langle g_{n, r}^{(\kappa, q)}(t) | g_{m, s}^{(\kappa, q)}(t) \rangle (2.18a)  \langle \hat{\mathcal{H}} \rangle^{(\kappa, q)}_{rs} (t) \hspace{-0.15cm} = \hspace{-0.15cm} \sum\limits_n \sum\limits_m \langle g_{n, r}^{(\kappa, q)}(t) | \langle \hat{H} \rangle^{(\kappa)}_{nm} (t) | g_{m, s}^{(\kappa, q)}(t) \rangle (2.18b) The two-layer ML-MCTDH, Eqs. (2.15a), (2.15b), and (2.16), can be recursively generalized to include an arbitrary number of layers by expanding the overall wave function via a hierarchical tensor contraction  |\Psi (t) \rangle\hspace{-0.15cm} = \hspace{-0.15cm} \sum\limits_{j_1} \sum\limits_{j_2} ... \sum\limits_{j_p} A_{j_1j_2...j_p}(t) \prod\limits_{\kappa=1}^{p} |\varphi_{j_\kappa}^{(\kappa)} (t) \rangle (2.19a)  |\varphi_{j_\kappa}^{(\kappa)}(t)\rangle \hspace{-0.15cm} = \hspace{-0.15cm} \sum\limits_{i_1} \sum\limits_{i_2} ... \sum\limits_{i_{Q(\kappa)}} \Big( B_{i_1i_2...i_{Q(\kappa)}}^{\kappa, j_\kappa}(t)\cdot \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \prod\limits_{q=1}^{Q(\kappa)} |v_{i_q}^{(\kappa, q)}(t) \rangle \Big) (2.19b)  |v_{i_q}^{(\kappa, q)}(t)\rangle = \sum\limits_{\alpha_1} \sum\limits_{\alpha_2} ... \sum\limits_{\alpha_{M(\kappa, q)}} \Big( C_{\alpha_1\alpha_2...\alpha_{M(\kappa, q)}}^{\kappa, q, i_q}(t)\cdot \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \prod\limits_{\gamma=1}^{M(\kappa, q)} |\xi_{\alpha_\gamma}^{(\kappa, q, \gamma)}(t) \rangle\Big)\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdots (2.19c) and the equations of motion can be derived by the same way as above for the two layer case using mathematical induction and they all have the same form [15-17]. The multilayer hierarchy is terminated at a particular level by expanding the SPFs in the deepest layer in terms of time-independent configurations. The introduction of the recursive, dynamically optimized layering scheme in the ML-MCTDH wave function provides a great deal of flexibility in the trial wave function, which results in a tremendous gain in one's ability to study large many-body quantum systems. This is demonstrated by many applications on simulating quantum dynamics of ultrafast electron transfer reactions in condensed phases [33-49]. The ML-MCTDH work of Manthe has introduced an even more adaptive formulation based on a layered correlation discrete variable representation [21, 50]. The form of the ML-MCTDH wave function in Eq.(2.19) has also received much attention recently in applied mathematics [51]. Ⅲ. CALCULATING EQUILIBRIUM REDUCED DENSITY MATRIX FOR THE SYSTEM-BATH MODEL USING ML-MCTDH The applicability of ML-MCTDH theory is not limited to real time quantum dynamics. By a trivial extension to imaginary time \tau =- it, ML-MCTDH can be used to evaluate the Boltzmann operator [35], thermal rate constant [52, 53], or eigenstates of the Hamiltonian [25]. Here we describe a procedure for calculating equilibrium reduced density matrix (ERDM) for the system-bath model. Starting from the generic Hamiltonian  \begin{eqnarray} H = H_{\rm S} + H_{\rm B} + H_{{\rm SB}} \end{eqnarray} (3.1) where H_{\rm S} is the Hamiltonian for the system, H_{\rm B} for the bath, and H_{{\rm SB}} for the system-bath coupling. ERDM \hat{\rho}^\beta is defined by tracing out the bath from the Boltzmann operator  \rho^\beta_{nm} \hspace{-0.15cm} = \hspace{-0.15cm} \frac{1}{Z_{\rm B}} {\rm tr} [{\rm e}^{-\beta \hat{H}} |n\rangle\langle m|] (3.2a)  Z_{\rm B} \hspace{-0.15cm} = \hspace{-0.15cm} {\rm tr} [{\rm e}^{-\beta \hat{H}_{\rm B}}] (3.2b) Here \beta =1/k_{\rm B}T, and |n\rangle denotes the basis states for the system, i.e.  \hat{\rho}^\beta\hspace{-0.15cm} \equiv \hspace{-0.15cm} \sum\limits_n \sum\limits_m |n\rangle \rho^\beta_{nm}\langle m| (3.3a)  \hat{H}_{\rm S}\hspace{-0.15cm} \equiv \hspace{-0.15cm} \sum\limits_n \sum\limits_m |n\rangle \epsilon_{nm}\langle m| (3.3b) Once ERDM is obtained, other relevant quantities can also be derived. For example, the system partition function  \begin{eqnarray} Z_{\rm s} = {\rm tr} [\hat{\rho}^\beta] = \sum\limits_n \rho^\beta_{nn} = \frac{Z}{Z_{\rm B}} \end{eqnarray} (3.4) with the total partition function  \begin{eqnarray} Z = {\rm tr} [{\rm e}^{-\beta \hat{H}}] \end{eqnarray} (3.5) the system internal energy  \begin{eqnarray}\label{uint} U_{\rm s} = -\frac{\partial {\rm ln} Z_{\rm s}}{\partial\beta} = \langle H \rangle_\beta - \langle H_{\rm B} \rangle_\beta \end{eqnarray} (3.6) with the internal energies  \langle H \rangle_\beta = \frac{1}{Z} {\rm tr} [{\rm e}^{-\beta \hat{H}} \hat{H}] (3.7a)  \langle H_{\rm B} \rangle_\beta = \frac{1}{Z_{\rm B}} {\rm tr} [{\rm e}^{-\beta \hat{H}_B} \hat{H}_{\rm B}] (3.7b) It should be emphasized that as long as the interaction between the system and the bath does not vanish, it is not possible to separate the "system" unambiguously from the "bath" thermodynamically. The above definition is a subjective classification that is often used in qualitative analysis. Practical evaluation of ERDM and other properties employs ML-MCTDH for the imaginary time propagation and an importance sampling procedure for Eq.(3.2). To do so we first introduce a zeroth-order Hamiltonian H_0 that can be diagonalized either analytically or numerically  \begin{eqnarray} {\rm e}^{-\beta \hat{H}_0} = \sum\limits_\alpha {\rm e}^{-\beta E_\alpha} |\chi_\alpha\rangle \langle \chi_\alpha| \end{eqnarray} (3.8) and insert it into the ERDM expression Eq.(3.2)  \begin{eqnarray} \label{prac-erdm} \rho^\beta_{nm} \hspace{-0.15cm}& = &\hspace{-0.15cm} \frac{1}{Z_{\rm B}} {\rm tr} [{\rm e}^{-\beta \hat{H}} |n\rangle\langle m|] \nonumber \\ &=&\hspace{-0.15cm}\frac{1}{Z_{\rm B}} {\rm tr} [{\rm e}^{-\beta \hat{H}_0} {\rm e}^{\beta \hat{H}_0/2} {\rm e}^{-\beta \hat{H}/2} |n\rangle\langle m| {\rm e}^{-\beta \hat{H}/2} {\rm e}^{\beta \hat{H}_0/2}] \nonumber \\ &= &\hspace{-0.15cm} \frac{Z_0}{Z_{\rm B}} \frac{1}{Z_0} \sum\limits_\alpha \Big({\rm e}^{-\beta E_\alpha} \langle \chi_\alpha| {\rm e}^{-\beta (\hat{H}-E_\alpha)/2}\nonumber\\ &&|n\rangle\langle m| {\rm e}^{-\beta (\hat{H}-E_\alpha)/2} |\chi_\alpha\rangle\Big) \end{eqnarray} (3.9) where we have defined the partition function  \begin{eqnarray} Z_0 = {\rm tr} [{\rm e}^{-\beta \hat{H}_0}] \end{eqnarray} (3.10) The quantum mechanical trace in Eq.(3.2) is obtained from importance sampling according to the weighting function {\rm e}^{-\beta E_\alpha}/Z_0 which is defined by a zeroth-order Hamiltonian H_0. The choice of H_0 is flexible as long as it is not too far away from H and can be diagonalized. For the system-bath Hamiltonian, an often good choice is H_0=H_{\rm S}+H_{\rm B}. Sometimes it is also useful to define H_0 in a more localized region to enhance the importance sampling [52]. Once the initial state |\chi_\alpha\rangle is selected, it is propagated in the imaginary time via ML-MCTDH to give {\rm e}^{-\beta (\hat{H}-E_\alpha)/2} |n\rangle. The final ERDM is calculated according to Eq.(3.9). Ⅳ. MULTI-ELECTRONIC-STATE PATH INTEGRAL MOLECULAR DYNAMICS The recently developed multi-electronic-state path integral molecular dynamics (MES-PIMD) approach [20] in principle offers a practical tool in either of the diabatic or adiabatic representations for studying exact quantum statistics of general MES systems when the Born-Oppenheimer approximation, Condon approximation, and harmonic bath approximation break down. The MES-PIMD approach employs a unified efficient thermostat scheme (the "middle" scheme) for PIMD (or MD) which applies to either stochastic or deterministic thermostats [53-57]. A. Imaginary time path integral formulation for multi-electronic-state systems Consider a Hamiltonian with N electronic states in the diabatic representation \hat{\bf{H}}=\hat{\bf{T}}+\hat{\bf{V}}, where \hat{\bf{V}}=\bf{V}(\hat{{\bf{R}}}) is the N$$\times$$N symmetric potential energy matrix as a function of the nuclear coordinate {{\bf{R}}} and \hat{\bf{T}} the (diagonal) kinetic energy matrix. The canonical partition function is  \begin{eqnarray} \label{partition} Z = {\rm tr}_{ne} [{\rm e}^{-\beta \hat{\bf{H}}}] \end{eqnarray} (4.1) and a specific physical property of interest is  \begin{eqnarray} \label{property} \langle \hat{B} \rangle = \frac{1}{Z} { {\rm tr}}_{ne} [{\rm e}^{-\beta \hat{\bf{H}}} \hat{B}] \end{eqnarray} (4.2) In Eq.(4.1) and Eq.(4.2) the trace is over both the nuclear and electronic degrees of freedom. That is,  \begin{eqnarray} \label{Etrace} { {\rm tr}}_{e} [\hat{B}]= \sum\limits_{n=1}^N \langle n|\hat{B}|n \rangle \end{eqnarray} (4.3) and  \begin{eqnarray} \label{Ntrace} { {\rm tr}}_{n} [\hat{B}]= \int {\rm d}{{\bf{R}}} \langle {{\bf{R}}}|\hat{B}|{{\bf{R}}} \rangle \end{eqnarray} (4.4) Inserting the resolution of the identity into Eq.(4.1) yields  \begin{eqnarray} \label{Pbeads} Z\hspace{-0.15cm}&=&\hspace{-0.15cm}\lim\limits_{P \to \infty} \int {{\rm d}{{\bf{R}}_1}...{\rm d}{{\bf{R}}_P}} \nonumber\\ &&\sum\limits_{n_1, \cdots, n_P=1}^N \prod\limits_{i=1}^P \langle {{\bf{R}}}_i, n_i|{\rm e}^{-\beta \hat{\bf{H}}/P}|{{\bf{R}}}_{i+1}, n_{i+1} \rangle \hspace{1cm} \end{eqnarray} (4.5) Eq.(4.5) then becomes  \begin{eqnarray} \label{Zexpress} Z\hspace{-0.15cm}&=&\hspace{-0.15cm}\lim\limits_{P \to \infty} \left | \frac{P \bf M}{2 \pi \beta \hbar^2} \right | ^{P/2} \int {{\rm d}{{\bf{R}}_1}...{\rm d}{{\bf{R}}_P}} \nonumber\\ &&{\rm exp} \left [-\frac{\beta}{2}\omega_P^2 \sum\limits_{i=1}^P ({\bf R}_i-{\bf R}_{i+1})^T {\bf M} ({\bf R}_i-{\bf R}_{i+1}) \right ] \nonumber\\ &&\times { {\rm Tr}}_e \left [ \prod\limits_{i=1}^P {\bf O}({\bf R}_i)^T {\bf O}({\bf R}_i) \right ] \hspace{0.5cm} \end{eqnarray} (4.6) where {\bf O}({\bf R}_i) depends on the splitting scheme used for expressing the element \langle {{\bf{R}}}_i, n_i|{\rm e}^{-\beta \hat{\bf{H}}/P}|{{\bf{R}}}_{i+1}, n_{i+1} \rangle in Eq.(4.5). We have proposed three splitting schemes, namely the "diagonalization", "first-order expansion", and "hyperbolic function" methods, the details of which are described in Ref.[20]. Because { {\rm Tr}}_e \left [ \prod\limits_{i=1}^P {\bf O}({\bf R}_i)^T {\bf O}({\bf R}_i) \right ] is often not positive-definite for general MES systems, regardless of which splitting scheme is employed. It is then not numerically favorable to use either { {\rm Tr}}_e \left [ \prod\limits_{i=1}^P {\bf O}({\bf R}_i)^T {\bf O}({\bf R}_i) \right ] or its absolute value to define an effective (real-valued) potential function \phi ({\bf R}_1, \cdots, {\bf R}_P) for performing PIMD. A reasonable way to define the effective potential [20] is  \begin{eqnarray} \label{EffPot} {\rm e}^{-\beta \phi^{{\rm (dia)}} ({\bf R}_1, \cdots, {\bf R}_P)} = {{\rm Tr}}_e \left [ {\prod\limits_{i=1}^P} {\rm e}^{-\beta {\bf V}_{{\rm diag}} ({\bf R}_i)/P} \right ] \end{eqnarray} (4.7) because the right-hand side (RHS) of Eq.(4.7) is always positive-definite. Here {\bf V}_{{\rm diag}} ({\bf R}) is a diagonal matrix, the diagonal elements of which are the same as those of the potential energy matrix {\bf V} ({\bf R}). The partition function (Eq.(4.6)) may then be evaluated by  \begin{eqnarray} \label{Zexp2} Z\hspace{-0.15cm}&=&\hspace{-0.15cm}\lim\limits_{P \to \infty} \left | \frac{P \bf M}{2 \pi \beta \hbar^2} \right | ^{P/2} \int {{\rm d}{{\bf{R}}_1}...{\rm d}{{\bf{R}}_P}} \nonumber\\ && {\rm exp} \left [-\beta U_{{\rm eff}}^{{\rm (dia)}} ({\bf R}_1, \cdots, {\bf R}_P) \right ] {\tilde Z}^{{\rm (dia)}} ({\bf R}_1, \cdots, {\bf R}_P)\nonumber\\ \end{eqnarray} (4.8) with the estimator for the partition function  \begin{eqnarray} \label{Zest} {\tilde Z}^{{\rm {\rm (dia)}}} ({\bf R}_1, \cdots, {\bf R}_P)=\frac{{{\rm Tr}}_e \left [ \displaystyle{{\prod\limits_{i=1}^P}} {\bf O}({\bf R}_i)^T {\bf O}({\bf R}_i) \right ]}{{{\rm Tr}}_e \left [ \displaystyle{{\prod\limits_{i=1}^P}} {\rm e}^{-\beta {\bf V}_{{\rm diag}} ({\bf R}_i)/P} \right ]} \end{eqnarray} (4.9) and  \begin{eqnarray} \label{Ueff} U_{{\rm eff}}^{{\rm (dia)}} ({\bf R}_1, \cdots, {\bf R}_P)=\frac{1}{2}\omega_P^2 \sum\limits_{i=1}^P ({\bf R}_i-{\bf R}_{i+1})^T\cdot \nonumber\\ {\bf M} ({\bf R}_i-{\bf R}_{i+1})+\phi^{{\rm (dia)}} ({\bf R}_1, \cdots, {\bf R}_P) \end{eqnarray} (4.10) Any physical property in Eq.(4.2) can be expressed as  \begin{eqnarray} \label{prop2} \langle \hat{B} \rangle \hspace{-0.15cm}&=&\hspace{-0.15cm} \lim\limits_{P \to \infty} \frac{ {\cal A }} {{\cal B}}\\ {\cal A}\hspace{-0.15cm}&=&\hspace{-0.15cm} \int {{\rm d}{{\bf{R}}_1}...{\rm d}{{\bf{R}}_P}} {\rm exp} \left [-\beta U_{{\rm eff}}^{{\rm (dia)}} ({\bf R}_1, \cdots, {\bf R}_P) \right ]\nonumber\\ && {\tilde B}^{{\rm (dia)}} ({\bf R}_1, \cdots, {\bf R}_P) \nonumber\\ {\cal B}\hspace{-0.15cm}&=&\hspace{-0.15cm}\int {{\rm d}{{\bf{R}}_1}...{\rm d}{{\bf{R}}_P}} {\rm exp} \left [-\beta U_{{\rm eff}}^{({\rm dia})} ({\bf R}_1, \cdots, {\bf R}_P) \right ]\nonumber\\ && {\tilde Z}^{{\rm (dia)}} ({\bf R}_1, \cdots, {\bf R}_P)\nonumber \end{eqnarray} (4.11) with {\tilde B}^{{\rm (dia)}} ({\bf R}_1, \cdots, {\bf R}_P) as the estimator for operator \hat B in the diabatic representation  \begin{array}{l} {{\tilde B}^{({\rm{dia}})}}({{{\bf{R}}}_1}, \cdots ,{{{\bf{R}}}_P}) = \frac{{\cal C}}{{\rm{T}}{{\rm{r}}_e}{\left[ {\mathop \prod \limits_{i = 1}^P {{\rm{e}}^{ - \beta {{\bf{V}}_{{\rm{diag}}}}({{{\bf{R}}}_i})/P}}} \right]}}\\ \;\;\;\;{\cal C} = \frac{1}{P}{\rm{T}}{{\rm{r}}_e}\left\{ {\sum\limits_{k = 1}^P {} } \right.\left[ {\left( {\mathop \prod \limits_{i = 1}^{k - 1} {\bf{O}}{{({{{\bf{R}}}_i})}^T}{\bf{O}}({{{\bf{R}}}_i})} \right)} \right.{\bf{O}}{({{{\bf{R}}}_k})^T}{\bf{B}}({{{\bf{R}}}_k})\\ \;\;\;\;\;\;\;\;{\bf{O}}({{{\bf{R}}}_k})\left. {\left. {\left( {\mathop \prod \limits_{i = k + 1}^P {\bf{O}}{{({{{\bf{R}}}_i})}^T}{\bf{O}}({{{\bf{R}}}_i})} \right)} \right]} \right\} \end{array} (4.12) For instance, the estimator for the potential energy is  \begin{array}{l} {{\tilde B}^{({\rm{dia}})}}({{{\bf{R}}}_1}, \cdots ,{{{\bf{R}}}_P}) = \frac{{\cal D}}{{\rm{T}}{{\rm{r}}_e}{\left[ {\mathop \prod \limits_{i = 1}^P {{\rm{e}}^{ - \beta {{\bf{V}}_{{\rm{diag}}}}({{{\bf{R}}}_i})/P}}} \right]}}\\ \;\;\;\;{\cal D} = \frac{1}{P}{\rm{T}}{{\rm{r}}_e}\left\{ {\sum\limits_{k = 1}^P {} } \right.\left[ {\left( {\mathop \prod \limits_{i = 1}^{k - 1} {\bf{O}}{{({{{\bf{R}}}_i})}^T}{\bf{O}}({{{\bf{R}}}_i})} \right)} \right.{\bf{O}}{({{{\bf{R}}}_k})^T}{\bf{V}}({{{\bf{R}}}_k})\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\bf{O}}({{{\bf{R}}}_k})\left. {\left. {\left( {\mathop \prod \limits_{i = k + 1}^P {\bf{O}}{{({{{\bf{R}}}_i})}^T}{\bf{O}}({{{\bf{R}}}_i})} \right)} \right]} \right\} \end{array} (4.13) and that for the electronic state density matrix element in the diabatic representation \rho_{ij}^{{\rm ds}} is  \begin{array}{l} \tilde \rho _{ij}^{({\rm{dia - ds}})}({{{\bf{R}}}_1}, \cdots ,{{{\bf{R}}}_P}) = \frac{{\cal E}}{{{\rm{T}}{{\rm{r}}_e}\left[ {\mathop \prod \limits_{i = 1}^P {{\rm{e}}^{ - \beta {{\bf{V}}_{{\rm{diag}}}}({{{\bf{R}}}_i})/P}}} \right]}}\\ \;\;{\cal E} = \frac{1}{P}{\rm{T}}{{\rm{r}}_e}\left\{ {\sum\limits_{k = 1}^P {} } \right.\left[ {\left( {\mathop \prod \limits_{i = 1}^{k - 1} {\bf{O}}{{({{{\bf{R}}}_i})}^T}{\bf{O}}({{{\bf{R}}}_i})} \right)} \right.{\bf{O}}{({{{\bf{R}}}_k})^T}\mathit{\boldsymbol{\rho }}_{ij}^{({\rm{ds}})}\\ \;\;\;\;\;\;\;\;\;{\bf{O}}({{{\bf{R}}}_k})\left. {\left. {\left( {\mathop \prod \limits_{i = k + 1}^P {\bf{O}}{{({{{\bf{R}}}_i})}^T}{\bf{O}}({{{\bf{R}}}_i})} \right)} \right]} \right\} \end{array} (4.14) where  \begin{eqnarray} \label{rhoij} {\boldsymbol{\rho}}_{ij}^{{\rm (ds)}}=\frac{1}{2}(|i \rangle \langle j|+|j \rangle \langle i|) \end{eqnarray} (4.15) with |i \rangle and |j \rangle as the i-th and j-th electronically diabatic states in the representation for the Hamiltonian \hat{\bf{H}}. B. Multi-electronic-state path integral molecular dynamics One may employ the MD or Monte Carlo scheme to perform the integrals in Eq.(4.11). Inserting fictitious momenta ({\bf p}_1, \cdots, {\bf p}_P) into each integral of Eq.(4.11) produces  \begin{array}{l} \langle \hat B\rangle = \mathop {\lim }\limits_{P \to \infty } \frac{{\cal F}}{{\cal G}}\\ {\cal F} = \int {{\rm{d}}{{{\bf{R}}}_1}...{\rm{d}}{{{\bf{R}}}_P}} \int {{\rm{d}}{{{\bf{p}}}_1}...{\rm{d}}{{{\bf{p}}}_P}} {\rm{exp}}\left[ - \right.\beta H_{{\rm{eff}}}^{({\rm{dia}})}\\ \;\;\;\;\;\;\;\;\;\;({{{\bf{R}}}_1}, \cdots ,{{{\bf{R}}}_P};{{{\bf{p}}}_1}, \cdots ,{{{\bf{p}}}_P})]{{\tilde B}^{({\rm{dia}})}}({{{\bf{R}}}_1}, \cdots ,{{{\bf{R}}}_P})\\ {\cal G} = \int {{\rm{d}}{{{\bf{R}}}_1}...{\rm{d}}{{{\bf{R}}}_P}} \int {{\rm{d}}{{{\bf{p}}}_1}...{\rm{d}}{{{\bf{p}}}_P}} {\rm{exp[}} - \beta H_{{\rm{eff}}}^{({\rm{dia}})}\\ \;\;\;\;\;\;\;\;\;\;({{{\bf{R}}}_1}, \cdots ,{{{\bf{R}}}_P};{{{\bf{p}}}_1}, \cdots ,{{{\bf{p}}}_P})]{{\tilde Z}^{({\rm{dia}})}}({{{\bf{R}}}_1}, \cdots ,{{{\bf{R}}}_P}) \end{array} (4.16) where the effective Hamiltonian for Eq.(4.16) is  \begin{eqnarray} \label{Heff} H_{{\rm eff}}^{{\rm (dia)}}&&\hspace{-0.3cm} ({\bf R}_1, \cdots, {\bf R}_P;{\bf p}_1, \cdots, {\bf p}_P)=\nonumber\\ && \frac{1}{2} \sum\limits_{j=1}^P {\bf p}_j^T {\tilde {\bf M}}_j^{-1} {\bf p}_j + U_{{\rm eff}}^{({\rm dia})} ({\bf R}_1, \cdots, {\bf R}_P) \end{eqnarray} (4.17) with the fictitious masses {\tilde {\bf M}}_j. It is trivial to verify that Eq.(4.16) is identical to Eq.(4.11) by performing the Gaussian integral over the fictitious momenta ({\bf p}_1, \cdots, {\bf p}_P) in the former. Thus the MES-PIMD equations of motion are given by  {\dot {\bf R}}_j = {\tilde {\bf M}}_j^{-1} {\bf p}_j \\ {\dot {\bf p}}_j = -\frac{\partial U_{{\rm eff}}^{{\rm (dia)}} ({\bf R}_1, \cdots, {\bf R}_P)}{\partial {\bf R}_j} \quad (j=1, \cdots, P) (4.18) Note that Eq.(4.18) must be coupled to a thermostating method to ensure a proper canonical distribution. For many thermostats, the integration in one time step \Delta t can be split into three parts [20, 53-57], the steps for updating coordinates, momenta, and thermostat, denoted as "R", "p", and "T", respectively. In this case the "equations of motion" may be expressed as  \begin{eqnarray} \label{EOMdt} \begin{bmatrix} {\rm d}{\bf{R}}\\ {\rm d}{\bf{p}} \end{bmatrix}\hspace{-0.1cm} &=&\hspace{-0.1cm} \underbrace{ \begin{bmatrix} \tilde{\bf{M}}^{-1}{\bf{p}}\\ 0 \end{bmatrix}{\rm d}t }_{\text{R}} + \underbrace{ \begin{bmatrix} 0 \\ -\nabla_{{\bf{R}}}U_{{\rm eff}}^{({\rm dia})}({\bf{R}}) \end{bmatrix}{\rm d}t }_{\text{p}} \nonumber\\ &&+\underbrace{ \begin{bmatrix} \text{thermostat} \\ \end{bmatrix} }_{\text{T}} \end{eqnarray} (4.19) For convenience we denote {\bf{R}}$$ \equiv$$({\bf R}_1, \cdots, {\bf R}_P) and {\bf{p}}$$\equiv({\bf p}_1, \cdots, {\bf p}_P) in Eq.(4.19) and thereafter in this section. A useful approach is to employ the forward Kolmogorov equation to express the evolution of the density distribution in the phase space [20, 53-57] \mathcal{P} ({\bf{R}}, {\bf{p}}).  \begin{eqnarray} \begin{aligned} \frac{\partial}{\partial t}\mathcal{P} =& \mathcal{L}\mathcal{P}\\ =& (\mathcal{L}_{\text{R}} + \mathcal{L}_{\text{p}} + \mathcal{L}_{\text{T}})\mathcal{P} \end{aligned} \end{eqnarray} (4.20) The relevant Kolmogorov operators for the 1st and 2nd terms of the RHS are  \mathcal{L}_{\text{R}}\mathcal{P} \hspace{-0.15cm} = \hspace{-0.15cm} -{\bf{p}}^T \tilde{\bf{M}}^{-1}\nabla_{{\bf{R}}}\mathcal{P} (4.21)  \mathcal{L}_{\text{p}}\mathcal{P} \hspace{-0.15cm} = \hspace{-0.15cm} \nabla{}_{{\bf{R}}} U_{{\rm eff}}^{{\rm (dia)}}({\bf{R}})\cdot\nabla{}_{\bf{p}}\mathcal{P} (4.22) respectively. The definition of \mathcal{L}_{\text{T}} depends on the specific thermostat. The phase space propagators for a time interval \Delta t for the three parts are {\rm e}^{\mathcal{L}_{\text{R}}\Delta t}, {\rm e}^{\mathcal{L}_{\text{p}}\Delta t}, and {\rm e}^{\mathcal{L}_{\text{T}}\Delta t}, respectively. The propagation in each time step with the velocity Verlet (VV) algorithm in the "middle" scheme [20, 53-57] reads  \begin{eqnarray} {\rm e}^{\mathcal{L}\Delta t} \approx {\rm e}^{\mathcal{L}_{\text{p}}\Delta t/2} {\rm e}^{\mathcal{L}_{\text{R}}\Delta t/2} {\rm e}^{\mathcal{L}_{\text{T}}\Delta t} {\rm e}^{\mathcal{L}_{\text{R}}\Delta t/2} {\rm e}^{\mathcal{L}_{\text{p}}\Delta t/2} \label{middle:eq2} \end{eqnarray} (4.23) The phase space propagator for the thermostat part {\rm e}^{\mathcal{L}_{\text{T}}\Delta t} is designed in the middle of the propagation. The numerical algorithm for Eq.(4.23) is  \begin{eqnarray} && \text{Update momenta for half a step:}\nonumber\\ &&\quad {\bf{p}} \leftarrow {\bf{p}} - \frac{\partial U_{{\rm eff}}^{({\rm dia})}}{\partial {\bf{R}}} \frac{\Delta t}{2} \nonumber\\ && \text{Update coordinates for half a step:} \nonumber\\ &&\quad {\bf{R}} \leftarrow {\bf{R}} + \tilde{\bf{M}}^{-1}{\bf{p}} \frac{\Delta t}{2}\nonumber \end{eqnarray} Thermostat for a full time step: thermostat_step  \begin{eqnarray}&& \text{Update coordinates for another half step:} \nonumber\\ && \quad {\bf{R}} \leftarrow {\bf{R}} + \tilde{\bf{M}}^{-1}{\bf{p}} \frac{\Delta t}{2}\nonumber \\ && \text{Update momenta for another half step:} \nonumber\\ &&\quad {\bf{p}} \leftarrow {\bf{p}} - \frac{\partial U_{{\rm eff}}^{{\rm (dia)}}}{\partial {\bf{R}}} \frac{\Delta t}{2}\nonumber \end{eqnarray} where thermostat_step represents the thermostat process for a time interval \Delta t, which is determined according to the stochastic or deterministic thermostat method of choice [20, 53-57]. Note that PIMD with the "middle" thermostat scheme is often more efficiently performed by employing the staging transformation of path integral beads [20, 53, 55] or the normal-mode transformation [59]. In the present paper we use the MES-PIMD approach [20] in the way that takes no advantage of the specific form of the system-bath model. It should be stressed that one could design a much more efficient method for performing MES-PIMD when the Condon approximation and harmonic-bath approximation are valid. Ⅴ. NUMERICAL EXAMPLE: THE SPIN-BOSON MODEL A. Model We consider practical application of our method to the spin-boson model that has been widely used in the context of electron transfer theory [60, 61]. This is a relatively simple system-bath type Hamiltonian where the system comprises only two electronic states (|\phi_1\rangle and |\phi_2\rangle) that are linearly coupled to a bath of harmonic oscillators. Using mass-weighted coordinates the Hamiltonian reads  \begin{eqnarray} \label{sbh} H =- \epsilon \sigma_z + \Delta\sigma_x + \frac{1}{2}\sum\limits_{j=1}^N (p^2_j+\omega_j^2 q^2_j) + \sigma_z \sum\limits_{j=1}^N c_j q_j \end{eqnarray} (5.1) where \sigma_x and \sigma_z are Pauli matrices  \sigma_x \hspace{-0.15cm} = \hspace{-0.15cm} |\phi_1\rangle \langle \phi_2 | + |\phi_2\rangle \langle \phi_1 | (5.2a)  \sigma_z \hspace{-0.15cm} = \hspace{-0.15cm} |\phi_1\rangle \langle \phi_1 | - |\phi_2\rangle \langle \phi_2 | (5.2b) The properties of the bath that influence the dynamics of the two-state subsystem are specified by the spectral density function [60, 61]  \begin{eqnarray} \label{Jw} J(\omega) = \frac{\pi}{2} \sum\limits_j \frac{c_j^2}{\omega_j}\delta(\omega - \omega_j) \end{eqnarray} (5.3) Different forms of J(\omega) provide different models of the phonon bath. In the examples below we use two forms: an Ohmic (linear) spectral density with an exponential cutoff  J(\omega) = \eta \omega \; {\rm e}^{-\omega/ \omega_c} (5.4a) and an Ohmic spectral density with a Lorentzian cutoff  \begin{eqnarray} \label{Jdebye} J(\omega) = \eta \omega \frac{1}{1+\left(\displaystyle{\frac{\omega}{\omega_{\rm c}}}\right)^2} \end{eqnarray} (5.4b) where \eta is the dimensionless system-bath coupling strength and \omega_{\rm c} is the cutoff frequency of the bath. The continuous bath spectral density of Eq.(5.4) can be discretized to the form of Eq.(5.3) via the relation [15, 62]  \begin{eqnarray} \label{coupl} c_j^2 = \frac{2}{\pi} \omega_j \frac{J(\omega_j)} {\rho(\omega_j)} \end{eqnarray} (5.5) in which the density of frequencies \rho(\omega) is defined from the integral relation  \begin{eqnarray}\label{den1} \int_0^{\omega_j} {\rm d}\omega \; \rho(\omega) = j, \hspace{.5in} j = 1, ..., N \end{eqnarray} (5.6a) In this work, \rho(\omega) is chosen as  \begin{eqnarray} \rho(\omega) = s \frac{J(\omega)}{\omega} \end{eqnarray} (5.6b) with Eq.(5.6a) enforced to obtain the scaling factor s. To model the condensed phase environment in the ERDM simulation, the number of bath modes N in the discretization is a convergence parameter and needs to be sufficiently large to represent the continuum for the property of interest. In this work N has been chosen in the range of a few tens to a few hundred. This number is somewhat smaller than a normal quantum dynamics simulation on a similar spin-boson model, suggesting ERDM is easier to converge than the real time quantum dynamics. B. Results We first consider a model for an Ohmic spectral density with an exponential cutoff, Eq.(5.4a). The parameters are given as: \epsilon =\Delta=60 cm^{-1} in Hamiltonian (Eq.(5.1)), \omega_c=100 cm^{-1} and \eta=0.3 in spectral density (Eq.(5.4a)). FIG. 1 shows all the elements of ERDM over the temperature range 20-300 K. Despite using completely different numerical strategies, ML-MCTDH and MES-PIMD produce results in agreement over the entire temperature range, for both diagonal and off-diagonal elements of ERDM. This is encouraging because both methods are seen to be capable of capturing the high temperature (classical) and low temperature (quantum) limits. It also provides benchmark results to test other methodologies. As mentioned in the introduction, a method that is accurate in principle may not be so in practice. The comparison here offers a reliable and unbiased test for both methods.  FIG. 1 Equilibrium reduced density matrix for a spin-boson model. An Ohmic spectral density with an exponential cutoff is employed. The parameters are: \epsilon=\Delta=60 cm^{-1}, \omega_c=100 cm^{-1}, and \eta=0.3. The lines are intended as guide to the eye. The example above is in the relatively weak coupling regime where both methods are easy to converge. For ML-MCTDH simulations only two layers are needed, with a few thousand configurations per layer to obtain accurate results. For MES-PIMD calculations the number of beads ranges from 32 at 300 K to 256 at 20 K, which also has moderate computational cost. The error bar for each method is smaller than the size of the symbols in FIG. 1. Next we consider a model with the spectral density function (Eq.(5.4b)). A symmetric two-level system is studied here, \epsilon=0. The parameters and physical variables are scaled with respect to the tunneling splitting of the bare two-level system 2\Delta, i.e. 2\Delta\rightarrow $$2\Delta/(2\Delta) =1, \omega_c$$\rightarrow $$\omega_c / (2\Delta), \omega$$\rightarrow $$\omega / (2\Delta), \beta$$\rightarrow$(2\Delta) \beta$. The parameters are given as: $\Delta$=0.5 in Hamiltonian (Eq.(5.1)), $\omega_c$=0.5, and $\eta$=1 in spectral density (Eq.(5.4b)). We note that the Hamiltonian

 $\begin{eqnarray} H = - \epsilon \sigma _x + \Delta \sigma _z + \frac{1}{2}\mathop \sum \limits_{j = 1}^N (p_j^2 + \omega _j^2 q_j^2 ) + \sigma _x \mathop \sum \limits_{j = 1}^N c_j q_j \\ \end{eqnarray}$ (5.7)

and that in Eq.(5.1) share the same thermodynamic properties of the total system such as the internal energy.

The results of ERDM agree between ML-MCTDH and MES-PIMD simulations. In addition, we consider a different property, the internal energy of the system as defined in Eq.(3.6). As shown in FIG. 2, the results from the two methods again agree well with each other. Since this example has a stronger coupling strength as well as a quantity (internal energy) more difficult to converge, both methods require more computational effort. ML-MCTDH requires 2 - 3 layers, with ~10$^5$ configurations for the top layer; whereas MES-PIMD requires 128 beads ( $\beta$=0.5) to 2048 beads ($\beta$=5). We note that both methods employ different sets of convergence parameters. The estimate of errors is about a few percent within each method's self check. This confirms again the reliability of both methods.

 FIG. 2 System internal energy, Eq.(3.6), for a spin-boson model. An Ohmic spectral density with a Lorentzian cutoff is employed. The parameters are: $\epsilon$=0, $\Delta$ =0.5, $\omega_c$=0.5, and $\eta$=1. The lines are intended as guide to the eye.
Ⅵ. CONCLUDING REMARKS

The ML-MCTDH and the MES-PIMD methods were both derived from exact principles. Both have been applied in other contexts and verified many times against available benchmark results. In this paper, we extended the methods and presented our implementations for computing reduced density matrices of system-bath models at equilibrium. We tested our approaches on the spin-boson model with two different spectral densities. Excellent agreement is achieved between the two methods. We have also found agreement between our methods for other parameter regimes. All this demonstrates the reliability of the two methods.

One advantage of both methods is that they are not restricted to treating a harmonic bath. It would be interesting to extend our implementation to treat anharmonic potential functions. This will provide a useful way to go beyond linear response model and investigate new physics and chemistry therein.

Ⅶ. ACKNOWLEDGMENTS

The work on the ML-MCTDH approach in the paper has been supported by the U.S. National Science Foundation CHE-1500285, and used resources from the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No.DE-AC02-05CH11231.

The work on the MES-PIMD approach in the paper has been supported by the Ministry of Science and Technology of China (No.2017YFA0204901 and No.2016YFC0202803, the National Natural Science Foundation of China (No.21373018 and No.21573007), the Recruitment Program of Global Experts, and Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase) under grant No.U1501501. We acknowledge the Beijing and Tianjin supercomputer centers and the High-performance Computing Platform of Peking University for providing computational resources.

Reference
 [1] R. Egger, and U. Weiss, Z. Phys. B 89 , 97 (1992). [2] R. Egger, and C. H. Mak, Phys. Rev. B 50 , 15210 (1994). [3] N. Makri, and D. E. Makarov, J. Chem. Phys. 102 , 4600 (1995). [4] E. Sim, and N. Makri, J. Phys. Chem. B 101 , 5446 (1997). [5] L. Mühlbacher, and R. Egger, J. Chem. Phys. 118 , 179 (2003). [6] R. Bulla, N. G. Tong, and M. Vojta, Phys. Rev. Lett. 91 , 170601 (2003). [7] F. B. Anders, and A. Schiller, Phys. Rev. B 74 , 245113 (2006). DOI:10.1103/PhysRevB.74.245113 [8] F. B. Anders, R. Bulla, and M. Vojta, Phys. Rev. Lett. 98 , 210402 (2007). DOI:10.1103/PhysRevLett.98.210402 [9] Y. Tanimura, and R. Kubo, J. Phys. Soc. Jpn. 58 , 101 (1989). DOI:10.1143/JPSJ.58.101 [10] Y. Yan, F. Yang, Y. Liu, and J. Shao, Chem. Phys. Lett. 395 , 216 (2004). [11] R. Xu, P. Cui, X. Li, Y. Mo, and Y. Yan, J. Chem. Phys. 122 , 041103 (2005). DOI:10.1063/1.1850899 [12] X. Zheng, J. Jin, S. Welack, M. Luo, and Y. Yan, J. Chem. Phys. 130 , 164708 (2009). [13] Y. Tanimura, J. Chem. Phys. 141 , 044114 (2014). [14] Y. Tanimura, J. Chem. Phys. 142 , 144110 (2015). DOI:10.1063/1.4916647 [15] H. Wang, and M. Thoss, J. Chem. Phys. 119 , 1289 (2003). DOI:10.1063/1.1580111 [16] H. Wang, and M. Thoss, J. Chem. Phys. 131 , 024114 (2009). DOI:10.1063/1.3173823 [17] H. Wang, J. Phys. Chem. A 119 , 7951 (2015). [18] A. Suarez, and R. Silbey, J. Chem. Phys. 95 , 9115 (1991). [19] N. Makri, J. Phys. Chem. B 103 , 2823 (1999). DOI:10.1021/jp9847540 [20] X. Liu, and J. Liu, J. Chem. Phys. 148 , 102319 (2018). DOI:10.1063/1.5005059 [21] U. Manthe, J. Chem. Phys. 130 , 054109 (2009). DOI:10.1063/1.3069655 [22] T. Hammer, and U. Manthe, J. Chem. Phys. 134 , 224305 (2011). DOI:10.1063/1.3598110 [23] T. Hammer, and U. Manthe, J. Chem. Phys. 136 , 054105 (2012). DOI:10.1063/1.3681166 [24] R. Wodraszka, J. Palma, and U. Manthe, J. Phys. Chem. A 116 , 11249 (2012). DOI:10.1021/jp3052642 [25] H. Wang, J. Phys. Chem. A 117 , 9253 (2014). [26] J. Frenkel, Wave Mechanics. Oxford: Clarendon Press (1934). [27] H. D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165 , 73 (1990). DOI:10.1016/0009-2614(90)87014-I [28] U. Manthe, H. D. Meyer, and L. S. Cederbaum, J. Chem. Phys. 97 , 3199 (1992). [29] M. H. Beck, A. Jäckle, G. A. Worth, and H. D. Meyer, Phys. Rep. 324 , 1 (2000). DOI:10.1016/S0370-1573(99)00047-2 [30] H. D. Meyer, and G. A. Worth, Theor. Chem. Acc. 109 , 251 (2003). DOI:10.1007/s00214-003-0439-1 [31] L. R. Tucker, Psychometrika 31 , 279 (1966). DOI:10.1007/BF02289464 [32] T. G. Kolda, and B. W. Bader, SIAM Rev. 51 , 455 (2008). [33] M. Thoss, I. Kondov, and H. Wang, Chem. Phys. 304 , 169 (2004). [34] M. Thoss, and H. Wang, Chem. Phys. 322 , 210 (2006). DOI:10.1016/j.chemphys.2005.07.011 [35] H. Wang, and M. Thoss, J. Chem. Phys. 124 , 034114 (2006). [36] I. Kondov, H. Wang, and M. Thoss, J. Phys. Chem. A 110 , 1364 (2006). [37] H. Wang, and M. Thoss, J. Phys. Chem. A 111 , 10369 (2007). DOI:10.1021/jp072367x [38] I. Kondov, M. Cizek, C. Benesch, H. Wang, and M. Thoss, J. Phys. Chem. C 111 , 11970 (2007). [39] M. Thoss, I. Kondov, and H. Wang, Phys. Rev. B 76 , 153313 (2007). [40] I. R. Craig, M. Thoss, and H. Wang, J. Chem. Phys. 127 , 144503 (2007). DOI:10.1063/1.2772265 [41] H. Wang, and M. Thoss, Chem. Phys. 347 , 139 (2008). DOI:10.1016/j.chemphys.2007.12.004 [42] H. Wang, and M. Thoss, New J. Phys. 10 , 115005 (2008). [43] D. Egorova, M. F. Gelin, M. Thoss, H. Wang, and W. Domcke, J. Chem. Phys. 129 , 214303 (2008). [44] K. A. Velizhanin, H. Wang, and M. Thoss, Chem. Phys. Lett. 460 , 325 (2008). DOI:10.1016/j.cplett.2008.05.065 [45] K. A. Velizhanin, and H. Wang, J. Chem. Phys. 131 , 094109 (2009). [46] H. Wang, and M. Thoss, Chem. Phys. 370 , 78 (2010). DOI:10.1016/j.chemphys.2010.02.027 [47] O. Vendrell, and H. D. Meyer, J. Chem. Phys. 134 , 044135 (2011). DOI:10.1063/1.3535541 [48] Y. Zhou, J. Shao, and H. Wang, Mol. Phys. 110 , 581 (2012). [49] H. Wang, and S. Shao, J. Chem. Phys. 137 , 22A504 (2012). [50] U. Manthe, J. Chem. Phys. 128 , 164116 (2008). DOI:10.1063/1.2902982 [51] L. Grasedyck, SIAM. J. Matrix Anal. & Appl. 31 , 2029 (2010). [52] H. Wang, D. E. Skinner, and M. Thoss, J. Chem. Phys. 125 , 174502 (2006). DOI:10.1063/1.2363195 [53] I. R. Craig, M. Thoss, and H. Wang, J. Chem. Phys. 135 , 064504 (2011). DOI:10.1063/1.3624342 [54] J. Liu, D. Li, and X. Liu, J. Chem. Phys. 145 , 024103 (2016). DOI:10.1063/1.4954990 [55] D. Li, X. Han, Y. Chai, C. Wang, Z. Zhang, Z. Chen, J. Liu, and J. Shao, J. Chem. Phys. 147 , 184104 (2017). DOI:10.1063/1.4996204 [56] Z. Zhang, X. Liu, Z. Chen, H. Zheng, K. Yan, and J. Liu, J. Chem. Phys. 147 , 034109 (2017). [57] D. Li, Z. Chen, Z. Zhang, and J. Liu, Chin. J. Chem. Phys. 30 , 735 (2017). DOI:10.1063/1674-0068/30/cjcp1711223 [58] J. Liu, Z. Zhang, K. Yan, Z. Chen, and X. Liu, Amber 2018 Reference Manual , 329–336 (2018). [59] J. Liu, D. Li, and X. Liu, J. Chem. Phys. 145 , 024103 (2016). [60] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59 , 1 (1987). DOI:10.1103/RevModPhys.59.1 [61] U. Weiss, Quantum Dissipative Systems, 2nd Edn. Singapore: World Scientific, (1999). [62] H. Wang, M. Thoss, and W. H. Miller, J. Chem. Phys. 115 , 2979 (2001). DOI:10.1063/1.1385561

a. 美国科罗拉多大学丹佛分校化学系, 科罗拉多州 80217-3364;
b. 北京大学化学与分子工程学院, 理论与计算化学研究所, 北京分子科学国家研究中心, 北京 100871