We introduce a scheme for approximating quantum time correlation functions numerically within the Feynman path integral formulation. Starting with the symmetrized version of the correlation function expressed as a discretized path integral, we introduce a change of integration variables often used in the derivation of trajectory-based semiclassical methods. In particular, we transform to sum and difference variables between forward and backward complex-time propagation paths. Once the transformation is performed, the potential energy is expanded in powers of the difference variables, which allows us to perform the integrals over these variables analytically. The manner in which this procedure is carried out results in an open-chain path integral (in the remaining sum variables) with a modified potential that is evaluated using imaginary-time path-integral sampling rather than requiring the generation of a large ensemble of trajectories. Consequently, any number of path integral sampling schemes can be employed to compute the remaining path integral, including Monte Carlo, path-integral molecular dynamics, or enhanced path-integral molecular dynamics. We believe that this approach constitutes a different perspective in semiclassical-type approximations to quantum time correlation functions. Importantly, we argue that our approximation can be systematically improved within a cumulant expansion formalism. We test this approximation on a set of one-dimensional problems that are commonly used to benchmark approximate quantum dynamical schemes. We show that the method is at least as accurate as the popular ring-polymer molecular dynamics technique and linearized semiclassical initial value representation for correlation functions of linear operators in most of these examples and improves the accuracy of correlation functions of nonlinear operators.

## I. INTRODUCTION

The numerical computation of equilibrium quantum time correlation functions remains one of the outstanding challenges in the theoretical treatment of nuclear quantum effects. This challenge constitutes a serious stumbling block to the accurate determination of both transport properties, which are determined from time correlation functions via the Green-Kubo theory, and vibrational spectra, which are fundamentally linked to integral transforms of time correlation functions. In principle, if the eigenvalues and eigenvectors of a quantum system could be determined, then the calculation of time correlation functions, which requires the time propagation of one of the operators in the correlation function, would be rendered straightforward. However, for condensed-phase systems, this is generally intractable. As an alternative to the direct calculation of eigenstates, the Feynman path integral formulations of quantum mechanics^{1–4} and quantum statistical mechanics are most often the approach of choice, and for static quantum properties, path integral techniques are both exact within the choice of the Hamiltonian and computationally efficient. Because static properties are derived from the trace of the Boltzmann operator exp(−*βĤ*), where *Ĥ* is the Hamiltonian and *β* = 1/*k*_{B}*T*, the path integral formulation is tantamount to sampling the configuration space of a set of ring-polymers composed of pseudoparticles or “beads” indexed by the discretization of an imaginary time parameter *τ* = *βℏ* determined by the temperature.^{5} Unfortunately, when it comes to quantum time correlation functions, the path integral formulation does not yield a computationally tractable problem, as the time propagators needed to evolve an operator in the Heisenberg picture, exp (±*iĤt*/*ℏ*), give rise to highly oscillatory phase factors in the path integral that are extremely difficult to converge.

Among the various approximations that are made in the calculation of quantum time correlation functions, the simplest method is to replace the quantum correlation function with its classical counterpart, which clearly constitutes a complete neglect of quantum effects. For vibrational spectra, quantum corrections can be incorporated via frequency-dependent prefactors based on an exact relation between classical and quantum time correlation functions for a harmonic oscillator. Not surprisingly, this approximation is satisfactory if the correlation function in question does not deviate strongly from either the harmonic or classical limits. Another class of highly popular techniques that exploit the connection between quantum and classical correlation functions are the centroid molecular dynamics^{6–9} (CMD) and ring-polymer molecular dynamics^{10} (RPMD) methods. These methods exploit the path integral molecular dynamics (PIMD) methodology,^{2,5,11–13} which is a widely used approach to evaluate imaginary time path integrals for static properties, ascribing a quantum dynamical character to the molecular dynamics evolution of the ring-polymer centroid or ring-polymer bead positions and using this dynamics to compute classical-like dynamical correlations of functions of these centroid or bead positions. These produce approximations to the Kubo-transformed quantum time correlation function. The primary difference between CMD and RPMD is the choice of fictitious dynamical masses of the beads and the manner in which thermostats are employed to produce Boltzmann sampling. However, because these are essentially imaginary time methods, they do not explicitly include any phase information. Consequently, systems with strong quantum coherence are treated rather poorly within the CMD and RPMD approximations. Furthermore, while various methods have been proposed,^{7,14–18} CMD and RPMD methods continue to exhibit difficulties in evaluating nonlinear correlation functions. Other approaches, including the thermal Gaussian molecular dynamics (TGMD) methods,^{19,20} which are conceptually similar to CMD, provide an alternative mapping of a quantum system to a classical system using an effective potential derived from the variational Gaussian wave packet approximation.^{21} While computationally efficient, TGMD provides slightly improved results when compared to RPMD and CMD for some correlation functions.

Beyond CMD and RPMD, it is also possible to employ semiclassical approaches that do treat phase factors albeit approximately via linearization techniques.^{22–36} These methods exist in a hierarchy, with some requiring generation of classical trajectories in very large numbers in order to converge the semi-classical expressions with an associated high computational overhead but possessing the ability to handle quantum coherence,^{22} and others requiring fewer trajectories but being more likely to break down when strong quantum coherence is present. The basic idea of most semiclassical trajectory-based approaches is the exploitation of a transformation to sum and difference variables between forward and backward time propagation paths, and the linearization is performed with respect to the difference variables. This type of transformation can be applied to the various formulations of the time correlation function including the standard^{24–26,28–30,33–35} and symmetrized versions^{31–33} and can even be carried out to second order^{31} or be used to create trajectory-based Monte Carlo algorithms.^{32} Finally, in certain cases, it has been shown that the phase factors can be computed via quadrature techniques using points sampled from the corresponding positive-definite quantum Boltzmann distribution.^{37–40} The putative scaling of such an approach, however, would be expected to grow exponentially with system size and is, therefore, difficult to apply to more than a small number of quantum degrees of freedom, although Makri and co-workers^{37–40} have shown that optimal Monte Carlo sampling of the grid points improves the scaling considerably.

