We apply the Feynman-Kleinert Quasi-Classical Wigner (FK-QCW) method developed in our previous work [Smith *et al.*, J. Chem. Phys. **142**, 244112 (2015)] for the determination of the dynamic structure factor of liquid *para*-hydrogen and *ortho*-deuterium at state points of (*T* = 20.0 K, *n* = 21.24 nm^{−3}) and (*T* = 23.0 K, *n* = 24.61 nm^{−3}), respectively. When applied to this challenging system, it is shown that this new FK-QCW method consistently reproduces the experimental dynamic structure factor reported by Smith *et al.* [J. Chem. Phys. **140**, 034501 (2014)] for all momentum transfers considered. This shows that FK-QCW provides a substantial improvement over the Feynman-Kleinert linearized path-integral method, in which purely classical dynamics are used. Furthermore, for small momentum transfers, it is shown that FK-QCW provides nearly the same results as ring-polymer molecular dynamics (RPMD), thus suggesting that FK-QCW provides a potentially more appealing algorithm than RPMD since it is not formally limited to correlation functions involving linear operators.

## I. INTRODUCTION

Low-temperature liquid *para*-hydrogen and *ortho*-deuterium have become benchmarks in the development of approximate quantum dynamics methods^{1–7} that allow for the practical evaluation of a general quantum time correlation function of the form

