1 Introduction↩︎

In 1928, van Vleck derived the semi-classical real-time propagator (the van Vleck propagator) for a quantum system in the limit as \(\hbar\rightarrow0\) [1]. The phase of the van Vleck propagator, with respect to a given initial and final time and position, is determined by the action of the corresponding classical trajectory. The van Vleck propagator reveals the intrinsic connection between quantum and classical mechanics, having had a significant impact on quantum physics [2], [3].

However, deriving the van Vleck propagator requires solving the classical Hamiltonian equations with fixed initial and final positions. This is not convenient for numerical calculations of realistic systems due to the root-finding problem. Consequently, the van Vleck propagator is not well-suited for numerical studies of atom or molecule systems. To overcome this problem, many authors [4]–[10] have attempted to develop initial-value representations (IVRs) for semi-classical real-time propagators. In the IVRs, these propagators are expressed as functionals of classical trajectories, determined by the given initial position and momentum, which can be easily derived numerically via standard algorithms, such as Runge-Kutta. A highly influential IVR for the semi-classical real-time propagator was derived by Herman and Kluk (HK) in 1984 [10], based on coherent states. Over the past forty years, the HK representation has been widely used in studies of chemical dynamics, and has proven to be a valuable tool for exploring the quantum effects of nuclear motion [11], [12]. Moreover, in 2006, Kay provided a rigorous derivation of the HK representation from the Schrödinger equation, demonstrating that this representation is the leading term in an asymptotic expansion of the quantum propagator in powers of \(\hbar\) [13].

In addition to the real-time propagator, the imaginary-time propagator, which is the matrix element of the Boltzmann operator \(e^{-\hat{H}/(k_{B}T)}\), also play a crucial role in quantum physics and chemistry. Here, \(k_B\) and \(T\) are the Boltzmann constant and temperature, respectively. In 1971, Miller derived the semi-classical imaginary-time propagator [14], through a direct \(t \rightarrow -i\tau\) transformation of the van Vleck propagator, where \(t\) is the real time and \(\tau = \hbar/(k_B T)\). Similar to the van Vleck propagator, the result obtained by Miller in Ref. [14] is determined by the classical trajectory with an inverted potential, with respect to fixed initial and final positions.

An semi-classical IVR for the imaginary-time propagator, which is similar to the HK representation for the real-time one, would clearly be very useful for numerical studies in atom physics, molecule physics and chemical physics. Intuitively speaking, one might expect to obtain such an IVR through the aforementioned \(t \rightarrow -i\tau\) transformation from the HK representation of the real-time propagator [15]. However, as pointed out by Yan, Liu, and Shao in Ref. Yan?, and?, Shao?, this is not the case, as the result of this transformation diverges in the limit \(T \rightarrow \infty\). In Appendix 6, we demonstrate this again with detailed calculations. This problem is crucial, as the semi-classical approximation should be applicable in the high-temperature limit. Although several alternative IVRs for the imaginary-time propagator have been derived [15], Zhao?, and?, Yan?, and?, Shao?, an IVR based on coherent-state-type wave functions, similar to the HK IVR, has yet to be found.

In this work we solve the above problem for the systems with the absolute value of potential-energy gradient (force intensity) having a finite upper bound. For these systems we derive a reasonable HK-like semiclassical IVR for the Boltzmann operator, via an approach generalized from the one of Kay in Ref. [13]. Our IVR is exact for free particles and harmonic oscillators, although the potential-energy gradient of the latter system is un-bounded. The applicability for other systems are illustrated via numerical examples. Our results are helpful for studies of the properties of atomic systems, molecular systems, and chemical reactions at finite temperatures.

The remainder of this paper is organized as follows. For the convenience of the readers, we first display our HK-like IVR for the Boltzmann operator in Sec. 2, and then show the derivation of this IVR as well as the condition for the semi-classical approximation in Sec. 3. The applicability of our HK-like IVR is illustrated with some examples in Sec. 4. In Sec. 5 there is a summary. Some details of the derivations and calculations are given in the appendixes.

2 Central Result↩︎

We consider a general multi-particle quantum system, with Cartesian components of the coordinates and momenta being denoted as \(({p}_{1},...,{p}_{N})\) and \(({q}_{1},...,{q}_{N})\), respectively, and the Hamiltonian being given by \[H=\sum_{j=1}^{N}\frac{{p}_{j}^{2}}{2m_{j}}+V({\boldsymbol{q}}).\label{h}\tag{1}\] Here \({\boldsymbol{q}}=({q}_{1},...,{q}_{N})\), and \(m_{j}\) (\(j=1,...,N\)) is mass of the particle to which the coordinate \(q_{j}\) belongs. Furthermore, as mentioned above, we assume the norm of potential energy gradient, i.e., \(|\nabla_{\boldsymbol{q}}V({\boldsymbol{q}})|\), having a finite upper bound in the \({\boldsymbol{q}}\)-space.

The Boltzmann operator of our system is \(e^{-\hat{H}/(k_{B}T)}\), where \(\hat{H}\) is the quantum operator of the Hamiltonian in Eq. (1 ). The matrix element of this operator in the coordinate representation (imaginary-time propagator) can be expressed as \[\begin{align} K_{\tau}(\tilde{\boldsymbol{x}},{\boldsymbol{x}}):= \langle\tilde{\boldsymbol{x}}|e^{-\hat{H}/(k_{B}T)}|\boldsymbol{x}\rangle= \langle\tilde{\boldsymbol{x}}|e^{-\hat{H}\tau/\hbar}|\boldsymbol{x}\rangle, \end{align}\] where \(|\boldsymbol{x}\rangle\) and \(|\tilde{\boldsymbol{x}}\rangle\) are eigen-states of the coordinate operators, with eigen-values \(\boldsymbol{x}=(x_{1},...,x_{N})\) and \(\tilde{\boldsymbol{x}}=(\tilde{x}_{1},...,\tilde{x}_{N})\), respectively, and \[\begin{align} \tau=\frac{\hbar}{k_{B}T}.\label{tau} \end{align}\tag{2}\]

The HK-like IVR for the Boltzmann operator under the semiclassical approximation, which we have derived in this work, can be expressed as: \[\begin{align} K_{\tau}(\tilde{\boldsymbol{x}},{\boldsymbol{x}}) & = &A\int d\boldsymbol{p}d\boldsymbol{q}\left[D_\tau e^{-\frac{1}{\hbar}(S_\tau+B_\tau+C_\tau)}\right]. \label{asu2} \end{align}\tag{3}\] Here, \({\boldsymbol{q}}=({q}_{1},...,{q}_{N})\) (as defined above), \(\boldsymbol{p}=(p_1,...,p_N)\), and \(\int d\boldsymbol{p},d\boldsymbol{q}=\int dp_{1}\cdots dp_{N},dq_{1}\cdots dq_{N}\). Moreover, the factors \(S_\tau\), \(A\), \(B_\tau\), \(C_\tau\) and \(D_\tau\) are all real and independent of \(\hbar\), with the definitions being introduced in the following.

Factor \(S_\tau\)↩︎

The factor \(S_\tau\) of Eq. (3 ) is a function of \(\boldsymbol{p}\) and \(\boldsymbol{q}\), and is independent of \({\boldsymbol{x}}\) and \(\tilde{\boldsymbol{x}}\). It is defined as \[S_\tau(\boldsymbol{q},\boldsymbol{p})=\int_{0}^{\tau}d\eta\bigg[\sum_{j=1}^{N}\frac{{{p}}_{\eta,j}^{2}}{2m_{j}}+V({q}_{\eta,1},...,{q}_{\eta,N})\bigg].\label{s}\tag{4}\] Here \(q_{\eta,1},...,q_{\eta,N}\) and \(p_{\eta,1},...,p_{\eta,N}\) satisfy the Hamilton’s equations with inverted potential \(-V\), initial position \({\boldsymbol{q}}\) and initial momenta \({\boldsymbol{p}}\), i.e., the equations \[\begin{align} \frac{d}{d\eta}{q}_{\eta,j} & =\frac{{p}_{\eta,j}}{m_{j}};\tag{5}\\ \frac{d}{d\eta}{p}_{\eta,j} & =\frac{\partial V(z_{1},...,z_{N})}{\partial z_{j}} \bigg\vert_{z_1=q_{\eta, 1};...;z_N=q_{\eta, N}},\tag{6}\\ &(j=1,...,N),\nonumber \end{align}\] and the initial conditions \[\begin{align} q_{\eta=0,j}=q_{j};\;\;\;p_{\eta=0,j}={p_{j}},(j=1,...,N). \label{ice} \end{align}\tag{7}\]

