text:
# Random Series and Discrete Path Integral methods: The Lévy-Ciesielski implementation

## I Introduction

Ever since their introduction over fifty years ago, path integral methods have been an intense research field for physicists, chemists, and mathematicians alike, even if these researchers have sometimes used arguments of a rather different nature. The field began with Feynman’s observation Fey48  that the time propagator of the Schrödinger equation can be represented as a “sum over histories,” effectively giving a formula for the propagator as a limit of integrals over spaces of increasing dimension Fey65 . In mathematical terms, the existence of this limit is problematic though several research directions are known Alb76 ; Joh00 ; Car95 ; deF91 .

In a significant development, Kaç noticed that the “imaginary time” version of the formula utilized by Feynman has a definite probabilistic sense and could be interpreted as an integral of a functional of the seemingly ubiquitous Brownian motion Kac51 . Such a formula could represent the Green’s function for a certain class of diffusion processes, as for instance the density matrix for the Bloch equation. The end product of their work is beautifully summarized by what is now called the Feynman-Kaç formula Sim79

$$\frac{\rho (x,x^{};\beta )}{\rho _{fp}(x,x^{};\beta )}=𝔼\mathrm{exp}\left\{\beta _0^1V\left[x_r(u)+\sigma B_u^0\right]du\right\},$$
(1)
where $`\rho (x,x^{};\beta )`$ is the density matrix for a monodimensional canonical system characterized by the inverse temperature $`\beta =1/(k_BT)`$ and made up of identical particles of mass $`m_0`$ moving in the potential $`V(x)`$. The stochastic element that appears in Eq. (1), $`\{B_u^0,u0\}`$, is a so-called standard Brownian bridge defined as follows: if $`\{B_u,u0\}`$ is a standard Brownian motion starting at zero, then the Brownian bridge is the stochastic process $`\{B_u|B_1=0,\mathrm{\hspace{0.17em}0}u1\}`$ i.e., a Brownian motion conditioned on $`B_1=0`$ Dur96 . In this paper, we shall reserve the symbol $`𝔼`$ to denote the expected value (average value) of a certain random variable against the underlying probability measure of the Brownian bridge $`B_u^0`$. To complete the description of Eq (1), we set $`x_r(u)=x+(x^{}x)u`$ (called the reference path), $`\sigma =(\mathrm{}^2\beta /m_0)^{1/2}`$, and let $`\rho _{fp}(x,x^{};\beta )`$ denote the density matrix for a similar free particle.

Rather than directly employing Eq. (1), chemical physicist’s arguments are usually constructed around the Trotter composition rule Tro59  that exploits the fact that $`\{e^{\beta H};\beta >0\}`$ is a semigroup of operators on $`L^2()`$, so that

$$e^{(\beta _1+\beta _2)H}=e^{\beta _1H}e^{\beta _2H}$$
(2)
or, in coordinate representation,

$$x|e^{(\beta _1+\beta _2)H}|x^{}=_{}dzx|e^{\beta _1H}|zz|e^{\beta _2H}|x^{}.$$
(3)
By writing $`\beta =_{k=1}^n\beta _k`$, repeatedly applying the Trotter rule and choosing an adequate “short-time” approximation, one ends up with a sequence of integrals on spaces of increasing dimension, converging to the density matrix as $`\mathrm{max}_{1kn}\beta _k0`$. Of course, this is much in the spirit of the original Feynman path integral approach. The methods deduced by this technique are usually called Discrete Path Integral (DPI) methods (see Ref. Mie01,  and the cited bibliography).

It has become apparent that the Kaç interpretation of the Feynman’s formula may, in fact, offer a valuable starting point for the general construction finite dimensional approximations to the density matrix. This is so because the Brownian motion is a mathematically well understood object, for which various constructions are known. For example, the use of the Ito-Nisio theorem Kwa92  has lead the present authors to the development of the random series path integral methods in a surprisingly general fashion Pre02 . This generality of the theory allows one to identify optimality criteria and eventually to answer questions as what the best representation is or how to modify the approach in order to improve the convergence.

In this paper, we shall look at the relation between the Kaç interpretation of the Feynman formula and discrete path integral methods. In the first part, we consider what we call the standard DPI methods. We shall again show the strength of the Kaç approach, at least in terms of generality and mathematical interpretation of the formulae. In the second part, we explore the connection between the random series technique, as particularized for the Lévy-Ciesielski series representation of the Brownian bridge, and certain DPI results from the chemical literature. While not the primary goal, we do obtain in an effortless manner a better way of implementing the latter results that features a “built in” fast sine-Fourier transform. We shall also derive the two basic modifications of the Lévy-Ciesielski path integral method (LCPI), the partial averaging Dol85 ; Pre02  and the reweighted techniques Pre02 , and establish their asymptotic law of convergence. We suggest that these results again emphasize the power of the Kaç interpretation of the Feynman formula. By providing a central framework for discussion and analysis, the Kaç approach significantly aids in characterizing the various methods and in establishing their interconnections, links that otherwise would be obscured by the multitude of possible representations.

## II The standard Discrete Path Integral method

### II.1 Trotter-Suzuki approach.

Our definition of the standard DPI method has to do with the particular short-time approximation that replaces the exact one in the Trotter product

$$e^{\beta H}=\underset{n+1\text{terms}}{\underset{}{e^{\frac{\beta }{n+1}H}\mathrm{}e^{\frac{\beta }{n+1}H}}}.$$
(4)
We follow closely the arguments of Suzuki Suz91 . The Hamiltonian of the system can be written as a sum between the kinetic operator and the potential operator in the form $`H=K+V`$. The coordinate representations for the two operators are analytically known to be

$$x|e^{\beta K}|x^{}=\rho _{fp}(x,x^{};\beta )$$

and

$$x|e^{\beta V}|x^{}=e^{\beta V(x)}\delta (x^{}x),$$

respectively. It is therefore natural to consider short-time approximations that can be expressed by a finite composition of the above density matrices. The simplest example is the $`2`$-term splitting formula

$$e^{\beta (K+V)}=e^{\beta K}e^{\beta V}[1+𝒪(\beta ^2)],$$
(5)
which is of order $`1`$. More generally, the order of a splitting formula is said to be $`k`$ if the relative error is $`𝒪(\beta ^{k+1})`$. The motivation for this is that if

$$e^{\beta (K+V)}=f_k(K,V;\beta )[1+𝒪(\beta ^{k+1})],$$

then

$$e^{\beta (K+V)}=\left[f_k(K,V;\frac{\beta }{n})\right]^n\left[1+𝒪\left(\frac{\beta ^{k+1}}{n^k}\right)\right]$$
(6)
i.e., the error of the final $`n`$-term Trotter product formula decays as fast as $`1/n^k`$. The relation (6) was actually proved by Suzuki Suz85  in terms of operator norms for bounded operators $`A`$ and $`B`$, but such an estimation also holds for $`K`$ and $`V`$ (which generally are unbounded operators; however, $`e^{\beta (K+V)}`$ is bounded for most of the potentials of physical interest and for all positive $`\beta `$).

A better splitting is offered by the $`3`$-term formula

$$e^{\beta (K+V)}=e^{\frac{1}{2}\beta V}e^{\beta K}e^{\frac{1}{2}\beta V}[1+𝒪(\beta ^3)],$$
(7)
or the one obtained by permuting $`V`$ with $`K`$. These are of order $`2`$ and go by the name of symmetrical trapezoidal Trotter short-time approximations Rae83 ; Suz91 . More generally, let us define a $`2l+1`$-term splitting formula of order $`k`$ by the expression

$`e^{\beta (K+V)}=`$ $`e^{a_0\beta V}e^{b_1\beta K}e^{a_1\beta V}\mathrm{}`$ (8)
$`\mathrm{}e^{b_l\beta K}e^{a_l\beta V}[1+𝒪(\beta ^{k+1})].`$
Symmetry arguments suggest that for the optimum $`2l+1`$-term splitting formula the sequences $`a_0,\mathrm{},a_l`$ and $`b_1,\mathrm{},b_l`$ should be *palindromic* i.e., if the coefficients are read left to right, they form the same numerical sequence as if they are read right to left. A look at the trapezoidal Trotter formula shows that this condition is natural, as one has little reason to believe that anything new can be achieved by considering some arbitrary $`e^{a\beta V}e^{\beta K}e^{(1a)\beta V}`$ decomposition. In fact, with the help of the Campbell-Baker-Haussdorf-Dynkin formula Var88 , it can be shown that this more general expression is an order $`2`$ splitting only if $`a=1/2`$ and that it is an order $`1`$ splitting otherwise. More generally, since the operator $`e^{\beta H}`$ is Hermitian, it is natural and, as argued by De Raedt and De Raedt Rae83 , optimal to approximate it by a sequence of Hermitian operators. It is straightforward to see that the $`n`$-order Trotter product (4) is Hermitian if and only if the short-time approximation (8) is Hermitian. In turn, this requires the palindromicity of the $`\{a_i\}`$ and $`\{b_i\}`$ sequences.

It is not difficult to see that if $`_ia_i=p`$ and $`_ib_i=q`$, then the $`n`$-order Trotter formula (4) converges to $`\mathrm{exp}[\beta (qK+pV)]`$. On the other hand, the equality

$$e^{\beta (qK+pV)}=e^{\beta (K+V)}$$

holds for arbitrary potentials $`V(x)`$ if and only if $`q=1`$ and $`p=1`$. Therefore, the additional constraints $`_ia_i=1`$ and $`_ib_i=1`$ must be enforced upon the sequences $`\{a_i\}`$ and $`\{b_i\}`$.

We considered this more general problem with the hope that by using a more advanced splitting one may improve the asymptotic order of the Trotter product formula. Now, there is one more restriction that we have to place on the sequence $`a_0,b_1,a_1,\mathrm{}`$, namely it should be made up of *real and positive* numbers only. Otherwise, the short-time approximations are either ill-defined or, by Trotter composition, generate algorithms that are numerically unstable at low temperatures. Unfortunately, the following theorem of Suzuki (see Theorem 3 of Ref. Suz91, ) says that

###### Theorem 1 (Suzuki nonexistance theorem)

There are no finite length splitting formulae (8) of order $`3`$ or more such that the coefficients $`a_0,b_1,a_1,\mathrm{}`$ are all real and positive.

This negative result shows that more general splitting formulas do not produce short-time approximations capable of improving upon the trapezoidal Trotter result, at least as far as the asymptotic order of the Trotter product rule is concerned. However, the product rule (4), which uses equally spaced time slices, does not provide the most general standard DPI expression. In the next section, we shall argue that this most general expression is of the form given by Eq. (8), for which the Suzuki nonexistance theorem does not apply.

### II.2 Direct quadrature of the Feynman-Kaç formula.

Let us notice that the form of the equation (8) is invariant under the Trotter composition (4) and so it can be regarded as the most general standard DPI approximation to the density matrix provided that we can give a recipe for choosing the sequences $`a_0,b_1,a_1,\mathrm{},b_l,a_l`$ in such a way that the correct result is recovered in the limit $`l\mathrm{}`$. While in the Trotter-Suzuki approach this may seem a daunting task, the problem has an easy solution by means of the Kaç interpretation of the Feynman formula. In this section, we shall derive a more general expression for the standard Discrete Path Integral method simply by replacing the monodimensional integral over $`u`$ in Eq. (1) with an approximate quadrature sum and then using the definition of the Brownian bridge to compute the expectation of the resulting functional. Given the Suzuki nonexistence theorem, it is hard to believe that one may eventually devise a standard DPI method with asymptotic convergence $`𝒪(1/n^3)`$ or better. However, before one starts to investigate the validity of this conjecture, one needs a more general statement of the standard DPI method.

For obvious reasons, the random process $`W_{x,x^{}}^\sigma (u)=x_r(u)+\sigma B_u^0`$ is called a Brownian bridge of variance $`\sigma ^2`$ and end points $`(x,x^{})`$. The Feynman-Kaç formula can be expressed in terms of the new process in the form

$$\frac{\rho (x,x^{};\beta )}{\rho _{fp}(x,x^{};\beta )}=𝔼\mathrm{exp}\left\{\beta _0^1V\left[W_{x,x^{}}^\sigma (u)\right]du\right\}.$$
(9)
A quite important property of the Brownian bridge $`W_{x,x^{}}^\sigma (u)`$ is the joint distribution of the variables $`W_{x,x^{}}^\sigma (u_1),\mathrm{},W_{x,x^{}}^\sigma (u_n)`$ for a given partitioning $`0<u_1<\mathrm{}<u_n<1`$ of the interval $`[0,1]`$. Let us set

$$p_t(x)=\frac{1}{\sqrt{2\pi t}}e^{x^2/2t}$$

and notice that $`\rho _{fp}(x,x^{};\beta )=p_{\sigma ^2}(x^{}x).`$ From the very definition of the Brownian motion Dur96a , the aforementioned joint distribution can be straightforwardly shown to be