In this paper, we consider a somewhat different perspective on the calculation of quantum time correlation functions, which, though inspired by the aforementioned semiclassical ideas, nevertheless, aims at a purely imaginary-time sampling scheme capable of capturing approximate phase information explicitly. The starting point is the symmetrized time correlation function alluded to above, which uses complex-time forward and backward propagators to evolve the quantum operators and which can be written as a path integral that possesses a positive-definite density and an oscillatory phase factor, both of which depend on the same path integral bead variables.^{41,42} The bead variables in this formulation of the symmetrized time correlation function are then transformed to sum and difference variables along the discretized complex-time paths, which, as previously noted, is not new.^{24–26,28–35} However, the manner and context in which it is applied here are. Once the transformation is made, the potential is then expanded as a power series in the difference variables, and the resulting integrals over these variables are then performed analytically. Following this procedure in the present context gives rise to a purely imaginary time path integral in the remaining sum variables that are subject to a modified potential. In this way, we calculate the approximate quantum time correlation function as a standard path-integral sampling problem rather than a costly trajectory-sampling problem. Nevertheless, phase information is included up to second order in the path difference variables. The resulting imaginary time path integral can be evaluated by any numerical scheme commonly used for path integrals, including Monte Carlo,^{43} molecular dynamics,^{13} hybrid Monte Carlo,^{13,44} or any path-integral enhanced sampling scheme.^{45,46} Moreover, the approximation can be systematically improved via a cumulant expansion carried out to any desired order. Within such a cumulant expansion, the integrals over the path difference variables can be performed analytically at each order, and the remaining path integral can be shown to be real and positive definite at any order and can, in principle, be evaluated using the same imaginary-time path integral sampling schemes mentioned above. Even at second order, the resulting correlation functions are found, in various cases, to be of comparable accuracy to those produced by the methods discussed above for linear operators and are generally of improved accuracy for nonlinear operators. Thus, it is possible to compute approximate quantum time correlation functions in which both the approximation made and a path to systematic improvement of the approximation are manifestly clear.

This paper is organized as follows: In Sec. II, we lay out the theory behind our new approach. Following this, a number of simple numerical examples are presented. Computational details of these examples are presented in Sec. III, and the results are shown in Sec. IV. Conclusions, further discussion, and prospects for future developments are given in Sec. V. Finally, in Appendixes A and B, we derive the first cumulant correction at third order in the path difference variables and the Suzuki-Chin^{47,48} formulation of our scheme.

## II. BACKGROUND AND THEORY

The theoretical development begins with the symmetrized quantum time correlation function *G*_{AB}(*t*), which is defined as

In this representation, *Ĥ* is the Hamiltonian, *Z* is the canonical partition function, *Â* and $B^$ are the operators of interest, and *τ*_{c} is a complex time variable such that *τ*_{c} = *t* − *iβℏ*/2. Thus, in the symmetrized time correlation function framework, the operator $B^$ is evolved in the complex time variable *τ*_{c} that mixes temperature and time. The symmetrized correlation function *G*_{AB}(*t*) and the true time correlation function, *C*_{AB}(*t*), are related through their Fourier transforms via the relation

In the analysis to be presented here, we will consider a single particle in one dimension with Hamiltonian $\u0124=p^2/2\u2009m+U(x^)$. Assuming *Â* and $B^$ are functions only of the position operator $x^$, we proceed by evaluating the trace in the coordinate basis, which gives

Unlike *C*_{AB}(*t*), *G*_{AB}(*t*) is a real object since the two matrix elements are complex conjugates of each other, and *a*(*x*) and *b*(*x*′), which are the eigenvalues of *Â* and $B^$ at $x$ and $x\u2032$, respectively, are real. The integrand in Eq. (3) can be interpreted as follows: start at the point *x* and measure the operator *Â* at that point; the point *x* is then propagated forward in complex time *τ*_{c} to *x*′; at *x*′, the operator $B^$ is measured; finally, the point *x*′ is propagated backward in complex time back to *x*. Using the Trotter factorization method, Krilov and Berne^{41} demonstrated that *G*_{AB}(*t*) can be written in a path integral formulation by expressing each matrix element as a discrete path integral with *P* beads. Thus, the path-integral representation of *G*_{AB}(*t*) contains 2*P* beads and takes the form

where *ρ*(*x*_{1}…*x*_{2P}) is a positive-definite distribution given by

and the phase factor *ϕ*(*x*_{1}…*x*_{2P}) is defined to be

The condition *x*_{2P+1} = *x*_{1} arises from the trace condition. Consequently, the path integral we obtain involves a ring polymer consisting of 2*P* beads. As Eqs. (4)–(6) show, both the positive definite distribution and the phase factor depend on $\beta $ and $t$; they also depend on the same set of path-integral variables *x*_{1}, …, *x*_{2P}, which has some advantages for numerical convergence over the standard correlation function *C*_{AB}(*t*), as it allows the phase factor to be treated as a weight in the configurational average. However, the rapid fluctuations of the phase factor caused by small changes in the complex-time paths still cause severe convergence issues.

Since the positive-definite distribution and phase factor in the symmetrized correlation function involve the same set of path-integral bead variables, a useful way to tame the phase factor is to introduce a change of the integration variables *x*_{1}, …, *x*_{2P} to new variables that represent an average of the forward and backward complex-time paths and the difference between these two paths. One possibility for such a transformation is the variable set *r*_{i} = (*x*_{i} + *x*_{P+i})/2, *s*_{i} = *x*_{i} − *x*_{P+i}, where *i* = 1, …, *P*. In this transformation, each point along the forward path is combined with the point along the reverse path having the same imaginary time index. This transformation has a few appealing properties, in particular, since *x*_{1} = *x*_{2P+1}, it follows that *r*_{P+1} = *r*_{1} and *s*_{P+1} = −*s*_{1}. Thus, if we expand the potential appearing in *ρ* and *ϕ* to second order in *s*_{1}, …, *s*_{P} and then integrate analytically over *s*_{1}, …, *s*_{P}, the remaining path integral over *r*_{1}, …, *r*_{P} will involve a ring polymer of *P* beads in an effective potential with contributions arising from the analytical integration over the path difference variables. However, as the impetus for expanding the potential in powers of *s*_{i} is that these variables are small, the unfortunate fact remains that the points *x*_{i} and *x*_{P+i} are separated by *P* intervening beads, which gives them a maximal separation in the original 2*P*-bead ring polymer so that the assumption under which we expand the potential in powers of *s*_{i} is not valid. Consider, however, the following transformation:^{24–26,28–35}

The variables *r*_{i} and *s*_{i} are also path average and difference variables, respectively; however, in this transformation, early points along the forward path are matched with corresponding *late* points along the reverse path and vice versa, as illustrated in Fig. 1. This transformation has the following inverse:

In many condensed-phase problems, there is significant decoherence, suggesting that the path difference variables *s*_{i} can be expected to be small. Thus, the transformation in Eq. (7), which makes these variables explicit, allows the resulting path integral to be expanded in powers of *s*_{i}, an expansion that would be expected to converge quickly as these powers increase. For practical implementation, if the expansion is truncated at second order, then the resulting integrals over *s*_{i} can be performed analytically. When this is done, the resulting path integral over *r*_{i} will involve an *open* or worm-like chain. Under this transformation, the harmonic bead-bead coupling terms in the distribution function *ρ* can be easily worked out, and we find the following:

Similarly, in the phase factor, we obtain

For the potential energy, we perform a power series expansion about *s*_{i} = 0 on the two types of terms that result from the transformation,

