The calculation of transport properties in quantum liquids using the maximum entropy numerical analytic continuation method: Application to liquid parahydrogen
See allHide authors and affiliations

Contributed by Bruce J. Berne
Abstract
We present a method based on augmenting an exact relation between a frequencydependent diffusion constant and the imaginary time velocity autocorrelation function, combined with the maximum entropy numerical analytic continuation approach to study transport properties in quantum liquids. The method is applied to the case of liquid parahydrogen at two thermodynamic state points: a liquid near the triple point and a hightemperature liquid. Good agreement for the selfdiffusion constant and for the realtime velocity autocorrelation function is obtained in comparison to experimental measurements and other theoretical predictions. Improvement of the methodology and future applications are discussed.
One of the major goals of, and perhaps the most challenging problem in, computational statistical mechanics is the simulation of quantum dynamics in condensed phases. In principle, the density matrix formalism provides all the tools necessary to study equilibrium and timedependent properties of any chemical system. In practice, however, the exact solution of the timedependent quantum Wigner–Liouville equation is possible for a very limited class of simple systems, and the numerical solution for a general manybody system is not possible because of the well known phase cancellation problem (the sign problem).
This problem has led to a variety of different techniques to include the effects of quantum fluctuations on the dynamic response of the system. One of the viable alternatives to the exact quantum mechanical solution is the use of techniques that are “semiclassical” in nature; namely, the dynamic response is obtained with the aid of classical trajectories (1). Although such techniques appear promising, technical issues have prevented their use in describing dynamics in realistic quantum liquids.
Another class of methods that has been used with success in a variety of problems involves sophisticated numerical analytical continuation of exact imaginarytime pathintegral Monte Carlo (PIMC) data (2, 3). These methods have been applied to a variety of condensed phase problems, including the dynamics of an excess electron solvated in water (4), helium and xenon (5), vibrational relaxation (6, 7), optical spectroscopy (6–8), adiabatic reaction dynamics (9, 10), dynamics in various quantum lattice models (11, 12), and density fluctuations in superfluid helium (13). However, the application of these approaches to study density fluctuations (13) and transport properties (4) in quantum liquids has not been completely successful.
In this paper, we show that analytic continuation methods can be used successfully to study the transport properties of a “realistic” liquid. We express the imaginary time velocity autocorrelation function, which is obtained from a suitable PIMC method (14), in terms of a frequencydependent diffusion constant and use the maximum entropy method to analytically continue the imaginary time data to real time and thus obtain the selfdiffusion constant and the velocity autocorrelation function. We use the method to study transport properties in fluid parahydrogen at two thermodynamic state points: a liquid near the triple point at T = 14K and ρ = 0.0235Å^{−3}, and a hightemperature liquid at T = 25K and ρ = 0.0190Å^{−3}. The results are compared with experimental observations (15) and with results obtained from a quantum modecoupling theory (16, 17).
The paper is structured as follows: In Section I we express the velocity autocorrelation function and the selfdiffusion constant in a form suitable for the analytic continuation method. We also provide a working expression for the imaginary time velocity autocorrelation function amenable for PIMC simulation techniques. In Section II we describe the maximum entropy (MaxEnt) method (11, 18) used to perform the analytic continuation. In Section III we apply the method to study the transport of liquid parahydrogen. Section IV concludes.
I. Analytic Continuation of the Velocity Autocorrelation Function
In this section, we outline a convenient approach to obtain the velocity autocorrelation function and the selfdiffusion constant that is suitable for the analytic continuation method. We start from the well known Green–Kubo relation: 1 where the realtime velocity autocorrelation function is given by 2 (from now on, we set ℏ = 1), where Z = Tre^{−βH} is the partition function, β = 1/k_{b}T is the inverse temperature, and v is the velocity vector of a tagged particle in the liquid.
For the analytic continuation of the velocity autocorrelation function, it is convenient to define a frequencydependent diffusion constant, D(ω), which is given in terms of the power spectrum of C_{υ}(t) 3 The selfdiffusion constant is the zero frequency value of D(ω) divided by 6. The velocity autocorrelation function can be obtained by inverting the Fourier relation in Eq. 3: 4 The frequencydependent diffusion constant is analogous to the spectral density used in the analytic continuation of spectral line shapes (4, 5) and to the frequencydependent rate constant used in analytic continuation of the flux–flux correlation function (9). By performing the replacement t → −iτ and by using the detailed balance relation D(−ω) = e^{−βω}D(ω), we obtain 5 where t, τ ≥ 0, and 6 The reason for introducing the imaginary time velocity autocorrelation function is that, unlike its realtime counterpart, it is straightforward to obtain by using an appropriate PIMC simulation technique (19, 20). The approach we adopt here to obtain G_{υ}(τ) is based on the method developed in ref. 14. The result to lowest order in ɛ = β/P, where P is the number of Trotter slices, is 7 where N is the total number of liquid particles, r_{j} is a shorthand notation for the position vector of all liquid particles associated with bead j, r is the position vector of liquid particle α of bead j, and P(r_{1}, ⋯ , r_{P}) is the regular sampling function used in the standard cyclic PIMC method (with r_{0} = r_{P}).
To obtain the frequencydependent diffusion constant and the realtime velocity autocorrelation function, one has to invert the integral Eq. 5 by using an analytic continuation method. A specific choice is the MaxEnt method, which is described in the next section.
II. MaxEnt Numerical Analytic Continuation
The MaxEnt analytic continuation method described below is identical to the one used in our work on the reactiveflux analytic continuation method and is outlined here for completeness. We seek to analytically continue the imaginary time velocity autocorrelation function given in Eq. 5. Because G_{υ}(τ) is analytic for 0 < τ < β, the analytic continuation is accomplished by inverting the integral Eq. 5 to obtain a solution for D(ω). The zero mode value D(0) would then correspond to the experimentally observable selfdiffusion constant. Because of the singular nature of the integration kernel, the inversion of Eq. 5 is an illposed problem. As a consequence, a direct approach to the inversion would lead to an uncontrollable amplification of the statistical noise in the data for G_{υ}(τ), resulting in an infinite number of solutions that satisfy Eq. 5. Clearly, in this case, little can be said about the realtime dynamics and the corresponding diffusion constants.
Recently, Bayesian ideas have been used to deal with the illposed nature of continuing the noisy imaginary time Monte Carlo data to real time (2, 3). One of the most widely used approaches is the MaxEnt method (3, 18). The method requires only that the transformation that relates the data and the solution be known. Furthermore, MaxEnt allows the inclusion of prior knowledge about the solution in a logically consistent fashion. As such, the method is well suited for solving illposed mathematical problems.
For the purpose of the MaxEnt approach, we rewrite the integral Eq. 5 8 In this notation, G(τ) ≡ G_{υ}(τ) is the data (in this case the imaginary time velocity autocorrelation function), K(τ, ω) = e^{−ωτ} + e^{(τ−β)ω} is the singular kernel, and A(ω) is the solution, referred to as the map, corresponding to D(ω). MaxEnt principles provide a way to choose the most probable solution consistent with the data through the methods of Bayesian inference. Typically, the data are known only at a discrete set of points {τ_{j}}, and we likewise seek the solution at a discrete set of points {ω_{k}}. The MaxEnt method selects the solution that maximizes the posterior probability or the probability of the solution A given a data set G. Using the Bayes theorem, one can show that (11, 18) the posterior probability is given by 9 Here, χ^{2} is the standard mean squared deviation from the data 10 where C_{jk} is the covariance matrix 11 with M being the number of measurements.
In Eq. 9, S is the information entropy, the form of which is axiomatically chosen to be 12 In this formulation, the entropy is measured relative to a default model m(ω), which can contain prior information about the solution, and α is a positive regularization parameter.
Obtaining the MaxEnt solution then involves finding a map A that maximizes the posterior probability and is therefore a maximization problem in N variables, where N is the number of points {ω_{k}} at which the solution is evaluated. The solution obtained in this way is still conditional on the arbitrary parameter α, which can be interpreted as a regularization parameter controlling the smoothness of the map. Large values of α lead to a result primarily determined by the entropy function and hence the default model. Small α, in turn, lead to a map determined mostly by the χ^{2} and thus to a closer fitting of the data. The principal drawback is that, along with the data, the errors would be fit as well.
In this study, we use a flat default map [m(ω)], which satisfies a known sum rule, such as the integral over D(ω), and α is selected according to the Lcurve method (21). In this context, we regard α as a regularization parameter controlling the degree of smoothness of the solution and entropy as the regularizing function. The value of α is selected by constructing a plot of log[−S(A)] vs. log χ^{2}. This curve has a characteristic Lshape, and the corner of the L, or the point of maximum curvature, corresponds to the value of α that is the best compromise between fitting the data and obtaining a smooth solution.
We use a maximization algorithm due to Bryan (22), which reduces the space in which the search for the solution is performed. The kernel is first factored by using singular value decomposition K = VΣU^{T}. The singular nature of the kernel ensures only a small number of eigenvalues of Σ will be nonsingular. Because the space spanned by the rows of K is the same as that spanned by the columns of U associated with nonsingular eigenvalues, the search for the solution can be performed in this singular space of dimensionality N_{s}, where N_{s} is the number of nonsingular eigenvalues. The solution in singular space is expressed in terms of the vector u, which is related to the N dimensional map space via 13 This exponential transformation is useful, because it ensures the positivity of the solution.
III. Application to SelfDiffusion of Liquid ParaHydrogen
In this section, we study the transport of liquid parahydrogen by using the above outlined analytic continuation method. Although it is known that liquid parahydrogen may be treated as a Boltzmann particle near its triple point (23), it still exhibits some of the hallmarks of a highly quantum liquid. In fact, recent neutron scattering experiments of Bermejo et al. (24) have uncovered the existence of collective excitations that are absent in the classical fluid. These quantum excitations are a precursor of some of the collective modes that exist in the superfluid state. The model potential we use to study liquid parahydrogen is based on the Silvera–Goldman potential (25, 26), where the entire H_{2} molecule is described as a spherical particle, so the potential depends only on the radial distance between particles. This potential has been used to study thermodynamic properties and phase equilibrium of fluid hydrogen (27, 28) and has also been used to study transport (17, 29, 30) and density fluctuations (16, 31, 32) for liquid parahydrogen. The Silvera–Goldman potential is given by 14 where the first term on the righthand side (RHS) accounts for shortrange repulsive interactions, the second set of terms on the RHS account for longrange attractive dispersion interactions, and the last term on the RHS is an effective threebody correction (25). The last two terms are multiplied by a damping function, which turns off these interactions at short distances and is given by 15 where θ(r) is the Heaviside function (step function). The parameters for the potential are given in Table 1.
To obtain the imaginary time velocity autocorrelation function required for the analytic continuation to real time, we have performed PIMC simulations at two thermodynamic state points: a liquid near the triple point at T = 14K and ρ = 0.0235Å^{−3}, and a hightemperature liquid at T = 25K and ρ = 0.0190Å^{−3}. The density for both state points was chosen to be the average density under zero pressure (27). The PIMC simulations were done by using the Canonical ensemble with N = 108 particles interacting via the Silvera–Goldman potential defined in Eq. 14. The staging algorithm (33) for Monte Carlo chain moves was used to compute the imaginary time velocity autocorrelation function. The number of Trotter slices was P = 50 and P = 28 for the low and hightemperatures, respectively. This choice ensures that structural properties are well converged (34). The Monte Carlo moves (1 × 10^{6}) were made, with an acceptance ratio of approximately 0.35 for both state points.
The results for the imaginary time velocity autocorrelation function are shown in Fig. 1 for the aforementioned two state points. We find that the statistical error is quite small, indicating the high precision of the PIMC simulation method used in this work. The high quality of the imaginary time data is necessary for the numerical analytic continuation to the realtime axis, as will be illustrated below. The integral over the imaginary time velocity autocorrelation function in the range [0:βℏ] divided by βℏ yields the Kubo transform (35) of the initial time velocity autocorrelation function. It can be shown that the Kubo transform of the velocity autocorrelation function at the origin of time evaluates to k_{B}T/m, i.e., the classical result. We have computed the Kubo transform of the initial time velocity autocorrelation function for both state points. Our numerical results are within 0.001% of the exact analytic result, again indicating the high quality of the imaginary time data.
The imaginary time velocity autocorrelation functions obtained from the PIMC simulation data were then used as input data for the MaxEnt numerical analytic continuation procedure. The covariance matrices required by the MaxEnt procedure were computed by block averaging the Monte Carlo data. Unlike the case found in analytic continuation of the reactiveflux (9), the covariance matrices for the velocity autocorrelation function are not block diagonal, and thus the proper procedure must be used to decorrelate the statistical noise of each data point (2). The MaxEnt procedure was then used to determine the frequencydependent diffusion constant corresponding to each state point, by inverting Eq. 5. As mentioned in Section II, we have used the Lcurve method to determine the regularization parameter, α. Given the high quality of imaginary time data, a plot of log[−S(A)] vs. log χ^{2} results in a Lshape curve with a very sharp corner, corresponding to the choice of the regularization parameter α. We note that the results shown below are not sensitive to the exact choice of α over a relatively wide range.
In Fig. 2, we plot the frequencydependent diffusion constant for both state points studied in this work. The results obtained from the MaxEnt analytic continuation method are compared with the recent results obtained from a selfconsistent quantum modecoupling theory (16, 17). The agreement between the MaxEnt analytic continuation method for the lowtemperature state point is remarkable over a very wide range of frequencies. In particular, the position and width of the peak in D(ω) and the selfdiffusion constant given by the zero value of D(ω) are in excellent agreement at T = 14K, ρ = 0.0235Å^{−3}. The agreement between the MaxEnt method and the quantum modecoupling theory at the higher temperature point is good, although the position of the peak obtained from the MaxEnt method is slightly shifted to higher frequencies. Given that the quantum modecoupling theory is best suited for dense liquids, such as those near the triple point (cf. T = 14K, ρ = 0.0235Å^{−3} for liquid parahydrogen), the discrepancies between the two methods might indicate the breakdown of the quantum hydrodynamic approach at the higher temperature state point.
As mentioned above, the selfdiffusion constant of liquid parahydrogen can be obtained from the zero frequency value of D(ω). The values of the selfdiffusion constants obtained from the MaxEnt analytic continuation method are 0.28Å^{2}ps^{−1} and 1.47Å^{2}ps^{−1} for T = 14K, ρ = 0.0235Å^{−3} and T = 25K, ρ = 0.0190Å^{−3}, respectively. These results are in good agreement with the experimental results (0.4Å^{2}ps^{−1} and 1.6Å^{2}ps^{−1}) (15), with the results obtained from the quantum modecoupling theory (0.30Å^{2}ps^{−1} and 1.69Å^{2}ps^{−1}) (17), and from the full Constrained Molecular Dynamics (CMD) method (0.32Å^{2}ps^{−1} and 1.54Å^{2}ps^{−1}) (36). The agreement obtained for the selfdiffusion constant between the different methods suggests that the observed deviations from the experimental results is mainly affected by the accuracy of the Silvera–Goldman interaction potential (Eq. 14) and not by the MaxEnt analytic continuation method.
The realtime velocity autocorrelation function period can be obtained directly from the frequencydependent diffusion constant by using the relation defined in Eq. 4. In Fig. 3, we compare the normalized velocity autocorrelation function obtained from the MaxEnt analytic continuation method to the velocity autocorrelation function obtained from a quantum modecoupling theory (17). In view of the above discussion for the frequencydependent diffusion constant, it is hardly surprising that the agreement between the two methods is better at the lower temperature state point, where we obtain quantitative agreement for the shorttime decay, for the position of the first minimum in C_{υ}(t) and for the overall decay rate. The best agreement between the two methods is obtained at short times, which is expected because the quantum modecoupling theory is exact to order t^{6}, and the statistical errors in the MaxEnt analytic continuation method are small at short times. The small deviations between the two methods at longer times may result from increasing statistical errors in the MaxEnt method or from the approximations introduced in the quantum modecoupling theory. However, the overall good agreement between the two methods is a strong indication for the robustness and accuracy of both methods.
IV. Conclusion
In this paper, we have presented a method to study transport properties in highly quantum liquids. We expressed the imaginary time velocity autocorrelation function, which was obtained from a suitable PIMC method (14) in terms of a frequencydependent diffusion constant and used the MaxEnt method to obtain the selfdiffusion constant and the velocity autocorrelation function by analytically continuing the imaginary time data to real time.
The accuracy of the method was tested for liquid parahydrogen at two thermodynamic state points. As far as we know, this is the first successful application of an analytic continuation method to the study of transport in a realistic liquid. We find that the selfdiffusion constants are in good agreement with the experimental results for both thermodynamic state points. Furthermore, the agreement with transport coefficients obtained from the quantum modecoupling theory (17) is excellent, indicating that any discrepancy found between the MaxEnt method and the experiments results from the approximated pairpotential used and not from the dynamical method itself. We have also calculated the realtime velocity autocorrelation function and obtained excellent agreement with the results obtained from the quantum modecoupling theory.
Although detailed comparison between the MaxEnt analytic continuation approach and other methods is not the major goal of the present work, it should be noted that the approach taken here has some very attractive advantages. First, the method presented in this paper is systematic and accurate within the noise level of the numerical imaginary time input, and uncontrolled approximations that are typically made in other methods on the basis of various semiclassical and mixed quantum/classical techniques are not necessary here. Second, situations where the static distribution is not described by Boltzmann statistics can easily be handled within the present framework, because the additional complication of proper particle statistics may be absorbed into the PIMC calculation of the imaginary time input. Last, there are numerous possible improvements for our method. On the one hand, more efficient sampling techniques are needed to reduce the noise level in the PIMC simulations. On the other hand, the numerical analytic continuation method used in this work is a very basic implementation of the maximum entropy method. We believe that considerable improvement to the numerical analytic continuation procedure can be achieved through a better utilization of the MaxEnt procedure. For example, rather than using a flat default model, one could use a more informative one. Such a model could be obtained from approximate methods, such as a quantum modecoupling approach (16, 17). Evaluating the sensitivity of the solution to the default model would lead to an increased confidence in its validity. In addition, recently it was shown that combining short realtime dynamical information with the imaginary time data can significantly improve the quality of the analytically continued results (10, 37, 38). All of these issues may be the subject of future investigation.
Acknowledgments
This work was supported by The Israel Science Foundation, founded by the Israel Academy of Sciences and Humanities (Grants 9060/99 and 34/00 to E.R.), by the Research Corporation Innovation Award (Grant RIO642 to D.R.R.), and by a grant to B.J.B. from the National Science Foundation. E.R. acknowledges financial support from the Israeli Council of Higher Education (Alon Fellowship).
Abbreviations
 PIMC,
 pathintegral Monte Carlo;
 MaxEnt,
 maximum entropy
 Accepted October 11, 2001.
 Copyright © 2002, The National Academy of Sciences
References
 ↵
 Miller W H
 ↵
 Jarrell M,
 Gubernatis J E
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Sim E,
 Krilov G,
 Berne B J
 ↵
 ↵
 ↵
 ↵
 Rabani E,
 Reichman D R
 ↵
 Esel'son B N,
 Blagoi Y P,
 Grigor'ev V V,
 Manzhelii V G,
 Mikhailenko S A,
 Neklyudov N P
 ↵
Rabani, E. & Reichman, D. R. (2002) Phys. Rev. E65, in press.
 ↵
 ↵
 Skilling J
 ↵
 Berne B J
 ↵
 ↵
 Lawson C L,
 Hanson R J
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Scharf D,
 Martyna G,
 Klein M L
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Kubo R,
 Toda M,
 Hashitsume N
 ↵
 ↵
 ↵