According to the above definitions and equations, the factor \(S_\tau\) of Eq. (4 ) is the action of the classical trajectory \(q_{\eta,1},...,q_{\eta,N}\) and \(p_{\eta,1},...,p_{\eta,N}\). Moreover, in the following we will consider \(q_{\eta,1},...,q_{\eta,N}\) and \(p_{\eta,1},...,p_{\eta,N}\) as functions of \(\eta\) and \(\{{\boldsymbol{q}},{\boldsymbol{p}}\}\).

Factors \(A\), \(B_\tau\) and \(C_\tau\)↩︎

The factors \(B_\tau\) and \(C_\tau\) of Eq. (3 ) are functions of both \({\boldsymbol{p}},\boldsymbol{q}\) and \(\boldsymbol{x}, \tilde{\boldsymbol{x}}\). They are defined as \[\begin{align} B_\tau ({\boldsymbol{p}},\boldsymbol{q};\boldsymbol{x}, \tilde{\boldsymbol{x}}) & =\sum_{j=1}^{N}\bigg[\gamma_j\left({x}_{j}-{q}_{j}\right)^{2}+\gamma_j\left(\tilde{{x}}_{j}-{{q}}_{\tau,j}\right)^{2}\bigg.\nonumber\\ & \;\;\;\;\;\;\;\;\;\;\bigg.-{p}_{j}\left({x}_{j}-{q}_{j}\right)+{{p}}_{\tau,j}\left(\tilde{{x}}_{j}-{{q}}_{\tau,j}\right)\bigg];\tag{8}\\ C_\tau ({\boldsymbol{p}},\boldsymbol{q};\boldsymbol{x}, \tilde{\boldsymbol{x}})& =\frac{1}{\tau}\sum_{j=1}^{N}m_{j}\bigg[\left(\tilde{{x}}_{j}-{{q}}_{\tau,j}\right)-\left({x}_{j}-{q}_{j}\right)\bigg]^{2},\tag{9} \end{align}\] where \(\gamma_{1,...,N}\) are arbitrary positive parameters. Additionally, the factor \(A\) is defined as: \[\begin{align} A=\prod_{j=1}^N\left(\frac{\gamma_j}{2\hbar^{3}\pi^{3}}\right)^{1/2}.\label{cast} \end{align}\tag{10}\] Notice that the factor \(C_\tau\), as well as the signs of the terms \({p}_{j}\left({x}_{j}-{q}_{j}\right)\) and \({{p}}_{\tau,j}\left(\tilde{{x}}_{j}-{{q}}_{\tau,j}\right)\) of the factor \(B_\tau\), cannot be obtained from direct \(t\rightarrow -i\tau\) transformation on the HK representation of real-time propagator.

Factor \(D_\tau\)↩︎

Similar to \(S_\tau\), the factor \(D_\tau\) in Eq. (3 ) is also a function of \({\boldsymbol{p}},\boldsymbol{q}\), and is independent of \(\boldsymbol{x}\) and \(\tilde{\boldsymbol{x}}\). To show the definition of \(D_\tau\), we first introduce four \(N\times N\) matrices \({\rm R}^{qu}(\eta)\), \({\rm R}^{qv}(\eta)\), \({\rm R}^{pu}(\eta)\), and \({\rm R}^{pv}(\eta)\), with elements: \[\begin{align} {\rm R}_{ij}^{qu}(\eta) & = & -\frac{2}{\eta}m_{i}\delta_{ij}+\left(2\gamma_j+\frac{2}{\eta}m_{j}\right)\frac{\partial q_{\eta,j}}{\partial q_{i}}-\frac{\partial p_{\eta,j}}{\partial q_{i}};\nonumber \\ \tag{11}\\ {\rm R}_{ij}^{qv}(\eta) & = & 2\gamma_j\delta_{ij}-\frac{2}{\eta}m_{j}\left(\frac{\partial q_{\eta,j}}{\partial q_{i}}-\delta_{ij}\right);\\ {\rm R}_{ij}^{pu}(\eta) & = & \left(2\gamma_j+\frac{2}{\eta}m_{j}\right)\frac{\partial q_{\eta,j}}{\partial p_{i}}-\frac{\partial p_{\eta,j}}{\partial p_{i}};\\ {\rm R}_{ij}^{pv}(\eta) & = & \delta_{ij}-\frac{2}{\eta}m_{j}\frac{\partial q_{\eta,j}}{\partial p_{i}},\tag{12}\\ &&(i,j=1,...,N), \end{align}\] with \(\delta_{ij}\) being the Kronecker symbol. We further define other four \(N\times N\) matrices \({\rm T}^{uq}(\eta)\), \({\rm T}^{vq}(\eta)\), \({\rm T}^{up}(\eta)\), and \({\rm T}^{vp}(\eta)\), which relate to the these R-matrices via \[\left( \begin{array}{cc} {\rm T}^{uq}(\eta) & {\rm T}^{up}(\eta)\\ {\rm T}^{vq}(\eta) & {\rm T}^{vp}(\eta) \end{array} \right) = \left( \begin{array}{cc} {\rm R}^{qu}(\eta) & {\rm R}^{qv}(\eta)\\ {\rm R}^{pu}(\eta) & {\rm R}^{pv}(\eta) \end{array} \right)^{-1}.\label{dt}\tag{13}\]

The factor \(D_\tau\) of Eq. (3 ) can be expressed in terms of the elements of the above T-matrices, which are denoted as \({\rm T}^{to}_{ij}(\eta)\) (\(t=u,v\); \(o=p,q\); \(i,j=1,...,N\)). Specifically, we have \[\begin{align} D_\tau(\boldsymbol{q},\boldsymbol{p})= e^{-\int_{0}^{\tau}g_{\eta}^ {}d\eta}, \label{bigc} \end{align}\tag{14}\] with \(g_{\eta}\) being a function of \((\boldsymbol{q},\boldsymbol{p})\): \[\begin{align} &g_{\eta}(\boldsymbol{q},\boldsymbol{p})\nonumber \\ = & -\frac{1}{2}\sum_{i,j=1}^{N}{\rm W}_{ij}(\eta)V_{ij}(\eta)\nonumber \\ & +\sum_{j=1}^{N}\bigg\{\left(\frac{1}{\eta}+\frac{\gamma_j}{m_{j}}\right)+\left(\frac{2\gamma_j^{2}}{m_{j}}+\frac{m_{j}}{\eta^{2}}+\frac{4\gamma_j}{\eta}\right){\rm W}_{jj}(\eta)\nonumber \\ & \;\;\;\;\;\;\;\;\;\;+\left(\frac{2m_{j}}{\eta^2}+\frac{4\gamma_j}{\eta}\right){\rm T}_{jj}^{uq}(\eta)-\frac{m_{j}}{\eta^2}{\rm T}_{jj}^{vq}(\eta)\bigg\},\nonumber \\ \label{ftau} \end{align}\tag{15}\] where \[\begin{align} V_{ij}(\eta)&=\left.\frac{\partial^{2}}{\partial z_{i}\partial z_{j}}V(z_{1},...,z_{N})\right\vert_{z_1=q_{\eta,1};...;z_N=q_{\eta,N}},\\ {\rm W}_{ij}(\eta) & =-\sum_{s=1}^{N}\left[{\rm T}_{is}^{uq}(\eta)\frac{\partial q_{\eta,j}}{\partial q_{s}}+{\rm T}_{is}^{up}(\eta)\frac{\partial q_{\eta,j}}{\partial p_{s}}\right]. \end{align}\] According to this definition, calculating \(D_\tau\) requires computing the derivatives \(\frac{\partial q_{\eta,i}}{\partial q_j}\), \(\frac{\partial q_{\eta,i}}{\partial p_j}\), \(\frac{\partial p_{\eta,i}}{\partial q_j}\) and \(\frac{\partial p_{\eta,i}}{\partial p_j}\) (\(i,j=1,...,N\)) for \(0\leq\eta\leq \tau\). As shown in Appendix 7, one can derive these function via solving the Hamilton’s equations (5 , 6 ) together with another group of ordinary differential equations.

3 Derivation of Eq. (3 )↩︎

Now we show our approach for deriving the HK-like IVR for the imaginary-time propagator, i.e., Eq. (3 ). For simplicity, we consider the system of a single particle in one-dimensional space (\(N=1\)), and the derivation can be directly generalized to the cases with arbitrary \(N\). Consequently, we will omit the subscript denoting the particle index in the following. In the following we present the main ideas and framework of the derivation. The details of the derivations are provided in Appendix 8.