A similar second-order expansion of the potential around the difference variables has been implemented in the local harmonic approximation/local Gaussian approximation (LHA/LGA) for the linearized semiclassical initial value representation (LSC-IVR) methods.^{49,50} The justification for the second-order expansion of the potential is to evaluate appropriate integrals analytically as shown below. Thus, in the positive-definite distribution, the potential energy becomes

and in the phase factor,

Putting the potential and harmonic coupling expressions together, *ρ* transforms to

and the phase factor becomes

Note that the distribution can be separated into a purely *r*_{i}-dependent term and a term that depends on both *r*_{i} and *s*_{i} according to

where

and

We now define a square tridiagonal (*P* − 1) × (*P* − 1) matrix *M* as

where *r* ≡ *r*_{2}, *…*, *r*_{P}, and

and a (*P* − 1)-dimensional vector *S* = (*s*_{2}, …, *s*_{P}). In terms of these, the distribution *ρ*_{2} can be rewritten as

The phase factor can be rewritten in a similar manner by introducing a (*P* − 1)-dimensional vector *K* defined as

where *r*′ ≡ *r*_{1}, …, *r*_{P+1}, and

In terms of *K* and *S*, the phase becomes

Thus, the overall expression for *G*_{AB}(*t*) within our approximation becomes

Within this approximation, the integration over *S* is simply a Gaussian integral which can be performed analytically, and while the second-order expansion and resulting integration are not new, the motivation of the expansion in this derivation allows for the elimination of the underlying imaginary terms in the integrals under the assumption that the eigenvalues of *M* are all positive definite. This last condition should be obeyed if *P* is sufficiently large, which contrasts with the conditions obtained in trajectory-based methods.^{31} The result of the integration is

and the symmetrized correlation function that results involves an imaginary-time path integral over open chains in the variables *r*_{1}, …, *r*_{P+1} is

Eq. (27) indicates that the open chain is subject to an effective or modified potential determined not only by the original potential *U* but also by the matrix *M*(*r*) and the vector *K*(*r*′). However, as Eq. (27) is also just an imaginary-time path integral, any of the usual sampling schemes typically applied for path integrals can be employed, including Monte Carlo,^{43} molecular dynamics,^{13} hybrid Monte Carlo,^{13,44} or any path-integral enhanced sampling scheme.^{45,46} Here, although it might not be optimal, for simplicity, we have employed a standard PIMD scheme. We thus introduce fictitious momenta and corresponding fictitious masses and write Eq. (27) as

where *p* ≡ *p*_{1}, …, *p*_{P+1}, and where the bead Hamiltonian is given by

Every term in this Hamiltonian, with the exception of the last term and $\omega p=P/|\tau c|$, is identical to standard open PIMD.^{51–53} The last term is an additional potential resulting from the Gaussian integration and is defined as

Note that, unlike *U*, the potential *Ṽ*(*r*_{1}, …, *r*_{P+1}) couples all of the beads, i.e., it cannot be written as a sum over terms that are functions of only one bead variable. It also depends on both $\beta $ and $t$.

At this point, a few comments are in order. First, since the approach is based on an open-chain path integral formulation of the symmetrized correlation function, we assign the acronym OPSCF to the method, which stands for Open Path Symmetrized Correlation Function. Second, although we have truncated the expansion in the path difference variables at second order, the scheme can be *systematically* improved by going to higher orders in *s*_{i} and employing a cumulant expansion. In Appendix A, we derive an example, computing to the first cumulant coming from an expansion up to third order in *s*_{i}. Importantly, the expansion and subsequent analytical integration that lead to Eq. (30) involve an expansion of the potential only about *s*_{i} = 0 rather than a global expansion of the potential as in Ref. 25. Thus, the full potential dependence on *r*_{i} is retained along with the extra terms in the effective potential appearing in Eq. (30). The quality of this effective potential determines the accuracy of the approximate correlation functions produced by OPSCF. Third, and following on the previous point, it is possible to improve the accuracy with the choice of *P* by, for example, designing colored-noise thermostatting approaches^{54,55} or by using a higher-order Trotter expansion, for example, the Suzuki-Chin approximation.^{47,48} Note that these schemes cannot be used to accelerate RPMD or CMD calculations, indicating a possible advantage of sampling schemes over trajectory-based methods. In Appendix B, we derive the corresponding Suzuki-Chin version of Eq. (27). The result, which is manifestly real, is the appearance of new, additional terms in the effective potential in Eq. (30). Fourth, the operators in Eq. (27) need not be linear functions of position. Since the transformation does not affect the *a*(*x*_{1}) = *a*(*r*_{1}) and *b*(*x*_{P+1}) = *b*(*r*_{P+1}) terms in the integral, the approximation should be of similar quality for both linear and nonlinear operators. Fifth, as previously mentioned, the path integral in Eq. (27) can be evaluated by any type of sampling algorithms, including those designed for Monte Carlo calculations.^{56} Sixth, the value of *t* is treated as an initial parameter in the simulation, and thus each value of *t* requires a separate simulation regardless of the sampling scheme. Finally, for general potential energy functions, the canonical ensemble is not preserved within the OPSCF scheme, and we will discuss this issue shortly.

The OPSCF method approximates a 2*P* ring polymer by a *P* + 1 open chain in an effective potential. This should not be surprising as the symmetrized correlation function, in contrast to the standard correlation function *C*_{AB}(*t*), is constructed such that the same bead variables define both the positive-definite distribution and the phase function. Consequently, any approximation made to tame the phase factor must be made to the positive-definite (canonical) distribution as well. Thus, the validity of the method relies on the accuracy of the second-order expansion of the potential made in both contributions so the correlation function. Note that the canonical ensemble distribution is generated only at *t* = 0 in the symmetrized correlation function, and since the phase factor vanishes at *t* = 0, it is always possible to treat *t* = 0 as a special case for which an exact 2*P*-bead cyclic path integral calculation is performed without making the OPSCF transformation and subsequent power-series approximation. If this is done, then the canonical ensemble will be preserved at *t* = 0. In the applications to be presented below, however, we will not employ this trick, as it is important to assess the degree to which OPSCF approximates the true canonical distribution at *t* = 0. The OPSCF method does preserve detailed balance, however, which is also true for CMD and RPMD but is not true for some LSC-IVR implementations. The only changes to the derivation needed to prove preservation of detailed balance are the sign of the imaginary term in Eq. (25) and an exchange of the bead indices at which operators are evaluated, i.e., [*b*(*r*_{1}) and *a*(*r*_{P+1})]. Clearly, the change of sign does not affect the result of the analytical integration of Eq. (25) over *s*_{2}, …, *s*_{P}. Moreover, the symmetry of the open chain allows the evaluation of the operators at either endpoint as long as both endpoints sample the entire phase space. It must be noted, however, that the argument made above assumes the validity of the positive-definite distribution *ρ*(*r*_{1}, …, *r*_{P+1}) of the OPSCF approximation. However, as this distribution is only exact for the harmonic oscillator, preservation of detailed balance in anharmonic systems will only be valid with respect to this (approximate) distribution.

