{'fraction_non_alphanumeric': 0.04545808359928388, 'fraction_numerical': 0.010352611504631432, 'mean_word_length': 4.530779164873009, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the approximation of Feynman-Kac path integrals for quantum statistical mechanics

Discretizations of the Feynman-Kac path integral representation of the quantum mechanical density matrix are investigated. Each infinite-dimensional path integral is approximated by a Riemann integral over a finite-dimensional function space, by restricting the integration to a subspace of all admissible paths. Using this process, a wide class of methods can be derived, with each method corresponding to a different choice for the approximating subspace. The traditional ``short-time'' approximation and ``Fourier discretization'' can be recovered from this approach, using linear and spectral basis functions respectively. As an illustration, a novel method is formulated using cubic elements and is shown to have improved convergence properties when applied to a simple model problem.

The symbol D [x (τ )] indicates that the integration is performed over the set of all differentiable curves, x : [0, βh] → R d , with x (0) = a and x (βh) = b. The integer d reflects the dimensionality, with d = 3N for a system of N -particles in 3-dimensional space. The functional Φ can be derived from the classical action by introducing a relationship between temperature and imaginary time (it = βh) [1]. In this Letter, we will restrict our attention to the quantum many-body system, for which Φ takes the following form: Calculating the path integral in (1) is a challenging task, which in general cannot be performed analytically. It is only for simple model problems, such as quadratic potentials that an exact solution can be obtained. For more challenging systems, the path integral has traditionally been estimated using either the "short-time" approximation (STA) [3] or "Fourier discretization" (FD) [4,5]. Many authors have proposed improvements to the standard STA and FD, using techniques such as improved estimators [6], partial averaging [7][8][9], higher-order exponential splittings [10], advanced reference potentials [11], semi-classical expansions [12], and extrapolation [13]. The fundamental approach is the same in all of these methods: the path integral is reduced to a high (but finite) dimensional Riemann integral, which is approximated using either a Monte Carlo or Molecular Dynamics.
The aim of this Letter is to provide a framework for the formulation of a wide class of methods for the discretization of quantum mechanical path integrals. The idea of approximating path integrals using a finite subset of basis functions has been suggested before in the literature. Davison was one of the first to consider the use of orthogonal function expansions in the representation of Feynman path integrals [14], although he did not explore truncating the expansion. In a related article on Wiener integration, Cameron proposed using a finite set of orthogonal basis functions, and investigated the convergence of Fourier (spectral) elements [15]. In this Letter, we do not require that the basis functions are orthogonal, allowing for the direct comparison of the STA and FD methods. Although other authors have explored fundamental connections between the STA and FD [16] methods, we are unaware of any comparison using the approach investigated here. In addition, our approach allows for the construction of new methods using general classes of orthogonal polynomials or finite elements. To illustrate the flexibility of this approach, we derive a new method, using compactly supported (Hermite) cubic splines (HCS), which is shown to exhibit improved efficiency when applied to model problems.
To illustrate how one can use a subspace approxima-tion to discretize the quantum density matrix in (1), we start by introducing a change of variables to simplify the boundary conditions and temperature dependence for each path integral: x (τ ) = a + (b − a) τ /βh + y (τ /βh).
Since the admissible paths, x, satisfy the boundary conditions x (0) = a and x (βh) = b, the reduced paths given by y, will satisfy Dirichlet boundary conditions, y (0) = y (1) = 0, independent of a, b, and β. Introducing this change of variables into (1), results in the following: Note that the ith component of each reduced path y, denoted by y i , is a real-valued function on the interval [0, 1], satisfying Dirichlet boundary conditions. For the systems considered in this article, we also require that the derivative of each y i is measurable (i.e., square-integrable).
Functions of this form are members of an infinite dimensional Sobolev space [17], defined by S 1 We proceed in the following manner to discretize (3): Consider a sequence of subspaces of increasing dimen- For convenience, let each subspace be defined as the span of a particular set of basis functions: V P = span {ψ 1 , · · · , ψ P }. Now, given a component function y i ∈ S 1 0 [0, 1], we can define its projection on V P uniquely by Using the projection y (P ) i (ξ) as an approximation of y i (ξ) reduces the infinite-dimensional path integral in (3) to a finite-dimensional Riemann integral over the coefficients, α i,k : Here, we have used simplified notation for the multidimensional integral, with dα ≡ k,i dα k,i . The constant J reflects the particular choice of variables, and can be readily calculated (as discussed later). The reader should note that (4) does not depend on the basis functions chosen to represent the approximating subspace. If both {ψ 1 , · · · , ψ P } and {ψ 1 , · · · ,ψ P } span V P , then there is an invertible linear transformation (i.e., change of variables) U such thatα = Uα.
To show in detail how subspace methods can be applied in practice, we consider the case of an N -body Hamiltonian system: Here, the coordinate and momentum operators are denoted by x i andp i respectively. For this system the functional Φ is given by (2), which when applied to the projected path, After expanding the τ -integrals, introducing a change of variables ξ = τ /βh, and using the boundary conditions of each ψ k , we have where α i ≡ [α 1,i · · · α P,i ] T . The "stiffness matrix", K ∈ R P ×P , has entries given by the inner-product K j,k = 1 0ψ j (ξ)ψ k (ξ) dξ . Substituting (5) into (4), we obtain a simplified expression for the approximate density matrix: For the Fourier case, one typically calculates J by requiring that the discretization be exact when applied to an ideal gas (i.e., V ≡ 0) [4,5,14]. Applying this same technique to a generic subspace method, and assuming that K is positive definite, one can solve for J in a straightforward manner: Before discussing particular choices for basis functions, we should mention that, in general, the one-dimensional ξ-integral in (6) cannot be performed analytically. This problem has been traditionally circumvented by using a discrete approximation, such as Gaussian quadrature [5,16]. For example, one can view the primitive STA as using the trapezoidal rule. If the quadrature scheme is of sufficiently high order its use will not reduce the asymptotic rate of convergence of the overall method. An optimal scheme must be efficient, since for nonlinear N -body systems evaluating V may be computationally expensive.
As mentioned above, the real benefit of using a general subspace approach is the flexibility afforded through the choice of basis functions. By considering a general class of pseudo-spectral or finite-element basis functions, a diverse group of discretizations can be constructed. Direct comparisons can be made between basis functions of varying smoothness and support. However, for brevity, we restrict our attention in this Letter to three different types of basis functions: linear, spectral, and cubic elements. Representative basis functions from each of these discretizations are shown in Figure 1. The traditional STA method can be constructed by considering polygonal paths, which can be represented by piecewise linear basis functions [17]. For a given number of linear segments, P +1, we can define an approximating subspace V P as the span of basis functions {ψ 1 , · · · , ψ P }, where each ψ k is defined by the following formula For this discretization, it is routine to show that the elements of the "stiffness" matrix, K ∈ R P ×P , can be determined by K i,j = −1δ i−1,j +2δ i,j −1δ i+1,j . In a similar manner, the FD method can be derived using the subspace approach by considering spectral basis functions of the form ψ k (ξ) = 1/k sin (k πξ). K is diagonal for this basis, with entries given by K i,j = π 2 /2 δ i,j .
A new method can be constructed by approximating the space of paths using piecewise (Hermite) cubic splines (HCS) [17]. Each spline is defined on an interval of width 2/P , with its shape uniquely determined by its function value and derivative at the ends of the interval. It is assumed here that P is an even integer. Each piecewise cubic path has a continuous derivative, and is described by linear combinations of the basis functions One can verify that the reduced path y (P ) (ξ) = α k ψ k (ξ) satisfies Dirichlet boundary conditions, and interpolates the interior grid points (2j/P, α j ) for integers 1 ≤ j < P/2. The derivative of the path at all the grid points is determined by the remaining P/2 + 1 coefficients, α k . Due to the compact support of the basis functions, the stiffness matrix is banded, with block structure: where the blocks are given by Note that the blocks are not all the same size, with K 3 of dimension (P/2 − 1) × (P/2 + 1). The determinant of K hcs may be calculated exactly, but for most purposes it is enough to know that it is a constant, which will cancel out when (6) is used to calculate averages. As a numerical experiment, we apply each path integral discretization to the problem of calculating the average energy of a particle in a one-dimensional double-well. We have chosen the same double-well potential considered in [5], which is as follows: V (x) = m ω 2 x 2 /2 + A/((x/a) 2 + 1). The parameter values are all in atomic units, with ω = 0.006, A = 0.009, a = 0.09, and m = 1836. At low temperatures, the energy is just above 0.006, which is below the barrier height of 0.009.
To measure the accuracy of each method, we compute the energy at a fixed temperature of T = 0.1hω/k, using Metropolis Monte Carlo to generate the canonically distributed configurations. The one-dimensional lineintegrals of the potential are approximated using Simpson's rule for the FD and HCS methods, and the traditional trapezoidal rule for the STA method. The number of integration nodes is set equal to the number of basis functions, P , resulting in the same number of potential evaluations for each method. For the STA method this results the potential is evaluated at the end points of each polygonal segment (consistent with its traditional implementation).
It has been previously observed that averaged quantities (such as energy) converge at different rates, depending on the system, reference potential, and the form of the estimator [5,6,18]. We use a virial estimator of the energy [6], E = V (x) + xV ′ (x)/2 , which is known to exhibit improved convergence properties in many problems. The accuracy of each average is determined by comparing with the "exact" solution, computed by summing over the 15 lowest energy levels as calculated with Numerov's method [19]. In Fig. 2, the error in the computed energy is shown as a function of (a) the number of basis functions and (b) normalized CPU time. When the number of basis functions (or potential evaluations) is used as a measure of the work, we find that the FD and HCS methods are comparable, and both are more efficient than the STA method. However, when compared on the basis of CPU time, the HCS method is dramatically more efficient than both other methods. The inefficiency of the FD method for low-dimensional problems can be explained by considering the work required to compute P points on the path. This work scales like O(P 2 ) for the FD method, since the spectral basis functions are not compactly supported. On the other hand, for the STA and HCS methods this cost scales linearly with P . Although for very high dimensional problems, the cost of evaluating the potential should dominate, and we expect that the differences in computational cost would not be as pronounced.
In summary, the problem of approximating Feynman-Kac path integrals can be addressed using the finite-dimensional subspace approach. This technique allows for the formulation of new methods through the choice of a suitable set of basis functions. In addition, traditional methods such as the short-time approximation and Fourier discretization methods can be compared using this framework. As an example, by considering (Hermite) cubic splines, a new method can be constructed which exhibits improved efficiency when applied to a onedimensional double-well problem.