*Z* being the partition function and *β* the inverse temperature 1/*k _{B}T*. The reason that these systems provide such good testing grounds for the development of these methods is because pronounced nuclear quantum effects are exhibited by their dynamical properties. This is due to their low molecular masses, which in turn cause their thermal de Broglie wavelength to be relatively large at low temperatures. However, these quantum effects are not significant enough that one must worry about the quantum statistics of molecular indistinguishability,

^{8}and in addition, there exists a relatively simple pair potential,

^{9}which provides a very accurate description of the molecular interactions.

^{9–11}Hence, low-temperature liquid

*para*-hydrogen and

*ortho*-deuterium are quantum liquids which are relatively easy to model, thus being ideal testing grounds for the development of novel approximate quantum dynamics methods.

Presently, the most successful approximate quantum dynamics methods are the Classical Wigner (CW) approximation,^{12–14} Centroid Molecular Dynamics (CMD),^{15} and Ring-Polymer Molecular Dynamics (RPMD).^{16} All of them have been shown to provide relatively accurate and practical approximations to Eq. (1), and, in addition, to become exact in the harmonic, high temperature and short time limits.^{12,13,15–17} However, each of these methods has its own downfalls. For example, CMD and RPMD begin to break down for correlation functions involving non-linear operators,^{15,16,18} while the CW approximation is equally valid for non-linear operators, but, in general, it does not produce time invariant thermodynamic properties for systems at thermal equilibrium. Explicitly, for $ A \u02c6 =1$, the exact quantum expression in Eq. (1) has the property that

while the purely classical propagation of the initially quantized phase space distribution in the CW approximation does not ensure this property. As shown in Ref. 19, this downfall can have a significant impact for slow processes like diffusion due to zero point energy leakage from intramolecular to intermolecular modes as the system is propagated.

Recently, Liu and Miller^{17,20–22} proposed a route to remedy this downfall of the CW approximation by replacing the classical propagation of the initial quantum phase space distribution with a form of dynamics that ensures the equality in Eq. (2). Similarly, we proposed in Ref. 23 a different form of dynamics that also ensures this property by requiring the dynamics to conserve the initial quantum ensemble within the Feynman-Kleinert (FK) approximation of the density operator.^{12,24,25}

In the present paper, version 2 of the Feynman-Kleinert Quasi-Classical Wigner method introduced in Ref. 23 will be used and referred to as FK-QCW. It was shown to greatly extend the accuracy of the Feynman-Kleinert linearized path-integral (FK-LPI) implementation of the CW approximation in the challenging model problems of both the quartic and double well potentials, in which numerically exact solutions are obtainable. Furthermore, we were able to show that one can introduce an arbitrary frequency function into the dynamics of this method, resulting in an entire class of ensemble conserving dynamics that preserve the equality in Eq. (2). The FK-QCW method was shown to recover the exact classical and high temperature limits of the quantum time correlation function and may be applied to systems with barriers.^{23} Furthermore, since this method is developed within the framework of the CW approximation, no problems arise for non-linear operators.

The purpose of this work is to test how well the FK-QCW method performs when applied to the standard benchmark systems of low-temperature liquid *para*-hydrogen and *ortho*-deuterium. We accomplish this task by computing the dynamic structure factor, which is experimentally accessible by inelastic X-ray scattering (IXS). We then provide a comparison between the present calculations and the experimental determinations, as well as with the ones obtained by RPMD and FK-LPI as previously published in Ref. 7. Specifically, it was found in Ref. 7 that for a momentum transfer of *k* = 20 nm^{−1}, the FK-LPI method fails to correctly reproduce the IXS spectrum for the *para*-hydrogen system, in which quantum effects are more prevalent. In addition, due to the increased non-linearity of the correlation function at high momentum transfers, RPMD has been shown in the case of *para*-hydrogen to be limited to a maximal momentum transfer value of *k* = 15 nm^{−1}, see Ref. 1. For larger values of *k*, RPMD becomes inaccurate. Hence, a challenging test case has been established for the development of improved methods, and it is therefore interesting to check how the FK-QCW method performs where these leading methods fail.

This paper is organized as follows: In Sec. II, we provide an introduction to the CW approximation to quantum time correlation functions, as well as the multidimensional generalization of the FK-QCW method. In addition, we also give a brief introduction to the theory of inelastic scattering. In Sec. III, we begin by discussing the computational details of our simulations, followed by a comparison of the FK-QCW method with the experimental dynamic structure factor and that obtained by RPMD and FK-LPI previously published in Ref. 7. The conclusions are presented in Sec. IV.

## II. THEORY AND METHODOLOGY

### A. Classical Wigner

The CW^{12–14} expression for a general quantum time correlation function of a many-body system is given by

where $ q ( t ) , p ( t ) $ are the classically evolved quantum phase space variables propagated from the initial quantum distribution $ q , p $. Here, the Wigner transform of a general operator $ C \u02c6 $ is given by

where $ q $ is the direct product of the single particle position kets.

Although the CW approximation has been shown to perform relatively well,^{3,7,12,13,26–28} as we previously noted, the classical evolution of the quantum phase space results in thermodynamic properties of equilibrium systems being incorrectly time dependent. Our newly developed FK-QCW method corrects this inconsistency by replacing the purely classical dynamics used within Eq. (3) with a time evolution that ensures that the initial quantum ensemble is conserved such that $\u3008 B \u02c6 ( t ) \u3009=\u3008 B \u02c6 ( 0 ) \u3009$, in accord with the exact quantum time correlation function. However, since these dynamics were developed within the FK approximation to the density operator,^{12,24,25} before we present their multidimensional generalization, we first provide an introduction to the FK density operator, which allows for a practical evaluation of the Wigner function $ [ e \u2212 \beta H \u02c6 A \u02c6 ] W q , p $ appearing in Eq. (3).

### B. Many-body Feynman-Kleinert density operator

The most difficult part in evaluating the CW expression of Eq. (3) is obtaining the Wigner transform of $ e \u2212 \beta H \u02c6 A \u02c6 $, since knowledge of the many-body density matrix is required. As in the FK-LPI method, we accomplish this within the FK-QCW method by combining the effective frequency variational theory of Feynman^{24} and Kleinert^{25} with the quasidensity operator (QDO) formalism of Jang and Voth.^{29} This Feynman-Kleinert approximation to the density operator allows for an efficient evaluation of the Wigner transform of $ e \u2212 \beta H \u02c6 A \u02c6 $ and has been shown to be very accurate when applied to realistic many-body systems.^{3,7,12,26–28} In addition, the FK approximation to the density operator gives the best local harmonic approximation to the systems free energy^{24,25} and becomes exact in the harmonic and high temperature limits.^{12}

For a many-body system, the FK approximation to the density operator is explicitly given by

where $ x c , p c $ are the 3*N* dimensional vectors of centroid positions and momenta and $ \delta \u02c6 FK ( x c , p c ) $ is the effective frequency QDO. The FK approximation to the centroid phase space density, for a system of *N* particles, is given by

$ W 1 x c $ being the FK approximation to the centroid potential. In Eq. (6), ** M** is the diagonal matrix of particle masses.

The FK variational effective frequency matrix is determined from the local curvature of the system’s Gaussian smeared potential by

** H**(

**) being the 3**

*q**N*× 3

*N*classical Hessian matrix and

**(**

*A*

*x*_{c}) the smearing width matrix which measures the importance of quantum fluctuations around the classical-like centroid positions. Defining

**(**

*U*

*x*_{c}) as the orthonormal matrix containing the eigenvectors of the effective frequency matrix, then

where *ω*^{2}(*x*_{c}) is the 3*N* dimensional vector of eigenvalues and *I* is the identity matrix. Using this, one can define the mass-weighted normal modes as

and the smearing width matrix can be diagonalized through

where

In terms of the eigenvalues of the effective frequency matrix, the centroid potential in Eq. (6) is explicitly written as

where

is the FK smeared potential. In 3*N* dimensions, the effective frequency QDO written in terms of the mass-weighted normal modes is simply a direct product of 1-dimensional QDOs and is given by

with

For all but the simplest potentials, determination of the FK effective frequency matrix is the main computational load in applying the FK approximation to the density operator since Eqs. (7) and (10) must be solved iteratively. For a discussion of efficient ways to determine the FK effective frequency matrix using different numerical schemes, the interested reader is referred to Ref. 27.

Using the FK approximation to the density operator, the Wigner transform of $ e \u2212 \beta H \u02c6 A \u02c6 $ becomes

Due to the Gaussian form of the QDO in Eq. (14), an analytical expression for the Wigner transform of $ \delta \u02c6 FK ( x c , p c ) A \u02c6 $ is readily obtained for any operator $ A \u02c6 $ depending only on position or momentum.

### C. FK-QCW in many dimensions

The multi-dimensional generalization of the FK-QCW dynamics derived in Eq. (67) of Ref. 23 simply becomes

where

are the *k*th elements of the dimensionless normal mode coordinates $ ( q \u0303 ( t ) , p \u0303 ( t ) ) $, and *f _{k}*(

*x*_{c}(

*t*)) is an arbitrary frequency function. As shown in Ref. 23, the 1-dimensional version of these dynamics gives the exact real, but not imaginary, part of the position autocorrelation function in the harmonic limit if the FK effective frequency is chosen for the frequency function. Similarly, one can show that this exact limit is also obtained in the multi-dimensional case if we choose

In Eq. (18), the time evolved normal modes are given by

where ** U**(

*x*_{c}(

*t*)) diagonalizes the effective frequency matrix evaluated at

*x*_{c}(

*t*), and the centroid dynamics are governed by the classical-like equations

the gradient with respect to *x*_{c}(*t*) being taken while holding the smearing width matrix constant.

Once the centroid and dimensionless normal mode coordinates have been propagated through Eqs. (17) and (21), the instantaneous quantum phase space variables can be obtained from

where

and in terms of the dimensionless normal mode coordinates, they have the elements

The FK-QCW approximation to the quantum time correlation function takes the same form as the CW expression in Eq. (3),

However, within FK-QCW, $ [ e \u2212 \beta H \u02c6 A \u02c6 ] W q , p $ is obtained through Eq. (16), and $ q ( t ) , p ( t ) $ are propagated from the initial quantum distribution $ q , p $ using the ensemble conserving dynamics in Eq. (17). It should be noted that the Wigner transform of $ B \u02c6 $ in this expression is easily obtainable when $ B \u02c6 $ is a function of only position or momentum operators; since in this case, one obtains the corresponding classical expression

or

For completeness, we note that the multi-dimensional generalization of the dynamics within the FK-QCW method that was derived in Eq. (43) of Ref. 23 is very similar to Eq. (17), the only difference is that $ p \u0303 k ( t ) $ in Eq. (18) takes the form

and Eq. (22) is changed accordingly. As shown in Ref. 23 (see Appendix C therein), these dynamics reduce to CMD when used within the expression for the Kubo transformed time correlation function and also become exact in the harmonic limit, if one uses Eq. (19) for the frequency function. However, they are not very practical when used within Eq. (25) for non-linear operators since they are ill-defined when *ω _{k}*(

*x*_{c}(

*t*)) becomes imaginary. This then makes the FK-QCW dynamics in Eq. (17) preferable for condensed phase systems since one can show that these dynamics are well defined in this case as long as |

*ω*(

_{k}

*x*_{c}(

*t*))| <

*π*/