All of these characteristics need to be contrasted to the CMD and RPMD approaches, in which it is not possible to perform the calculations by any means other than molecular dynamics, higher-order factorizations cannot be used, and which degrade when applied to nonlinear operators. The fact that the OPSCF scheme presented here can be improved in a systematic and practically computable way also stands in contrast to CMD and RPMD for which it is not straightforward to see how to obtain systematically better results than what the basic schemes produce. Furthermore, the OPSCF transformation and subsequent power-series expansion are not explicitly limited to the symmetrized correlation function but could also be applied to the Kubo-transformed correlation function. The OPSCF method should also be contrasted with trajectory-based schemes, which also require the generation of dynamical trajectories and carry a computational overhead when attempting to accurately capture quantum coherence. Finally, we point out that the OPSCF scheme is exact for harmonic potentials, for which the expansion in *s*_{i} is exact, and converges to the correct classical limit, where the difference between the forward and backward path vanishes. The obvious downside is the need to compute the second derivative of the potential.

In Secs. III and IV, we will examine a number of simple one-dimensional examples that are commonly used to study the performance of the CMD, RPMD, and semiclassical approaches. We first describe the computational details of our simulations followed by a presentation of the results for each of the examples.

## III. COMPUTATIONAL DETAILS

Two correlation functions on which we will focus in this work are the position autocorrelation function $(\xc2(x^)=B^(x^)=x^)$ and the position-squared autocorrelation function $(\xc2(x^)=B^(x^)=x^2)$. We explored four different potentials using the Open Path Symmetrized Correlation Function (OPSCF) method. These are a harmonic oscillator

a mildly anharmonic potential

a quartic potential

and a symmetric double well potential

These various potentials have been extensively used to benchmark approximation schemes for computing quantum time correlation functions.^{7,10,17,18,32,33,57–61} In all examples, we use natural units in which *ℏ* = 1. Two different temperatures were explored for each potential. Specifically, *β* = 1 and *β* = 8, except for the symmetric double well, in which the low temperature was *β* = 4. The position-squared correlation function, $Gx2x2(t)$, was excluded from the main body of the text for the mildly anharmonic potential and quartic potential. We found that $Gx2x2(t)$ for these potentials at the lowest temperatures considered here were so small as to be within the sampling error that no useful information concerning the performance of the method could be discerned.^{62,63} At the higher temperatures, however, we were able to obtain quality $Gx2x2(t)$ curves that we could benchmark against exact results. Section IV will discuss the particular issues encountered in generating $Gx2x2(t)$ at low temperatures, specifically for the harmonic oscillator and double well potentials.

We implemented the stochastic isokinetic thermostat^{64} scheme solely to ensure proper canonical sampling. The OPSCF simulations were performed using staging variables^{13} for open chains.^{53} The time step of the simulations ranged from 0.005 to 0.05 natural units of time depending on the temperature, the value of *t*, and the bead number *P*. We chose the value of *P* such that |*τ*_{c}|/*P* remains approximately constant as *t* is increased, and as *P* is increased, the time step must typically be decreased, or multiple time-step methods can be employed.^{16,81} Figure S1 of the supplementary material demonstrates convergence of the correlation functions with respect to the bead number. Simulations were run for a minimum of 10^{7} steps for *β* = 1 and 2 × 10^{7} steps for *β* = 4, 8 until proper convergence is reached. Standard errors of the data points of each simulation are depicted as error bars in various figures. Finally, the Suzuki-Chin factorization for open path-integral molecular dynamics (OPIMD)^{47,48,65,66} (see Appendix B) was implemented for the mildly anharmonic and quartic potentials (Figs. S2–S5 of the supplementary material) in order to determine if a higher-order factorization scheme improved the efficiency by allowing smaller values of *P*.

We benchmarked the results of the OPSCF method against RPMD correlation functions as well as correlation functions obtained from the linearized semiclassical initial value representation (LSC-IVR). For the RPMD simulations, 100 000 configurations were accumulated from a long thermostatted PIMD simulation of 10^{8} steps for the various potentials and temperatures. From these configurations, we sampled the initial velocities from the Maxwell-Boltzmann distribution. Trotter numbers of 8 and 32 were used for the high and low temperatures, respectively. From each of these initial phase-space points, we performed the RPMD simulations in primitive coordinates with a time step of 0.0005 natural units of time to converge the dynamics fully, following the protocol of Ref. 57. Each RPMD simulation was run for 80 000 steps, which is the equivalent to 40 natural units of time. Since RPMD produces the Kubo transformed correlation functions, *K*_{AB}(*t*), the symmetrized correlation function for RPMD was computed by performing a Fourier transform of *K*_{AB}(*t*) and generating the Fourier transform of *G*_{AB}(*t*) via

followed, finally, by an inverse Fourier transform to give the RPMD approximation to *G*_{AB}(*t*). Figures S6–S9 in the supplementary material provide the Kubo transformed correlation functions directly from the RPMD simulations prior to their transformation to *G*_{AB}(*t*).

For the LSC-IVR calculations, the thermal Gaussian approximation (TGA-LSC-IVR) was used following the scheme devised by Liu and Miller.^{27} Specifically, 41 imaginary time trajectories were generated for the potentials at the desired temperatures. For these 41 trajectories, the differential equations were solved using a fourth-order Runge-Kutta algorithm. Subsequently, approximately 10^{7} real time trajectories were sampled from each imaginary time trajectory and the dynamics were evolved using the velocity Verlet integrator with a time step of 0.02. TGA-LSC-IVR generates the symmetrized autocorrelation function, *G*_{AA}(*t*), directly. While TGA-LSC-IVR might not be the most robust semiclassical method, it has been shown to perform well in condensed-phase systems^{28} and is an appropriate semiclassical method to compare to OPSCF. Exact quantum time correlation functions were calculated by solving the eigenvalues and eigenvectors of each system using a harmonic basis set, which was performed in Mathematica.^{67} In these calculations, 100 basis functions were sufficient to converge all exact quantum time correlation functions.

## IV. RESULTS AND DISCUSSION

### A. Harmonic oscillator

The OPSCF method is exact for harmonic systems, as is also the case for RPMD, CMD, and TGA-LSC-IVR. Therefore, simulations of the harmonic oscillator provide a benchmark for understanding the strengths and weaknesses associated with numerical evaluation of Eq. (27). Figure 2(a) depicts the position autocorrelation function *G*_{xx}(*t*) at *β* = 1. RPMD, TGA-LSC-IVR, and OPSCF reproduce the exact result as expected. Figure S1 of the supplementary material shows the convergence of the OPSCF correlation functions as a function of the number of beads. In Sec. III, we discussed the convergence criterion for OPSCF calculations, specifically, that the value |*τ*_{c}|/*P* should be kept approximately constant. In this way, the proposed criterion resembles the convergence criteria often used in imaginary-time path integral simulations (i.e., *β*/*P*) remain constant as *β* is varied. In OPSCF, however, |*τ*_{c}| depends on both time and temperature. Thus, each of the data points in the correlation function plots to be presented in this section corresponds to the minimum value of *P* required to reach convergence.

