We develop two classes of quasi-classical dynamics that are shown to conserve the initial quantum ensemble when used in combination with the Feynman-Kleinert approximation of the density operator. These dynamics are used to improve the Feynman-Kleinert implementation of the classical Wigner approximation for the evaluation of quantum time correlation functions known as Feynman-Kleinert linearized path-integral. As shown, both classes of dynamics are able to recover the exact classical and high temperature limits of the quantum time correlation function, while a subset is able to recover the exact harmonic limit. A comparison of the approximate quantum time correlation functions obtained from both classes of dynamics is made with the exact results for the challenging model problems of the quartic and double-well potentials. It is found that these dynamics provide a great improvement over the classical Wigner approximation, in which purely classical dynamics are used. In a special case, our first method becomes identical to centroid molecular dynamics.

## I. INTRODUCTION

One of the main goals in quantum statistical mechanics is the calculation of quantum time correlation functions, which in the canonical ensemble takes the general form

*Z* being the partition function and *β* the inverse temperature 1/*k _{B}T*. The time correlation functions defined by Eq. (1) are of pivotal interest since they are accessible by most experimental spectroscopic techniques. For example, the dynamic structure factor, which is measured by either inelastic x ray or neutron scattering, the diffusion constant, and absorption spectra are just some of the quantities that can be related to quantum time correlation functions of the form of Eq. (1). Unfortunately, to compute these quantities exactly, we must either compute a full real time path-integral or solve the corresponding time-dependent Schrödinger equation, which is not practical for most realistic systems due to the intense computational load. In order to overcome this difficulty, it is desirable to develop approximate quantum methods which are as accurate and efficient as possible, and therefore enable the practical evaluation of Eq. (1). To date, there exists multiple approximation schemes including the Classical Wigner (CW) approximation,

^{1–3}Centroid Molecular Dynamics (CMD),

^{4}as well as Ring-Polymer Molecular Dynamics (RPMD)

^{5}in which semi-classical methods are used to calculate quantum time correlation functions approximately.

While all of these methods provide relatively accurate and practical approximations to Eq. (1), each has their own downfalls. For example, while both CMD and RPMD have been shown to be exact in the harmonic, high temperature, and short time limits, both of these methods begin to break down for correlation functions involving non-linear operators.^{4–6} In addition, the intrinsic dynamics of the ring-polymer in RPMD can introduce artificial frequencies in the spectrum of the correlation function, which become especially problematic when simulating absorption spectra.^{7} On the other hand, the CW approximation is also exact in the harmonic, high temperature, and short time limits, even for correlation functions involving non-linear operators.^{1,8} However, the CW approximation does not, in general, correctly produce time invariant thermodynamic properties of thermal equilibrium systems. Explicitly, for $ A \u02c6 =1$, the exact quantum expression in Eq. (1) has the property that

This is however generally not true for the CW approximation due to the classical dynamics used. While both CMD and RPMD also use a form of classical dynamics, these methods are able to retain the desirable property in Eq. (2) of the exact quantum time correlation function because the form of classical dynamics used in these methods conserves the initial quantum ensemble, while the purely classical dynamics used in the CW approximation do not. Recently, Liu and Miller^{8–10} have proposed a route to remedy this downfall of the CW approximation by replacing the classical propagation of the initial phase space distribution with a form of dynamics that conserves the initial quantum ensemble. Similarly, we set out to present a simple form of dynamics which conserve the ensemble and can be used to improve upon the Feynman-Kleinert (FK) implementation of the CW approximation known as the Feynman-Kleinert linearized path-integral (FK-LPI)^{1} method.

This paper is organized as follows: An introduction to the CW approximation and the FK approximation of the density operator is presented in Secs. II and III. In Sec. IV, we present the main idea of how to extend the FK-LPI dynamics scheme so as to conserve the ensemble. We then present two different classes of ensemble conserving dynamics in detail in Secs. V and VI. In Sec. VII, we apply these dynamics within the CW approximation for both the quartic and double-well potentials, and a comparison with the exact quantum time correlation function, as well as RPMD, CMD, and FK-LPI is made. The conclusions are presented in Sec. VIII.

## II. CLASSICAL WIGNER APPROXIMATION

The CW^{1–3} expression for a general quantum time correlation function of a system in one dimension is given by

where (*q _{t}*,

*p*) are coordinates in the classically evolved quantum phase space, (

_{t}*q*,

*p*) being the coordinates of the initial quantum distribution, and the Wigner transform of a general operator $ C \u02c6 $ is given by

The CW approximation can be shown to follow not only from a path linearization approximation implemented within the semiclassical initial value representation (SC-IVR) of the propagator^{2} but also from a linearization of the action difference between the forward-backward time paths in the corresponding exact path-integral expression for a general time correlation function.^{1,11} Furthermore, the CW approximation is exact in the limit that *t* → 0 since in this limit Eq. (3) reduces to the Wigner representation of $\u3008 A \u02c6 B \u02c6 \u3009$.

In fact, the CW approximation can be conjectured from the exact Wigner representation of the quantum time correlation function given by

since this approximation simply amounts to replacing the exact quantum evolution of $ B \u02c6 ( t ) = e i H \u02c6 t / \u0127 \u2009 B \u02c6 \u2009 e \u2212 i H \u02c6 t / \u0127 $ with the classical evolution of the quantum phase space such that

However, the classical evolution of the quantum phase space does not, in general, conserve the initial quantum phase space distribution, which results in thermodynamic properties of equilibrium systems being, incorrectly, time dependent. Liu and Miller^{8–10} have proposed a route to remedy this downfall of the CW approximation by replacing the classical propagation of the initial phase space distribution with a form of dynamics that preserves the time invariance of thermodynamic properties, such that $\u3008 B \u02c6 ( t ) \u3009=\u3008 B \u02c6 ( 0 ) \u3009$ in accord with the exact correlation function. As we show in Secs. V and VI, the CW expression in Eq. (3) can be made to share this property, if we use the FK approximation to the density operator and tie the dynamics of the phase-space points (*q _{t}*,

*p*) to the evolving centroid variables in the FK model. We present two classes of dynamics in Secs. V and VI that preserve the ensemble and can therefore be used to improve the CW approximation.

_{t}## III. FEYNMAN-KLEINERT DENSITY OPERATOR

Regardless of the dynamics used, the first step in invoking the CW approximation is obtaining the corresponding matrix elements of the Boltzmann operator since the Wigner transform of $ e \u2212 \beta H \u02c6 A \u02c6 $ is required. As in the FK-LPI method, we accomplish this by combining the effective frequency variational theory of Feynman^{12} and Kleinert^{13} with the quasidensity operator formalism of Jang and Voth.^{14} This FK approximation to the Boltzmann operator enables one to obtain an analytical expression for the Wigner transform of $ e \u2212 \beta H \u02c6 A \u02c6 $ which in turn allows for an efficient numerical evaluation of the CW approximation. The FK approximation of the density operator is exact in the harmonic and high temperature limits. Furthermore, the FK approximation to $ e \u2212 \beta H \u02c6 $ gives the best local harmonic approximation to the systems’ free energy and has been shown to be very accurate.^{1,15–19}

The FK approximation for $ e \u2212 \beta H \u02c6 $ in one dimension is given by

where (*x _{c}*,

*p*) are the classical centroid phase space variables describing the average position and momentum of a particle during thermal time

_{c}*β*ħ and are defined as

The FK approximation to the centroid phase space density *ρ _{FK}*(

*x*,

_{c}*p*) is given by

_{c}*W*_{1}(*x _{c}*) being the FK approximation to the centroid potential. The effective frequency quasidensity operator $ \delta \u02c6 F K ( x c , p c ) $ is

where the centroid dependent variational effective frequency Ω(*x _{c}*) is determined from the local curvature of the system’s Gaussian smeared potential and is given by

The smeared potential, *V*_{a2}(*x _{c}*), which accounts for quantum-statistical path fluctuations is given by

the smearing width *a*^{2}(*x _{c}*) being

which measures the importance of quantum fluctuations around the classical-like position *x _{c}*. In Eq. (10),

*α*(

*x*) is related to the smearing width

_{c}*a*

^{2}(

*x*) by

_{c} Note that when using Eq. (11) to determine the effective frequency Ω(*x _{c}*), the derivative is taken while treating

*a*

^{2}(

*x*) as a constant. Furthermore, by using Eqs. (11) and (12), one can write the explicit form that Ω

_{c}^{2}(

*x*) takes in terms of

_{c}*a*

^{2}(

*x*), and it is given by

_{c}Once Ω(*x _{c}*) is determined by solving Eqs. (13) and (15) iteratively, the FK approximation to the centroid potential

*W*

_{1}(

*x*) is given by

_{c}Using this FK approximation of the density operator, the Wigner distribution function, $ [ e \u2212 \beta H \u02c6 ] W ( q , p ) $, takes the form

where the Wigner transform of the QDO is written explicitly as

## IV. A NEW PRINCIPLE THAT PRESERVES THE CANONICAL ENSEMBLE UNDER PROPAGATION—AN INTUITIVE PICTURE

As it turns out, there are multiple ways of generating dynamics that conserve the initial quantum ensemble and can be used to evaluate the CW approximation. For example, Liu and Miller^{9} proposed a clever way of generating ensemble conserving dynamics by making an analogy to Liouville’s theorem in classical mechanics and set

and then proceed to generate the dynamics from this relation by choosing

from which it will follow that the effective force will be given by

However, if one takes this equilibrium Liouville dynamics^{9} route for the FK approximation of the Wigner distribution function then, as shown in Ref. 21, the dynamics involve an integration over the centroid distribution at each time step, which adds a significant computational load.

The main result of this paper is two new ensemble conserving extensions of the original FK-LPI method that indeed *are* computationally efficient. Before considering the detailed derivations of them, which involve some algebra, we present here the simplest of them, the Feynman-Kleinert Quasi-Classical Wigner method 1: FK-QCW(1) and leave the finer mathematical details for later.

We first consider why the original FK-LPI method conserves the canonical ensemble under time propagation in the case of a globally harmonic potential but fails to do so more generally. Understanding this problem guides us in formulating a corrected form of quantum dynamics that has the property that it always conserves the ensemble.

In the original formulation, the FK-LPI propagation of the Wigner distribution of the Boltzmann operator proceeds as follows. The centroids (*x _{c}*,

*p*) are sampled as classical-like variables and used for constructing a set of “local” quasi-density operators (QDOs). The Wigner transform of these operators $ [ \delta \u02c6 F K ( x c , p c ) ] W ( q , p ) $ is local Gaussian functions and they sum up to the global FK Boltzmann Wigner distribution (see Eq. (17)). Each of the $ [ \delta \u02c6 F K ( x c , p c ) ] W ( q , p ) $ provides a means to sample locally to obtain a set of phase-space points (

_{c}*x*,

*p*), and each point is then driven forward independently by classical dynamics. In Figure 1, we illustrate how four such points, sampled from a local QDO Wigner distribution $ [ \delta \u02c6 F K ( x \u2032 , p \u2032 ) ] W ( q , p ) $ centered at (

*x*′,

*p*′), are propagated.

By sampling and propagating a sufficient number of (*x*, *p*) points, the whole QDO Wigner distribution can be time-evolved (see Figure 1). From the evolved positions of the (*x*, *p*) points, a new centroid position (*x*″, *p*″) can be calculated. For harmonic systems, where the classical propagation is exact, the propagated QDO Wigner distribution is easily shown to retain its shape around the centroid position as time evolves. We note that for a globally harmonic potential, in the set of initial FK QDO Wigner distributions $ { [ \delta \u02c6 F K ( x c , p c ) ] W ( q , p ) } $, all members have the same shape. It therefore follows that a time-evolved QDO Wigner distribution, which is now centred at (*x*″, *p*″) is simply equal to $ [ \delta \u02c6 F K ( x \u2033 , p \u2033 ) ] W ( q , p ) $, i.e., it is equal to the original FK QDO Wigner function that was centred at (*x*″, *p*″) at the initial time. Moreover, the statistical weight associated with the centroid position, i.e., $exp { \u2212 \beta ( p c 2 2 m + W 1 ( x c ) ) } $ (see Eq. (9)) is identical when evaluated for (*x _{c}*,

*p*) = (

_{c}*x*′,

*p*′) and (

*x*,

_{c}*p*) = (

_{c}*x*″,

*p*″). The propagated QDO Wigner function is, therefore, not only still a FK QDO Wigner function from the original set but also of equal weight as the one that at the initial time was associated with (

*x*″,

*p*″). Hence, the ensemble is conserved. We will refer to this property—that a propagated QDO Wigner distribution replaces one from the initial time with a QDO Wigner distribution with identical shape and weight—as a replacement mapping. This property holds for classical dynamics and globally harmonic potentials but need not for other potentials (see Figure 2). Therefore, in general, the FK approximation to the canonical ensemble is not conserved under classical dynamics. Our new algorithm FK-QCW(1) is designed specifically so that replacement mapping holds for a general potential.

In FK-QCW(1), the centroids, rather than the individual (*x*, *p*) phase space points, are driven forward in time and this is done by centroid molecular dynamics using the Feynman-Kleinert centroid potential *W*_{1}(*x _{c}*). The motion of the phase-space points (

*x*,

*p*) is determined

*relative*to their centroid. Note that the dynamics of the centroid point do not depend on the dynamics of the (

*x*,

*p*) phase-space points that evolve around it. The relative motion is chosen as harmonic, and as frequency, we pick the effective FK frequency Ω(

*x*). Nevertheless, as in the original FK-LPI method, the direct application of this idea will not lead to conservation of the ensemble since the FK effective harmonic frequencies depend on the centroid position (

_{c}*x*″,

*p*″). Figure 2 still applies. However, it turns out that by performing the relative dynamics in dimensionless coordinates defined in terms of the local effective frequencies, all FK QDO Wigner functions will have an identical shape, and (see below) the equations of motion solved in these scaled coordinates introduce the replacement mapping of the FK QDOs. Hence, with a position-dependent scaling of the coordinates, we can recover ensemble conservation.

## V. THE FEYNMAN-KLEINERT QUASI-CLASSICAL WIGNER METHOD: FK-QCW(1)

The FK approximation to the canonical ensemble average of an operator $ B \u02c6 $ is

The phase-space points (*q*, *p*) may be used as starting points for trajectories, propagated by, e.g., classical dynamics, as in the classical Wigner model or by some other form of dynamics. (*q*, *p*) is sampled “locally” from a QDO Wigner function $ [ \delta \u02c6 F K ( x c , p c ) ] W ( q , p ) $. Depending on the dynamics, the ensemble average $ B \u02c6 ( t ) $ may or may not be time-dependent, and it becomes

In passing, we notice that it is important to always keep (*q _{t}*,

*p*) in the inner integral, as we will in a moment let (

_{t}*q*,

_{t}*p*) depend on the centroids (

_{t}*x*,

_{c}*p*) of the QDO Wigner function it was sampled from. We will find a criterium for (

_{c}*q*,

_{t}*p*) so that $ B \u02c6 ( t ) = B \u02c6 $ for an arbitrary operator $ B \u02c6 $, whereby the ensemble is conserved. As discussed in Sec. IV, we will let the centroids (

_{t}*x*,

_{c}*p*) of the QDO Wigner functions evolve by centroid molecular dynamics, i.e.,

_{c} and formulate the motion of (*q _{t}*,

*p*) relative to (

_{t}*x*(

_{c}*t*),

*p*(

_{c}*t*)). In the following, we abbreviate $ \delta \u02c6 F K ( t ) \u2261 \delta \u02c6 F K ( x c ( t ) , p c ( t ) ) $. Let us establish a requirement on the dynamics of the unknown (

*q*,

_{t}*p*) trajectories if they are to conserve the ensemble. As shown in Appendix A, if the trajectories (

_{t}*q*,

_{t}*p*) do not cross, then Eq. (22) can be written equivalently by virtue of the multi-dimensional change-of-variables theorem

_{t}^{20}as

where *J*_{(qt,pt)} is the Jacobian pertaining to the coordinate transform from old (*q*, *p*) to propagated coordinates (*q _{t}*,

*p*),

_{t} see Appendix A. Taking the time-derivative of *ρ _{FK}*(

*x*(

_{c}*t*),

*p*(

_{c}*t*)) and using Eqs. (24) and (25), it follows that

*ρ*(

_{FK}*x*(

_{c}*t*),

*p*(

_{c}*t*)) =

*ρ*(

_{FK}*x*,

_{c}*p*). Then, Eq. (26) may be written as

_{c} If Eq. (29) is satisfied, then the FK ensemble will be invariant under time-propagation of the (*q _{t}*,

*p*) trajectories. The Jacobian in Eq. (27) is explicitly given by

_{t}^{9}

As we now show, we are able to find a set of dynamics that fulfil this relation and therefore is guaranteed to conserve the ensemble. This is accomplished by working in terms of the dimensionless relative coordinates

Thus, using these relations in Eq. (29), we see that the ensemble will be conserved as long as the dynamics fulfil

Therefore, if we require the dynamics to satisfy

or equivalently

and furthermore require

then the ensemble will be conserved.

This then gives us two conditions that must be met. However, it does not give us a way to generate dynamics without first assuming a form for either $ q \u0303 \u0307 ( t ) $ or $ p \u0303 \u0307 ( t ) $. Therefore, to guide us in this matter, we look to the harmonic limit in which Ω is constant and the exact dynamics of the centroid and quantum phase space variables are given by their corresponding classical equations. In this limit, it is straightforward to show that the exact relation for $ q \u0303 \u0307 ( t ) $ is

Thus, assuming this harmonic form and furthermore introducing an arbitrary frequency function *f* such that for a general potential, $ q \u0303 \u0307 ( t ) $ takes the form

then from the condition in Eq. (36), we have that

Furthermore, the second condition in Eq. (37) that must be met for these dynamics to conserve the ensemble is

However, by using Eq. (31), we can make a change of variables for the partial derivatives and this condition can be equivalently written as

which we recognize is always satisfied as long as our arbitrary function *f* does not depend on *q _{t}* or

*p*. Therefore, the dynamics of the general form

_{t}Therefore, since we are able to introduce an arbitrary function *f*(*x _{c}*(

*t*),

*p*(

_{c}*t*)) into the dynamics, we have found an entire class of dynamics that conserves the ensemble and can be used within the CW expression for the time correlation function. Furthermore, the direct propagation of the (

*q*,

_{t}*p*) variables through the dynamics in Eq. (44) is not necessary due to the fact that once the dimensionless and centroid variables have been propagated according to Eqs. (43), (24), and (25), respectively, we can obtain the instantaneous quantum phase space variables through

_{t}There is one clear drawback to the entire class of dynamics since the momentum becomes imaginary when the FK effective frequency takes on imaginary values. Explicitly, when $ \Omega t 2 <0$, the Ω_{t}*α _{t}* term in Eq. (46) becomes negative, which does not allow us to transform from $ p \u0303 ( t ) $ to

*p*. Thus, there exists a limitation in these dynamics; since in order to use them for barriers, one must adopt the convention to set

_{t}*p*=

_{t}*p*(

_{c}*t*), which is the limiting value of

*p*for Ω

_{t}_{t}→ 0, as seen from Eq. (46). While this convention does not present a problem when used within the CW expression for the time correlation of linear operators, it would be desirable to have a similar method that conserves the ensemble and is also able to handle imaginary frequencies in a robust way. This is the motivation for an additional class of dynamics which are presented in Sec. VI. We also note that, as we show in Appendix C, the entire class of dynamics described above when used within the CW approximation for the Kubo-transformed time-correlation function reduces to CMD for linear operators.

We next consider which form of the *f* function that is appropriate.

### A. Harmonic limit

Although a whole class of dynamics exists which conserve the initial quantum ensemble, one can easily conjecture that only a subset of these would actually be useful and provide a reasonable approximation to the real quantum dynamics. Therefore, we now provide a way to tailor the general frequency function introduced into the dynamics such that we recover exact quantum dynamics in the harmonic limit.

For a harmonic potential, Ω_{t} and *α _{t}* are simply constants. Furthermore, as we alluded to in Eq. (38), the exact dynamics of the dimensionless variables in Eq. (31) in the harmonic limit take the form

then these dynamics will be exact in the harmonic limit since in this limit Ω(*x _{c}*(

*t*)) → Ω. Furthermore, since the CW approximation is exact in this limit, the time correlation function will also be exact.

### B. Classical and high temperature limits

In both the classical (ħ → 0) and high temperature (*β* → 0) limits, it is easy to show from Eq. (46) that *q _{t}* and

*p*become fixed at

_{t}*x*(

_{c}*t*) and

*p*(

_{c}*t*), respectively. Furthermore, the centroid equations of motion reduce to the corresponding classical equations of motion since the smearing width

*a*

^{2}(

*x*(

_{c}*t*)) → 0 in both of these limits. Using these properties, it is straightforward to show that when the FK approximation to the density operator is used, the CW expression for the time correlation function in Eq. (3) reduces to the corresponding classical expression. Thus, the entire class of dynamics in Eq. (43) is able to recover both the classical and high temperature limits of the quantum time correlation function exactly.

### C. Low temperature limit

At low temperatures, the centroid distribution in Eq. (9) singles out *p _{c}* = 0 as well as the

*x*value corresponding to the global minimum of the centroid potential

_{c}*W*

_{1}(

*x*), denoted $ x c m $. Therefore, from Eqs. (24) and (25), we have that (

_{c}*x*(

_{c}*t*),

*p*(

_{c}*t*)) will be fixed at $ ( x c m , 0 ) $, and the dynamics of Eq. (44) reduce to

*f*= Ω(

*x*(

_{c}*t*)), we have in the low temperature limit that $f\u2192\Omega ( x c m ) $ and the solution to Eq. (49) is

Therefore, since in the low temperature limit, the exact dynamics are dominated by the ground and first excited states, then, similar to CMD,^{22} this method produces coherent quantum dynamics in this limit when the FK effective frequency is used in the dynamics.

## VI. AN ADDITIONAL FEYNMAN-KLEINERT QUASI-CLASSICAL WIGNER METHOD: FK-QCW(2)

The imaginary frequency problem encountered in the dynamics of Sec. V occurs because the Ω_{t}*α _{t}* factor in Eq. (31) becomes negative for imaginary frequencies. Looking back at how we defined the dimensionless variables, we see that the root of this problem is the exp{ − (

*p*−

_{t}*p*(

_{c}*t*))

^{2}/

*m*ħΩ

_{t}

*α*} factor contained within the FK Wigner distribution function in Eq. (18), which in addition becomes divergent for imaginary frequencies. However, the FK Wigner distribution function is Gaussian in the centroid momentum, and after integrating this variable out, it takes the form

_{t}where we have defined

and

One can show that this form is well defined for imaginary frequencies as long as |Ω_{0}| < *π*/ħ*β*.