*β*ħ. In addition, by adopting the same convention as FK-LPI

^{12}in which we set

*ν*(

_{k}*t*) = 0 for frequencies outside of this range, we are then able to apply these FK-QCW dynamics for imaginary frequencies as long as |

*ω*(

_{k}

*x*_{c}(

*t*))| < 2

*π*/ħ

*β*, which is the entire range where the FK approximation to the density operator is well defined.

### D. Inelastic scattering

The quantum time correlation function we apply the FK-QCW method to in this paper is the intermediate scattering function given by

where $ \cdots $ denotes a canonical ensemble average, and $ x \u02c6 i ( t ) $ is the time-dependent position operator of the *i*th particle, given in the Heisenberg picture by $ x \u02c6 i ( t ) = e i H \u02c6 t / \u0127 \u2009 x \u02c6 i \u2009 e \u2212 i H \u02c6 t / \u0127 $. This correlation function is related to the dynamic structure factor, *S*(** k**,

*ω*), through

The dynamic structure factor gives the spectrum of density fluctuations and Van Hove^{30} showed that within the first Born approximation, this quantity is proportional to the inelastic scattering cross section. Measured by either inelastic neutron or X-ray scattering, the inelastic scattering cross section measures the probability that a neutron or photon transfers momentum ħ** k** = ħ(

*k*_{f}−

*k*_{i}) and energy ħ