Figure 2(b) demonstrates the increased difficulty of obtaining correlation functions for non-linear operators over linear operators. This problem has challenged schemes such as CMD and RPMD, and various approaches have been introduced to address it.^{7,14–18} Specifically, Fig. 2(b) shows the position-squared autocorrelation function $Gx2x2(t)$ at *β* = 1. RPMD misses some of the extrema, particularly the minima, of this correlation function. The OPSCF method and TGA-LSC-IVR, by contrast, are seen to reproduce the exact result within statistical error bars. Equation (27) is agnostic to the functional form of operators *Â* and $B^$ as long as the operators are solely functions of positions. The accuracy of $Gx2x2(t)$ suggests that the OPSCF method should reproduce any linear and nonlinear correlation function for the harmonic oscillator as long as the configuration space is sufficiently sampled. This highlights the strength of the OPSCF approximation in representing correlation functions for nonlinear operators.

Reducing the temperature of the harmonic oscillator ensures that higher energy states are less accessible, thus creating a ground-state dominance, thereby increasing the quantum character of the system. In principle, since the OPSCF method is exact for the harmonic oscillator, we should be able to reproduce the exact correlation functions for both linear and nonlinear operators. *G*_{xx}(*t*) [Fig. 3(a)] is reproduced very well with OPSCF, RPMD, and TGA-LSC-IVR although a larger value of *P* is needed for the OPSCF and RPMD compared to *β* = 1. Some of the data points lie just outside the exact curve but are within statistical accuracy denoted by the error bars.

Finally, we discuss the efficiency of the OPSCF method in reproducing $Gx2x2(t)$ for the harmonic oscillator at *β* = 8 [Fig. 3(b)]. As stated, the OPSCF method is exact in the harmonic oscillator and should reproduce any position-dependent correlation function. However, Fig. 3(b) demonstrates that while many points are within statistical error of the exact correlation function, the data points do not produce the correct periodicity compared to the exact curve. Given that the oscillations for the $Gx2x2exact(t)$ range within an order of 10^{−4}, the level of sampling is not able to resolve such small differences. Note that while the periodicity of RPMD is quite close to the exact result, the curve is shifted considerably from the exact curve. Interestingly, however, this numerical resolution issue is less severe for the Kubo-transformed correlation function, $Kx2x2(t)$ (Fig. S6 of the supplementary material), which RPMD approximates directly. As Fig. S6(d) of the supplementary material shows, the amplitude of the oscillations is larger than that of $Gx2x2(t)$. Furthermore, the results from TGA-LSC-IVR are excluded, as the evaluation of the symmetrized correlation function is numerically unfavorable for this method. However, it has been shown that LSC-IVR is able to recover the correlation function $Cx2x2(t)$ exactly for the harmonic oscillator at low temperatures.^{62,63} This issue regarding $Gx2x2(t)$ at low temperatures is present in both the mildly anharmonic potential and the quartic potentials. Therefore, we have excluded the data for $Gx2x2(t)$ at *β* = 8 from this present work owing to the difficulty of converging the curve to a sufficiently small numerical error that a meaningful comparison could be made. We will, however, continue to present $Gx2x2(t)$ for *β* = 1.

### B. Mildly anharmonic potential

The mildly anharmonic potential, Eq. (32), and its associated time correlation functions have been extensively studied as a nontrivial example to test approximation schemes.^{10,17,18,32,57–59} The slight anharmonicity in the potential produces small deviations from the simple harmonic oscillator. We can, therefore, test the validity of the second-order truncation of the potential energy within the OPSCF method using this anharmonic potential. Once again, at *β* = 1, Fig. 4(a), all three methods reproduce the correlation function, *G*_{xx}(*t*). All methods exhibit a very slight shift in frequency at longer times. When comparing to the harmonic oscillator, the OPSCF requires a larger value of *P* to accurately capture the autocorrelation function.

In Fig. 4(b), we see that the OPSCF does remarkably well for the nonlinear case of $Gx2x2(t)$ compared to RPMD, and OPSCF compares well against TGA-LSC-IVR. In particular, RPMD accurately captures some of the extrema of the curve but fails at the minima, while OPSCF outperforms RPMD at many of the extrema. So far, the OPSCF method has proven to be an effective tool at the higher temperature *β* = 1 even for $Gx2x2(t)$. The Suzuki-Chin factorization implemented for this potential at *β* = 1 [see Fig. S2(b) of the supplementary material] improves the performance of $Gx2x2(t)$ compared to the traditional Trotter factorization. Specifically, smaller values of *P* are needed to converge some of the individual data points.

At *β* = 8, we begin to see the limitations of all three methods. Interestingly, Fig. 5 demonstrates that OPSCF produces the correct frequency in the curve but with a shift in the amplitude, while for RPMD and TGA-LSC-IVR, there is a noticeable shift in the frequency, but the amplitude is more accurate. Furthermore, one may notice that the RPMD curve does not begin at the correct value of *G*_{xx}(0). This is a result of the Fourier transform described in Eq. (35). Figure S7(b) of the supplementary material demonstrates the accuracy of the RPMD simulations compared to the exact Kubo-transformed correlation functions, *K*_{xx}(*t*), for the mildly anharmonic potential. Even TGA-LSC-IVR does not have the correct *t* = 0 limit in this case. However, the semiclassical method is able to recover the correct *t* = 0 limit for anharmonic potential using different correlation functions,^{62,63} thus demonstrating the difficulty in converging the symmetrized correlation function at low temperatures. The OPSCF method is able to capture the correct *t* = 0 limit at this temperature, but it begins to degrade soon after. Even a slight anharmonicity in the potential will be affected by the truncation of the potential to second order in Eq. (11). In Appendix A, we propose an expansion of the potential up to third order in *s*_{i} allowing us to calculate the first cumulant. In theory, this would improve the results, specifically for *t* > 0. As stated in Appendix A, this cumulant will be incorporated into the time-dependent phase factor, and thus *t* = 0 would not be affected by the first cumulant. This is ideal for the mildly anharmonic potential since the short time limit is accurately obtained for the OPSCF method. The Suzuki-Chin implementation of the OPSCF method for *β* = 8, Fig. S3 of the supplementary material, provides moderate improvement to *G*_{xx}(*t*), especially within the first few data points.

### C. Quartic potential