To derive Eq. (3 ), we express the imaginary-time propagator as the integration \[\begin{align} K_\tau(\tilde{x},x)=A\int dpdq\left[ F_D e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}\right], \label{test} \end{align}\tag{16}\] with \(A\), \(S_\tau\), \(B_\tau\) and \(C_\tau\) being defined in Eqs. (10 ), (4 ), (8 ) and (9 ), respectively, and \(F_D\) being a to-be-determined function of \(\tau\) and \((p,q)\), which is independent of \(x\) and \(\tilde{x}\).

In the following, we first prove that the function \(F_D\) satisfies the “initial condition" \(\lim_{\tau\rightarrow 0}F_{D}=1\). Then we derive the differential equation of \(F_D\), and solve this equation together with this”initial condition", using the semi-classical approximation. We will find that the solution is just \(F_D=D_\tau\), with \(D_\tau\) being given by Eq. (14 ) for \(N=1\). Substituting this result into Eq. (16 ), we obtain the result of (3 ) for \(N=1\).

3.1 “Initial Condition" of \(F_{D}\)↩︎

To prove \(\lim_{\tau\rightarrow 0}F_{D}=1\), we first calculate the limitation \(\lim_{\tau\rightarrow 0}A\int dpdq\; e^{-\frac{(S_\tau+B_\tau+D_\tau)}{\hbar}}\). Notice that since the potential-energy gradient intensity \(|dV(q)/dq|\) has a finite upper bound in the total \(q\)-space, in the short-\(\tau\) limit the solutions of the Hamiltonian equations (5 ) and (6 ) are just those for free motion. Specifically, in this limit we have \(q_\eta=q+\eta p/m\) and \(p_\eta=p\). Substituting this solution into the definitions (4 , 8 , 9 ) of \(S_\tau\), \(B_\tau\), and \(C_\tau\), we find that (Appendix 8.1) \[\begin{align} &&\lim_{\tau\rightarrow 0}A\int dpdq \; e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}\nonumber\\ &=&\lim_{\tau\rightarrow 0}\left({\frac{m}{2\pi\hbar\tau}}\right)^{1/2} e^{-\frac{m(x-\tilde{x})^2}{2\hbar\tau}}\tag{17}\\ \nonumber\\ &=&\delta(x-\tilde{x}).\tag{18} \end{align}\] On the other hand, it is clear that \(\lim_{\tau\rightarrow 0} K_\tau(\tilde{x},x)=\delta(x-\tilde{x})\). Thus, we have \[\begin{align} \lim_{\tau\rightarrow 0} K_\tau(\tilde{x},x)=\lim_{\tau\rightarrow 0}A\int dpdq\;e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}. \end{align}\] Comparing this result and Eq. (16 ), we find that the function \(F_D\) satisfies \[\begin{align} \lim_{\tau\rightarrow 0}F_{D}=1.\label{ic} \end{align}\tag{19}\]

3.2 Equation of \(F_D\) and Semi-Classical Approximation↩︎

Now we derive the differential equation for the function \(F_D\). To this end, we introduce the correction operator [16]–[18] for our system, which is defined as \[\begin{align} \hat{\Lambda}_{\tilde{x}}:= \hbar\frac{\partial}{\partial\tau}-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial\tilde{x}^{2}}+V(\tilde{x}) \label{lam2}. \end{align}\tag{20}\] It is clear that the imaginary-time propagator \(K_\tau(\tilde{x}, x)\) satisfies \(\hat{\Lambda} \left[ K_\tau(\tilde{x}, x) \right] = 0\). By substituting Eq. (16 ) into this equation, we find that the function \(F_D\) satisfies \[\begin{align} \hat{\Lambda}\bigg[\int dpdq\left( F_D e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}\right)\bigg]=0.\label{seqa1} \end{align}\tag{21}\] As detailed in Appendix 8.2, using the method generalized from Ref. [13], we find that the sufficient condition for Eq. (21 ) can be expressed as a differential equation for \(F_D\): \[\begin{align} \bigg(\hat{\tilde{L}}_0+\hbar \hat{\tilde{L}}_1+\hbar^2\hat{\tilde{L}}_2+...\bigg)F_D=0,\label{eqf} \end{align}\tag{22}\] where \(\hat{\tilde{L}}_{0,1,2,...}\) are \(\hbar\)-independent operators. Specially, the operator \(\hat{\tilde{L}}_0\) is given by \[\begin{align} \hat{\tilde{L}}_0=\frac{\partial }{\partial \tau}+g_\tau(q,p), \end{align}\] where \(g_\tau(q,p)\) is just the function given by Eq. (15 ) with \(N=1\).

Note that the operator in the l.h.s of Eq. (22 ) is expanded as a power series of \(\hbar\). Under the semi-classical approximation, we further ignore the terms proportional to \(\hbar^{n}\) (\(n \geq 1\)) in this series, and approximate Eq. (22 ) as \[\begin{align} \frac{\partial }{\partial \tau}F_D+g_\tau F_D=0.\label{eqf2} \end{align}\tag{23}\] Furthermore, it is clear that under the initial condition (19 ), the solution of Eq. (23 ) is just \[\begin{align} F_D= e^{-\int_{0}^{\tau}g_{\eta}^ {}d\eta}.\label{fcr} \end{align}\tag{24}\]

3.3 Final Derivation↩︎

Eq. (24 ) yields that \(F_D=D_\tau,\) with \(D_\tau\) being the one defined in Eq. (14 ) with \(N=1\). Substituting this result into Eq. (16 ), we finally obtain Eq. (3 ) for \(N=1\).

3.4 Condition of the Semi-Classical Approximation↩︎

At the end of this section, we discuss the condition underlying the semiclassical approximation employed in deriving our HK-like IVR, namely, the condition under which our HK-like IVR provides a good approximation to the exact matrix element of the Boltzmann operator. For comparison, recall that the semiclassical approximation used in the derivation of the van Vleck propagator requires the characteristic length scale \(l_V\) of the potential energy variation to be much larger than the de Broglie wavelength [19]. The condition adopted here is similar: specifically, for a system with \(N\) cooridnates, \(l_V\) should be much larger than both \(\sqrt{\hbar/\gamma_j}\) and \(\sqrt{\hbar\tau/m_j}=\sqrt{\hbar/(k_BTm_j)}\) (\(j=1,...,N\)). Thus, the semi-classical approximation works well in the high temperature systems. Additionally, in practical calculations for realistic systems, appropriate values of \(\gamma_j\) (\(j=1,...,N\)) should be chosen to ensure that this condition is satisfied.

4 Examples↩︎

In the previous two sections, we presented our HK-like IVR for the imaginary-time propagator and its derivation. In this section, we demonstrate its applicability. Specifically, we analytically prove that the IVR is exact for free particles and harmonic oscillators, and numerically illustrate its performance with an anharmonic system.

4.1 Free Particles↩︎

We first consider the system of \(N\) free particles, i.e., \(V=0\). As mentioned above, in this case the solutions of the Hamilton’s equations (5 , 6 ) are just \(q_{\eta,j}=q_j+\eta p_j/m_j\) and \(p_{\eta,j}=p_j\) (\(j=1,...,N\)). Consequently, the integration in the r. h. s. of Eq. (3 ) can be performed analytically. With straightforward calculations, we find that the right-hand side (r. h. s.) of Eq. (3 ) is \[\begin{align} \prod _{j=1}^N \left({\frac{m_j}{2\pi\hbar\tau}}\right)^{1/2} e^{-\frac{m_j}{2\hbar\tau}\left(x_j-\tilde{x}_j\right)^2}, \end{align}\] which is same as the exact imaginary-time propagator for this system. Thus, our HK-like IVR is exact for the free particles.

4.2 Harmonic Oscillators↩︎

Now we consider \(N\) harmonic oscillators with potential energy \(V=\sum_{j=1}^N m_j\omega_j^2 x_j^2/2\). The Hamilton’s equations (5 , 6 ) for this system can also be solved analytically, and we have \(q_{\eta,j}=c_je^{\omega_j\eta}+d_je^{-\omega_j\eta}\) and \(p_{\eta,j}=m_j\omega_j(c_je^{\omega\eta}-d_je^{-\omega\eta})\), where \(c_j=(p_j +m_j \omega_j q_j)/(2m_j\omega_j)\) and \(d_j=(m_j \omega_j q_j-p_j)/(2m_j\omega_j)\). Using these results, we can also analytically perform the integration in the r. h. s. of Eq. (3 ), and find that the r. h. s. of Eq. (3 ) is just (Appendix 9) \[\begin{align} \prod _{j=1}^N\left(\frac{m_j\omega_j}{2\pi\hbar\sinh(\omega_j \tau)}\right)^{1/2}e^{-\frac{m_j\omega_j\left[ \left(x_j^2+\tilde{x}_j^2\right)\cosh(\omega_j \tau)-2x_j\tilde{x}_j \right]}{2\hbar\sinh(\omega_j \tau)}},\nonumber\\ \label{ho3} \end{align}\tag{25}\] which is same as the exact imaginary-time propagator of these oscillators. Therefore, as for the free particles, our IVR is also exact for harmonic oscillators.