*ω*= ħ(

*ω*−

_{f}*ω*) to the sample.

_{i}The shape of the dynamic structure factor is defined by its *n*th spectral moments,

The comparison between measured and computed values of the spectral moments can be used as a metric to judge the quality of a theoretical simulation. The zeroth moment of the dynamic structure factor is referred to as the static structure factor,^{31} *S*(** k**), which is related to the spatial Fourier transform of the pair distribution function by

The first moment of *S*(** k**,

*ω*) is of special interest since for a system which interacts through a momentum-independent potential,

^{32}this quantity is exactly given by

where *m* is the molecular mass. This relation, in principle valid for monatomic systems, can be extended to molecular ones provided that the rotational and vibrational motions of the molecule can be neglected in the probed dynamic range.

Since the dynamic structure factor is related to the temporal Fourier transform of a quantum time correlation function, it obeys the principle of detailed balance,

which is straightforward to show by working in the basis of energy eigenstates. Using this principle of detailed balance along with the symmetry of Eq. (29), *F*(** k**,

*t*) =

*F*(

**, −**

*k**t*)

^{∗}, one can show that the dynamic structure factor can be equivalently written using only the real part of the intermediate scattering function as

For an isotropic system such as a liquid, the intermediate scattering function depends only on the magnitude of ** k**, such that

*F*(

**,**

*k**t*) =

*F*(

*k*,

*t*). Hence when computing the intermediate scattering function of an isotropic system, we are free to choose the direction of

**.**

*k*The FK-QCW approximation to the intermediate scattering function takes the form

where, after choosing ** k** to be parallel to the x-axis,

and in terms of the FK approximation to the density operator,

where

## III. RESULTS

### A. Computational details

To obtain the dynamic structure factor for liquid *para*-hydrogen and *ortho*-deuterium at the state points (*T* = 20.0 K, *n* = 21.24 nm^{−3}) and (*T* = 23.0 K, *n* = 24.61 nm^{−3}), respectively, the FK-QCW approximation to the intermediate scattering function in Eq. (36) was evaluated by using the Silvera-Goldman (SG) potential.^{9} The SG potential has been used in a number of previous studies^{1,2,4–6,17} and has been shown to provide very accurate descriptions of the fluid and solid thermodynamics, except at extremely high pressure.^{10,11} This semi-empirical isotropic pair potential, applicable to both *para*-hydrogen and *ortho*-deuterium, treats each molecule as a spherical particle which is justifiable at low temperatures since only the *J* = 0 rotational state is populated in each isotopoloque. To expedite the determination of the FK centroid potential, we represented the SG potential as a sum over four Gaussian functions whose parameters can be found in Table II of Ref. 3.