As we now show, we are able to use this form to find another class of ensemble conserving dynamics that are well defined for this range of imaginary frequencies. This is accomplished by multiplying Eq. (51) by

after which we obtain

We may use this expression for calculating the time-independent ensemble average of an operator $ B \u02c6 $,

Similarly as in Sec. V, we may perform a change of integration variables in this equation, by appealing to the multi-dimensional change of variables theorem, see Appendix A. Equation (56) then becomes

where the Jacobian is again given by Eq. (27). As noticed before, *ρ _{FK}*(

*x*(

_{c}*t*),

*p*(

_{c}*t*)) =

*ρ*(

_{FK}*x*,

_{c}*p*) and so we arrive at

_{c} We may also write the possibly time-dependent ensemble average of $ B \u02c6 $ using the propagated trajectories (*q _{t}*,

*p*) as

_{t}If Eq. (59) is identical to Eq. (58), then the canonical ensemble is “conserved” under time-propagation. Comparing these equations, we see that this is the case, if

However, similar to our previous analysis in Sec. V, this condition for the dynamics to conserve the ensemble can be greatly simplified by working in new dimensionless coordinates defined as

Furthermore, one can show, by using an argument similar to that in Appendix B, that the absolute value of the Jacobian in terms of these new dimensionless coordinates takes the form