4.3 Anharmonic System↩︎

Finally, we consider a single particle in an anharmonic potential, with Hamiltonian \(H = p^2/(2m) + V(q)\). Here \(m\), \(p\) and \(q\) are the mass, momentum and coordinate of this particle, respectively. The anharmonic potential \(V(q)\) is given by \[V(q)=-\frac{m\omega^2L^2}{1-\frac{q^2}{2L^2}+ \frac{q^4}{4L^2q_0^2}}+V_0,\label{vq}\tag{26}\] where \(\omega\), \(q_0\), and \(L\) are positive parameters, and \(V_0\) is a \(q\)-independent constant: \[V_0=\frac{m\omega^2 L^2}{1-\frac{q_0^2}{4L^2}},\label{v0}\tag{27}\] which is chosen such that the minimum value of \(V(q)\) is zero. Specifically, \(\omega\) has the dimension of frequency, and both \(q_0\) and \(L\) have the dimension of length, and satisfy \(q_0<2L\).

In Fig. 1 (a) and Fig. 2 (a) we illustrate the potential \(V(q)\) for the following two groups of parameters:

• : \(L=10l_0\), \(q_0=6l_0\);

• : \(L=10l_0\), \(q_0=3l_0\),

where \(l_0\) is defined as \(l_0=\sqrt{\hbar/(m\omega)}\). It is clear that the gradient of this potential is bounded in the \(q\)-space. Furthermore, as shown in Fig. 1 (b) and Fig. 2 (b), in the region around \(q=0\), \(V(q)\) is a double-well with minimum points being localized at \(q=\pm q_0\).

We calculate the matrix elements of the Boltzmann operator for this particle, \(K(\tilde{x}, x) \equiv \langle \tilde{x} | e^{-\hat{H}/(k_B T)} | x \rangle,\) using the parameters (i) and (ii) for temperatures \(T = {10\hbar \omega}/{k_B}\), \(T = {\hbar \omega}/{k_B}\), and \(T = {0.5\hbar \omega}/{k_B}.\) Specifically, we derive the results from our HK-like IVR method with \(\gamma =5 {\hbar}/{l_0^2},\) as well as those from exact numerical diagonalization of the Hamiltonian operator \(\hat{H}\). These are denoted as \(K_{\rm IVR}(\tilde{x}, x)\) and \(K_{\rm Exact}(\tilde{x}, x)\), respectively. In Figs. 1 (c-h) and 2 (c-h), we compare \(K_{\rm IVR}(\tilde{x}, x)\) with \(K_{\rm Exact}(\tilde{x}, x)\) for \(\tilde{x} = \pm x\). It is shown that they are in perfect quantitative agreement.

In Fig. 3, we further compare \(K_{\rm IVR}(\tilde{x}, x)\) with \(K_{\rm Exact}(\tilde{x}, x)\) for additional values of \(\tilde{x}\) and \(x\), and present the relative error of our HK-like IVR approach, defined as \({|K_{\rm Exact}(\tilde{x}, x) - K_{\rm IVR}(\tilde{x}, x)|}/{|K_{\rm Exact}(\tilde{x}, x)|}.\) It is shown that for \(T = 10\hbar\omega / k_B\) and \(T = \hbar\omega / k_B\), the relative error remains below \(10^{-2}\) for both parameters (i) and (ii). Moreover, for \(T = 0.5\hbar\omega / k_B\), the relative error stays below \(10^{-2}\) for parameter (i), while it can reach up to \(8 \times 10^{-2}\) for parameter (ii). Specifically, for parameter (ii) with \(T = 0.5\hbar\omega / k_B\), the relative error is below \(4 \times 10^{-2}\) in most of the \((\tilde{x}, x)\) region, except at some places (the squares enclosed by black lines in Fig. 3 (r)) where \(K_{\rm Exact}(\tilde{x}, x)\) is very small compared to its maximum value in the entire \((\tilde{x}, x)\) domain.

These results clearly demonstrate the applicability of our HK-like IVR method. Furthermore, the relatively large relative error for parameter (ii) with \(T = 0.5\hbar\omega / k_B\) is consistent with the fact that the semi-classical approximation performs well at high temperatures, when the characteristic length scale \(q_0\) of the potential energy is much larger than \(\sqrt{\hbar / \gamma} = l_0/\sqrt{5}\), as discussed in Sec. 3.4.

5 Summary↩︎

In this work, we derive an HK-like semiclassical IVR for the Boltzmann operator, applicable to systems in which the gradient of the potential energy (i.e., the force) is bounded in real space, as shown in Eq. (3 ). Unlike the direct analytical continuation of the HK representation for the real-time propagator, our IVR converges in the high-temperature limit \(T \rightarrow \infty\) (\(\tau \rightarrow 0\)). Our IVR is exact for free particles and harmonic oscillators, and its applicability to other systems is demonstrated through examples. Our HK-like IVR method is useful for calculating the partition function and various physical properties of molecular systems in thermal equilibrium states in future research.

Acknowledgments↩︎

We express our deep gratitude to Prof. Jiushu Shao for the insightful discussions and for providing numerous important suggestions. This work is supported by the Innovation Program for Quantum Science and Technology (Grant No. 2023ZD0300700) and the National Key Research and Development Program of China (Grant No. 2022YFA1405300).

6 Divergence of Direct Extensions of HK Representation↩︎

In this appendix, we demonstrate that directly applying the \(t \rightarrow -i\tau\) transformation to the HK representation [10] results in divergence in the high-temperature limit. Without loss of generality, we consider the system of a single particle in one-dimensional (1D) space as an example. Our analysis can be easily extended to systems in arbitrary dimensions and with an arbitrary particle number.

We first consider the free-particle case with Hamiltonian \(\hat{H}=\hat{p}^2/(2m)\). In this case the HK representation of the real-time propagator \(\langle \tilde{x}|e^{-i\hat{H}t/\hbar}|x\rangle\) of this particle is proportional to \[\begin{align} \int dpdq e^{-\frac{\gamma(x-q)^2}{\hbar}-i\frac{p(x-q)}{\hbar}}e^{-\frac{\gamma(\tilde{x}-q_t)^2}{\hbar}+i\frac{p_t(\tilde{x}-q_t)}{\hbar}}e^{-i\frac{p^2}{2m\hbar}t},\label{HKf} \end{align}\tag{28}\] where \(\gamma\) is an arbitrary positive number, and \(q_{t}\) and \(p_{t}\) satisfy the classical Hamiltonian equations \[\begin{align} \frac{dq_{t}}{d{t}}=\frac{p_t}{m},\frac{dp_{t}}{d{t}}=0.\label{hef} \end{align}\tag{29}\] Now we apply the transformation \(t\rightarrow -i\tau\) to Eq. (28 ). Due to Eq. (29 ), this transformation leads to another transformations \(q_{t}\rightarrow q_{\tau}\) and \(p_{t}\rightarrow ip_{\tau}\), with \(q_{\tau}\) and \(p_{\tau}\) satisfying \(\frac{dp_{\tau}}{d{\tau}}=0\) and \(p_{\tau}=m\frac{dq_{\tau}}{d{\tau}}\), i.e., \[\begin{align} q_{\tau}=q+p{\tau}/m,\;\;\;\;\;p_{\tau}=p. \label{freea} \end{align}\tag{30}\] Substituting these results into Eq. (28 ), we find that the result of the transformation is \[\begin{align} I&=&\int dpdq e^{-\frac{\gamma(x-q)^2}{\hbar}-\frac{p(x-q)}{\hbar}}e^{-\frac{\gamma(\tilde{x}-q-p\tau/m)^2}{\hbar}-\frac{p(\tilde{x}-q-p\tau/m)}{\hbar}} e^{-\frac{p^2}{2m\hbar}\tau}\nonumber\\ &\propto& \int dp e^{\frac{(m\tau-\tau^2\gamma)}{2\hbar m^2}p^2+\frac{(x-\tilde{x})(m-\tau\gamma)}{\hbar m}p-\frac{(x-\tilde{x})^2\gamma}{2\hbar }}.\label{THK} \end{align}\tag{31}\] Clearly, \[\begin{align} I=\infty,{\rm for}\;\;\tau<m/\gamma. \end{align}\] Thus, the result of the transformation \(t\rightarrow -i\tau\) on HK representation diverges in the limit \(\tau\rightarrow 0\), i.e., the high-temperature limit \(T\rightarrow \infty\).