Generating accurate time correlation functions for the quartic potential, Eq. (33), has been difficult for most approximation schemes.^{6,10,32,33,58–60} None of the methods presented replicate the long time behavior of the correlation function. In fact, RPMD and OPSCF methods are intriguingly similar, while TGA-LSC-IVR falls short. *G*_{xx}(*t*) at *β* = 1, Fig. 6(a), for RPMD and OPSCF decay soon after the first oscillation. However, neither extrema of the first oscillation are perfectly captured. Most current methods fail to reproduce multiple oscillations. At a slightly lower temperature, *β* = 2, Matsubara dynamics^{68,69} are able to reproduce more oscillations compared to RPMD and CMD for the quartic potential. While Matsubara dynamics are computationally expensive, its basis in semiclassical theory may provide integral information in improving the OPSCF method to produce multiple oscillations for the quartic potential.

The OPSCF method outperforms both RPMD and TGA-LSC-IVR for $Gx2x2(t)$ (Fig. 6) by locating the correct extrema in the first oscillation. This has not been achieved either with RPMD or with CMD nor have improvements to CMD, such as the use of CMD in conjunction with maximum entropy analytical continuation^{16} (CMD-MEAC) succeeded in this regard. Thus far, one of the strengths of OPSCF is the accuracy of generating nonlinear quantum correlation functions for *β* = 1. Capturing the early extrema in $Gx2x2(t)$ is a positive outcome for the OPSCF method even if the OPSCF method cannot extend past the first oscillation; the short time behavior of the correlation function is well reproduced.

The low-temperature regime still proves to be a challenge for the OPSCF method as currently implemented (Fig. 7). The low temperature of *β* = 8 demonstrates the first instance in this paper when $GxxOPSCF(0)\u2260Gxxexact(0)$. *t* = 0 is a special case within the OPSCF scheme when the phase factor in Eq. (4) disappears resulting in a 2*P* imaginary-time path integral. Thus, by performing the transformation of Eq. (7), we are effectively approximating a closed-chain imaginary time path integral with an open-chain imaginary time path integral containing a modified form of the additional potential of Eq. (30), $V\u0303=12\beta ln(det[M(r)])$, which includes the approximate second-order expansion of the potential energy around *s*_{i}. In other words, the *t* = 0 statistics are approximated as they are at other times, and this stands in contrast to other semiclassical approaches that correctly reproduce the equilibrium statistics. This is a consequence of the use of the symmetrized correlation function, as the same bead variables appear in both the positive-definite density and the phase factor, and hence, approximation made to the latter also affects the former. The data at *t* = 0 confirm that the truncation is insufficient for the quartic potential at low temperatures, while at *β* = 1, quantum effects are less pronounced and we find the truncation to be sufficient in the short time limit. Note, however, that the *t* = 0 value is still a slight improvement over RPMD, which also fails to produce the correct *t* = 0 value of the correlation function. TGA-LSC-IVR interestingly correctly captures the *t* = 0 limit even though it failed to do so for the anharmonic potential. In a similar manner to the mildly anharmonic potential, a cumulant expansion provides a systematic improvement to the method. However, it is clear that a higher-order cumulant would be necessary (not shown in Appendix A) in order to bring in a correction term at *t* = 0. Moreover, as previously noted, $t=0$ can be treated as a special case, and the exact $2P$-path integral could be performed without approximation in order to obtain the correct value at this point.

For *t* > 0, the OPSCF method also compares poorly to the exact correlation function, *G*_{xx}(*t*), in the low temperature regime. The data points fail to describe both the amplitude and frequency of the correlation function. This system would need to be evaluated more extensively to determine if and how a systematic improvement could be developed. One such approach could be the incorporation of the first- and second-order cumulant expansions. As stated above, at *t* = 0, the second-order expansion of the potential energy is insufficient in describing the quartic potential at low temperatures, and this would hold true for *t* > 0. Even though, TGA-LSC-IVR correctly obtains the short time limit, this method decays too rapidly compared to RPMD. However, a more robust semiclassical method would presumably perform better at longer times. RPMD does a much better job in the low temperature regime over a long time compared to *β* = 1. It has been argued by Voth and co-workers^{6,58,60} that the low temperature regime of the quartic oscillator creates a two-state system that effectively evolves in a quasi-harmonic environment. Both RPMD and CMD are exact in the harmonic limit and thus will perform reasonably well in this regime. However, the dynamics of the higher frequency modes in the RPMD scheme cause the dampening effect of the correlation function further out in time. Although the OPSCF method is also exact in the harmonic limit, it appears to be more sensitive to the deviations produced by the quartic potential. Unlike the mildly anharmonic potential, the Suzuki-Chin factorization (Figs. S4 and S5 of the supplementary material) does not appear to affect the convergence of the correlation function with respect to *P*.

### D. Symmetric double well

The symmetric double well potential, Eq. (34), is the final potential used to test the OPSCF method and has been explored for a variety of the methods currently discussed.^{58,61–63} Not only does this potential deviate greatly from the harmonic oscillator but also it has been studied in both enhanced sampling methods^{45} and quantum transition state theory,^{70} specifically in understanding quantum tunneling.

Figure 8(a) shows the position correlation function in the short time regime for *β* = 1. All methods perform similarly for *G*_{xx}(*t*). However, OPSCF improves when examining $Gx2x2(t)$ [Fig. 8(b)] confirming the OPSCF method’s ability to determine accurately the nonlinear correlation function for the higher temperature (*β* = 1) compared to RPMD and TGA-LSC-IVR. While the higher temperature regime favors OPSCF, the method continues to struggle with the lower temperatures. Figure 9 illustrates the limitations of the method at *β* = 4. An intermediate temperature was used instead of *β* = 8 to ensure that sufficient sampling over the barrier could be achieved. The OPSCF method reproduces the correlation function *G*_{xx}(*t*) at short time only [Fig. 9(a)], while $Gx2x2(t)$ [Fig. 9(b)] is rather poor at all times. By contrast, RPMD captures the short time limit of both correlation functions rather well. TGA-LSC-IVR is excluded for this temperature since convergence of the symmetrized correlation function was poor, meaning that we could not run enough imaginary-time or real-time trajectories to yield converged results. However, TGA-LSC-IVR has performed well in the short time limit for the double well potential at low temperatures for the Kubo and standard correlation functions.^{62,63} Systematic improvements to the OPSCF method, specifically first-and second-order cumulants, could lead to more accurate results.

Another issue arises with the OPSCF method in regards to the symmetric double well potential. Specifically, the OPSCF simulations for the double well potential begin to fail for times *t* > 3. A preliminary analysis of the forces produced by the simulations revealed that the forces spike to large values (many orders of magnitude) throughout the simulation. This alters the dynamics and the simulations no longer follow the equations of motions defined by the Hamiltonian of Eq. (29). We are currently investigating the nature of this behavior and if either an alternative sampling scheme or other improvements can be made to alleviate this problem. Monte Carlo would provide an alternative sampling mechanism if the evaluation of forces proves to be the cause of the failed simulations.