Starting from an equilibrated centroid configuration, the real part of the intermediate scattering function was evaluated using the FK-QCW approximation by performing a molecular dynamics simulation of the centroid variables, which are propagated according to Eq. (21). For each centroid trajectory, 100 sets of initial dimensionless normal mode coordinates $ ( q \u0303 ( 0 ) , p \u0303 ( 0 ) ) $ were sampled according to their Gaussian distribution

which appears in the FK approximation of $ [ e \u2212 \beta H \u02c6 A \u02c6 ] W q , p $ (see, for example, Eq. (38)). Using the frequency function in Eq. (19), the dimensionless normal mode coordinates were then propagated through Eq. (17) and the instantaneous (** q**(

*t*),

**(**

*p**t*)) values were obtained from the relations in Eq. (22). The real part of the intermediate scattering function was constructed by averaging over 1000 consecutive 3 ps centroid trajectories, in which the centroid momentum was resampled at the beginning of each trajectory and a time step of 1 fs was used. In order to obtain statistical uncertainties, this entire process was repeated six times. For a more detailed outline of the algorithm used for the integration of the centroid and dimensionless normal mode coordinates over one time step Δ

*t*, which allows for the propagation of (

**(**

*q**t*),

**(**

*p**t*)) → (

**(**

*q**t*+ Δ

*t*),

**(**

*p**t*+ Δ

*t*)), the interested reader is referred to the Appendix.

Within the simulation, we employed cubic periodic boundary conditions with the minimum image convention and a spherical cutoff at half box length. Furthermore, due to the isotropic nature of the SG potential, the momentum transfer ** k** was chosen, without loss of generality, to be in the

*x*-direction. In order to fulfill the Laue condition

^{33}

*k*= 2

*πn*/

*l*, where

*l*is the length of the simulation cell and

*n*is an integer, the number of particles

*N*treated in the simulation had to be varied for each momentum transfer.

*N*,

*l*, and the density used for each

*k*are listed in Table I.

. | D_{2}
. | H_{2}
. | ||||
---|---|---|---|---|---|---|

k (nm^{−1})
. | N
. | l (bohr)
. | ρ (nm^{−3})
. | N
. | l (bohr)
. | ρ (nm^{−3})
. |

5.5 | 37 | 21.65 | 24.61 | 32 | 21.66 | 21.24 |

12.8 | 78 | 27.76 | 24.61 | 68 | 27.85 | 21.24 |

15.3 | 109 | 31.03 | 24.61 | 94 | 31.03 | 21.24 |

20.0 | 95 | 29.64 | 24.61 | 82 | 29.65 | 21.24 |

. | D_{2}
. | H_{2}
. | ||||
---|---|---|---|---|---|---|

k (nm^{−1})
. | N
. | l (bohr)
. | ρ (nm^{−3})
. | N
. | l (bohr)
. | ρ (nm^{−3})
. |

5.5 | 37 | 21.65 | 24.61 | 32 | 21.66 | 21.24 |

12.8 | 78 | 27.76 | 24.61 | 68 | 27.85 | 21.24 |

15.3 | 109 | 31.03 | 24.61 | 94 | 31.03 | 21.24 |

20.0 | 95 | 29.64 | 24.61 | 82 | 29.65 | 21.24 |

The dynamic structure factor was obtained from the real part of the FK-QCW approximation to the intermediate scattering function using the relation in Eq. (35), which ensures that this quantity fulfills the detailed balance condition. The imaginary part of the intermediate scattering function was then obtained from the inverse Fourier transform of this quantity. By processing the data in this manner, we ensure that the approximate quantum time correlation function fulfils the detailed balance condition. Furthermore, since the detailed balance condition in Eq. (34) holds for the spectrum of any quantum time correlation function, this procedure could be applied to any correlation function involving hermitian operators;^{34} since in this case, the general relation

### B. The experimental dynamic structure factor

In an inelastic x-ray scattering experiment on a liquid, one must cope with spurious scattering from the container which must be carefully subtracted from the raw intensity. Due to random temperature drifts in the analyzer crystal^{7} which introduce an uncertainty in the zero of energy transfer (ħ*ω* = 0) between the sample + container and empty container scattering measurements, such a subtraction was not completely straightforward for the experimental data used in our previous study.^{7} In order to deal with these spurious experimental effects, we developed in Ref. 7 a method that enables one to use the dynamic structure factor obtained from a simulation as an input to fix this unknown spectral shift.