Next, we consider the systems with non-zero potential energy \(V(x)\). The Hamiltonian operator of such a system is \(\hat{H}=\hat{p}^2/(2m)+V(\hat{x})\). After the \(t\rightarrow -i\tau\) transformation, \(q_{\tau}\) and \(p_{\tau}\) satisfy the classical Hamiltonian equations with inverted potential, i.e., \(dq_\tau/d\tau=p_\tau/m\) and \(dp_\tau/d\tau=\frac{dV(x)}{dx}\vert_{x=q_\tau}\). As mentioned in the main text, we consider systems in which \(|dV(x)/dx|\) has a finite upper bound over the entire \(q\)-space. For these systems, in the high-temperature limit (\(\tau\rightarrow 0\)), the behaviors of \(q_{\tau}\) and \(p_{\tau}\) is always same as the ones of a free particle, i.e., still satisfying Eq. (30 ). As a result, the above analysis for this limit is still applicable. Thus, the result of the \(t\rightarrow -i\tau\) transformation on HK representation always diverges for \(\tau\rightarrow 0\) or \(T\rightarrow \infty\). Moreover, if \(|dV(q)/dq|\) is un-bounded, then \(q_\tau\) and \(p_\tau\) may increase with \(p\) even faster than the ones in Eq. (30 ). Therefore, the representation given by the transformation may still diverges.

We also notice that in Ref. [15] the authors propose another two extensions of HK representation, which are also based on the \(t\rightarrow -i\tau\) transformation. For a single free particle, they are proportional to \(\int dpdq e^{-\frac{\gamma(x-q)^2}{\hbar}-\frac{p(x-q)}{\hbar}}e^{-\frac{\gamma(x-q-p\tau)^2}{\hbar}-\frac{p(x-q-p\tau)}{\hbar}} e^{-\frac{p^2}{2m\hbar}\tau}\) and \(\int dpdq e^{-\frac{\gamma(x-q)^2}{\hbar}+\frac{p(x-q)}{\hbar}}e^{-\frac{\gamma(x-q-p\tau)^2}{\hbar}+\frac{p(x-q-p\tau)}{\hbar}} e^{-\frac{p^2}{2m\hbar}\tau},\) respectively. Using the same method as above, we directly find that both of these two extensions also diverge in the high-temperature limit, for systems in arbitrary dimensions and with arbitrary potential energy and particle numbers.

7 Equations for \(\frac{\partial q_{\eta,i}}{\partial q_j}\), \(\frac{\partial q_{\eta,i}}{\partial p_j}\), \(\frac{\partial p_{\eta,i}}{\partial q_j}\) and \(\frac{\partial p_{\eta,i}}{\partial p_j}\)↩︎

In this appendix we present the ordinary differential equations satisfied by \(\frac{\partial q_{\eta,i}}{\partial q_j}\), \(\frac{\partial q_{\eta,i}}{\partial p_j}\), \(\frac{\partial p_{\eta,i}}{\partial q_j}\) and \(\frac{\partial p_{\eta,i}}{\partial p_j}\). We first consider the 1D single-particle case (\(N=1\)). In this case \(\{q_{\eta},p_{\eta}\}\) satisfy the Hamilton’s equations \(\frac{d}{d\eta}q_{\eta} = \frac{ p_{\eta}}{m}\), \(\frac{d}{d\tau} p_{\eta} = d V(z)/dz\vert_{z=q_\eta},\) with initial conditions \(q_{\eta=0} = q\); \(p_{\eta=0}=p.\) We consider \(\{q_{\eta},p_{\eta}\}\) as functions of \((\eta,q,p)\), and define \[\begin{align} \xi^{(1)}(\eta)=\frac{\partial}{\partial q} q_{\eta}; \xi^{(2)}(\eta)=\frac{\partial}{\partial q} p_{\eta}. \end{align}\] By calculating \(\partial/\partial q\) in both sides of the aforementioned Hamilton’s equations satisfied by \(\{q_{\eta},p_{\eta}\}\), we find that \(\xi^{(1,2)}(\eta)\) satisfy the equations \[\begin{align} \frac{d}{d \eta} \xi^{(1)}(\eta)=\frac{\xi^{(2)}(\eta)}{m}; \frac{d}{d \eta} \xi^{(2)}(\eta)=\frac{d^2V(z)}{dz^2}\bigg\vert_{z=q_\eta} \xi^{(1)}(\eta). \label{deq} \end{align}\tag{32}\] Moreover, it is clear that \(\xi^{(1,2)}(\eta)\) also satisfy the initial conditions \[\begin{align} \xi^{(1)}(\eta=0)=1; \xi^{(2)}(\eta=0)=0. \label{dic} \end{align}\tag{33}\] Thus, one can derive \(\partial q_{\tau}/\partial q\) and \(\partial p_{\tau}/\partial q\) (i.e., \(\xi^{(1)}(\tau)\) and \(\xi^{(2)}(\tau)\)) by solving Eq. (32 ) with initial condition Eq. (33 ). Similarly, one can calculate \(\partial q_{\tau}/\partial p\) and \(\partial p_{\tau}/\partial p\) solving Eq. (32 ) with another initial condition \(\{\xi^{(1)}(\eta=0)=0; \xi^{(2)}(\eta=0)=1\}\). The solutions just relate to \(\partial q_{\tau}/\partial p\) and \(\partial p_{\tau}/\partial p\) via \(\partial q_{\tau}/\partial p=\xi^{(1)}(\tau)\) and \(\partial p_{\tau}/\partial p=\xi^{(2)}(\tau)\).

The above approach can be directly generalized to the cases with arbitrary \(N\). Specifically, to derive \(\frac{\partial q_{\eta,i}}{\partial q_j}\) and \(\frac{\partial p_{\eta,i}}{\partial q_j}\) for each fixed \(j\), one can solve equaitons \[\begin{align} \frac{d}{d \eta} \xi^{(1)}_k(\eta)=\frac{\xi^{(2)}_k(\eta)}{m_k}; \frac{d}{d \eta} \xi^{(2)}_k(\eta)=\sum_{s=1}^N \frac{\partial^2V(z_1,...,z_N)}{\partial z_k\partial z_s}\bigg\vert_{z_1=q_{\eta,1};z_2=q_{\eta,2};...;z_N=q_{\eta,N}} \xi^{(1)}_s(\eta), (k=1,...,N), \label{deq2} \end{align}\tag{34}\] with initial condition \(\{\xi^{(1)}_s(\eta=0)=\delta_{sj}, \xi^{(2)}_s(\eta=0)=0\;(s=1,...,N)\}\). The solutions are related to \(\frac{\partial q_{\eta,i}}{\partial q_j}\) and \(\frac{\partial p_{\eta,i}}{\partial q_j}\) via \(\frac{\partial q_{\eta,i}}{\partial q_j}=\xi^{(1)}_i(\eta)\) and \(\frac{\partial p_{\eta,i}}{\partial q_j}=\xi^{(2)}_i(\eta)\) (\(i=1,...,N\)). Similarly, to derive \(\frac{\partial q_{\eta,i}}{\partial p_j}\) and \(\frac{\partial p_{\eta,i}}{\partial p_j}\) for each fixed \(j\), one can solve Eq. (34 ), with another initial condition \(\{\xi^{(1)}_s(\eta=0)=0, \xi^{(2)}_s(\eta=0)=\delta_{sj}\;(s=1,...,N)\}\). The solutions for this initial condition are related to \(\frac{\partial q_{\eta,i}}{\partial p_j}\) and \(\frac{\partial p_{\eta,i}}{\partial p_j}\) via \(\frac{\partial q_{\eta,i}}{\partial p_j}=\xi^{(1)}_i(\eta)\) and \(\frac{\partial p_{\eta,i}}{\partial p_j}=\xi^{(2)}_i(\eta)\), (\(i=1,...,N\)).

8 Details of the Derivation of Eq. (3 ) for \(N=1\)↩︎