$`P\left\{W_{x,x^{}}^\sigma (u_1)[x_1,x_1+\mathrm{d}x_1],\mathrm{},W_{x,x^{}}^\sigma (u_n)[x_n,x_n+\mathrm{d}x_n]\right\}`$
$`=p_{\sigma ^2u_1}(x_1x)p_{\sigma ^2(u_2u_1)}(x_2x_1)\mathrm{}p_{\sigma ^2(u_nu_{n1})}(x_nx_{n1})`$ (10)
$`\times p_{\sigma ^2(1u_n)}(x^{}x_n)/p_{\sigma ^2}(x^{}x)\mathrm{d}x_1\mathrm{}\mathrm{d}x_n.`$
One may use the above joint distribution density to compute the expectations of the functionals of the Brownian bridge which are of the form

$`p_{\sigma ^2}(x^{}x)𝔼\{f[W_{x,x^{}}^\sigma (u_1),\mathrm{},W_{x,x^{}}^\sigma (u_n)]\}={\displaystyle _{}}dx_1\mathrm{}{\displaystyle _{}}dx_nf(x_1,\mathrm{},x_n)`$
$`\times p_{\sigma ^2u_1}(x_1x)p_{\sigma ^2(u_2u_1)}(x_2x_1)\mathrm{}p_{\sigma ^2(u_nu_{n1})}(x_nx_{n1})p_{\sigma ^2(1u_n)}(x^{}x_n),`$ (11)
where $`f(x_1,\mathrm{},x_n)`$ is some integrable $`n`$-dimensional function. As a direct application of Eq. (II.2), consider a quadrature scheme on the interval $`[0,1]`$ specified by the points $`0=u_0<u_1<\mathrm{}<u_n<u_{n+1}=1`$ and the corresponding *nonnegative* weights $`w_0,w_1,\mathrm{},w_{n+1}`$. Replacing the monodimensional integral in the Feynman-Kaç formula (9) by its quadrature form, we obtain an approximation to the density matrix of the form

$$\frac{\rho _n^{\text{DPI}}(x,x^{};\beta )}{p_\sigma (x^{}x)}=𝔼\mathrm{exp}\left\{\beta \underset{i=0}{\overset{n+1}{}}w_iV\left[W_{x,x^{}}^\sigma (u_i)\right]\mathrm{d}u\right\}.$$
(12)
The expectation value of this formula can be *exactly* reduced to a finite dimensional integral with the help of the formula (II.2). We call Eq. (12) the standard Discrete Path Integral (DPI) method and we expect it to converge to the correct result for all continuous and bounded from below potentials $`V(x)`$. In this respect, remember that with probability one the Brownian paths are continuous and therefore so is $`V[W_{x,x^{}}^\sigma (u)]`$ as a function of $`u`$. Also remember that by definition a quadrature scheme on $`[0,1]`$ is constructed so that it eventually integrates all continuous functions on $`[0,1]`$.

Formula (12) can indeed be formally deduced starting with the Trotter composition rule (2) and a carefully chosen sequence of short-time approximations. More precisely, for $`i=0,1\mathrm{},n`$, define $`\theta _i=u_{i+1}u_i`$. Then the equation (12) is nothing else but the Trotter product

$`e^{w_0\beta V}e^{\theta _0\beta K}e^{w_1\beta V}e^{\theta _1\beta K}\mathrm{}`$
$`\mathrm{}e^{w_n\beta V}e^{\theta _n\beta K}e^{w_{n+1}\beta V},`$ (13)
which is of course of the type given by the formula (8). Finally, let us notice that we always have $`_i\theta _i=1`$. We also require that $`_iw_i=1`$. In fact, one is not interested in relaxing these equalities because they are the necessary and sufficient conditions to obtain the exact free particle density matrix at all levels of approximation. The reader can directly verify this fact by assuming that the potential $`V(x)`$ in Eq. (II.2) is constant but not zero. Moreover, the additional restriction of the integration schemes to those for which the sequences $`\theta _i`$ and $`w_i`$ are palindromic is justified by the requirement that the approximate density matrices be Hermitian.

To summarize, the advantage of the equation (12) is the novel interpretation for the sequences $`\theta _i=u_{i+1}u_i`$ and $`w_i`$, leading us to a more general convergence problem: what is the best convergence order for the standard DPI approach (12) and for what types of quadrature schemes is it attained? As we suggested in the beginning of this subsection, it is very plausible that the answer to the above question is two and is attained for almost all “sensible” quadrature schemes. We illustrate this by studying the convergence of the diagonal matrix element $`\rho (0;\beta )=0|e^{\beta H}|0`$ of an harmonic oscillator for the following quadrature techniques: the trapezoidal rule (TT) and the Gauss-Legendre method (GL) Pre92 . In both cases the condition $`_iw_i=1`$ and the palindromicity of the sequences $`\theta _i`$ and $`w_i`$ are respected. We leave it for the reader to show that if the trapezoidal rule is used for integration, then one recovers the classical trapezoidal Trotter formula.

If $`Mt`$ stands for any of the methods studied and if $`\alpha _{Mt}`$ represents the convergence order of the corresponding matrix element $`\rho _n^{Mt}(0;\beta )`$, then the convergence constant is defined by

$$c_{Mt}=\underset{n\mathrm{}}{lim}n^{\alpha _{Mt}}[\rho (0;\beta )\rho _n^{Mt}(0;\beta )].$$

The above relation can be cast in the more intuitive but equivalent form

$$\rho (0;\beta )\rho _n^{Mt}(0;\beta )+\frac{c_{Mt}}{n^{\alpha _{Mt}}},$$

with an appropriate definition of the symbol $``$. These convergence orders and convergence constants can be evaluated numerically as follows. For each method, we compute

$$\alpha _n^{Mt}=(n^21/4)\mathrm{log}\left[1+\frac{\rho _{4n2}^{Mt}(0;\beta )\rho _{4n+2}^{Mt}(0;\beta )}{\rho _{4n+2}^{Mt}(0;\beta )\rho (0;\beta )}\right],$$

where $`\rho _{4n+2}^{Mt}(0;\beta )`$ represents the DPI approximation of order $`4n+2`$ for the method $`Mt`$. The evaluation of the matrix elements $`\rho _{4n+2}^{Mt}(0;\beta )`$ is discussed in Appendix B. Then, as argued in Ref. (Pre02, ), $`\alpha _n^{Mt}`$ as a function of $`n`$ is asymptotically a straight line, whose slope gives the convergence order. Therefore, $`\alpha _{Mt}=lim_n\mathrm{}\alpha _{n+1}^{Mt}\alpha _n^{Mt}`$. As to the convergence constants, they can be evaluated by studying the asymptotic slopes of

$$c_n^{Mt}=(4n+2)^{\alpha _{Mt}}(n+1/2)\left[\rho _{4n+2}^{Mt}(0;\beta )\rho (0;\beta )\right],$$

once $`\alpha _{Mt}`$ is known. The computations were performed in atomic units for a particle of mass $`m_0=1`$ and for the harmonic oscillator $`V(x)=x^2/2`$. The inverse temperature was $`\beta =10`$.

As shown in Fig. 1, the asymptotic convergence order of both methods is 2. The convergence constants are found to be $`c_{TT}=0.103`$ and $`c_{GL}=0.127`$, respectively. One notices that at the temperature $`\beta =10`$, the trapezoidal Trotter method is slightly faster. However for $`\beta =1`$, one computes $`c_{TT}=0.033`$ and $`c_{GL}=0.005`$, which indicates that for this temperature the Gauss Legendre method is faster.

The conclusion we draw from this analysis is that “better” integration schemes do not necessarily improve upon the convergence of the standard DPI methods. Why is this so? Again the Kaç interpretation of the Feynman formula gives us an explanation which is not obvious from the Trotter composition rule. A famous theorem due to Paley, Wiener, and Zygmund Dur96b  says that with probability one the paths of the Brownian motion are continuous but not differentiable at every point. Therefore, $`V[x_r(u)+B_u^0]`$ as a function of $`u`$ is not differentiable even if the potential $`V`$ is. As emphasized by Press *et al.* Pre92 , higher order quadrature schemes do not automatically translate into better convergence, unless the integrand is well behaved. In our case, there is a limit upon the rate of convergence of the quadrature schemes which is set by the properties of the Brownian motion paths rather than by the properties of the potential, provided that the latter has a continuous first order derivative. In conclusion, one expects that there is an intrinsic limit for the convergence order of the standard DPI methods. Moreover, the Suzuki nonexistence theorem predicts that none of the classical quadrature formulas for equally spaced abscissas (e.g. Simpson’s rules, Bode’s rule, etc.) are going to improve upon the asymptotical convergence of the trapezoidal rule. This is strong evidence that the intrinsic limit for the convergence order of the standard DPI methods is 2.

### II.3 Kinetic energy diagonalization for the standard DPI technique

In practical applications, it is generally difficult to work directly with the formula (II.2) because this involves a correlated Gaussian multidimensional distribution. As shown by Butler and Friedman But55 , this correlated distribution can be replaced with an independent one by simple algebraic manipulations. Later, Coalson Coa86  used a similar technique in order to demonstrate, on an intuitive basis, the relation between the discrete and the Fourier path integral methods. The two approaches mentioned above are technically different and in fact there are an infinite number of such transformations. As we shall see in this section, they are related by simple orthogonal transformations and in Section III we shall propose a new approach, which allows for faster numerical implementations.

We begin by performing a coordinate transformation so as to diagonalize the positive definite quadratic form associated with the kinetic operator. More precisely, let us introduce the transformation of coordinates $`z_n=[x_nx_r(u_n)]/\sigma `$. By using the condition $`_i\theta _i=1`$, it is straightforward to show that the formula (II.2) becomes

$`𝔼\{f[W_{x,x^{}}^\sigma (u_1),\mathrm{},W_{x,x^{}}^\sigma (u_n)]\}={\displaystyle _{}}dz_1\mathrm{}{\displaystyle _{}}dz_nf[z_1\sigma +x_r(u_1),\mathrm{},z_n\sigma +x_r(u_n)]`$
$`\times p_{\theta _0}(z_1)p_{\theta _1}(z_2z_1)\mathrm{}p_{\theta _{n1}}(z_nz_{n1})p_{\theta _n}(z_n),`$
or, in an even more compact notation,

$`𝔼\{f[W_{x,x^{}}^\sigma (u_1),\mathrm{},W_{x,x^{}}^\sigma (u_n)]\}`$
$`={\displaystyle _{}}dz_1\mathrm{}{\displaystyle _{}}dz_n{\displaystyle \frac{1}{\sqrt{(2\pi )^ndet(A)^1}}}\mathrm{exp}\left({\displaystyle \frac{1}{2}}\overline{z}^TA\overline{z}\right)`$
$`\times f[z_1\sigma +x_r(u_1),\mathrm{},z_n\sigma +x_r(u_n)],`$
where the matrix $`A`$ is the $`n`$-dimensional tridiagonal matrix defined by $`A_{i,i}=1/\theta _i+1/\theta _{i1}`$ for $`1in`$ and $`A_{i,i+1}=A_{i+1,i}=1/\theta _i`$ for $`1in1`$.

By construction, the matrix $`A`$ is symmetric and positive definite \[otherwise, the integrability of $`\mathrm{exp}\left(\overline{z}^TA\overline{z}/2\right)`$ would be violated\] and can be diagonalized by an orthogonal matrix $`S`$. Defining the new coordinates $`\overline{y}=S^T\overline{z}`$, and letting $`\{\lambda _i;1in\}`$ be the set of the $`n`$ real and (strictly) positive eigenvalues of $`A`$, we have

$`𝔼\{f[W_{x,x^{}}^\sigma (u_1),\mathrm{},W_{x,x^{}}^\sigma (u_n)]\}=\left[{\displaystyle \underset{i=1}{\overset{n}{}}}\left({\displaystyle \frac{\lambda _i}{2\pi }}\right)\right]^{1/2}`$
$`\times {\displaystyle _{}}\mathrm{d}y_1\mathrm{}{\displaystyle _{}}\mathrm{d}y_n\mathrm{exp}({\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{n}{}}}\lambda _iy_i^2)`$
$`\times f[x_r(u_1)+\sigma {\displaystyle \underset{j=1}{\overset{n}{}}}S_{1,j}y_j,\mathrm{},x_r(u_n)+\sigma {\displaystyle \underset{j=1}{\overset{n}{}}}S_{n,j}y_j].`$
Finally, setting $`a_i=\lambda _i^{1/2}y_i`$, one ends up with