Denoting *I _{raw}*(

**,**

*k**ω*) and

*I*(

_{EC}**,**

*k**ω*) as the measured intensities scattered by the sample + container system and by the empty container, respectively, the unknown shift in the zero of energy transfer, $ \delta \u2212 \theta $, can be taken into account within the experimental dynamic structure factor by writing this quantity as

where *α*(** k**) is the proportionality factor that relates the dynamic structure factor to the inelastic scattering cross section, and

*T*(

**) is the transmission coefficient of the sample.**

*k*^{7}Due to the finite resolution of the experiment, the experimental dynamic structure factor,

*S*(

_{exp}**,**

*k**ω*), in Eq. (42) is related to the sample’s true dynamic structure factor,

*S*(

**,**

*k**ω*), through a convolution with the instrumental resolution function,

*R*(

*ω*), such that

Using the procedure developed in Ref. 7, the refined experimental dynamic structure factor is obtained by performing a least squares fit of Eq. (42) to the results of a theoretical simulation which have been convoluted with the instrument resolution function, *R*(*ω*), for the determination of the unknown spectral shifts (*δ*, *θ*). Within this fitting procedure, *α*(** k**) is determined for each (

*δ*,

*θ*) through

which ensures that the refined experimental quantity fulfills the first moment sum rule in Eq. (33). We note that *S*(** k**) in Eq. (44) is the static structure factor, which in the present work has been obtained from the theoretical simulation. However, if the experimental measurement of this quantity was available, one could in principal use it within the fitting procedure.

Within this study, we have reperformed this refinement process of the experimental quantity using the FK-QCW method as an input, and the resulting fitting parameters are shown in Table II. In Sec. III C, we provide a comparison of the new experimental dynamic structure factors obtained using the FK-QCW method with the experimental quantities previously published in Ref. 7, which were obtained using FK-LPI as an input. In addition, we also compare the theoretical dynamic structure factors obtained using the FK-QCW method with the results of FK-LPI and RPMD, also published in Ref. 7.

D_{2}
. | ||||||
---|---|---|---|---|---|---|

k (nm^{−1})
. | α()
. k | T()
. k | δ (meV)
. | θ (meV)
. | |δ − θ|(meV)
. | R^{2}
. |

5.5 | 7.51 × 10^{4} | 0.897 | −0.616 | −0.545 | 0.070 | 0.97 |

12.8 | 1.41 × 10^{5} | 0.896 | −0.090 | 0.143 | 0.233 | 0.99 |

15.3 | 2.12 × 10^{5} | 0.895 | −0.391 | −0.205 | 0.186 | 0.99 |

20.0 | 3.78 × 10^{5} | 0.893 | −0.678 | 0.228 | 0.906 | 0.99 |

H_{2} | ||||||

|δ − θ| | ||||||

k (nm^{−1}) | α() k | T() k | δ (meV) | θ (meV) | (meV) | R^{2} |

5.5 | 5.27 × 10^{4} | 0.731 | −0.134 | 0.233 | 0.367 | 0.97 |

12.8 | 1.79 × 10^{5} | 0.954 | 0.183 | 0.531 | 0.348 | 0.93 |

15.3 | 1.57 × 10^{5} | 0.740 | 0.188 | 0.929 | 0.742 | 0.97 |

20.0 | 4.05 × 10^{5} | 0.952 | −0.734 | −0.703 | 0.032 | 1.00 |

D_{2}
. | ||||||
---|---|---|---|---|---|---|

k (nm^{−1})
. | α()
. k | T()
. k | δ (meV)
. | θ (meV)
. | |δ − θ|(meV)
. | R^{2}
. |

5.5 | 7.51 × 10^{4} | 0.897 | −0.616 | −0.545 | 0.070 | 0.97 |

12.8 | 1.41 × 10^{5} | 0.896 | −0.090 | 0.143 | 0.233 | 0.99 |

15.3 | 2.12 × 10^{5} | 0.895 | −0.391 | −0.205 | 0.186 | 0.99 |