Therefore, using these relations in Eq. (60), we have that the ensemble will be conserved as long as these new dynamics fulfil

which we recognize is the same condition as in Eq. (34), only now written in terms of these new dimensionless variables. Therefore, as long as these new dynamics simultaneously fulfil

and

the ensemble will be conserved.

Thus, after assuming the same harmonic like form of $ q \u0303 \u0307 ( t ) $ as in Sec. V, one can show that these two conditions determine these new dynamics as

Thus, we have found another class of dynamics which conserve the ensemble and can be used within the CW approximation of the quantum time correlation function. Furthermore, the direct propagation of the (*q _{t}*,

*p*) variables is again not necessary for this new class of dynamics since the instantaneous quantum phase space variables can be obtained through

_{t}*p*is real as long as |Ω

_{t}_{t}| <

*π*/ħ

*β*and we are therefore able to use these dynamics in the presence of barriers. Furthermore, by adopting the same convention as FK-LPI

^{1}in which we set

*p*= 0 for frequencies outside of this range, we are able to apply this method for imaginary frequencies as long as |Ω

_{t}_{t}| < 2

*π*/ħ

*β*, which is the entire range where the FK approximation to the thermal density operator is well defined. Thus, this new class of dynamics allows us to handle imaginary frequencies in a robust way. Similar to the dynamics in Sec. V, the multi-dimensional generalization of this method is relatively straightforward since it only amounts to propagating the dimensionless coordinates written in terms of the mass-weighted normal mode coordinates defined in Ref. 1 through the dynamics in Eq. (67).