At this point, a comment is in order on how one might use OPSCF to compute a quantum rate since rates require the evaluation of a different nonlinear correlation function. The alternative methods of RPMD,^{71,72} CMD,^{73} and LSC-IVR^{50,74} have the ability to calculate rate constants by utilizing the flux-side correlation function. The symmetrized version of this correlation function is

where $F^=12\u2009mp^\delta (x^(0)\u2212x\u2021)+\delta (x^(0)\u2212x\u2021)p^$ is the flux operator and $\u0125=h(x^(t)\u2212x\u2021)$ is the Heaviside function, and *x*^{‡} is the location of the transition state (dividing surface) such that the reaction rate is

where *Q*_{r}(*T*) is the reactant partition function. In a similar fashion to Eq. (3), the symmetrized version of the flux-side correlation function can be written as

The incorporation of a momentum-dependent operator (flux operator) yields a few more terms in Eq. (38) compared to Eq. (3). Path integral representations of symmetrized correlation functions of momentum operators do not result in a closed ring polymer [see Eq. (3) and Fig. 1 for position operators], but instead an open chain is formed. Therefore, the transformation of Eq. (7) from an open chain to a different open chain and the second-order expansion of the potential become much less straightforward for momentum operators. However, calculating momentum operators can be accomplished by taking an appropriate time-derivative of position operators. In the case of the flux-side correlation function, a reaction rate can be expressed by a time derivative of the side-side correlation function

The symmetrized side-side correlation function,

is solely dependent on position operators and can be used within the OPSCF method resulting in an estimator for the side-side correlation function,

The side-side correlation function has not been explored extensively in the literature outside its relation to the flux-side correlation function. However, it has been shown that the symmetrized side-side correlation function can be analytically solved for a simple parabolic barrier.^{75} Since the transformation and second-order expansion of the potential described in Sec. II are exact for quadratic potentials, one can recover analytically the side-side symmetrized correlation function for the parabolic barrier within the OPSCF method. For more complex systems, Eq. (41) is a proposed sampling scheme to determine side-side correlation functions necessary for reaction rates. Written in this form, the estimator yields an interesting picture such that the estimator is non-zero when one end of the open chain is on the “product” side of the transition state and the other end is on the “reactant” side. While reaction rates can be calculated in principle, further exploration of the OPSCF method in relation to quantum rates is necessary in order to understand any limitations as well as improvements to the current derivation of the method.

We make one final point on the OPSCF method pertaining to the double well potential. An interesting advantage of the OPSCF method is the potential use of enhanced sampling methods to obtain quantum time correlation functions. As stated previously, Eq. (27) can use any preferred sampling technique to produce quantum time correlation functions. Since the observables of interest, *a*(*x*_{1}) and *b*(*x*_{P+1}), are functions of the chain endpoints, enhanced sampling methods can be used to ensure that the end points of the open chain are sampled sufficiently. Furthermore, we can propose an enhanced sampling technique over the entire configuration space in systems where the open chain cannot thermally access regions of phase space, not just the end to end distances. There is an emerging line of research focused on the development of enhanced sampling of “rare” events in quantum mechanics.^{45,46} While this has not been explored in the present work, quantum time correlation functions of “rare” events would be an interesting area of research to explore, which neither CMD or RPMD can achieve on their own.

## V. CONCLUSION

Given the challenge of computing quantum time correlation functions numerically, new perspectives on this problem are needed. In this paper, we have proposed an approach that draws inspiration from semiclassical methodology but that aims to work within an imaginary-time bead sampling paradigm rather than a trajectory-based one. In so doing, we have generated a scheme that can leverage all of the available sampling and enhanced sampling schemes applicable for imaginary-time path integrals, thereby offering the potential for considerable savings over trajectory-based formulations. The approach exploits a transformation that is often used in semiclassical schemes, namely, to sums and differences of the forward and backward complex-time paths of the symmetrized correlation function *G*_{AB}(*t*), together with an expansion of the potential in powers of the difference variables. The transformation is applied, in the present context, in such a way that, after analytical integration over the difference variables, the result is an open-chain imaginary-time path integral subject to a modified potential that couples all of the path-integral beads. Therefore, we have given the new method the acronym of the OPSCF (Open Path Symmetrized Correlation Function) method. Despite it being an imaginary time formulation, it, nevertheless, explicitly accounts for phase information up to second-order in the path difference variables. The new scheme is applied to a set of one-dimensional problems that are often used to benchmark new approximate quantum dynamical methods, as exact results are available for them. It is important to note that the new formulation can be derived within a higher-order Trotter factorization, such as the Suzuki-Chin^{47,48} technique (see Appendix B). In addition, higher-order terms in the path difference variables can be employed within a cumulant expansion approach, thus suggesting one possible path for systematic improvement of the approximation. The integrals over the path difference variables can still be done analytically, and the result is manifestly real and formulated as an imaginary-time sampling problem. An interesting consequence of this formulation as an imaginary-time sampling problem, as can be seen in the mathematical structure of Eq. (27), is that the behavior of a correlation function in time *t* can be rather easily anticipated qualitatively: One only needs to predict the rough behavior of the distribution of the endpoints of an open chain as the length of the chain grows, which necessarily occurs when *t* increases.

In the present paper, we employed open path-integral molecular dynamics (OPIMD) to sample the resulting expression for the approximate correlation function, however, any number of sampling methods could be applied, including Monte Carlo,^{43} molecular dynamics,^{13} hybrid Monte Carlo,^{13,44} or path-integral metadynamics.^{45,46} At the highest temperature studied (*T* = 1), the results are at least as accurate as the popular ring-polymer molecular dynamics (RPMD) approach and the TGA-LSC-IVR method for linear operators and show improved accuracy for nonlinear operators. At the lower temperatures studied, we start to see the breakdown of the second-order expansion. Interestingly, while RPMD and TGA-LSC-IVR also breaks down at these lower temperatures, the behavior of these methods is rather different, something that merits further investigation. At low temperatures, all methods reproduce the short-time dynamics reasonably well for linear operators. For nonlinear operators, the case of the harmonic oscillator shows that OPSCF is considerably better in the short-time limit at the low temperature studied.