$`𝔼\{f[W_{x,x^{}}^\sigma (u_1),\mathrm{},W_{x,x^{}}^\sigma (u_n)]\}={\displaystyle _{}}da_1\mathrm{}{\displaystyle _{}}da_n`$
$`\times \left(2\pi \right)^{n/2}\mathrm{exp}({\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{n}{}}}a_i^2)f[x_r(u_1)`$ (14)
$`+\sigma {\displaystyle \underset{j=1}{\overset{n}{}}}S_{1,j}{\displaystyle \frac{a_j}{\lambda _j^{1/2}}},\mathrm{},x_r(u_n)+\sigma {\displaystyle \underset{j=1}{\overset{n}{}}}S_{n,j}{\displaystyle \frac{a_j}{\lambda _j^{1/2}}}].`$
This formula is advantageous for numerical applications because the integration is performed over independent identically distributed Gaussian distributions. Given a quadrature scheme, one diagonalizes the tridiagonal matrix $`A`$ and tabulates the values of $`S_{i,j}`$ and $`\lambda _i`$. For the case of equally spaced time slices, the eigenvectors and eigenvalues of the matrix $`A`$ are known analytically:

$$S_{i,j}=\sqrt{\frac{2}{n+1}}\mathrm{sin}\left(\frac{ij\pi }{n+1}\right),1i,jn$$

and

$$\lambda _i=4(n+1)\mathrm{sin}^2\left[\frac{i\pi }{2(n+1)}\right],1in,$$

respectively.

Similar to the invariance of the Brownian bridge at a change of basis as shown by the Ito-Nisio theorem, the formula (II.3) is invariant to arbitrary orthogonal transformations of the vectors $`\overline{a}=(a_1,\mathrm{},a_n)`$. Indeed, let $`Q`$ be an arbitrary $`n`$-dimensional orthogonal matrix, and consider the coordinate transformation $`\overline{a}^{}=Q^T\overline{a}`$. Notice that $`_ia_i^2=_ia_{}^{}{}_{i}{}^{2}`$ and define the matrix

$$T_{i,j}=\underset{k=1}{\overset{n}{}}\frac{S_{i,k}Q_{k,j}}{\lambda _k^{1/2}}.$$
(15)
Then a little algebra shows that

$`𝔼\{f[W_{x,x^{}}^\sigma (u_1),\mathrm{},W_{x,x^{}}^\sigma (u_n)]\}={\displaystyle _{}}da_1\mathrm{}{\displaystyle _{}}da_n`$
$`\times \left(2\pi \right)^{n/2}\mathrm{exp}({\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{n}{}}}a_i^2)f[x_r(u_1)`$ (16)
$`+\sigma {\displaystyle \underset{j=1}{\overset{n}{}}}T_{1,j}a_j,\mathrm{},x_r(u_n)+\sigma {\displaystyle \underset{j=1}{\overset{n}{}}}T_{n,j}a_j].`$
Because of the additional degrees of freedom, the last formula is more useful in practical applications than the transformation (II.3). A good part of the computational time is spent with the evaluation of the current paths. For a monodimensional system, one usually needs a number of operations proportional to $`n^2`$ in order to compute the vector $`T\overline{a}`$ by matrix multiplication. However, if equally spaced time slices are used, the $`n`$ elements of the form $`_{j=1}^nS_{i,j}a_j/\lambda _j^{1/2}`$ from Eq. (II.3) can be computed by fast sine-Fourier transform in a number of operations proportional to $`n\mathrm{log}_2(n)`$, provided that $`n=2^k1`$ with $`k1`$ Mie01 ; Pre92a . Equivalently, one may say that there must be some orthogonal matrix $`Q`$ such that the associated matrix $`T`$ defined by the relation (15) is a sparse matrix with at most $`k`$ nonvanishing elements on any line. Therefore, the evaluation of the elements $`T\overline{a}`$ by direct matrix multiplication requires only $`𝒪(kn)`$ operations. In this paper, we shall directly find such a matrix $`T`$ by means of the Lévy-Ciesielski representation of the Brownian bridge, which is discussed in the next section \[see formula (III.2)\].

## III The Lévy-Ciesielski representation of the Feynman-Kaç formula

As we discussed in the previous section, the transformation (II.3) was utilized by Coalson in order to establish a connection between the discrete path integral methods and the Wiener-Fourier path integral technique Coa86 . However, strictly speaking the Wiener-Fourier sequence of approximations is not equivalent to any discretization scheme. That is, for any $`n`$, there is no sequence of short-time approximations which by Trotter composition would generate the $`n`$th order Wiener-Fourier approximation. A more precise statement of this assertion is given at the end of Section III.B. Then, a natural question arises: Is there any random series for which at least a particular subsequence can be thought of as a DPI method? The answer is positive and is furnished by the Lévy-Ciesielski random series construction of the Brownian bridge.

In this section, we shall specialize the general theory of the random series representation of the Feynman-Kaç formula Pre02  for the particular case of the Lévy-Ciesielski representation of the Brownian motion. The respective method will be designated by the acronym LCPI. We shall also derive the three associated methods: the primitive LCPI, the partial averaging LCPI, and the reweighted LCPI. Moreover, with the help of the Lévy-Ciesielski series representation, we shall prove the Trotter product rule for the case $`n=2^k1`$ and for this subsequence, we shall show that each of the above modifications of the LCPI method can be interpreted as the $`n`$-order Trotter product of some appropriate short-time approximations. In doing so, we establish a direct connection between the discrete and the random series path integral techniques. As a practical application, we shall obtain a sparse matrix $`T`$ of the form (15) which requires only $`𝒪(kn)`$ operations to compute the vector $`T\overline{a}`$ by matrix multiplication.

### III.1 The Lévy-Ciesielski path integral method.

Some of the arguments we use in the following introduction to the Lévy-Ciesielski representation of the Brownian bridge can be found in Ref. McK69, . For $`k=1,2,\mathrm{}`$ and $`j=1,2,\mathrm{},2^{k1}`$, the Haar function $`f_{k,j}`$ is defined by