In this appendix we show the details of the derivation of Eq. (3 ) for a 1D single-particle system. As mentioned in Sec. 3, we express the matrix element of the Boltzmann operator as \[\begin{align} K_\tau(\tilde{x},x)=A\int dpdq\left[ F_D e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}\right], \label{tapp} \end{align}\tag{35}\] where the factors \(A\), \(S_\tau\), \(B_\tau\) and \(C_\tau\) are given by Eqs. (10 ), (4 ), (8 ) and (9 ), respectively, for \(N=1\), i.e., we have \[\begin{align} A&=&\left(\frac{\gamma}{2\hbar^{3}\pi^{3}}\right)^{1/2};\tag{36}\\ S_\tau&=&\int_{0}^{\tau}d\eta \bigg[ \frac{{p}_{\eta}^{2}}{2m}+V(q_{\eta})\bigg];\tag{37}\\ B_\tau & = & \gamma\left(x-q\right)^{2}+\gamma\left(\tilde{x}-{q}_{\tau}\right)^{2}-p\left(x-q\right)+p_{\tau}\left(\tilde{x}-q_{\tau}\right); \tag{38} \\C_\tau & = & \frac{m}{\tau}\left[\left(\tilde{x}-q_{\tau}\right)-\left(x-q\right)\right]^{2}.\tag{39} \end{align}\] Here \(p_{\eta}\) and \(q_{\eta}\) (\(0\leq \eta\leq\tau\)) satisfy the Hamilton’s equations and initial condition: \[\begin{align} \frac{d}{d\eta}q_{\eta} & = & \frac{ p_{\eta}}{m}; \frac{d}{d\eta} p_{\eta} = \frac{d V(z)}{dz}\bigg\vert_{z=q_\eta},\tag{40}\\ q_{\eta=0} & = & q; p_{\eta=0}=p.\tag{41} \end{align}\] Furthermore, in Eq. (35 ) \(F_D\) is a to-be-determined function of \((\tau,p,q)\).

8.1 Proof of Eqs. (17 , 18 )↩︎

Now we prove Eq. (17 , 18 ). To this end, we need to calculate the integration \(A\int dpdq \; e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}\) in the limit \(\tau\rightarrow 0\). Let us first perform the integration for \(p\). As mentioned in Sec. 3, since \(|dV(q)/dq|\) has a finite upper bound in the total real space, in the limit \(\tau\rightarrow 0\) the solution of the Hamilton’s equations (40 41 ) uniformly converges to the one of free motion, i.e., \(q_\eta=q+\eta p/m\) and \(p_\eta=p\). We substitute this solution into the definitions (37 , 38 , 39 ) of \(S_\tau\), \(B_\tau\), and \(C_\tau\), and then obtain \[\begin{align} A\int _{-\infty}^{+\infty}dp \; e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}=F_a(\tau)e^{-\frac{F_b(\tau,q)}{\hbar}}, \end{align}\] where \[\begin{align} F_a(\tau)&=&\frac{m}{\hbar\pi}\sqrt{\frac{\gamma}{\tau(m+2\tau\gamma)}};\\ F_b(\tau,q)&=&\frac{m^2 (x - \tilde{x})^2 + 2 m \tau (2 q^2 - 4 q x + 3 x^2 - 2 x \tilde{x} + {\tilde{x}}^2) \gamma + 4 \tau^2 (q - x)^2 \gamma^2}{2 \tau (m + 2 \tau \gamma)} +\int_0^\tau V\left(q+\frac{\eta p}{m}\right)d\eta. \end{align}\] Thus, we have \[\begin{align} \lim_{\tau\rightarrow 0}A\int dpdq \; e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}= \lim_{\tau\rightarrow 0}\int_{-\infty}^\infty dq F_a(\tau)e^{-\frac{F_b(\tau,q)}{\hbar}}.\label{inta} \end{align}\tag{42}\] Furthermore, in the limit \(\tau\rightarrow 0\), we can expand \(F_a\) and \(F_b\) in the integrand of Eq. (42 ) as powers of \(\tau\), and ignore the terms proportional to \(\tau^n\) (\(n> 0\)). This approach leads to \[\begin{align} \lim_{\tau\rightarrow 0}A\int dpdq \; e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}} = \frac{1}{\hbar\pi}\sqrt{\frac{m\gamma}{\tau}}\lim_{\tau\rightarrow 0} \int_{-\infty}^\infty dq e^{-\frac{1}{\hbar}\left[\frac{m\left(x - \tilde{x}\right)^2}{2\tau} + 2\left(q^2 \gamma- 2 q x\gamma + x^2 \gamma\right)\right]}.\label{inta2} \end{align}\tag{43}\] Performing the integration in the r. h. s. of Eq. (43 ) directly, we obtain Eq. (17 ) of Sec. 3, i.e., \[\begin{align} \lim_{\tau\rightarrow 0}A\int dpdq \; e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}} &=&\lim_{\tau\rightarrow 0}\left({\frac{m}{2\pi\hbar\tau}}\right)^{1/2} e^{-\frac{m(x-\tilde{x})^2}{2\hbar\tau}}.\label{inta1} \end{align}\tag{44}\] Additionally, substituting the result \(\left({\frac{m}{2\pi\hbar\tau}}\right)^{1/2} e^{-\frac{m(x-\tilde{x})^2}{2\hbar\tau}}=\langle \tilde{x} | e^{ -\frac{{\hat{p}}^2}{2m}\tau }|x\rangle\) into Eq. (44 ), we further obtain \[\begin{align} \lim_{\tau\rightarrow 0}A\int dpdq \; e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}&=&\lim_{\tau\rightarrow 0} \langle \tilde{x} | e^{ -\frac{{\hat{p}}^2}{2m}\tau }|x\rangle =\delta(x-\tilde{x}),\label{int2a} \end{align}\tag{45}\] which is just Eq. (18 ) of Sec. 3.

8.2 Derivation of Eq. (22 )↩︎

1. Preliminary Calculations↩︎

In the following we derive the equation of \(F_D\), i.e., Eq. (22 ) of Sec. 3. We begin from Eq. (21 ), i.e., \[\begin{align} \hat{\Lambda}\Big[\int dpdq\left( F_D e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}\right)\Big]=0.\label{seqa} \end{align}\tag{46}\] For our system with \(N=1\), the factors \(S_\tau\), \(B_\tau\) and \(C_\tau\) are given by Eqs. (37 ), (38 ) and Eq. (39 ), respectively. Additionally, the correction operator \(\hat{\Lambda}\) is given by Eq. (20 ), i.e., \[\begin{align} \hat{\Lambda}= \hbar\frac{\partial}{\partial\tau}-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial\tilde{x}^{2}}+V(\tilde{x}) \label{lam2a}. \end{align}\tag{47}\] Direct calculation yields \[\begin{align} \hat{\Lambda}\Big[\int dpdq\left( F_D e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}\right)\Big]=\int dpdq\left\{ F_D e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}\left[J_\tau+\hbar\frac{\dot{F_D}}{F_D}-\dot{S}_\tau+V(\tilde{{x}})\right]\right\} ,\label{c17} \end{align}\tag{48}\] where the dot above symbols means \[\begin{align} \dot{(...)}\equiv \frac{\partial (...)}{\partial\tau}, \end{align}\] and \[\begin{align} J_\tau & = & \hbar\left(\frac{1}{\tau}+\frac{\gamma}{m}\right)-\frac{ p_{\tau}^{2}}{2m}+\frac{2}{m} p_{\tau}\left\{ \frac{m\left[(x-q)-(\tilde{x}-q_{\tau})\right]}{\tau}-\gamma(\tilde{x}-q_{\tau})\right\} \nonumber \\ & & -\frac{2}{m}\left\{ \frac{m\left[(x-q)-(\tilde{x}-q_{\tau})\right]}{\tau}-\gamma(\tilde{x}-q_{\tau})\right\} ^{2}+\frac{m}{\tau^{2}}\left[(x-q)-(\tilde{x}-q_{\tau})\right]^{2}-(\tilde{x}-q_{\tau})\dot{ p}_{\tau}\nonumber \\ & & + p_{\tau}\dot{q}_{\tau}+2\gamma(\tilde{x}-q_{\tau})\dot{q}_{\tau}-\frac{2m}{\tau}\left[(x-q)-(\tilde{x}-q_{\tau})\right]\dot{q}_{\tau}.\label{bigj} \end{align}\tag{49}\] For the convenience of the following calculations, we define \[\begin{align} u & = & \tilde{x}-q_{\tau};\tag{50}\\ v & = & x-q,\tag{51} \end{align}\] and \[\begin{align} \Phi_\tau&=&-(B_\tau+C_\tau+S_\tau);\label{phitau}\\ V^{(n)}(z)&=&\frac{d^{n}}{dz^{n}}V(z). \end{align}\tag{52}\] Substituting Eqs. (40 , 37 ) into Eqs. (48 , 49 ), and using the above definitions and the fact \[V(\tilde{x})=V(q_\tau)+\sum_{n=1}^{\infty}u^n V^{(n)}(q_\tau),\] we find that Eq. (48 ) can be re-expressed as \[\begin{align} & & \hat{\Lambda}\Big[\int dpdq\left( F_D e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}\right)\Big] = \int dpdq\left\{ F_D \left[G_\tau+\hbar\frac{\dot{F_D }}{F_D }+\sum_{s=2}^{\infty}\frac{1}{s!}u^{s}V^{(s)}(q_{\tau})\right]e^{\Phi_\tau/\hbar}\right\} ,\label{res1} \end{align}\tag{53}\] where \[G_\tau=\hbar\left(\frac{1}{\tau}+\frac{\gamma}{m}\right)-\frac{2}{m}\left[\frac{m\left(v-u\right)}{\tau}-\gamma u\right]^{2}+\frac{m}{\tau^{2}}\left(u-v\right)^{2}.\]