20.0 | 3.78 × 10^{5} | 0.893 | −0.678 | 0.228 | 0.906 | 0.99 |

H_{2} | ||||||

|δ − θ| | ||||||

k (nm^{−1}) | α() k | T() k | δ (meV) | θ (meV) | (meV) | R^{2} |

5.5 | 5.27 × 10^{4} | 0.731 | −0.134 | 0.233 | 0.367 | 0.97 |

12.8 | 1.79 × 10^{5} | 0.954 | 0.183 | 0.531 | 0.348 | 0.93 |

15.3 | 1.57 × 10^{5} | 0.740 | 0.188 | 0.929 | 0.742 | 0.97 |

20.0 | 4.05 × 10^{5} | 0.952 | −0.734 | −0.703 | 0.032 | 1.00 |

### C. Dynamic structure factor

The FK-QCW results for the intermediate scattering function are shown in Fig. 1. Also included in Fig. 1 are the results obtained by FK-LPI and RPMD previously published in Ref. 7. Interestingly, for low momentum transfers where RPMD is least troubled by the non-linearity of the correlation function in Eq. (29), the FK-QCW method predicts nearly the same results as RPMD. This suggests that FK-QCW provides a method with an accuracy comparable to RPMD. Yet, FK-QCW does not suffer from the non-linear operator problem, since it was developed within the framework of the CW approximation and can therefore be used to evaluate correlation functions involving non-linear operators.

When comparing the intermediate scattering functions of FK-QCW and FK-LPI in Fig. 1, one notices that a significantly longer decay time is predicted by FK-QCW for both *para*-hydrogen and *ortho*-deuterium at *k* = 20.0 nm^{−1}. In Ref. 7, we attributed the failure of FK-LPI for *para*-hydrogen at this momentum transfer to the purely classical propagation within this method not being able to correctly account for the long-time behavior of the correlation function. The fact that the dynamics within FK-QCW predicts a significantly longer decay time suggests that this method may be able to stand up to this challenging test case where FK-LPI failed and RPMD was inaccurate.

The hypothesis just proposed is confirmed in Fig. 2, where for the case of *para*-hydrogen at *k* = 20.0 nm^{−1}, FK-QCW is in almost exact agreement with the experimental dynamic structure factor obtained using this method as an input. This then shows that the ensemble conserving dynamics of FK-QCW extends the accuracy of the FK-LPI method to longer times, where the classical propagation within the CW approximation fails. In addition, one sees that the experimental quantities obtained using FK-QCW or FK-LPI as an input are in relatively good agreement for all of the momentum transfers considered. The fact that this is true for *para*-hydrogen at *k* = 20.0 nm^{−1}, where FK-LPI fails, shows that the process we developed in Ref. 7 to refine the experimental dynamic structure factor provides a robust method that allows one to correct for spurious experimental effects, even when the theoretical input is relatively inaccurate.

For the case of *ortho*-deuterium at *k* = 20.0 nm^{−1}, one sees in Fig. 3 that, similar to the case of *para*-hydrogen, the FK-QCW method reproduces the experimental quantity almost exactly. While FK-LPI is in relatively good agreement with the experimental quantity for this momentum transfer, the fact that FK-QCW is more accurate shows that the longer decay time predicted in the correlation function by the ensemble conserving dynamics of this method, once again, better reflects the true dynamics of the system.

As seen in Figs. 2 and 3, the similarity of the FK-QCW and RPMD intermediate scattering functions for both *para*-hydrogen and *ortho*-deuterium at low momentum transfers translates to similar predictions for their corresponding dynamic structure factors. While not completely systematic due to the case of *para*-hydrogen at *k* = 12.8 nm^{−1}, overall one sees a better agreement between the predictions of these two methods and the experimental quantity for both systems, as compared to FK-LPI. This then suggests that the more pronounced oscillations in the intermediate scattering functions of FK-QCW and RPMD more accurately describe the large wavelength density fluctuations of both the *para*-hydrogen and *ortho*-deuterium systems, consistent with the ensemble conserving dynamics used within FK-QCW being comparable to RPMD and more accurate than the purely classical dynamics of FK-LPI.