### A. Harmonic limit

In general, there does not exist a frequency function, *f*(*x _{c}*(

*t*)), such that we are able to recover the harmonic limit exactly for these dynamics. This hindrance can be traced back to the way that we defined $ p \u0303 ( t ) $, in which

*p*(

_{c}*t*) is absent. However, if one chooses for the frequency function the FK effective frequency

*f*= Ω

_{t}, then one can show that in the harmonic limit these dynamics used within the CW approximation give for the position autocorrelation function

while the exact expression is

Thus, one sees that these dynamics correctly reproduce the real part of the position autocorrelation function, but not the imaginary part. In addition, one can show that

and

such that Eq. (71) reduces to the exact expression in both the high and low temperature limits.

### B. Classical and high temperature limits

Using an argument similar to that presented in Sec. V B, one can show that this entire class of dynamics produces the correct classical and high temperature limits of the quantum time correlation function for position dependent operators when used within the CW approximation.

### C. Low temperature limit

For low temperatures, the dynamics in Eq. (67) with *f* = Ω(*x _{c}*(

*t*)) reduce to the same low temperature limit as Eq. (50) in Sec. V C. Thus, it follows that when the FK effective frequency is used in the dynamics of Eq. (67), one also obtains coherent quantum dynamics in the low temperature limit that approximate the exact dynamics.