2. Eliminating \(u\) and \(v\)↩︎

Now we eliminate the factors \(u\) and \(v\) in Eq. (53 ). Using the fact \[\begin{align} \frac{\partial S_\tau}{\partial q} & = & -p+\frac{\partial q_{\tau}}{\partial q} p_{\tau}; \frac{\partial S_\tau}{\partial p} = \frac{\partial q_{\tau}}{\partial p} p_{\tau}, \end{align}\] we find that that the factor \(\Phi_\tau\) defined in Eq. (52 ) satisfies \[\left(\begin{array}{c} \frac{\partial\Phi_\tau}{\partial q}\\ \\ \frac{\partial\Phi_\tau}{\partial p} \end{array}\right) =\left(\begin{array}{cc} {\rm R}^{qu} & {\rm R}^{qv}\\ {\rm R}^{pu} & {\rm R}^{pv} \end{array}\right) \left(\begin{array}{c} u\\ \\ v \end{array}\right). \label{pru}\tag{54}\] Here the parameters \({\rm R}^{qu}\), \({\rm R}^{qv}\), \({\rm R}^{pu}\) and \({\rm R}^{pv}\) are given by \[\begin{align} {\rm R}^{qu} & = & -\frac{2}{\tau}m+\left(2\gamma+\frac{2}{\tau}m\right)\frac{\partial q_{\tau}}{\partial q}-\frac{\partial p_{\tau}}{\partial q};\tag{55}\\ {\rm R}^{qv} & = & 2\gamma-\frac{2}{\tau}m\left(\frac{\partial q_{\tau}}{\partial q}-1\right);\\ {\rm R}^{pu} & = & \left(2\gamma+\frac{2}{\tau}m\right)\frac{\partial q_{\tau}}{\partial p_{}}-\frac{\partial p_{\tau}}{\partial p};\\ {\rm R}^{pv} & = & 1-\frac{2}{\tau}m\frac{\partial q_{\tau}}{\partial p_{}}.\tag{56} \end{align}\] We further define another four parameters \({\rm T}^{qu}\), \({\rm T}^{qv}\), \({\rm T}^{pu}\) and \({\rm T}^{pv}\) via the relation: \[\left(\begin{array}{cc} { \rm T}^{uq} & { \rm T}^{up}\\ { \rm T}^{vq} & { \rm T}^{vp} \end{array}\right) = \left(\begin{array}{cc} {\rm R}^{qu} & {\rm R}^{qv}\\ {\rm R}^{pu} & {\rm R}^{pv} \end{array}\right)^{-1}.\] Thus, Eq. (54 ) yields that \[\left(\begin{array}{c} u\\ \\ v \end{array}\right)= \left(\begin{array}{cc} { \rm T}^{uq} & { \rm T}^{up}\\ { \rm T}^{vq} & { \rm T}^{vp} \end{array}\right) \left(\begin{array}{c} \frac{\partial\Phi_\tau}{\partial q}\\ \\ \frac{\partial\Phi_\tau}{\partial p} \end{array}\right). \label{uve}\tag{57}\] Note that these R- and T-parameters are simply the R-matrices and T-matrices defined in Eqs. (11 13 ) of Sec. (2), with \(N=1\), respectively.

For convenience of the following calculations, we introduce differential operators \(\mathbb{D}_{u}\) and \(\mathbb{D}_{v}\): \[\begin{align} \mathbb{D}_{u} & = & { \rm T}^{uq}\frac{\partial}{\partial q}+{ \rm T}^{up}\frac{\partial}{\partial p}; \mathbb{D}_{v} = { \rm T}^{vq}\frac{\partial}{\partial q}+{ \rm T}^{vp}\frac{\partial}{\partial p}, \end{align}\] as well as the parameters \[\begin{align} Q_{\alpha\beta} & = &\mathbb{D}_{\alpha}\left[\beta\right];(\alpha,\beta=u,v).\label{dq} \end{align}\tag{58}\] Substituting the definitions (50 , 51 ) of \(u\) and \(v\) into Eq. (58 ), we further obtain the expressions of the \(Q\)-parameters: \[\begin{align} Q_{uu} & = & -{ \rm T}^{uq}\frac{\partial q_{\tau}}{\partial q}-{ \rm T}^{up}\frac{\partial q_{\tau}}{\partial p}; Q_{uv} = -{ \rm T}^{uq};\tag{59}\\ Q_{vu} & = & -{ \rm T}^{vq}\frac{\partial q_{\tau}}{\partial q}-{ \rm T}^{vp}\frac{\partial q_{\tau}}{\partial p}; Q_{vv} = -{ \rm T}^{vq}.\tag{60} \end{align}\]

Using the differential operators \(\mathbb{D}_{u}\) and \(\mathbb{D}_{v}\), we can re-express Eq. (57 ) as: \[\begin{align} \alpha e^{\Phi_\tau/\hbar} & = & \hbar\mathbb{D}_{\alpha}\left[e^{\Phi_\tau/\hbar}\right];(\alpha=u,v).\label{ruve} \end{align}\tag{61}\] Moreover, Eq. (61 ) leads to \[\begin{align} \alpha\beta e^{\Phi_\tau/\hbar}&=&\hbar\alpha\mathbb{D}_{\beta}\left[e^{\Phi_\tau/\hbar}\right]\nonumber\\ &=& \hbar\mathbb{D}_{\beta}\left[\alpha e^{\Phi_\tau/\hbar}\right]-\hbar e^{\Phi_\tau/\hbar}\mathbb{D}_{\beta}\left[\alpha \right],(\alpha,\beta=u,v).\label{ee1} \end{align}\tag{62}\] Substituting Eqs. (61 ) and (58 ) into Eq. (62 ), we further obtain \[\begin{align} \alpha\beta e^{\Phi_\tau/\hbar}=\hbar^2\mathbb{D}_{\beta}\left\{\mathbb{D}_{\alpha}\left[ e^{\Phi_\tau/\hbar}\right]\right\}-\hbar Q_{\beta\alpha}e^{\Phi_\tau/\hbar},(\alpha,\beta=u,v).\label{ee2} \end{align}\tag{63}\] i.e., \[\begin{align} u^2e^{\Phi_\tau/\hbar}&=&\hbar^{2}\mathbb{D}_{u}^2\left[e^{\Phi_\tau/\hbar}\right]-\hbar Q_{uu}e^{\Phi_\tau/\hbar},\\ v^2e^{\Phi_\tau/\hbar}&=&\hbar^{2}\mathbb{D}_{v}^2\left[e^{\Phi_\tau/\hbar}\right]-\hbar Q_{vv}e^{\Phi_\tau/\hbar},\\ uve^{\Phi_\tau/\hbar}&=&\hbar^{2}\mathbb{D}_{v}\mathbb{D}_{u}\left[e^{\Phi_\tau/\hbar}\right]-\hbar Q_{vu}e^{\Phi_\tau/\hbar}. \end{align}\] Repeating this technique, we can express any term of the form \(u^m v^n e^{\Phi_\tau/\hbar}\) (\(m,n=1,2,...\)) as a series in \(\hbar\), with each coefficient taking the form \(C_1 C_2 \dots [e^{\Phi_\tau/\hbar}]\), where each \(C_1, C_2, \dots\) is either a \(Q\)-factor or a \(\mathbb{D}\)-operator. Using this approach and Eqs. (59 , 60 ), we can re-express Eq. (53 ) as \[\begin{align} \hat{\Lambda}\Big[\int dpdq\left( F_D e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}\right)\Big] = \hbar\int dpdq\left\{ F_D \left[\hat{L}_{\tau}+\frac{\dot{F_D }}{F_D }\right]e^{\Phi_\tau/\hbar}\right\} ,\label{res1a} \end{align}\tag{64}\] where \[\hat{L}_{\tau}=\hat{L}_{\tau}^{(0)}+\hbar\hat{L}_{\tau}^{(1)}+\hbar^{2}\hat{L}_{\tau}^{(2)}+...\label{bigla}\tag{65}\] Here \(\hat{L}_{\tau}^{(0,1,2,...)}\) are \(\hbar\)-independent operators. Specially, we have \[\begin{align} \hat{L}_{\tau}^{(0)}&=&\left(\frac{1}{\tau}+\frac{\gamma}{m}\right)+\left(\frac{2\gamma^{2}}{m}+\frac{m}{\tau^{2}}+\frac{4\gamma}{\tau}\right)Q_{uu}+\left(\frac{2m}{\tau^{2}}+\frac{4\gamma}{\tau}\right){\rm T}^{uq}-\frac{m}{\tau^{2}}{\rm T}^{vq}-\frac{1}{2}Q_{uu}V^{(2)}(q_{\tau}),\nonumber\\ &=&g_\tau(q,p), \end{align}\] and \[\begin{align} \hat{L}_{\tau}^{(1)}&=&-\left(\frac{2\gamma^{2}}{m}+\frac{m}{\tau^{2}}+\frac{4\gamma}{\tau} -\frac{1}{2}V^{(2)}(q_{\tau}) \right)\mathbb{D}_u^2 +\left(\frac{2m}{\tau^{2}}+\frac{4\gamma}{\tau}\right)\mathbb{D}_u \mathbb{D}_v -\frac{m}{\tau^{2}}\mathbb{D}_v^2 \nonumber\\ && -\frac{1}{3!}V^{(3)}(q_{\tau}) \bigg( \mathbb{D}_{u}Q_{uu} +2Q_{uu} \mathbb{D}_{u} \bigg)+\frac{3}{4!}V^{(4)}(q_{\tau})Q_{uu}^2. \end{align}\] Notice that \(g_\tau(q,p)\) is just the function given by Eq. (15 ) of Sec. 2, with \(N=1\),