We note that for *para*-hydrogen, particularly at *k* = 12.8 nm^{−1}, there is disagreement between the FK-QCW and the experimental dynamic structure factors. At first glance, this appears to be due to the area underneath the FK-QCW dynamic structure factor being too small, i.e., a shortcoming of the method. However, this area is precisely equal to the static structure factor *S*(** k**) =

*F*(

**,**

*k**t*= 0), for which FK-QCW and FK-LPI are equivalent, to within statistical accuracy. The FK-LPI static structure factor was published in Ref. 7, where it was found to be in agreement with the experimental measurements of Ref. 35. It is reasonable, then, to look for alternative explanations. One possibility is an uncertainty in the absolute scale of the empty container scattering contribution. Adjustment of that scale was found to play a role for the cases of

*para*-hydrogen at

*k*= 5.5 nm

^{−1}and

*k*= 15.3 nm

^{−1}.

The first moment of the dynamic structure factor obtained by FK-QCW, FK-LPI, and RPMD is shown in Fig. 4. The different dynamics of these three methods become relevant for the first moment since this quantity depends on the time derivative of the intermediate scattering function through Eq. (31). However, as seen in Fig. 4, FK-QCW reproduces the exact values of the first moment of the dynamic structure factor to within the uncertainty of the simulation for all of the momentum transfers considered. While FK-LPI is also able to reproduce this quantity for all of the momentum transfers considered, this is not true for RPMD. The fact that these new ensemble conserving dynamics are shown to obey this exact quantum mechanical sum rule testifies even more to the accuracy of this new method.

## IV. CONCLUSIONS

We have shown that the ensemble conserving dynamics of the FK-QCW method can be easily applied to realistic condensed phase systems, such as low temperature *para*-hydrogen and *ortho*-deuterium. Furthermore, it was found that the dynamics of the FK-QCW method greatly extend the accuracy of the FK-LPI approximation when the long-time behavior of the quantum time correlation function becomes important. This was evidenced by the FK-QCW method nearly exactly reproducing the experimental dynamic structure factor of *para*-hydrogen at *k* = 20.0 nm^{−1}, where the purely classical dynamics used within the FK-LPI approximation fail and RPMD is too inaccurate.

In addition, where RPMD was applicable due to the approximate linearity of Eq. (29), we found that FK-QCW provides an accuracy comparable to RPMD. This then suggests that FK-QCW may be a top contender in the realm of approximate quantum dynamics methods, since FK-QCW is able to assess correlation functions involving non-linear operators, and will not encounter artificial frequencies in the simulated absorption spectra arising from unphysical high frequency oscillations within the dynamics.^{19,36,37}

Readers wanting to use these results can obtain the data set containing the refined experimental dynamic structure factors obtained using FK-QCW as an input in Figs. 2 and 3 as well as the non-convoluted FK-QCW dynamic structure factors of *para*-hydrogen and *ortho*-deuterium for all of the momentum transfers considered from the supplementary material.^{38}

## 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: FK-QCW MOLECULAR DYNAMICS ALGORITHM

Defining the 3*N* dimensional vector that contains the FK smeared force as

where the gradient is taken with respect to *x*_{c}(*t*), the algorithm for integrating the FK-QCW method dynamics one time step Δ*t* using the velocity Verlet algorithm goes as follows:

It is important to note intermediate steps (A6) and (A7). These are performed to ensure that the dimensionless normal modes $ q \u0303 ( t + \Delta t ) $ and $ p \u0303 ( t + \Delta t ) $ are formed properly for the use in steps (A8) and (A9) since the eigenvector within the *j*th column of the orthonormal matrices ** U**(

*x*_{c}(

*t*)) and

**(**

*U*

*x*_{c}(

*t*+ Δ

*t*)) will typically not correspond to the same normal mode due to the numerical method used for the diagonalization of

**Ω**

^{2}(

*x*_{c}(

*t*)) and

**Ω**

^{2}(

*x*_{c}(

*t*+ Δ

*t*)).

## REFERENCES

We note that even though the operators in Eq. (29) are not hermitian, the relation *F*(** k**,

*t*) =

*F*(

**, −**

*k**t*)

^{∗}still holds; since in this case,

*A*

^{∗}=

*B*, and for non-hermitian operators, one can show that

*C*(−

_{AB}*t*)

^{∗}=

*C*

_{B∗A∗}(

*t*).