## VII. APPLICATION TO MODEL PROBLEMS

As we have shown, the CW approximation in Eq. (3) evaluated with the dynamics of Secs. V and VI, which we term FK-QCW(1) and FK-QCW(2), respectively, is able to recover the exact classical and high temperature limits of the quantum time correlation function. In addition, the FK-QCW(1) method is able to recover the exact harmonic limit when the function *f* = Ω_{t} is used in the dynamics. To find out how well these dynamics perform outside of these limits, we compute the quantum time correlation function of two challenging systems, the quartic potential and the double-well potential, in which quantum coherence effects are very important. However, since the CW approximation in Eq. (3) lacks the phase information necessary to capture these effects, the main goal of studying these model problems is to see how FK-QCW(1) and FK-QCW(2) compare to exact results and other approximate methods over relatively short times. This is because ultimately we would like to apply these dynamics to realistic condensed phase systems, in which the coupling to the many degrees of freedom quenches long time quantum coherence effects,^{1,23} and the behavior of the correlation function over relatively short times is most important. To this end, we compute both linear and non-linear correlation functions for these two model systems using FK-QCW(1) and FK-QCW(2) and then provide a comparison to the exact results, as well as to those obtained by RPMD, CMD, and FK-LPI. It should be noted that the results we present here can also be compared with those previously obtained by Jang and Voth,^{4} Craig and Manolopoulos,^{5} and Liu.^{10}