We feel that the OPSCF scheme offers an interesting perspective on the formulation of quantum time correlation functions as an imaginary-time sampling problem and that it shows considerable promise. Future work on the OPSCF method will need to address several issues. First, different sampling techniques need to be explored to determine which can best treat the modified potential *Ṽ*(*r*_{1}, …, *r*_{P+1}) appearing in Eq. (27). Given that this potential involves the Hessian, it is expected that either Monte Carlo or reweighted molecular dynamics will be the most efficient. However, this raises a second issue that must be addressed, namely, how best to incorporate the Hessian. Is the exact, analytical Hessian truly needed or could it be approximated in some way? Third, we will need to see how much the higher-order terms in the path difference variables and the cumulant expansion approach improve the accuracy of the results. Fourth, it could be useful to use the OPSCF results as input to a Bayesian^{76} or other types of machine learning schemes as part of an iterative improvement approach. Fifth, instead of performing a second-order expansion and integrating over the difference variables, it would be interesting to consider a complex Langevin equation approach^{77–80} applied to both *r* and *s* variables, allowing the latter to be expanded up to a few higher orders. Complex Langevin dynamics may improve the OPSCF method for position correlation functions and also allow for direct calculation of momentum correlation functions. Sixth, the OPSCF transformation could be applied to the Kubo-transformed correlation function to see if this function is more amenable to the second-order expansion in the path-difference variables. Towards the end of Sec. IV D, we discussed the modification of the OPSCF method needed for the calculation of momentum correlation functions. As noted, the equivalent of Eq. (4) for momentum operators does not result in a cyclic ring polymer but rather, an open chain, which makes the transformation of Eq. (7) and the second-order expansion less straightforward. While time derivatives of position operators can determine some momentum-dependent operators, it is necessary to devise alternative schemes, like complex Langevin dynamics, that may be beneficial in the computation of momentum correlation functions and general correlation functions for systems when the second-order expansion fails. Finally, it will be important to test the approach on a true condensed-phase system, such as liquid *para*-H_{2} or water in order to compare to other methods that have been used to obtain spectral or transport properties from correlation functions for these systems.

## SUPPLEMENTARY MATERIAL

See supplementary material for additional figures reporting (1) convergence of the OPSCF correlation functions of a harmonic oscillator at *β* = 1 and *β* = 8 as a function of the number of path-integral beads; (2) OPSCF correlation functions of the weakly anharmonic and quartic potentials computed using the Suzuki-Chin factorization at *β* = 1 and *β* = 8; and (3) Kubo-transformed correlation functions from RPMD simulations for all four potentials studied.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge Nancy Makri for informative and stimulating discussions. This work was supported by the National Science Foundation Award No. CHE-1566085 and the Margaret Strauss Kramer Fellowship.

### APPENDIX A: FIRST CUMULANT EXPANSION

As discussed throughout the text, cumulant expansions provide a route to the systematic improvement to the OPSCF method by rendering the integrals generated when including higher-order terms in *s*_{i} amenable to analytical evaluation. Depending on the order of expansion, different contributions to the method would be improved. Specifically, the positive definite distribution, Eq. (14), involves only even orders of *s*_{i} while the phase factor, Eq. (15), involves only odd orders of *s*_{i}. As an illustration of this idea, this section of the appendix will examine the first cumulant corresponding to an expansion of the potential energy up to third order in *s*_{i}, which occurs through the phase factor and, therefore, can only affect *t* ≠ 0 results. We begin with the expansion of the potential energy,

This results in a modified phase factor different from Eq. (15) while the distribution *ρ* remains the same,

If we assume that the third-order term is merely a correction, we can write the phase factor as

where

The integration of *S* = (*s*_{2}, …, *s*_{P}) can now be written as

where the average is taken with respect to the second-order expansion. We approximate the average of the phase correction,

From the above equation, we can derive a correction term to the OPSCF method. Specifically, we can define vectors *V*_{3,i} = *U*‴(*r*_{i}) and $S3,i=si3$, where *i* = 2, *…*, *P* such that

In order to perform this integral, we must first redefine *M* = *O*^{T}*DO* where *D* is a diagonal matrix, O is the orthogonal matrix which diagonalizes M, *s* = *O*^{T}*σ* where, since det*O* = 1, *ds* = *dσ*, and *K*^{T} = *λ*^{T}*O*,

Expanding the cubic function accordingly results in a sum of Gaussian integrals

By completing the square,

we can define an appropriate substitution, $yi=\sigma i\u2212i\lambda iDii$, whose inverse is $\sigma i=yi+i\lambda iDii$ to solve the set of integrals,

Expanding the terms yields

We recognize that the first four terms have odd exponents in at least one of the *y*_{m} terms, which, when multiplied by the Gaussian, will integrate to 0. Furthermore, for the next three terms, the integral is non-zero when the indices of *y*_{m} are the same, and the final term is a constant. Thus, the only remaining terms, which interestingly are all imaginary, will be

When the Gaussian integrals are performed, the result we obtain is

If necessary, Eq. (A15) can be transformed back to depend explicitly on the original vector *K*(*r*′) and matrix *M*(*r*). Specifically,

and

where

With this short hand, we can write the average of the correction term as

from which it can be seen that exp(*i*⟨*δϕ*⟩) is purely real. While this functional form is rather complex and would be cumbersome to implement directly into a molecular dynamics scheme, this clearly shows a way to improve the OPSCF method. In principle, higher-order cumulants can be derived as well to better improve the methodology, and the resulting path integral will remain positive definite.

### APPENDIX B: SUZUKI-CHIN FACTORIZATION

The derivation of the OPSCF method using the Suzuki-Chin^{47,48,65,66} factorization follows many of the similar steps as in Sec. II. However, we will point out some of the differences in this appendix. Unlike the Trotter factorization, the Suzuki-Chin factorization decomposes the complex-time propagator as follows:

where *ϵ* = *iτ*_{c}/2*ℏ* and $\u03f5\u2009=\u2009\u2212i\tau c*/2\u210f$ for the forward and backward complex-time propagators, respectively, and $Ue^$ and $Um^$, in one-dimensions for a single particle, are defined by

where *α* is a constant in the range of [0,1] and specifically chosen to be zero in our implementation. Following closely the derivation of the path integral representation of forward and backward complex-time propagators of Krilov,^{76} we obtain new expressions for the phase factor and positive-definite distribution that contain first-order derivatives of the potential,

After performing the transformation described in Eq. (7), we obtain a separable distribution function, *ρ*(*r*_{1}…*r*_{P+1}, *s*_{2}…*s*_{P}) = *ρ*_{1}(*r*_{1}…*r*_{P+1})*ρ*_{2}(*r*_{2}…*r*_{P}, *s*_{2}…*s*_{P}),

The phase factor, under the transformation, is now defined as

With the new definition of *ρ*_{2}, the matrix *M* with the same (*P* − 1) × (*P* − 1) dimensions is

where we introduce a potential energy *U*_{M} that is dependent on the parity of the bead number,

A similar procedure is done for the phase factor and the vector *K* in which we introduce the potential *U*_{K},

With these definitions of *M*(*r*) and *K*(*r*′), the rest of the derivation follows exactly as in Sec. II. The resulting Hamiltonian for a path integral scheme is as follows:

In theory, one could introduce a weighting function to mimic a more traditional OPSCF scheme,^{66} but we implemented the above Hamiltonian for this current project. As noted in Sec. IV, this higher-order factorization improves the correlation functions for the anharmonic potential but is insufficient for the quartic oscillator.

## REFERENCES

^{4}He

^{N}as real probabilities on C

^{N}and an application to path integrals