3. Equation of \(F_D\)↩︎

As in Ref. [13], by repeated integration by parts and using the fact that the integrated terms tend to zero for \(q\rightarrow \pm\infty\) or \(p\rightarrow \pm\infty\), we can further re-express Eq. (53 ) as \[\begin{align} \hat{\Lambda}\Big[\int dpdq\left( F_D e^{-\frac{(S_\tau+B_\tau+C_\tau)}{\hbar}}\right)\Big] = \hbar\int dpdq\left( e^{\Phi_\tau/\hbar} \hat{\tilde{L}}_{\tau}[F_D ]\right) ,\label{fntbkqiw} \end{align}\tag{66}\] where \[\hat{\tilde{L}}_{\tau}=\hat{\tilde{L}}_{\tau}^{(0)}+\hbar\hat{\tilde{L}}_{\tau}^{(1)}+\hbar^{2}\hat{\tilde{L}}_{\tau}^{(2)}+....\label{biglta}\tag{67}\] Here \(\hat{\tilde{L}}_{\tau}^{(0,1,2,...)}\) are also a group of \(\hbar\)-independent operators. For instance, we have \[\begin{align} \hat{\tilde{L}}_{\tau}^{(0)}&=&\frac{\partial}{\partial \tau}+\hat{L}_{\tau}^{(0)}=\frac{\partial}{\partial \tau}+g_\tau(q,p), \end{align}\] and \[\begin{align} \hat{\tilde{L}}_{\tau}^{(1)}&=&-\tilde{\mathbb{D}}_u^2\left(\frac{2\gamma^{2}}{m}+\frac{m}{\tau^{2}}+\frac{4\gamma}{\tau} -\frac{1}{2}V^{(2)}(q_{\tau}) \right) +\tilde{\mathbb{D}}_v \tilde{\mathbb{D}}_u \left(\frac{2m}{\tau^{2}}+\frac{4\gamma}{\tau}\right) -\tilde{\mathbb{D}}_v^2 \frac{m}{\tau^{2}} \nonumber\\ && -\frac{1}{3!} \bigg( Q_{uu} \tilde{\mathbb{D}}_{u} +2 \tilde{\mathbb{D}}_{u}Q_{uu} \bigg)V^{(3)}(q_{\tau}) +\frac{3}{4!}V^{(4)}(q_{\tau})Q_{uu}^2, \end{align}\] with \(\tilde{\mathbb{D}}_{u}\) and \(\tilde{\mathbb{D}}_{v}\) being defined as \[\begin{align} \tilde{\mathbb{D}}_{u}[...] & = & \frac{\partial}{\partial q}[{ \rm T}^{uq}...]+\frac{\partial}{\partial p}[{ \rm T}^{up}...]; \tilde{\mathbb{D}}_{v}[...] = \frac{\partial}{\partial q}[{ \rm T}^{vq}...]+\frac{\partial}{\partial p}[{ \rm T}^{vp}...]. \end{align}\] Substituting Eq. (64 ) and (67 ) into Eq. (46 ), we finally obtain Eq. (22 ) of of Sec. 3, i.e., \[\begin{align} \bigg(\hat{\tilde{L}}_0+\hbar \hat{\tilde{L}}_1+\hbar^2\hat{\tilde{L}}_2+...\bigg)F_D =0.\label{eqfa} \end{align}\tag{68}\]

9 Harmonic Oscillators↩︎

In this appendix we show that our HK-like IVR for the Boltzmann operator is exact for harmonic oscillators. We first consider a single harmonic oscillator with frequency \(\omega\), and use the natural unit \(\hbar=m=\gamma=1\). As shown in Sec. 4, for this system we have \(V=\omega^2 x^2/2\), and thus \(q_{\eta}=c e^{\omega\eta}+de^{-\omega\eta}\) and \(p_{\eta}=\omega(ce^{\omega\eta}-de^{-\omega\eta})\), where \(c=(p +\omega q)/(2\omega)\) and \(d=(\omega q-p)/(2\omega)\). Substituting these results into the definition of \(g_\tau\), we find that \[\begin{align} g_\tau=-\frac{-8 e^{\omega\tau} \omega + e^{2\omega\tau} (-2 + \omega) \big[4 - 4 \omega\tau + \tau^2 (-2 + \omega)\omega\big] + (2 + \omega) \big[4 + 4 \tau \omega + \tau^2\omega(2 + \omega)\big]}{2 \tau[-2 + e^{\tau \omega} (-2 + \omega) - \omega] \big\{4 + e^{\tau \omega} \big[-4 + \tau (-2 + \omega)\big]+ \tau(2 + \omega)\big\}}. \end{align}\] Thus, the function \(D_\tau=e^{-\int_0^\tau g_\eta d\eta}\) can be expressed as \[\begin{align} D_\tau=\frac{1}{2\sqrt{2}}\sqrt{ \frac{e^{-\tau \omega} \big[-2 + e^{\tau \omega} (-2 + \omega) - \omega\big] \big\{4 + e^{\tau \omega} \big[-4 + \tau (-2 + \omega)\big] + \tau (2 + \omega)\big\}}{\omega\tau} }.\label{hoe} \end{align}\tag{69}\] One can verify this result by substituting Eq. (69 ) into the equation \(\frac{d}{d\tau}D_\tau=-g_\tau D_\tau\) and the condition \(D_{\tau=0}=1\).

Substituting Eq. (69 ) and the above expressions of \(q_\eta\) and \(p_\eta\) into Eq. (3 ), we find that the integration \(\int dpdq...\) is just a Gaussian integration. Perform this integration analytically, we find that in the SI, the r. h. s. of Eq. (3 ) is just \(\left(\frac{m\omega}{2\pi\hbar\sinh(\omega \tau)}\right)^{1/2}e^{-\frac{m\omega\left[ \left(x^2+\tilde{x}^2\right)\cosh(\omega \tau)-2x\tilde{x} \right]}{2\hbar\sinh(\omega \tau)}}\), i.e., the exact imaginary-time propagator of this harmonic oscillator. Furthermore, the above calculation can be directly generalized to the general cases with \(N\) harmonic oscillators, and we can find that for these cases the r. h. s. of Eq. (3 ) is Eq. (25 ) of Sec. 4.

References↩︎