In all of the following simulations, we have used natural units in which *m* = ħ = *k _{B}* = 1. The quantum time correlation function was evaluated using the FK-QCW(1) and FK-QCW(2) methods by performing a molecular dynamics simulation of the centroid variables, which are propagated according to Eqs. (24) and (25). For each centroid trajectory, 100 $ ( q \u0303 ( 0 ) , p \u0303 ( 0 ) ) $ values were sampled according to the Gaussian distribution $exp ( \u2212 q \u0303 ( 0 ) 2 \u2212 p \u0303 ( 0 ) 2 ) $. For FK-QCW(1) (FK-QCW(2)), these dimensionless coordinates were then propagated through Eq. (43) (Eq. (67)) and then by using Eq. (46) (Eq. (70)) the instantaneous (

*q*,

_{t}*p*) values were obtained. For all of the time correlation functions presented in this section, we constructed the CW approximation of $\u3008 A \u02c6 ( 0 ) B \u02c6 ( t ) \u3009$ in Eq. (3) by averaging over 50 000 consecutive centroid trajectories, in which the centroid momentum was resampled at the beginning of each trajectory and a time step of 0.005 was used. To illustrate the applicability of multiple sets of dynamics, we have performed the simulations using two different frequency functions,

_{t}and

where Ω_{t} = Ω(*x _{c}*(

*t*)) is the FK effective frequency and

*f*

_{2}is a temperature dependent effective frequency such that

*f*

_{2}≥ Ω

_{t}, and the equality is met in the limit that

*T*→ 0.

The evaluation of CMD using the FK approximation of the density operator was accomplished by also performing a molecular dynamics simulation of the centroid variables in which the centroid momentum was resampled at the beginning of each trajectory and a time step of 0.005 was used. The Kubo transformed time correlation function was constructed by averaging over 50 000 consecutive centroid trajectories. For the implementation of RPMD, a molecular dynamics simulation in the extended ring-polymer phase space was performed in which we used *n* = 48 beads. At the beginning of each trajectory, the momentum was resampled from the Boltzmann distribution at an inverse temperature of *β*/*n*, and a time step of 0.005 was used. The Kubo transformed time correlation function was also constructed by averaging over 50 000 consecutive trajectories. For both RPMD and CMD, the standard time correlation function was obtained from the Kubo transformed quantity by inverse Fourier transforming the relation in Eq. (C2).

### A. The double-well potential

The first model system is the double-well potential in which

This system presents a very challenging case due to the presence of a barrier, the maximum of the barrier being located at *x* = 0. For imaginary frequencies encountered within FK-QCW(1), we used the convention to set *p _{t}* =

*p*(

_{c}*t*). For imaginary frequencies |Ω| ≥

*π*/ħ

*β*encountered within FK-QCW(2), we set

*p*= 0 which is the limiting value for |Ω

_{t}_{t}| →

*π*/ħ

*β*in Eq. (70).

In Fig. 3, we present a representative trajectory that really shows the differences between these new ensemble conserving dynamics and the classical dynamics used within the normal CW approximation. As seen, the coupling of the quantum phase space to the time evolved centroid through the dynamics in Eq. (67) results in relative motion about the centroid. Furthermore, one sees that when the initial momentum is not enough to make it over the barrier for the classical dynamics, a particle undergoing relative motion about the centroid is essentially pulled over the barrier by the centroid, which evolves on a smeared potential.