$$f_{k,j}(t)=\{\begin{array}{cc}2^{(k1)/2},& t[(l1)/2^k,l/2^k]\\ 2^{(k1)/2},& t[l/2^k,(l+1)/2^k]\\ 0,& \text{elsewhere,}\end{array}$$
(17)
where $`l=2j1`$. Together with $`f_01`$, these functions make up a complete orthonormal basis in $`L^2([0,1])`$. Their primitives

$$F_{k,j}(t)=\{\begin{array}{cc}2^{(k1)/2}[t(l1)/2^k],& t[(l1)/2^k,l/2^k]\\ 2^{(k1)/2}[(l+1)/2^kt],& t[l/2^k,(l+1)/2^k]\\ 0,& \text{elsewhere}\end{array}$$
(18)
are called the *Schauder functions*. As McKean puts it McK69 , the Schauder functions are “little tents,” which can be obtained one from the other by dilatations and translations. In modern terminology, this has to do with the fact that the original Haar wavelet basis is a multiresolution analysis of $`L^2([0,1])`$ organized in “layers” indexed by $`k`$ Mal99 . If we disregard the factor $`2^{(k1)/2}`$, the Schauder functions make up a pyramidal structure as shown in Fig. 2.

Let $`\{a_{k,j};k=1,2,\mathrm{};j=1,2,\mathrm{},2^{k1}\}`$ be i.i.d. standard normal variables and define $`Y_0(u,\overline{a})0`$ and

$$Y_k(u,\overline{a})=\underset{j=1}{\overset{2^{k1}}{}}a_{k,j}F_{k,j}(u).$$

Then by Ito-Nisio theorem,

$$B_u^0(\overline{a})=\underset{k=1}{\overset{\mathrm{}}{}}Y_k(u,\overline{a})$$
(19)
is equal in distribution to a standard Brownian bridge and the convergence of the right hand side random series is uniform almost surely.

Let us now define the primitive, partial averaged, and reweighted LCPI methods, which are the standard techniques that can be derived from a series representation Pre02 . They will be denoted in the following discussion by the acronyms Pr-LCPI, PA-LCPI, and RW-LCPI, respectively. The $`n`$th order Pr-LCPI term is obtained by approximating the Brownian bridge by the $`n`$-dimensional process

$$S_u^n(\overline{a})=\underset{l=1}{\overset{k}{}}Y_l(u,\overline{a})+\underset{l=1}{\overset{j}{}}a_{k+1,l}F_{k+1,l}(u),$$

where $`k`$ and $`j`$ are the unique numbers such that $`n=2^k+j1`$, with $`k0`$ and $`1j2^k`$. However, it appears natural to utilize only the subsequence of the form $`n=2^k1`$ with $`k0`$, corresponding to $`k`$ complete layers and from now on we shall restrict our attention to this subsequence, for which

$$S_u^n(\overline{a})=\underset{l=1}{\overset{k}{}}Y_l(u,\overline{a}).$$
(20)
Using the notation introduced in Ref. Pre02, , we denote the tail of the series (19) by

$$B_u^n(\overline{a})=\underset{l=k+1}{\overset{\mathrm{}}{}}Y_l(u,\overline{a}).$$

To define the PA-LCPI method, besides the sum (20), we need to evaluate

$$\mathrm{\Gamma }_n^2(u)=\sigma ^2𝔼[B_u^n(\overline{a})^2]=\sigma ^2\underset{l=k+1}{\overset{\mathrm{}}{}}\underset{j=1}{\overset{2^{l1}}{}}F_{l,j}(u)^2.$$

This quantity must be computed explicitly because it enters the final PA-LCPI formula by means of the “effective” potential

$$\overline{V}_{u,n}(x)=_{}\frac{1}{\sqrt{2\pi \mathrm{\Gamma }_n^2(u)}}\mathrm{exp}\left[\frac{z^2}{2\mathrm{\Gamma }_n^2(u)}\right]V(x+z)dz.$$

For more information, the reader is referred to the Section III of Ref. Pre02, . For all $`lk+1`$ and $`1j2^l`$, the functions $`F_{l,j}`$ are zero on the points $`u_p=p/2^k`$ with $`p=0,1,\mathrm{}2^k`$. Let us define the *support* of the function $`F_{l,j}`$ as the set $`\text{supp}(F_{l,j})=\{u[0,1]:F_{l,j}(u)0\}`$. Moreover, for $`1p2^k`$, let $`I_p=\{(l,j):lk+1,\mathrm{\hspace{0.17em}1}j2^l,\text{supp}(F_{l,j})[u_{p1},u_p]\}`$ and define

$$W_p(u,\overline{a})=\underset{(l,j)I_p}{}a_{l,j}F_{l,j}(u).$$

Then a little thought and the use of the Ito-Nisio theorem show that $`W_p(u,\overline{a})`$ is a Brownian bridge on the interval $`[u_{p1},u_p]`$ of variance $`1/2^k`$. In addition, if $`p_1p_2`$, then the Brownian bridges $`W_{p_1}(u,\overline{a})`$ and $`W_{p_2}(u,\overline{a})`$ are independent because they are functions of the independent Gaussian random variables $`\{a_{l,j}\}`$ with $`(l,j)I_{p_1}`$ and $`(l,j)I_{p_2}`$, respectively and the sets of indexes $`I_{p_1}`$ and $`I_{p_2}`$ are disjoint. It is convenient to denote by $`𝔼_p`$ the conditional expectation over the random variables $`a_{l,j}`$ with $`(l,j)I_p`$. Then, we have

$$B_u^n(\overline{a})=\underset{p=1}{\overset{2^k}{}}W_p(u,\overline{a})$$

and

$$𝔼[B_u^n(\overline{a})^2]=\underset{p=1}{\overset{2^k}{}}𝔼[W_p(u,\overline{a})^2]=\underset{p=1}{\overset{2^k}{}}𝔼_p[W_p(u,\overline{a})^2].$$
(21)
However, one computes

$$\gamma _{n,1}^2(u)=𝔼_1[W_1(u,\overline{a})^2]=\{\begin{array}{cc}u(12^ku),\hfill & 0u<2^k,\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}$$
(22)
and then by translation

$$\gamma _{n,p}^2(u)𝔼_p[W_p(u,\overline{a})^2]=\gamma _{n,1}^2[u(p1)/2^k].$$
(23)
Clearly, the functions $`\gamma _{n,p}^2(u)`$ have *disjoint* support. Finally, Eq. (21) becomes

$$\mathrm{\Gamma }_n^2(u)=\sigma ^2\underset{p=1}{\overset{2^k}{}}\gamma _{n,p}^2(u),$$
(24)
which concludes the definition of the PA-LCPI method.

The reweighted technique is yet another way of improving the convergence of the primitive method. It has the advantage that it does not require the evaluation of the Gaussian transform of the potential. As discussed in Ref. Pre02, , the main idea is to simulate the effect of the partial averaging method by replacing the tail series $`B_u^n(\overline{a})`$ in the full series expansion by a collection of random variables $`\{R_u^n(b_1,\mathrm{},b_{n+q})\}_{0u1}`$ defined over an $`n+q`$ dimensional probability space ($`q`$ is a small integer which does not depend upon $`n`$). We ask that (i) the variance at the point $`u`$ of $`R_u^n(b_1,\mathrm{},b_{n+q})`$, denoted by $`\mathrm{\Gamma }_n^2(u)`$, be as close as possible to $`\mathrm{\Gamma }_n^2(u)`$ and (ii) the variables $`S_u^n(a_1,\mathrm{},a_n)`$ and $`R_u^n(b_1,\mathrm{},b_{n+q})`$ be independent and their sum have a joint distribution as close to a Brownian bridge as possible. One candidate for our approach is $`R_u^n(b_1,\mathrm{},b_n)=_{p=1}^{n+q}b_p\mathrm{\Omega }_p(u),`$ where $`b_1,\mathrm{},b_n`$ are i.i.d. standard normal random variables. Condition (ii) above is realized in the Ito-Nisio theorem by insuring that the collection $`\{F_{l,j}(u),\omega _p(u)\}`$ with $`1lk`$, $`1j2^{l1}`$, and $`1pn+q`$ is orthogonal and we shall look for such a collection. Here, $`\omega _p(u)`$ is the derivative of $`\mathrm{\Omega }_p(u)`$ and it is not required to be normalized.

As opposed to the Wiener-Fourier series, for the Lévy-Ciesielski series it is possible to enforce the condition (i) exactly. The analysis done for the partial averaging method showed that we can represent $`\mathrm{\Gamma }_n^2(u)`$ as the sum of $`2^k=n+1`$ functions of *disjoint* support and which can be obtained one from the other by translation. Intuitively, we must set $`q=1`$ and replace the $`2^k`$ Haar functions making up the $`k+1`$ layer by

$$\omega _p(u)=\frac{d}{du}\gamma _n[u(p1)/2^k].$$

It is easy to notice that the functions $`\omega _p(u)`$ are orthogonal among themselves because they have disjoint support. Moreover, it is not difficult to see that the Haar functions $`f_{l,j}(u)`$ are constant on the intervals $`[u_{p1},u_p]`$ for all $`lk`$ and therefore, they are orthogonal on the $`\omega _p(u)`$ functions because

$$_0^1\omega _p(u)du=_{x_{p1}}^{x_p}\omega _p(u)du=\gamma _{n,1}(1/2^k)\gamma _{n,1}(0)=0.$$

In consequence, the $`n`$-order RW-LCPI approximation uses the series

$$S_u^n(\overline{a})=\underset{l=1}{\overset{k}{}}Y_l(u,\overline{a})+\underset{p=1}{\overset{2^k}{}}a_{k+1,p}\gamma _{n,p}(u)$$
(25)
for its implementation. \[It is customary to define the approximation order by the dimensionality of the underlying probability space. We shall not apply this rule in the present paper in order to keep the unity of the exposition. The squares of the functions $`\gamma _{n,p}(u)`$ are given by the relation (23).\] A look at Fig 3 shows that the $`n`$-order RW-LCPI method is identical to the $`2n+1`$ order Pr-LCPI, except for the replacement of the last layer of functions $`F_{k+1,p}(u)`$ with $`\gamma _{n,p}(u)`$.

### III.2 Properties of the Lévy-Ciesielski path integral method.

As announced in the beginning of the section, the LCPI method for $`n=2^k1`$ is virtually a reformulation of the Discrete Path Integral method with appropriate short-time approximations. This is shown by the following result.

###### Theorem 2

If $`n=2^k1`$, then the following relations hold true:

1. The Trotter theorem

$`\rho (x,x^{};\beta )={\displaystyle _{}}dx_1\mathrm{}{\displaystyle _{}}dx_n\rho (x,x_1;{\displaystyle \frac{\beta }{n+1}})`$
$`\mathrm{}\rho (x_n,x^{};{\displaystyle \frac{\beta }{n+1}}).`$
2. If $`Mt`$ stands for any of the LCPI methods, then

$`\rho _n^{Mt}(x,x^{};\beta )={\displaystyle _{}}dx_1\mathrm{}{\displaystyle _{}}dx_n\rho _0^{Mt}(x,x_1;{\displaystyle \frac{\beta }{n+1}})`$
$`\mathrm{}\rho _0^{Mt}(x_n,x^{};{\displaystyle \frac{\beta }{n+1}}),`$
where the short-time approximations $`\rho _0^{Mt}(x_n,x^{};\beta )`$ are defined as follows:

$`{\displaystyle \frac{\rho _0^{Pr}(x,x^{};\beta )}{\rho _{fp}(x,x^{};\beta )}}=\mathrm{exp}\left\{\beta {\displaystyle _0^1}V[x+(x^{}x)u]du\right\},`$
$`{\displaystyle \frac{\rho _0^{PA}(x,x^{};\beta )}{\rho _{fp}(x,x^{};\beta )}}=\mathrm{exp}\left\{\beta {\displaystyle _0^1}\overline{V}_{u,0}[x+(x^{}x)u]du\right\},`$
and

$`{\displaystyle \frac{\rho _0^{RW}(x,x^{};\beta )}{\rho _{fp}(x,x^{};\beta )}}={\displaystyle _{}}dz\left(2\pi \right)^{1/2}e^{z^2/2}`$
$`\times \mathrm{exp}\left\{\beta {\displaystyle _0^1}V[x+(x^{}x)u+z\mathrm{\Gamma }_0(u)]du\right\},`$
respectively.

*Proof.* We only prove the first point of the theorem. It is not difficult to see that the Trotter theorem is in fact a part of the latter case with $`\rho _0^{Mt}(x,x^{};\beta )=\rho _n^{Mt}(x,x^{};\beta )=\rho (x,x^{};\beta )`$. As such, the second point follows by arguments similar to the first one and is left to the reader.

Let us remember that for all $`lk+1`$ and $`1j2^l`$, the functions $`F_{l,j}(u)`$ are zero on the points $`u_p=p/2^k`$ with $`p=0,\mathrm{}2^k`$. This means that the joint distribution of the Brownian bridge at these points is *uniquely* determined by the series (20). For this proof, it is important to notice that the inverse result is also true: knowledge of the joint distribution of the points $`u_p`$ with $`p=0,\mathrm{}2^k`$ uniquely determines the series (20) because the latter is *linear* on the intervals $`[u_{p1},u_p]`$. It follows that the variables $`B_{u_p}^0`$ are independent of the displaced and rescaled Brownian bridges $`W_p(u,\overline{a})`$. Using this information together with the joint distribution density for the random variables $`x_r(u_p)+\sigma B_{u_p}^0`$, which is given by the formula (II.2) as shown in the previous section, one computes

$`p_{\sigma ^2}(x^{}x)𝔼\mathrm{exp}\left\{\beta {\displaystyle _0^1}V\left[x_r(u)+\sigma B_u^0\right]du\right\}=p_{\sigma ^2}(x^{}x)𝔼\mathrm{exp}\left\{\beta {\displaystyle \underset{p=0}{\overset{n}{}}}{\displaystyle _{u_p}^{u_{p+1}}}V\left[x_r(u)+\sigma B_u^0\right]du\right\}`$
$`={\displaystyle _{}}dx_1\mathrm{}{\displaystyle _{}}dx_n{\displaystyle \underset{p=0}{\overset{n}{}}}\left(p_{\frac{\sigma ^2}{n+1}}(x_{p+1}x_p)𝔼_p\mathrm{exp}\left\{\beta {\displaystyle _{u_p}^{u_{p+1}}}V\left[x_p+{\displaystyle \frac{(uu_p)}{(u_{p+1}u_p)}}(x_{p+1}x_p)+\sigma W_p(u,\overline{a})\right]du\right\}\right)`$
$`={\displaystyle _{}}dx_1\mathrm{}{\displaystyle _{}}dx_n{\displaystyle \underset{p=0}{\overset{n}{}}}\left(p_{\frac{\sigma ^2}{n+1}}(x_{p+1}x_p)𝔼\mathrm{exp}\left\{{\displaystyle \frac{\beta }{n+1}}{\displaystyle _0^1}V\left[x_p+u(x_{p+1}x_p)+{\displaystyle \frac{\sigma }{\sqrt{n+1}}}B_u^0(\overline{a})\right]du\right\}\right).`$
which proves the first point of the theorem. The latter follows by a similar line of thought: for instance, the result for the primitive method is obtained by setting $`W_p(u,\overline{a})=0`$ in the previous formula. $`\mathrm{}`$

From the theoretical point of view, the importance of the Theorem 2 consists of the fact that it establishes a direct connection between the random series and the discrete path integral techniques, even if only for the $`n=2^k1`$ subsequence. As such, we notice that the primitive result was employed by Makri and Miller Mak88 ; Mak89  and by Mielke and Truhlar Mie01  as the ZOP-DPI method. The latter authors found that the asymptotic convergence of the method was $`O(1/n)`$. This result is in good agreement with the present analysis because the primitive Lévy-Ciesielski method cannot exceed the convergence rate of the most rapidly convergent series, the Wiener-Fourier series, which behaves asymptotically as $`O(1/n)`$ Pre02 .

The partial averaging result is not new either. The DPI formulation was used by Kole and De Raedt Kol01  to treat systems with negative coulombic singularities, for which the non-averaged methods are known to be ill-behaved. However, Kole and De Raedt were not aware of the fact that they were using the partial averaging method in a special setting and regarded their approach as a separate method. It has been shown in a mathematically rigorous way that the partial averaging method is convergent for such potentials at least as far as the pointwise density matrix, the partition function, and related integral expressions are concerned Pre02b , for all series representations of the Feynman-Kaç formula. Therefore, the Theorem 2 can be used to give a mathematically rigorous proof of the Kole and De Raedt result, which conversely can be thought of as an argument demonstrating the desirable properties of the partial averaging strategy.

In a related situation, the partial averaging method as specialized for the Wiener-Fourier series representation was used to treat the polaron problem by Alexandrou, Fleischer, and Rosenfelder Ale90 . Later, the DPI formulation of the partial averaging technique was applied by Titantah, Pierleoni, and Ciuchi Tit01  for the same polaron problem and regarded once again as a separate technique. We hope we have convinced the reader that given the multitude of series representations that may enter the Feynman-Kaç formula, there are an infinite number of ways in which the partial averaging idea can be implemented. The Wiener-Fourier and the Lévy-Ciesielski series representations as well as the related DPI implementation are only some instances (although perhaps the most important ones).

In Section II.C, we promised that we would find a quick way to compute the current paths for the standard DPI methods by means of the Lévy-Ciesielski series representation. For the LCPI formulation, it is straightforward to notice that the computational time necessary to compute the current path at a point $`u`$ is proportional to $`k=\mathrm{log}_2(n+1)`$ for the Pr-LCPI and PA-LCPI methods and $`1+\mathrm{log}_2(n+1)`$ for the RW-LCPI method, respectively. This is so because given a point $`u`$, the only Schauder function from the layer $`l`$ that is non-zero at the point $`u`$ is $`F_{l,j}(u)`$ with $`j=[2^{l1}u]+1`$, where $`[x]`$ denotes the integral part of $`x`$. For the RW-LCPI method, we have in addition that the only function $`\gamma _{n,j}(u)`$ which is non-zero at the point $`u`$ is the one with $`j=[2^ku]+1`$. In fact, going back to the proof of Theorem 2, we remember that the joint distribution of the points $`u_p`$ with $`p=0,\mathrm{}2^k`$ uniquely determines the series (20) because, in a more mathematical notation, we have

$$S_{u_p}^n(\overline{a})=B_{u_p}^0(\overline{a})=\underset{l=1}{\overset{k}{}}Y_l(u_p,\overline{a}),\mathrm{\hspace{0.33em}1}pn.$$
(26)
Equation (26) allows us to write the following special form for Eq. (II.3):

$`𝔼\{f[W_{x,x^{}}^\sigma (u_1),\mathrm{},W_{x,x^{}}^\sigma (u_n)]\}=𝔼\{f[x_r(u_1)+\sigma S_{u_1}^n(\overline{a}),\mathrm{},x_r(u_n)+\sigma S_{u_n}^n(\overline{a})]\}`$
$`={\displaystyle _{}}da_1\mathrm{}{\displaystyle _{}}da_n\left(2\pi \right)^{n/2}\mathrm{exp}\left({\displaystyle \frac{1}{2}}{\displaystyle \underset{l=1}{\overset{k}{}}}{\displaystyle \underset{i=1}{\overset{2^{l1}}{}}}a_{l,i}^2\right)`$ (27)
$`\times f[x_r(u_1)+\sigma {\displaystyle \underset{l=1}{\overset{k}{}}}F_{l,[2^{l1}u_1]+1}(u_1)a_{l,[2^{l1}u_1]+1},\mathrm{},x_r(u_n)+\sigma {\displaystyle \underset{l=1}{\overset{k}{}}}F_{l,[2^{l1}u_n]+1}(u_n)a_{l,[2^{l1}u_n]+1}].`$
This proves that the standard DPI method can be implemented so that the number of operations necessary to compute the current paths is $`O(kn)`$.

As we said at the beginning of the Section III, as opposed to the $`n=2^k1`$ subsequence of the Lévy-Ciesielski representation, no subsequence of the Wiener-Fourier representation can be rationalized as a DPI method. The precise meaning of this is that if

$$S_u^n(\overline{a})=\sqrt{\frac{2}{\pi ^2}}\underset{k=1}{\overset{n}{}}a_k\frac{\mathrm{sin}(k\pi u)}{k},$$

then there is no sequence $`0=u_0<u_1,\mathrm{},u_n<u_{n+1}=0`$ such that

$`𝔼\{f[W_{x,x^{}}^\sigma (u_1),\mathrm{},W_{x,x^{}}^\sigma (u_n)]\}`$ (28)
$`=`$ $`𝔼\{f[x_r(u_1)+\sigma S_{u_1}^n(\overline{a}),\mathrm{},x_r(u_n)+\sigma S_{u_n}^n(\overline{a})]\}`$
for all functions $`f(x_1,\mathrm{},x_n)`$. Indeed, remembering $`W_{x,x^{}}^\sigma (u)=x_r(u)+\sigma B_u^0`$ and choosing $`f(x)=[xx_r(u_p)]/\sigma ^2`$ for some interior point $`0<u_p<1`$, one computes

$$𝔼\{f[W_{x,x^{}}^\sigma (u_1),\mathrm{},W_{x,x^{}}^\sigma (u_n)]\}=𝔼(B_{u_p}^0)^2=u_p(1u_p).$$

On the other hand,

$`𝔼\{f[x_r(u_1)+\sigma S_{u_1}^n(\overline{a}),\mathrm{},x_r(u_n)+\sigma S_{u_n}^n(\overline{a})]\}`$
$`=`$ $`𝔼[S_{u_p}^n(\overline{a})]^2={\displaystyle \frac{2}{\pi ^2}}{\displaystyle \underset{k=1}{\overset{n}{}}}{\displaystyle \frac{\mathrm{sin}^2(k\pi u_p)}{k^2}}.`$
Clearly, the equality (28) cannot hold because

$$u_p(1u_p)\frac{2}{\pi ^2}\underset{k=1}{\overset{n}{}}\frac{\mathrm{sin}^2(k\pi u_p)}{k^2}=\frac{2}{\pi ^2}\underset{k=n+1}{\overset{\mathrm{}}{}}\frac{\mathrm{sin}^2(k\pi u_p)}{k^2}$$

does not vanish on the interval $`(0,1)`$. To prove this, it is enough to notice that the zeros of $`\mathrm{sin}^2[(n+1)\pi u_p]`$ and $`\mathrm{sin}^2[(n+2)\pi u_p]`$ are strictly interlaced.

The following theorem, whose proof is left to the reader, provides the necessary and sufficient conditions for an $`n`$-order term of an arbitrary series to admit a particular $`m`$-order DPI representation.

###### Theorem 3

Let

$$\mathrm{\Gamma }_n^2(u)=\sigma ^2\left[u(1u)\underset{k=1}{\overset{n}{}}\mathrm{\Lambda }_k(u)^2\right]$$

and let $`0=u_0<u_1<\mathrm{}<u_m<u_{m+1}=1`$. Then

$`𝔼\{f[W_{x,x^{}}^\sigma (u_1),\mathrm{},W_{x,x^{}}^\sigma (u_m)]\}`$
$`=`$ $`𝔼\{f[x_r(u_1)+\sigma S_{u_1}^n(\overline{a}),\mathrm{},x_r(u_m)+\sigma S_{u_m}^n(\overline{a})]\}`$
for all $`f:^m`$ if and only if $`\mathrm{\Gamma }_n^2(u_p)=0`$ for all $`1pm`$.

The fact that the Wiener-Fourier representation cannot be rationalized as a DPI method should not be surprising. Indeed, we presented enough evidence in Ref. Pre02,  to support the idea that the convergence of the partial averaging and the reweighted Wiener-Fourier path integral methods is $`𝒪(1/n^3)`$ for sufficiently smooth potentials. On the other hand, the analysis performed in Section II suggests that we cannot expect an asymptotic convergence of the DPI methods better than $`𝒪(1/n^2)`$. In fact, as we will show in the next subsection, the $`n=2^k1`$ subsequence of the PA-LCPI and RW-LCPI methods can have at most $`𝒪(1/n^2)`$ asymptotic convergence.

### III.3 Convergence of the PA-DPI and of the RW-DPI methods

In this subsection, we shall study the convergence of the Trotter product formulae having as short-time approximations the partial averaging and the reweighted zero order formulae given in Theorem 2. It is natural to call this methods the PA-DPI and the RW-DPI methods, respectively. In particular, by virtue of Theorem 2, we obtain the asymptotic rates of convergence for the subsequences $`n=2^k1`$ of the corresponding LCPI methods. To anticipate, the convergence of the partition function and of the density matrix will be shown to be $`𝒪(1/n^2)`$ for both methods if the potential is smooth enough. More precisely, we limit our discussion to the class of potentials introduced in Ref. Pre02b, , which are the Kato-class potentials Chu95  having finite Gaussian transform. In this section, a potential is called smooth if it lies in the local Sobolev space $`W_{\text{loc}}^{1,2}(^d)`$ and if the squares of the potentials and of the first order derivatives have finite Gaussian transform. We remind the reader that the local Sobolev space $`W_{\text{loc}}^{1,2}(^d)`$ is made up of all $`L_{\text{loc}}^2(^d)`$ functions whose first order distributional derivatives are also $`L_{\text{loc}}^2(^d)`$ functions i.e.,

$$_D\left[V(x)^2+\underset{i=1}{\overset{d}{}}|V(x)/x_i|^2\right]dx<\mathrm{}$$

for all bounded domains $`D^d`$. We warn the reader that the $`𝒪(1/n^2)`$ convergence of the density matrix and of the partition function for this class of potentials does not automatically imply similar convergence for the energy estimators, for which additional restrictions upon the class of potentials might be necessary.

To simplify the notation, we prove the convergence results for the monodimensional case and only state the multidimensional analogues. Let us start with the asymptotic convergence of the partial averaging method. If we set

$$U(x,x^{},\beta ;\overline{a})=_0^1V\left[x_r(u)+\sigma B_u^0(\overline{a})\right]du,$$
(29)
a little algebra shows that

$`\left|\rho (x,x^{};\beta )\rho _0^{PA}(x,x^{};\beta )\right|=\rho (x,x^{};\beta )\rho _0^{PA}(x,x^{};\beta )`$
$`=\rho _0^{PA}(x,x^{};\beta )𝔼\left\{e^{\beta [U(x,x^{},\beta ;\overline{a})𝔼U(x,x^{},\beta ;\overline{a})]}1\right\}.`$
The first equality follows from the fact that zero order PA density matrix is always smaller than the true density matrix, according to equation (18) of Ref. Pre02, . However, for $`\beta `$ small, we can expand the exponential in a Taylor series in order to establish the order of the short-time approximation. We have:

$`\rho (x,x^{};\beta )\rho _0^{PA}(x,x^{};\beta )=\rho _0^{PA}(x,x^{};\beta )`$
$`\times \{{\displaystyle \frac{\beta ^2}{2}}𝔼[U(x,x^{},\beta ;\overline{a})𝔼U(x,x^{},\beta ;\overline{a})]^2`$ (30)
$`+{\displaystyle \frac{\beta ^3}{3!}}𝔼[U(x,x^{},\beta ;\overline{a})𝔼U(x,x^{},\beta ;\overline{a})]^3+𝒪(\beta ^4)\}.`$
Notice that the term of order one in the Taylor expansion cancels, so the asymptotic behavior is dictated by the variance of the function $`U(x,x^{},\beta ;\overline{a})`$. However, looking at the expression (29), we see that this variance must also decay to zero as $`\beta 0`$, because $`\sigma 0`$. The same is true for the third order moment and a gradient expansion similar to the one performed in Appendix A for the variance of the function $`U(x,x^{},\beta ;\overline{a})`$ shows that

$$\frac{\beta ^3}{3!}𝔼\left[U(x,x^{},\beta ;\overline{a})𝔼U(x,x^{},\beta ;\overline{a})\right]^3$$

decays to zero as fast as $`𝒪(\beta ^{4.5})`$.

We shall be more careful in establishing a proper bound on the variance of the function $`U(x,x^{},\beta ;\overline{a})`$ because this will eventually dictate the asymptotic rate of convergence. As shown in Appendix A, we have

$`\beta T_1(x,x^{};\beta )`$ $``$ $`𝔼\left[U(x,x^{},\beta ;\overline{a})𝔼U(x,x^{},\beta ;\overline{a})\right]^2`$ (31)
$``$ $`\beta T_2(x,x^{};\beta ),`$
where the functions $`T_1(x,x^{};\beta )`$ and $`T_2(x,x^{};\beta )`$ satisfy the relation

$`T(x,x^{})=\underset{\beta 0}{lim}T_1(x,x^{};\beta )=\underset{\beta 0}{lim}T_2(x,x^{};\beta )`$
$`={\displaystyle \frac{\mathrm{}^2}{m_0}}{\displaystyle _0^1}du{\displaystyle _0^1}d\tau {\displaystyle \frac{u+\tau 2u\tau |u\tau |}{2}}`$ (32)
$`\times V^{(1)}[x_r(u)]V^{(1)}[x_r(\tau )].`$
In particular, the inequalities

$`\rho _0^{PA}(x,x^{};\beta )\rho (x,x^{};\beta )`$
$`\rho _0^{PA}(x,x^{};\beta )\left[1+{\displaystyle \frac{\beta ^3}{2}}T_2(x,x^{};\beta )+𝒪(\beta ^4)\right]`$
show that the zero order partial averaging formula is of convergence order 2. Therefore, the assertion of Makri and Miller Mak88  that $`\rho _0^{PA}(x,x^{};\beta )`$ is not an order 2 short-time approximation is wrong. Also, notice that

$$T(x,x)=\frac{\mathrm{}^2}{12m_0}V(x)^2,$$
(33)
because

$$_0^1du_0^1d\tau \frac{u+\tau 2u\tau |u\tau |}{2}=\frac{1}{12}.$$

Trotter composing the relation (31) $`n`$ times and noticing that $`𝒪(\beta ^4)`$ eventually contributes a term decaying as fast as $`1/n^3`$, it is but a simple task to establish the identity

$`{\displaystyle \frac{\beta ^3}{2(n+1)^3}}{\displaystyle \underset{j=0}{\overset{n}{}}}{\displaystyle _{}}dx_1{\displaystyle _{}}dx_2\rho _{j1}^{PA}(x,x_1;{\displaystyle \frac{j\beta }{n+1}})\rho _0^{PA}(x_1,x_2;{\displaystyle \frac{\beta }{n+1}})\rho _{nj1}^{PA}(x_2,x^{};{\displaystyle \frac{(nj)\beta }{n+1}})`$
$`\times T_1(x_1,x_2;{\displaystyle \frac{\beta }{n+1}})\rho (x,x^{};\beta )\rho _n^{PA}(x,x^{};\beta ){\displaystyle \frac{\beta ^3}{2(n+1)^3}}{\displaystyle \underset{j=0}{\overset{n}{}}}{\displaystyle _{}}\mathrm{d}x_1{\displaystyle _{}}\mathrm{d}x_2\rho _{j1}^{PA}(x,x_1;{\displaystyle \frac{j\beta }{n+1}})`$
$`\times \rho _0^{PA}(x_1,x_2;{\displaystyle \frac{\beta }{n+1}})\rho _{nj1}^{PA}(x_2,x^{};{\displaystyle \frac{(nj)\beta }{n+1}})T_2(x_1,x_2;{\displaystyle \frac{\beta }{n+1}}),`$ (34)
with the understanding that $`\rho _1^{PA}(x,x^{};0)=\delta (x^{}x)`$. The above inequality is valid to the order of $`𝒪(\beta ^4/n^3)`$. Now, notice that in the sense of distributions, we have

$`\underset{\beta 0}{lim}\rho _0^{PA}(x_1,x_2;\beta )T_1(x_1,x_2;\beta )=\underset{\beta 0}{lim}\rho _0^{PA}(x_1,x_2;\beta )`$
$`\times T_2(x_1,x_2;\beta )=\delta (x_1x_2)T(x_1,x_1).`$
Multiplying it by $`2(n+1)^2/\beta ^3`$ and using the previous observation, the formula (III.3) becomes

$`{\displaystyle \frac{1}{(n+1)}}{\displaystyle \underset{j=0}{\overset{n}{}}}{\displaystyle _{}}dx_1\rho _{j1}^{PA}(x,x_1;{\displaystyle \frac{j\beta }{n+1}})\rho _{nj1}^{PA}(x_1,x^{};{\displaystyle \frac{(nj)\beta }{n+1}})T(x_1,x_1)`$
$`{\displaystyle \frac{2(n+1)^2}{\beta ^3}}\left[\rho (x,x^{};\beta )\rho _n^{PA}(x,x^{};\beta )\right]{\displaystyle \frac{1}{(n+1)}}{\displaystyle \underset{j=0}{\overset{n}{}}}{\displaystyle _{}}dx_1\rho _{j1}^{PA}(x,x_1;{\displaystyle \frac{j\beta }{n+1}})`$
$`\times \rho _{nj1}^{PA}(x_1,x^{};{\displaystyle \frac{(nj)\beta }{n+1}})T(x_1,x_1),`$
in the limit that $`n`$ is large. Again in the same limit, the Riemann sum from the above expression transforms into an integral on the interval $`[0,1]`$ and combining everything we obtain the following theorem

###### Theorem 4

$`\underset{n\mathrm{}}{lim}{\displaystyle \frac{2(n+1)^2}{\beta ^3}}\left[\rho (x,x^{};\beta )\rho _n^{PA}(x,x^{};\beta )\right]`$
$`={\displaystyle \frac{\mathrm{}^2}{12m_0}}{\displaystyle _0^1}x|e^{\theta \beta H}V^2e^{(1\theta )\beta H}|x^{}d\theta `$ (35)
It is convenient to write Eq. (4) as

$`\rho (x,x^{};\beta )\rho _n^{PA}(x,x^{};\beta )+{\displaystyle \frac{\mathrm{}^2\beta ^3}{24m_0(n+1)^2}}`$
$`\times {\displaystyle _0^1}x|e^{\theta \beta H}V^2e^{(1\theta )\beta H}|x^{}\mathrm{d}\theta .`$ (36)
The $`d`$-dimensional version of Theorem 4 can be formally obtained by replacing $`V^2/m_0`$ with

$$\underset{i=1}{\overset{d}{}}\frac{1}{m_{0,i}}\left[\frac{}{x_i}V(x_1,\mathrm{},x_d)\right]^2.$$

Finally, we turn our attention to the convergence of the RW-DPI method. It was previously proved Pre02b  that for the class of potentials considered in this section the density matrix and the partition function of any partial averaging method is convergent to the correct result. However, this might not be true of the primitive and the reweighted methods, as well as of the standard DPI methods. Indeed, it is known that the non-averaged methods suffer from what is called “classical collapse” for potentials with negative coulombic singularities Law69 ; Kle95 ; Mus97 ; Kol01 , for which the partial averaging method is, however, convergent. For such systems it happens that the $`n`$ order partition functions of the primitive, reweighted, and standard DPI methods are always $`+\mathrm{}`$, yet the true quantum partition function is finite. This situation can be prevented by requiring that the *classical* partition function be finite. For instance, for the case of the primitive random series the Jensen’s inequality implies

$`Z_n^{Pr}(\beta )={\displaystyle \frac{1}{\sqrt{2\pi \sigma ^2}}}{\displaystyle _{}}\mathrm{d}x{\displaystyle _\mathrm{\Omega }}\mathrm{d}P(\overline{a})\mathrm{exp}\{\beta {\displaystyle _0^1}V[x`$
$`+\sigma {\displaystyle \underset{k=1}{\overset{n}{}}}a_k\mathrm{\Lambda }_k(u)]\mathrm{d}u\}{\displaystyle \frac{1}{\sqrt{2\pi \sigma ^2}}}{\displaystyle }_{}\mathrm{d}x{\displaystyle }_\mathrm{\Omega }\mathrm{d}P(\overline{a})`$
$`\times {\displaystyle _0^1}\mathrm{d}u\mathrm{exp}\{\beta V[x+\sigma {\displaystyle \underset{k=1}{\overset{n}{}}}a_k\mathrm{\Lambda }_k(u)]\}.`$
By changing the order of integration, one ends up with

$$Z_n^{Pr}(\beta )\frac{1}{\sqrt{2\pi \sigma ^2}}_{}e^{\beta V(x)}dx=Z_{\text{cl}}(\beta )<\mathrm{},$$
(37)
which proves our assertion. The inequality (37) holds for the reweighted methods and the standard DPI methods, too (for the latter techniques one uses the condition $`_iw_i=1`$ and the discrete analog of the Jensen’s inequality). In this paper, the condition $`Z_{\text{cl}}(\beta )<\mathrm{}`$ is assumed to hold any time one deals with the non-averaged methods.

Going back to the asymptotic convergence problem, we may follow the reasoning for the partial averaging method provided that we interpret $`𝔼^{}`$ to mean the average against the Gaussian measure

$$\mathrm{d}\mu (z)=\frac{1}{\sqrt{2\pi }}e^{z^2/2}\mathrm{d}z.$$

By Jensen’s inequality one proves the inequality

$`{\displaystyle \frac{\rho _0^{RW}(x,x^{};\beta )}{\rho _{fp}(x,x^{};\beta )}}={\displaystyle _{}}d\mu (z)\mathrm{exp}\left\{\beta {\displaystyle _0^1}V[x_r(u)+z\mathrm{\Gamma }_0(u)]du\right\}`$
$`\mathrm{exp}\left\{\beta {\displaystyle _0^1}{\displaystyle _{}}d\mu (z)V[x_r(u)+z\mathrm{\Gamma }_0(u)]du\right\}={\displaystyle \frac{\rho _0^{PA}(x,x^{};\beta )}{\rho _{fp}(x,x^{};\beta )}}.`$
Therefore, $`\rho _n^{RW}(x,x^{};\beta )\rho _n^{PA}(x,x^{};\beta )`$. Moreover, the following analog of Eq. (III.3) holds

$`\rho _0^{RW}(x,x^{};\beta )\rho _0^{PA}(x,x^{};\beta )=\rho _0^{PA}(x,x^{};\beta )`$
$`\times \{{\displaystyle \frac{\beta ^2}{2}}𝔼^{}[U^{}(x,x^{},\beta ;z)𝔼^{}U^{}(x,x^{},\beta ;z)]^2`$
$`+{\displaystyle \frac{\beta ^3}{3!}}𝔼^{}[U^{}(x,x^{},\beta ;z)𝔼^{}U^{}(x,x^{},\beta ;z)]^3+𝒪(\beta ^4)\},`$
where we now define $`U^{}(x,x^{},\beta ;z)=V[x_r(u)+z\mathrm{\Gamma }_0(u)]`$. As discussed in Appendix A, we have

$`\beta T_1^{}(x,x^{};\beta )`$ $``$ $`𝔼^{}\left[U^{}(x,x^{},\beta ;z)𝔼^{}U^{}(x,x^{},\beta ;z)\right]^2`$ (38)
$``$ $`\beta T_2^{}(x,x^{};\beta ),`$
where the functions $`T_1^{}(x,x^{};\beta )`$ and $`T_2^{}(x,x^{};\beta )`$ satisfy the relation

$`T^{}(x,x^{})=\underset{\beta 0}{lim}T_1^{}(x,x^{};\beta )=\underset{\beta 0}{lim}T_2^{}(x,x^{};\beta )`$
$`={\displaystyle \frac{\mathrm{}^2}{m_0}}{\displaystyle _0^1}du{\displaystyle _0^1}d\tau \sqrt{u(1u)\tau (1\tau )}`$ (39)
$`\times V^{(1)}[x_r(u)]V^{(1)}[x_r(\tau )].`$
We also have

$$T^{}(x,x)=\frac{\pi ^2\mathrm{}^2}{64m_0}V(x)^2$$

because

$$_0^1du_0^1d\tau \sqrt{u(1u)\tau (1\tau )}=\frac{\pi ^2}{64}.$$

We leave it for the reader to rework the previous arguments for the partial averaging case and show that for large $`n`$ we have

$`\underset{n\mathrm{}}{lim}{\displaystyle \frac{2(n+1)^2}{\beta ^3}}\left[\rho _n^{RW}(x,x^{};\beta )\rho _n^{PA}(x,x^{};\beta )\right]`$
$`={\displaystyle \frac{\pi ^2\mathrm{}^2}{64m_0}}{\displaystyle _0^1}x|e^{\theta \beta H}V^2e^{(1\theta )\beta H}|x^{}d\theta `$
Since $`\pi ^2/64>1/12`$, the previous result demonstrates that for $`n`$ large enough $`\rho _n^{RW}(x,x^{};\beta )\rho (x,x^{};\beta )`$, so that the convergence of the RW-DPI is eventually from above. Combining with Theorem 4, one obtains

###### Theorem 5

$`\underset{n\mathrm{}}{lim}{\displaystyle \frac{2(n+1)^2}{\beta ^3}}\left[\rho (x,x^{};\beta )\rho _n^{RW}(x,x^{};\beta )\right]={\displaystyle \frac{\mathrm{}^2}{4m_0}}`$
$`\times \left({\displaystyle \frac{\pi ^2}{16}}{\displaystyle \frac{1}{3}}\right){\displaystyle _0^1}x|e^{\theta \beta H}V^2e^{(1\theta )\beta H}|x^{}d\theta .`$
As for the partial averaging case, the statement of Theorem 5 can be written in the short form

$`\rho (x,x^{};\beta )\rho _n^{RW}(x,x^{};\beta ){\displaystyle \frac{\mathrm{}^2\beta ^3}{8m_0(n+1)^2}}\left({\displaystyle \frac{\pi ^2}{16}}{\displaystyle \frac{1}{3}}\right)`$
$`\times {\displaystyle _0^1}x|e^{\theta \beta H}V^2e^{(1\theta )\beta H}|x^{}\mathrm{d}\theta .`$ (40)
From the Theorems (4) and (5) and by using cyclic invariance, one easily proves the following relations:

###### Corollary 1

$`{\displaystyle \frac{Z(\beta )Z_n^{PA}(\beta )}{Z(\beta )}}{\displaystyle \frac{\mathrm{}^2\beta ^3}{24m_0(n+1)^2}}{\displaystyle \frac{_{}\rho (x;\beta )V(x)^2dx}{_{}\rho (x;\beta )dx}}`$
and

$`{\displaystyle \frac{Z_n^{RW}(\beta )Z(\beta )}{Z(\beta )}}{\displaystyle \frac{\mathrm{}^2\beta ^3}{8m_0(n+1)^2}}\left({\displaystyle \frac{\pi ^2}{16}}{\displaystyle \frac{1}{3}}\right)`$
$`\times {\displaystyle \frac{_{}\rho (x;\beta )V(x)^2dx}{_{}\rho (x;\beta )dx}}.`$
*Observation 1* We have $`\pi ^2/161/30.284<1/3`$, so one may be tempted to say that the reweighted technique converges at a faster rate than the partial averaging method. However, as previously mentioned in the text, both the $`n`$-order LCPI and DPI reweighted techniques actually uses $`2n+1`$ random variables to parameterize the paths. If the convention of denoting the order of an approximation by the number of variables used to parameterize the paths is obeyed, then the constant $`\pi ^2/161/3`$ should be increased four times. In this case, we have $`40.284=1.134`$ which means that the partial averaging is about $`1.134/(1/3)=3.4`$ times faster than the reweighted technique.

*Observation 2* The asymptotic relative errors for the partition functions shown in Corollary 1 can be evaluated during the Monte Carlo procedure if so desired. It is a fact established in several occasions Pre02 ; Pre02b  that the convergence of the partial averaging density matrix and partition functions for all series representations is monotonically from below. In particular, the PA-DPI subsequence $`n=2^k1`$ has the same property since it is identical to the respective subsequence of the PA-LCPI method. However, it might be possible that the partition function for the reweighted methods are monotonically decreasing from above for the $`n=2^k1`$ subsequence. In fact, Golden Gol65  and Thompson Tho65  have shown that the partition function for the subsequence $`n=2^k1`$ of the trapezoidal Trotter DPI method is monotonically decreasing and this might be true of the RW-DPI method, too.

Let us remember that there are potentials, as for instance the potentials with negative coulombic singularities, for which the non-averaged methods do not converge. Conversely, there are smooth and bounded from below potentials, as for instance $`V(x)=\mathrm{exp}(x^4)`$ for which the non-averaged methods are convergent to the correct result yet the partial averaging method is not convergent because $`V(x)=\mathrm{exp}(x^4)`$ does not have a finite Gaussian transform. For such potentials, it is expected that the Theorem 5 as well as the second part of the Corollary 1 are still true.

We shall reinforce the conclusions of this section by verifying the theorems (4) and (5) for the simple case of the quadratic potential $`V(x)=m_0\omega ^2x^2/2`$. Again we use atomic units and set $`m_0=1`$, $`\omega =1`$, and $`\beta =10`$. The evaluation of the $`n`$-order partial averaging and reweighted elements $`\rho _n^{\text{PA}}(0;\beta )`$ and $`\rho _n^{\text{RW}}(0;\beta )`$ is analyzed in Appendix B. As discussed in Ref. Pre02, , for each method $`Mt`$ the convergence constant

$$c_{Mt}=\underset{n\mathrm{}}{lim}\frac{\rho (0;\beta )\rho _n^{Mt}(0;\beta )}{(n+1)^2}$$

can be obtained numerically by analyzing the asymptotic slope of the equation

$$c_n^{Mt}=(4n+2)^2(n+1/2)\left[\rho (0;\beta )\rho _{4n+2}^{Mt}(0;\beta )\right],$$

as a function of $`n`$. More precisely, we have $`c_{Mt}=lim_n\mathrm{}c_n^{Mt}c_{n1}^{Mt}`$. On the other hand, with the help of the exact density matrix $`\rho (x,x^{};\beta )`$ of the quadratic potential Fey94  and of the relations (III.3) and (III.3), one computes

$`c_{PA}={\displaystyle \frac{\beta ^3}{24}}{\displaystyle _0^1}d\theta {\displaystyle _{}}dx\rho (0,x;\theta \beta )`$
$`\times \rho [x,0;(1\theta )\beta ]x^2=0.0713`$
and $`c_{RW}=[(3\pi ^2/16)1]c_{PA}=0.0606`$, respectively.

The plots in Fig. 4 show that indeed, the current slopes $`c_n^{Mt}c_{n1}^{Mt}`$ converge to the corresponding values predicted by the theorems (4) and (5).

To conclude this section, we analyze how smooth realistic three dimensional potentials must be to fit the hypothesis of the the theorems (4) and (5). A prototypical example is the tridimensional spherical potential

$$V(r)=\frac{1}{2}r^2+\frac{1}{r^\alpha },0<\alpha 1,$$
(41)
for which both the partial averaging and the reweighted DPI methods are convergent because the potential $`V(r)`$ is a Kato-class potential having a finite Gaussian transform and because $`Z_{cl}(\beta )<\mathrm{}`$. The reader may easily verify that $`1/r^\alpha ^2=\alpha ^2/r^{2\alpha +2}`$ is locally integrable if and only if $`\alpha <1/2`$. Therefore, if $`\alpha <1/2`$, the theorems (4) and (5) apply and the convergence of both methods is $`𝒪(1/n^2)`$.

On the other hand, if $`\alpha 1/2`$, the convergence cannot be $`𝒪(1/n^2)`$ because the convergence constants are $`+\mathrm{}`$. This can be proved by using the additional information that the density matrix for the Kato-class potentials is continuous and strictly positive. In particular, there is $`ϵ>0`$ and $`\eta >0`$ such that $`\rho (x,\beta )ϵ`$ for all $`r\eta 1`$. Therefore, looking at the bounds for the partition functions given by Corollary 1, we have

$`{\displaystyle _^3}\rho (\overline{x};\beta )V(\overline{x})^2d\overline{x}4\pi ϵ\alpha ^2{\displaystyle _0^\eta }r^{2\alpha }dr`$
$`4\pi ϵ\alpha ^2{\displaystyle _0^1}r^1dr=+\mathrm{}.`$
We have treated this problem explicitly in order to show that the nature of the singularities of the potential affects the rate of convergence even if the singularities are oriented “upward.” Therefore, in “pushing” the Monte Carlo simulation to the limits, the reader may want to actually remove these singularities if they are physically irrelevant. He/she can do this either by a simple truncation or by approximating the singularity with a better behaved one.

Other prototypical examples of potentials are those having negative singularities

$$V(r)=\frac{1}{2}r^2\frac{1}{r^\alpha },0<\alpha 1.$$
(42)
For such potentials, the classical partition function is not finite and the reweighted technique does not properly converge. However, the partial averaging method is convergent and, if $`\alpha <1/2`$, the asymptotic convergence is $`𝒪(1/n^2)`$. The findings of this section demonstrate that “smooth enough” potentials may actually be discontinuous in the three dimensional space.

## IV  Summary and Discussion

A central theme of the present paper has been the characterization of various path integral approaches and the exploration of their interconnections. The Kaç interpretation of the Feynman approach is a valuable tool for such an analysis. We notice that it is difficult and unnatural to introduce the random series representation by means of the Trotter product rule. Indeed, in order to show that the path integrals

$$_0^1V[x_r(u)+B_u^0(\overline{a})]du$$
(43)
are correctly defined, one utilizes the fact that, with probability one, the Brownian paths are continuous. This property of the Brownian motion is not readily available from the Trotter product rule. However, as we have shown in Section II, the discrete methods can be directly derived from the Feynman-Kaç formula by simply replacing the integrals given by Eq. (43) with appropriate quadrature sums.

We have explored at some length two particular implementations of path integral methods: the Lévy-Ciesielski approach and the associated DPI technique. We have considered primitive, partial averaged and reweighted versions of this methods. As discussed in Section III, the Lévy-Ciesielski approach is of particular importance because its $`n=2^k1`$ subsequence can be rationalized both as a series and as a discrete method. This dual character is valuable for several reasons. For example, it provides a convenient and rigorous reformulation of Coalson’s findings linking series and discrete path integral methods, and, as illustrated by Eq. (III.2), suggests a means for reducing the numerical overhead associated with path construction. Using the unified framework the Lévy-Ciesielski approach provides, we have shown that the methods introduced by Kole and De Raedt Kol01  for systems with negative coulombic singularities as well as those introduced by Titantah, Pierleoni and Ciuchi Tit01  for the polaron problem are discrete versions of the partial averaging approach. Furthermore, Theorem 2 of Section III suggests that these previous methods can be implemented in a more robust manner using the Lévy-Ciesielski series approach.

We have been able to characterize the convergence properties of the partial averaging and reweighted DPI approaches and, therefore, of the $`n=2^k1`$ subsequence of the corresponding LCPI techniques. In this respect, Theorems 4 and 5 of Section IV provide sharp estimates of the convergence constants for the calculation of density matrix elements for both the partial averaged and reweighted methods. To our knowledge, this is the first time that such exact convergence constants have been established. Beyond their intrinsic interest, knowledge of these convergence constants can be used to devise an improved numerical implementation of the Feynman-Kaç approach. In particular, the results of Section IV indicate that the convergence constants for the reweighted and partial averaged methods are related by the formula $`c_{RW}=[(3\pi ^2/16)1]c_{PA}`$ for all pairs of points $`(x,x^{})`$ and for all $`\beta >0`$. Because the leading terms in $`1/n^2`$ thus cancel, the approach defined by the equation

$$\rho _n^{}(x,x^{};\beta )=\frac{\rho _n^{\text{RW}}(x,x^{};\beta )+[(3\pi ^2/16)1]\rho _n^{\text{PA}}(x,x^{};\beta )}{3\pi ^2/16}$$

has an asymptotic convergence better than $`𝒪(1/n^2)`$, i.e.

$`\underset{n\mathrm{}}{lim}{\displaystyle \frac{2(n+1)^2}{\beta ^3}}\left[\rho (x,x^{};\beta )\rho _n^{}(x,x^{};\beta )\right]=0.`$
In fact, we believe that if the potential $`V(x)`$ has also a well behaved second derivative, the convergence order of the new method is $`𝒪(1/n^3)`$.

Finally, we note that with the help of Theorems 4 and 5, the asymptotic behavior of the so-called T-method and H-method energy estimators (c.f. Section IV of Ref. Pre02, ) can be examined. In particular, it should be possible to deduce the convergence constants for these estimators from those of the corresponding density matrix expressions. We leave a detailed analysis of such issues for future discussion.

###### Acknowledgements.

The authors acknowledge support from the National Science Foundation through awards CHE-0095053 and CHE-0131114.

## Appendix A

It is well known that if $`A,B>0`$, and $`\alpha =C/\sqrt{AB}`$ such that $`|\alpha |<1`$, then the following Mehler’s formula Sze75  holds for all $`f`$ and $`g`$ whose squares have finite Gaussian transforms

$`\overline{[fg]}_{ABC}(x_0,y_0)`$ $`=`$ $`{\displaystyle _{}}dx{\displaystyle _{}}dy{\displaystyle \frac{1}{2\pi }}{\displaystyle \frac{1}{\sqrt{ABC^2}}}\mathrm{exp}\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{x^2B+y^2A2xyC}{ABC^2}}\right)f(x_0+x)g(y_0+y)`$
$`=`$ $`{\displaystyle _{}}dx{\displaystyle _{}}dy{\displaystyle \frac{1}{2\pi }}{\displaystyle \frac{1}{\sqrt{1\alpha ^2}}}\mathrm{exp}\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{x^2+y^22xy\alpha }{1\alpha ^2}}\right)f(x_0+x\sqrt{A})g(y_0+y\sqrt{B})`$
$`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}\alpha ^k{\displaystyle _{}}dx{\displaystyle _{}}dye^{(x^2+y^2)/2}H_k(x)H_k(y)f(x_0+x\sqrt{A})g(y_0+y\sqrt{B}).`$
In the above, the functions $`H_k(x)`$ are the normalized Hermite polynomials corresponding to the Gaussian weight

$$\mathrm{d}\mu (x)=\frac{1}{\sqrt{2\pi }}e^{x^2/2}.$$

They form a complete orthonormal basis in the Hilbert space $`L_\mu ^2()`$, which is endowed with the scalar product

$$\psi |\varphi =_{}\psi (x)\varphi (x)d\mu (x).$$

Let us notice that according to our hypothesis, the functions $`f(x_0+x\sqrt{A})`$ and $`g(y_0+x\sqrt{B})`$ as functions of $`x`$ are square integrable against $`\mathrm{d}\mu (x)`$ and thus they lie in the Hilbert space $`L_\mu ^2()`$.

By repeated integration by parts, the formula (A) is shown to equal

$$\overline{[fg]}_{ABC}(x_0,y_0)=\underset{k=0}{\overset{\mathrm{}}{}}\frac{C^k}{k!}\overline{f}_A^{(k)}(x_0)\overline{g}_B^{(k)}(y_0),$$
(45)
where in general $`\overline{f}_A^{(k)}(x_0)`$ is the $`k`$-order derivative of

$$\overline{f}_A(x_0)=_{}\frac{1}{\sqrt{2\pi A}}e^{z^2/(2A)}f(x_0+z)dz.$$

Let us notice that the series (A) can be extended to the case $`\alpha =1`$, too. Indeed, the last series in Eq. (A) for the case $`\alpha =1`$ is nothing else but the Bessel series

$$\underset{k=0}{\overset{\mathrm{}}{}}H_k|f(x_0+\sqrt{A})H_k|g(y_0+\sqrt{B}),$$

which is convergent to

$`f(x_0+\sqrt{A})|g(y_0+\sqrt{B})`$
$`={\displaystyle _{}}f(x_0+x\sqrt{A})g(y_0+x\sqrt{B})d\mu (x).`$
Next, we proceed to establish the inequality (31) from section III.C. We start with the identity

$`𝔼\left[U(x,x^{},\beta ;\overline{a})𝔼U(x,x^{},\beta ;\overline{a})\right]^2`$
$`=𝔼U(x,x^{},\beta ;\overline{a})^2\left[𝔼U(x,x^{},\beta ;\overline{a})\right]^2.`$ (46)
Clearly, we have

$$𝔼U(x,x^{},\beta ;\overline{a})=\overline{V}_{u,0}[x_r(u)].$$
(47)
Moreover,

$`𝔼U(x,x^{},\beta ;\overline{a})^2`$
$`={\displaystyle _0^1}du{\displaystyle _0^1}d\tau 𝔼V[x_r(u)+\sigma B_u^0]V[x_r(\tau )+\sigma B_\tau ^0]`$ (48)
and the variables $`B_u^0`$ and $`B_\tau ^0`$ have a joint Gaussian distribution of covariances

$`𝔼(B_u^0)^2=u(1u),𝔼(B_\tau ^0)^2=\tau (1\tau )`$
$`\text{and}𝔼(B_u^0B_\tau ^0)={\displaystyle \frac{u+\tau 2u\tau |u\tau |}{2}}.`$
This covariance matrix is independent of any particular representation of the Brownian bridge and therefore can be computed with the help of any basis. For instance, using the Wiener-Fourier basis, the last term of the above formula reads

$$𝔼(B_u^0B_\tau ^0)=\frac{2}{\pi ^2}\underset{k=1}{\overset{\mathrm{}}{}}\frac{\mathrm{sin}(k\pi u)\mathrm{sin}(k\pi \tau )}{k^2}$$
(49)
and the sum of the above series can be shown to equal $`(u+\tau 2u\tau |u\tau |)/2`$. It is useful to define the quantities $`G_0(u,\tau )=\sigma ^2𝔼(B_u^0B_\tau ^0)`$ and

$$\mathrm{\Delta }_0^2(u,\tau )=\mathrm{\Gamma }_0^2(u)\mathrm{\Gamma }_0^2(\tau )G_0(u,\tau )^2.$$

Then,

$`𝔼U(x,x^{},\beta ;\overline{a})^2={\displaystyle _0^1}du{\displaystyle _0^1}d\tau {\displaystyle _{}}dx{\displaystyle _{}}dy{\displaystyle \frac{1}{2\pi \mathrm{\Delta }_0(u,\tau )}}`$
$`\times \mathrm{exp}\left\{{\displaystyle \frac{1}{2}}{\displaystyle \frac{x^2\mathrm{\Gamma }_0^2(\tau )+y^2\mathrm{\Gamma }_0^2(u)2xyG_0(u,\tau )}{\mathrm{\Delta }_0^2(u,\tau )}}\right\}`$
$`\times V[x_r(u)+x]V[x_r(\tau )+y].`$
Using the expansion (45), one may write the above integral as the sum of the series

$`𝔼U(x,x^{},\beta ;\overline{a})^2`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{k!}}{\displaystyle _0^1}du{\displaystyle _0^1}d\tau G_0(u,\tau )^k`$
$`\times `$ $`\overline{V}_{u,0}^{(k)}[x_r(u)]\overline{V}_{\tau ,0}^{(k)}[x_r(\tau )],`$
where $`\overline{V}_{u,0}^{(k)}(x)`$ is the $`k`$ order derivative of $`\overline{V}_{u,0}(x).`$ With the help of Eq. (47), one recognizes the first term of the above series to be $`[𝔼U(x,x^{},\beta ;\overline{a})]^2`$, so that Eq. (A) becomes

$`𝔼\left[U(x,x^{},\beta ;\overline{a})𝔼U(x,x^{},\beta ;\overline{a})\right]^2={\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{k!}}`$
$`\times {\displaystyle _0^1}\mathrm{d}u{\displaystyle _0^1}\mathrm{d}\tau G_0(u,\tau )^k\overline{V}_{u,0}^{(k)}[x_r(u)]\overline{V}_{\tau ,0}^{(k)}[x_r(\tau )].`$ (50)
Now, we make an important observation: as its eigenfunction expansion (49) shows, $`G_0(u,\tau )`$ is a positive definite integral kernel $`L^2([0,1])L^2([0,1])`$ and it is not difficult to verify that all $`G_0(u,\tau )^k`$ are positive definite. Therefore,

$$_0^1du_0^1d\tau G_0(u,\tau )^k\overline{V}_{u,0}^{(k)}[x_r(u)]\overline{V}_{\tau ,0}^{(k)}[x_r(\tau )]0$$

for all $`k1`$. Considering only the first term in the series (A), we obtain the inequality

$`𝔼\left[U(x,x^{},\beta ;\overline{a})𝔼U(x,x^{},\beta ;\overline{a})\right]^2\beta T_1(x,x^{};\beta )`$ (51)
where

$`T_1(x,x^{};\beta )={\displaystyle \frac{\mathrm{}^2}{m_0}}{\displaystyle _0^1}du{\displaystyle _0^1}d\tau 𝔼[B_u^0B_\tau ^0]`$
$`\times \overline{V}_{u,0}^{(1)}[x_r(u)]\overline{V}_{\tau ,0}^{(1)}[x_r(\tau )].`$
It is not difficult to see that as $`\beta 0`$ we have $`\mathrm{\Gamma }_0^2(u)0`$ and so,

$`T(x,x^{})=\underset{\beta 0}{lim}T_1(x,x^{};\beta )={\displaystyle \frac{\mathrm{}^2}{m_0}}`$
$`\times {\displaystyle _0^1}\mathrm{d}u{\displaystyle _0^1}\mathrm{d}\tau 𝔼[B_u^0B_\tau ^0]V^{(1)}[x_r(u)]V^{(1)}[x_r(\tau )].`$ (52)
To prove the second inequality in Eq. (31), one uses the inequality $`1/k!1/(k1)!`$ and the positivity of the terms of the series (A) to establish the inequality

$`𝔼\left[U(x,x^{},\beta ;\overline{a})𝔼U(x,x^{},\beta ;\overline{a})\right]^2{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{k!}}{\displaystyle _0^1}du{\displaystyle _0^1}d\tau G_0(u,\tau )^{k+1}\overline{V}_{u,0}^{(k+1)}[x_r(u)]\overline{V}_{\tau ,0}^{(k+1)}[x_r(\tau )]`$
$`={\displaystyle _0^1}du{\displaystyle _0^1}d\tau G_0(u,\tau ){\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{k!}}G_0(u,\tau )^k\overline{V}_{u,0}^{(k+1)}[x_r(u)]\overline{V}_{\tau ,0}^{(k+1)}[x_r(\tau )]={\displaystyle _0^1}du{\displaystyle _0^1}d\tau G_0(u,\tau ){\displaystyle _{}}dx{\displaystyle _{}}dy`$
$`\times {\displaystyle \frac{1}{2\pi \mathrm{\Delta }_0(u,\tau )}}\mathrm{exp}\left\{{\displaystyle \frac{1}{2}}{\displaystyle \frac{x^2\mathrm{\Gamma }_0^2(\tau )+y^2\mathrm{\Gamma }_0^2(u)2xyG_0(u,\tau )}{\mathrm{\Delta }_0^2(u,\tau )}}\right\}V^{(1)}[x_r(u)+x]V^{(1)}[x_r(\tau )+y].`$
Therefore,

$$𝔼\left[U(x,x^{},\beta ;\overline{a})𝔼U(x,x^{},\beta ;\overline{a})\right]^2\beta T_2(x,x^{};\beta ),$$
(53)
where

$`T_2(x,x^{};\beta )={\displaystyle \frac{\mathrm{}^2}{m_0}}{\displaystyle _0^1}du{\displaystyle _0^1}d\tau 𝔼[B_u^0B_\tau ^0]{\displaystyle _{}}dx{\displaystyle _{}}dy`$
$`\times {\displaystyle \frac{1}{2\pi \mathrm{\Delta }_0(u,\tau )}}\mathrm{exp}\left\{{\displaystyle \frac{1}{2}}{\displaystyle \frac{x^2\mathrm{\Gamma }_0^2(\tau )+y^2\mathrm{\Gamma }_0^2(u)2xyG_0(u,\tau )}{\mathrm{\Delta }_0^2(u,\tau )}}\right\}`$
$`\times V^{(1)}[x_r(u)+x]V^{(1)}[x_r(\tau )+y]`$
Again, as $`\beta 0`$ we have $`\mathrm{\Gamma }_0^2(u)0`$ and

$`T(x,x^{})=\underset{\beta 0}{lim}T_2(x,x^{};\beta )={\displaystyle \frac{\mathrm{}^2}{m_0}}`$
$`\times {\displaystyle _0^1}\mathrm{d}u{\displaystyle _0^1}\mathrm{d}\tau 𝔼[B_u^0B_\tau ^0]V^{(1)}[x_r(u)]V^{(1)}[x_r(\tau )].`$ (54)
The relations (51),(A),(53), and (A) combined prove the equations (31) and (III.3) from Section III.C. The relations (38) and (III.3) follow by a similar reasoning and their proof is left to the reader. We only mention that one starts with the fact that the series (A) is well defined and convergent for $`\alpha =1`$ too, as shown in the beginning of the present appendix.

## Appendix B

In this section we discuss the computation of the matrix element $`0|e^{\beta H}|0`$ for the quadratic potential $`V(x)=m_0\omega ^2x^2/2`$ by means of the standard DPI method and of the partial averaging and the reweighted DPI methods. The density matrix for the quadratic potential is known analytically (see Ref. Fey94, ) and we do not reproduce it here. For a standard DPI method specified by the quadrature points $`0=u_0<u_1<\mathrm{}<u_n<u_{n+1}=1`$, by the increments $`\theta _i=u_{i+1}u_i`$, and by the weights $`w_0,w_1,\mathrm{},w_{n+1}`$, the formula (12) becomes

$`\rho _n^{\text{DPI}}(0;\beta )={\displaystyle _{}}dx_1\mathrm{}{\displaystyle _{}}dx_np_{\sigma ^2\theta _1}(x_1)`$
$`\times p_{\sigma ^2\theta _2}(x_2x_1)\mathrm{}p_{\sigma ^2\theta _n}(x_nx_{n1})p_{\sigma ^2\theta _n}(x_n)`$ (55)
$`\times \mathrm{exp}\left\{{\displaystyle \frac{m_0\omega ^2\beta }{2}}{\displaystyle \underset{i=1}{\overset{n}{}}}w_ix_i^2\mathrm{d}u\right\}.`$
Remember that $`\sigma ^2=\beta \mathrm{}^2/m_0`$. If we set $`\overline{x}^T=(x_1,x_2,\mathrm{},x_n)`$, the above $`n`$-dimensional integral can be written in the compact form

$`\rho _n^{\text{DPI}}(0,0;\beta )=\left({\displaystyle \underset{i=0}{\overset{n}{}}}{\displaystyle \frac{1}{2\pi \sigma ^2\theta _i}}\right)^{1/2}{\displaystyle _^n}e^{\overline{x}^TA\overline{x}/2}d\overline{x}`$
$`=\left({\displaystyle \underset{i=0}{\overset{n}{}}}{\displaystyle \frac{1}{2\pi \sigma ^2\theta _i}}\right)^{1/2}\left[det\left({\displaystyle \frac{A}{2\pi }}\right)\right]^{1/2},`$ (56)
where the matrix $`A`$ is the tridiagonal matrix defined by

$$A_{i,i}=\frac{1}{\sigma ^2\theta _{i1}}+\frac{1}{\sigma ^2\theta _i}+m_0\omega ^2\beta w_i\text{for}\mathrm{\hspace{0.33em}1}in$$

and

$$A_{i,i+1}=A_{i+1,i}=\frac{1}{\sigma ^2\theta _i}\text{for}\mathrm{\hspace{0.33em}1}in1.$$

The values of the quadrature points and the corresponding weights for the trapezoidal rule are well known, while for the Gauss-Legendre quadrature scheme the reader may use the routine given in Ref. (Pre92, ).

The zero order partial averaging density matrix for the quadratic potential has the explicit expression

$`\rho _0^{\text{PA}}(x,x^{};\beta )=p_{\sigma ^2}(x^{}x)`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{m_0\omega ^2\beta }{6}}\left(x^2+x_{}^{}{}_{}{}^{2}+xx^{}+\sigma ^2/2\right)\right].`$ (57)
Using Eq. (B), the reader may easily deduce that the corresponding $`n`$-order PA-DPI density matrix is

$`\rho _n^{\text{PA}}(0;\beta )=\left({\displaystyle \frac{n+1}{2\pi \sigma ^2}}\right)^{(n+1)/2}`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{\beta ^2\mathrm{}^2\omega ^2}{12(n+1)}}\right]\left[det\left({\displaystyle \frac{A}{2\pi }}\right)\right]^{1/2},`$ (58)
where the tridiagonal matrix $`A`$ is defined by the relations

$$A_{i,i}=2\left[\frac{n+1}{\sigma ^2}+\frac{m_0\omega ^2\beta }{3(n+1)}\right]\text{for}\mathrm{\hspace{0.33em}1}in$$

and

$$A_{i,i+1}=A_{i+1,i}=\frac{n+1}{\sigma ^2}+\frac{m_0\omega ^2\beta }{6(n+1)}\text{for}\mathrm{\hspace{0.33em}1}in1.$$

Finally, the zero order reweighted density matrix has the form

$`\rho _0^{\text{RW}}(x,x^{};\beta )=p_{\sigma ^2}(x^{}x)(1+\beta ^2\mathrm{}^2\omega ^2/6)^{1/2}`$
$`\times \mathrm{exp}\{{\displaystyle \frac{m_0\omega ^2\beta }{2}}[{\displaystyle \frac{x^2+x_{}^{}{}_{}{}^{2}+xx^{}}{3}}`$ (59)
$`{\displaystyle \frac{\pi ^2}{16^2}}{\displaystyle \frac{\beta ^2\mathrm{}^2\omega ^2(x+x^{})^2}{1+\beta ^2\mathrm{}^2\omega ^2/6}}]\},`$
which can be deduced by direct integration. Let us set

$$\eta _n^2=1+\frac{\beta ^2\mathrm{}^2\omega ^2}{6(n+1)^2}.$$

Then,

$$\rho _n^{\text{RW}}(0;\beta )=\left(\frac{n+1}{2\pi \sigma ^2}\frac{1}{\eta _n^2}\right)^{\frac{n+1}{2}}\left[det\left(\frac{A}{2\pi }\right)\right]^{1/2},$$
(60)
where the tridiagonal matrix $`A`$ is defined by the relations

$$A_{i,i}=2\left[\frac{n+1}{\sigma ^2}+\frac{m_0\omega ^2\beta }{3(n+1)}\left(\frac{\pi }{16}\right)^2\frac{m_0\mathrm{}^2\beta ^3\omega ^4}{\eta _n^2(n+1)^3}\right]$$

for $`1in`$ and

$$A_{i,i+1}=A_{i+1,i}=\frac{n+1}{\sigma ^2}+\frac{m_0\omega ^2\beta }{6(n+1)}\left(\frac{\pi }{16}\right)^2\frac{m_0\mathrm{}^2\beta ^3\omega ^4}{\eta _n^2(n+1)^3}$$

for $`1in1`$.