As seen in Fig. 4, for $\u3008 x \u02c6 ( 0 ) x \u02c6 ( t ) \u3009$ at *β* = 8, the exact dynamics exhibit significant coherent tunneling, and none of the approximate methods are able to reproduce this behavior since they all show significant dephasing. While FK-QCW(1) and FK-QCW(2) give results comparable to RPMD, FK-CMD penetrates the potential barrier the deepest of the approximate methods. However, as we show in Appendix C, the entire class of dynamics in FK-QCW(1) reduces to FK-CMD when used in the expression for the Kubo transformed time correlation function for linear operators. This then suggests that there exists an advantage in obtaining the standard correlation function from the Kubo transformed version when possible (e.g., linear operators). In comparing FK-QCW(1) and FK-QCW(2) to FK-LPI, one sees that even in this challenging case where none of the approximate methods are able to correctly produce the behavior of the exact correlation function except at very short times, the ensemble conserving dynamics used within these two new methods improves the accuracy of the CW approximation to significantly longer times.

Shown in Fig. 5 are the results for $\u3008 x \u02c6 2 ( 0 ) x \u02c6 2 ( t ) \u3009$ at *β* = 8. Due to the non-linear operators, we did not apply RPMD or FK-CMD for this correlation function. As seen, the use of *f*_{1} = Ω_{t} is able to produce somewhat coherent oscillations which agree at least qualitatively with the exact results. However, the use of *f*_{2} in Eq. (76) results in a slightly higher frequency in the correlation function that is closer to the exact. Once again, in comparison to FK-LPI, it is evident that these new dynamics extend the accuracy of the CW approximation.

### B. The quartic potential

The next model system we look at is the quartic potential in which the potential energy takes the form

This model system presents a challenging test case due to the lack of a harmonic term. The results for $\u3008 x \u02c6 ( 0 ) x \u02c6 ( t ) \u3009$ at *β* = 8 are shown in Fig. 6. As seen, both FK-QCW(1) and FK-QCW(2) give essentially the same results for the linear correlation function. Furthermore, both of these methods are able to maintain coherent oscillations longer than RPMD, which in comparison show significant dephasing. Out of all of the approximate methods employed, FK-CMD seems to maintain coherent oscillations the longest for this correlation function. This again suggest that there is an advantage in working with the Kubo transformed correlation function, since FK-QCW(1) and FK-CMD are equivalent for this quantity. In comparison to FK-LPI, which is almost completely quenched after the first oscillation, it is evident that these new ensemble conserving dynamics extend the accuracy of the CW approximation to much longer times.

Presented in Fig. 7 is $\u3008 x \u02c6 2 ( 0 ) x \u02c6 2 ( t ) \u3009$ at *β* = 8. Once again, due to the limitations of RPMD and CMD for non-linear operators, neither of these methods were applied for this correlation function. One notices in Fig. 7 that the use of *f*_{1} = Ω_{t} in FK-QCW(1) and FK-QCW(2) results in a slight frequency shift as compared to the exact correlation function, while, similar to the results in Fig. 5, the use of *f*_{2} in Eq. (76) results in a slightly higher frequency which seems to corrects this. Furthermore, while FK-QCW(1), FK-QCW(2), and FK-LPI are all very accurate for short times, once again, the ensemble conserving dynamics used in FK-QCW(1) and FK-QCW(2) maintain coherent oscillations much longer than FK-LPI, which suffers from strong dephasing.

## VIII. CONCLUSIONS

We have developed two classes of quasi-classical dynamics that have been shown to conserve the initial quantum ensemble. When used within the CW approximation of the quantum time correlation function, both classes produce the exact classical and high temperature limits, and one set produces the exact harmonic limit (FK-QCW(1) with *f* = Ω_{t}). Although FK-QCW(1) and FK-QCW(2) were shown to fail at capturing the long time quantum coherence effects for the model potentials studied, overall, these dynamics maintain coherent oscillations much longer than the classical dynamics implemented within FK-LPI and therefore significantly extend the accuracy of the CW approximation. Furthermore, the fact that these new dynamics fail in this regard is not surprising since they lack the phase information necessary to capture these effects. However, this does not, in general, present a problem for realistic condensed phase systems since the coupling to the many degrees of freedom acts as a bath that quenches quantum coherence effects.^{1,23}

The fact that these new dynamics were shown to be comparable to both RPMD and CMD suggests that they provide a potentially more appealing algorithm than these methods, since one is not limited to correlation functions involving linear operators. In addition, one will not encounter artificial frequencies when using FK-QCW(1) or FK-QCW(2) to simulate absorption spectra, as compared to RPMD^{7} (although thermostatted RPMD has recently been shown to reduce this problem^{24}). Furthermore, the practical application of these two new methods to realistic condensed phase systems is not out of reach since they present only a minimum amount of additional computational load as compared to FK-LPI. This potential will be realized in a forthcoming paper^{26} in which we present the multidimensional generalization of these dynamics and apply FK-QCW(2) for the determination of the dynamic structure factor in low temperature liquid para-hydrogen and ortho-deuterium and compare the results with those presented in Ref. 19.

## Acknowledgments

P.J.R. acknowledges the support of this research by the U.S. National Science Foundation (No. CHE-1362381), with additional support provided by the R. A. Welch Foundation (No. F-0019). G.N. acknowledges support from the Swedish Research Council.

### APPENDIX A: THE MULTI-DIMENSIONAL CHANGE OF VARIABLES THEOREM

Without loss of generality, we consider integration in two dimensions. The multi-dimensional change of variables theorem states the equality of the two integrals,^{20}

where the absolute value of the Jacobian

has been introduced. The regions of integration, *A* and *B*, are related through the change-of-variables mapping (*x*(*t*, *v*), *y*(*t*, *v*)), so that when (*t*, *v*) goes through *B*, *A* is produced: *A* = {(*x*, *y*)|(*x*, *y*) = (*x*(*t*, *v*), *y*(*t*, *v*)), (*t*, *v*) ∈ *B*}. The theorem holds when the mapping (*x*(*t*, *v*), *y*(*t*, *v*)) is one-to-one so that $ J $ is non-zero. A well-known example is the mapping of polar coordinates to Cartesian. If we write (*t*, *v*) = (*r*, *θ*), we have *x*(*r*, *θ*) = *r*cos(*θ*), *y*(*r*, *θ*) = *r*sin(*θ*), and *J* = *r*. Thus, the replacement becomes

Equations (26) and (57) can both be derived from the above theorem. We only consider the derivation of Eq. (26) since the procedure is the same. So let us derive Eq. (26) from Eq. (22),

As will be evident, the proposed dynamics of (*x _{c}*(

*t*),

*p*(

_{c}*t*),

*q*,

_{t}*p*) are “reasonable,” meaning that trajectories do not cross. It thus follows that the dynamics evolve in a one-to-one manner. Then, the mapping defined by time-propagation

_{t} where $ x c \u2032 ( t ) = x c \u2032 ( x c \u2032 , p c \u2032 , q \u2032 , p \u2032 ; t ) $ and similarly for $ p c \u2032 ( t ) $, *q*′(*t*) and *p*′(*t*), may be used for changing variables in Eq. (A4). Thus, the replacement becomes

Substituting into Eq. (A4), we obtain

where

The last equality follows from the fact that the centroid variables $ ( x c \u2032 ( t ) , p c \u2032 ( t ) ) $ are taken to evolve independently of the $ ( q t \u2032 , p t \u2032 ) $. It is well-known that if *A*, *B*, and *C* are all *n* × *n* matrices, then we may deduce for the determinant of the larger 2*n* × 2*n* matrix,

Thus, we have

The first Jacobian, $ J ( x c \u2032 ( t ) , p c \u2032 ( t ) ) $, is unity since the centroid variables fulfill Liouville’s theorem, i.e., by following the centroid trajectories, the density of trajectories does not change in time.^{25} Notice that $ J ( q t \u2032 , p t \u2032 ) $ depends on the centroid coordinates since the phase-space points $ ( q t \u2032 , p t \u2032 ) $ will. Substituting Eq. (A11) into Eq. (A7), we obtain Eq. (26) after removing the prime on all integration variables.

### APPENDIX B: DETERMINANT OF THE JACOBIAN MATRIX IN TERMS OF THE DIMENSIONLESS COORDINATES

The dimensionless coordinates in Eq. (31) can be written as

*x*(

_{c}*t*),

*p*(

_{c}*t*)) is independent of (

*q*,

_{t}*p*) since the centroid variables evolve according to the classical like dynamics in Eqs. (24) and (25), then from Eq. (B1), we have that

_{t}such that $ J ( q t , p t ) $ in terms of these coordinates takes the form

Part of this expression can be explicitly integrated to give

which after using Eq. (B2) becomes

### APPENDIX C: EQUIVALENCE OF FK-QCW(1) AND FK-CMD FOR THE KUBO-TRANSFORMED TIME-CORRELATION FUNCTION OF LINEAR OPERATORS

The Kubo transformed time-correlation function

is related to the standard quantum time correlation function in Eq. (1) through

where *C _{AB}*(

*ω*) and $ C \u0303 A B ( \omega ) $ are the Fourier transforms of the standard and Kubo transformed time correlation functions, respectively. This quantity can be expressed in terms of the exact centroid distribution and QDO in which it takes the form

^{14}

where $ A c ( x c , p c ) \u2261Tr ( \delta \u02c6 ( x c , p c ) A \u02c6 ) $. This expression is exact so long as $ A \u02c6 $ is linear in position and/or momentum.^{14}

Using the Wigner representation for the trace of two general operators

we can rewrite Eq. (C3) as

Invoking the CW approximation, $ [ B \u02c6 ( t ) ] W ( q , p ) \u2192 [ B \u02c6 ] W ( q t , p t ) $, along with the FK approximation to the density operator, Eq. (C5) becomes

where a factor of 2*π*ħ has been absorbed into the definition of *ρ _{FK}*. Using Eq. (29), the above equation may be written as

If we apply the multi-dimensional change of variables theorem for the (*q*, *p*) variables, then the inner integral in Eq. (C7) may be rewritten, and the result is

By renaming the integration variables: (*q _{t}*,

*p*) → (

_{t}*q*′,

*p*′), Eq. (C8) can be equivalently expressed as

Furthermore, switching back from the Wigner representation to the trace, we have that

where we have defined $ B c ( x c ( t ) , p c ( t ) ) \u2261Tr ( \delta \u02c6 F K ( x c ( t ) , p c ( t ) ) \u2009 B \u02c6 ) $. Recognizing this expression as being equivalent to using the FK approximation of the density operator within the framework of CMD, we see that the entire class of quasi-classical dynamics in Sec. V reduces to CMD when used to evaluate the Kubo transformed correlation function of linear operators. Furthermore, the fact that one can alternatively arrive at the expression in Eq. (C10) by only invoking the CMD approximation^{4}

shows that the CW approximation is equivalent to making the CMD approximation for the Kubo transformed correlation function, so long as the dynamics retain the property in Eq. (29).