# Quantitative Convergence Analysis of Path Integral Representations for
Quantum Thermal Average

Xuda Ye Beijing International Center for Mathematical Research, Peking
University, Beijing, 100871, P. R. China. Email<EMAIL_ADDRESS>Zhennan
Zhou Beijing International Center for Mathematical Research, Peking
University, Beijing, 100871, P. R. China. Email<EMAIL_ADDRESS>
###### Abstract

The quantum thermal average is a central topic in quantum physics and can be
represented by the path integrals. For the computational perspective, the path
integral representation (PIR) needs to be approximated in a finite-dimensional
space, and the convergence of such approximation is termed as the convergence
of the PIR. In this paper, we establish the Trotter product formula in the
trace form, which connects the quantum thermal average and the Boltzmann
distribution of a continuous loop in a rigorous way. We prove the qualitative
convergence of the standard PIR, and obtain the explicit convergence rates of
the continuous loop PIR. These results showcase various approaches to
approximate the quantum thermal average, which provide theoretical guarantee
for the path integral approaches of quantum thermal equilibrium systems, such
as the path integral molecular dynamics.

Keywords quantum thermal average, path integral representation, Trotter
product formula
AMS subject classifications 82B31, 81S40

##  1 Introduction

The quantum thermal average stands as a pivotal concept within the realm of
quantum physics, serving not only to characterize the quantum canonical
ensemble comprehensively but also to find extensive utility in elucidating the
thermal properties exhibited by intricate quantum systems. These applications
encompass the ideal quantum gases [1], chemical reaction rates [2, 3], the
density of states in crystals [4] and the quantum phase transitions [5].
Nonetheless, closed-form expressions of such quantities are rarely available,
and the simulation cost of the direct discretization methods grow
exponentially with the spatial dimension. Therefore, the exact calculation of
the quantum thermal average in high dimensions can be difficult.

A transformative milestone arrived with the advent of Feynman’s path integral
[6], which provides a powerful approach to address the calculation in the
quantum physics. While its original form revolves around the real-time quantum
dynamics, Kac in 1947 made a breakthrough by conceiving the notion of the
imaginary-time path integral [7]. This innovation culminates in the
Feynman–Kac formula, an instrumental construct that captures the solution to
the parabolic and elliptic equations through the expectation of a stochastic
process. The success of the Feynman–Kac formula serves as a catalyst for the
evolution of the path integral representation (PIR) [8, 9, 10, 11], which
represents the quantum thermal average in the expectation of a continuous
loop.

Delving further, we introduce the PIR on a quantum Hamiltonian system in
$\mathbb{R}^{d}$:

$\hat{H}=\frac{\hat{p}^{2}}{2}+V(\hat{q}).$ (1.1)

Here, $\hat{q}$ and $\hat{p}$ denote position and momentum operators in
$\mathbb{R}^{d}$, and $V(q)$ is a real-valued potential function. When the
system exists at a constant temperature $T=1/\beta$, its state is described by
the canonical ensemble with the density operator $e^{-\beta\hat{H}}$, and thus
the partition function is given by
$\mathcal{Z}=\mathrm{Tr}[e^{-\beta\hat{H}}]$. Following the concept in [12,
13, 14], the quantum system is expressed as a continuous loop in the torus
$[0,\beta]$, denoted as $x(\tau)\in C([0,\beta];\mathbb{R}^{d})$. The central
to this approach is the energy function

$\mathcal{E}(x)=\int_{0}^{\beta}\bigg{[}\frac{1}{2}|x^{\prime}(\tau)|^{2}+V(x(\tau))\bigg{]}\mathrm{d}\tau,$
(1.2)

and the counterpart of the canonical ensemble $e^{-\beta\hat{H}}$ is the
formal Boltzmann distribution $\pi(x)\propto\exp(-\mathcal{E}(x))$. At the
core of the matter, the quantum thermal average is defined by the average of
the observable operator $O(\hat{q})$ in the canonical ensemble
$e^{-\beta\hat{H}}$, where $O(q)$ is real-valued function. The expression
takes shape as follows:

$\langle{O(\hat{q})}\rangle_{\beta}=\frac{\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]}{\mathrm{Tr}[e^{-\beta\hat{H}}]}.$
(1.3)

Meanwhile in the PIR, the quantum thermal average unravels through a different
form:

$\langle{O(\hat{q})}\rangle_{\beta}=\int\bigg{(}\frac{1}{\beta}\int_{0}^{\beta}O(x(\tau))\mathrm{d}\tau\bigg{)}\pi(x)\mathcal{D}[x],$
(1.4)

where $\mathcal{D}[x]$ embodies the formal Lebesgue measure in the function
space $C([0,\beta];\mathbb{R}^{d})$.

While the PIR stands as a potent tool for theoretical investigations, its
direct application to compute the quantum thermal average still presents
challenges. This arises from the intricate nature of the formal distribution
$\pi(x)$, which proves arduous to ascertain analytically or sample
numerically. Consequently, especially when confronted with high-dimensional
scenarios, the quest for approximation methods becomes imperative in the
pursuit of effectively computing the quantum thermal average. Within the scope
of this paper, we focus on two types of approximation methods.

  1. 1.

Standard path integral representation (std-PIR)
The std-PIR arises from the classical theory of the path integral Monte Carlo
(PIMC) [15, 16, 17] and the path integral molecular dynamics (PIMD) [18, 19,
20, 21], and stands as a prominent and widely adopted technique in
computational physics and theoretical chemistry for computing the quantum
thermal average. In the std-PIR, the continuous loop $x(\tau)$ undergoes the
approximation by utilizing its grid values, denoted as $x_{j}=x(j\beta_{D})$,
where the index $j=0,1,\cdots,D-1$. Here, $D\in\mathbb{N}$ signifies the
number of grid points, and $\beta_{D}=\beta/D$. Consequently, the energy
function $\mathcal{E}(x)$ in (1.2) can be discretized with the finite-
difference approximation. This gives rise to the the approximated quantum
thermal average of the std-PIR, designated as
$\langle{O(\hat{q})}\rangle_{\beta,D}^{\mathrm{std}}$.

  2. 2.

Continuous loop path integral representation (CL-PIR)
The CL-PIR is an innovation method introduced in our recent paper [22] to
calculate the quantum thermal average. In contrast to spatial coordinates, CL-
PIR embraces normal mode coordinates as its primary variables. By truncating
the number of normal modes to a finite integer $N\in\mathbb{N}$, we arrive at
the truncated CL-PIR, and the resulting approximated quantum thermal average
is denoted as $\langle{O(\hat{q})}\rangle_{\beta,N}$. Upon this foundation,
should we proceed to engage numerical integration on the grid values in the
truncated CL-PIR, we engender the discretized truncated CL-PIR. This fully
discretized representation yields an approximated quantum thermal average
denoted as $\langle{O(\hat{q})}\rangle_{\beta,N,D}$, where $D\in\mathbb{N}$
denotes the number of grid points.

Note that we have employed the symbol $D\in\mathbb{N}$ to denote the number of
grid points in the interval $[0,\beta]$ for both the std-PIR and the CL-PIR.
In the special case $N=D$, it can be shown that the CL-PIR is almost
equivalent to the std-PIR except for the coefficients in the normal mode
coordinates, see Section 2.3.

Notably, both the std-PIR and the CL-PIR pivot on the integer parameters, and
the accuracy in approximating the quantum thermal average hinges on the
progressive growth of these parameters. As such, a fundamental query surfaces:
how rapidly do the outcomes of these approximation methods converge towards
the true quantum thermal average? This paper’s principal contribution lies in
the establishment and evaluation of the convergence rates of the std-PIR and
the CL-PIR, which are stated as follows:

  1. 1.

$\displaystyle\lim_{D\rightarrow\infty}\langle{O(\hat{q})}\rangle_{\beta,D}^{\mathrm{std}}=\langle{O(\hat{q})}\rangle_{\beta}$
for the std-PIR;

  2. 2.

$\displaystyle\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N}\big{|}\lesssim\frac{1}{\sqrt{N}}$
for the truncated CL-PIR;

  3. 3.

$\displaystyle\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N,D}\big{|}\lesssim\frac{1}{\sqrt{N}}+\frac{1}{\sqrt{D}}$
for the discretized truncated CL-PIR.

These convergence results unveils a fundamental connection between the quantum
thermal average and the statistical average of the continuous loop in the PIR,
forming the bedrock of the mathematical underpinning for the PIMD. Even within
high-dimensional spaces, the promise remains that by judiciously selecting
sufficiently large integer parameters—$N$ and $D$—the accuracy of the quantum
thermal average approximation can be assured. In particular, the convergence
results of the CL-PIR quantitatively estimate the rates at which these
approximation methods converge towards the true quantum thermal average.

We provide a concise introduction to the mathematical tools in substantiating
the convergence results. Our proof framework begins by casting the continuous
loop $x(\tau)$ as a Gaussian stochastic process in the torus $[0,\beta]$.
Notably, we establish a significant link between the quantum thermal average
$\langle{O(\hat{q})}\rangle_{\beta}$ and the expectation calculated from this
continuous loop. Among the three convergence results described above, the
first result’s proof employs the Trotter product formula, a foundational
mechanism that sets forth the Feynman–Kac formula. Different from Kac’s
original form in [7], the Feynman–Kac formula in this paper is based on the
Section 3.2 of [23], which involves the expectation of the continuous path
with fixed endpoints. The second result’s validation centers around the
spectral structure intrinsic to the continuous loop. The third result’s proof
capitalizes on the Hölder continuity of the continuous loop.

The paper is organized as follows. In Section 2 we review the path integral
representations in this paper. In Section 2 we review the std-PIR and the CL-
PIR. In Section 3 we prove the convergence of these PIRs.

##  2 Review of the path integral representations

In this section we review the std-PIR and the CL-PIR aforementioned in the
introduction.

###  2.1 Standard path integral representation

In the std-PIR, a finite-difference approximation of the energy function
$\mathcal{E}(x)$ from (1.2) comes to the forefront. Let $D\in\mathbb{N}$ be
the number of grid points in $[0,\beta]$, and we employ the grid values
$x_{j}=x(j\beta_{D})$ of the continuous loop to represent the energy
$\mathcal{E}(x)$, where we presume $\beta_{D}=\beta/D$. Utilizing the finite-
difference approximation

$x^{\prime}(j\beta_{D})\approx\frac{x((j+1)\beta_{D})-x(j\beta_{D})}{\beta_{D}}=\frac{x_{j+1}-x_{j}}{\beta_{D}},~{}~{}~{}~{}j=0,1,\cdots,D-1,$
(2.1)

the energy function $\mathcal{E}_{D}(x)$ can be approximated as

$\mathcal{E}_{D}^{\mathrm{std}}(x)=\frac{1}{2\beta_{D}}\sum_{j=0}^{D-1}|x_{j}-x_{j+1}|^{2}+\beta_{D}\sum_{j=0}^{D-1}V(x_{j}),~{}~{}~{}~{}x\in\mathbb{R}^{dD}.$
(2.2)

Here, $\\{x_{j}\\}_{j=0}^{D-1}$ represent the grid values expected to align
with the continuous loop $x(\tau)$, characterized by the approximation:

$x_{j}\approx x(j\beta_{D}),~{}~{}~{}~{}j=0,1,\cdots,D-1.$ (2.3)

Consequently, the Boltzmann distribution linked to the energy function
$\mathcal{E}_{D}^{\mathrm{std}}(x)$ is

$\pi_{D}^{\mathrm{std}}(x)=\frac{1}{Z_{D}^{\mathrm{std}}}\exp\big(-\mathcal{E}_{D}^{\mathrm{std}}(x)\big{missing}),~{}~{}~{}~{}Z_{D}^{\mathrm{std}}=\int_{\mathbb{R}^{dD}}\exp\big(-\mathcal{E}^{\mathrm{std}}_{D}(x)\big{missing})\mathrm{d}x.$
(2.4)

And the quantum thermal average $\langle{O(\hat{q})}\rangle_{\beta}$ finds
approximation through the expression:

$\langle{O(\hat{q})}\rangle_{\beta}\approx\langle{O(\hat{q})}\rangle_{\beta,D}^{\mathrm{std}}=\int_{\mathbb{R}^{dD}}\bigg{(}\frac{1}{D}\sum_{j=0}^{D-1}O(x(j\beta_{D}))\bigg{)}\pi_{D}^{\mathrm{std}}(x)\mathrm{d}x.$
(2.5)

As $D\rightarrow\infty$, the formal continuum limit of the std-PIR (2.5) is
given by (1.4).

To compute the statistical average
$\langle{O(\hat{q})}\rangle_{\beta,D}^{\mathrm{std}}$ numerically, one usually
simulates an ergodic Langevin process whose invariant distribution is exactly
$\pi_{D}^{\mathrm{std}}(x)$, see [21, 22].

###  2.2 Continuous loop path integral representation

In the context of the CL-PIR, the continuous loop $x(\tau)$ finds
representation through normal mode coordinates. Beginning with the eigenvalue
problem associated with the second-order differential operator,

$-\ddot{c}_{k}(\tau)=\omega_{k}^{2}c_{k}(\tau),~{}~{}~{}~{}\tau\in[0,\beta],~{}~{}~{}~{}k=0,1,\cdots,$
(2.6)

the eigenvalues $\\{\omega_{k}\\}_{k=0}^{\infty}$ and the eigenfunctions
$\\{c_{k}(\tau)\\}_{k=0}^{\infty}$ manifest as

$\displaystyle\omega_{0}$ $\displaystyle=0,$ $\displaystyle c_{0}(\tau)$
$\displaystyle=\sqrt{\frac{1}{\beta}};$ (2.7) $\displaystyle\omega_{2k-1}$
$\displaystyle=\frac{2k\pi}{\beta},$ $\displaystyle c_{2k-1}(\tau)$
$\displaystyle=\sqrt{\frac{2}{\beta}}\sin\frac{2k\pi\tau}{\beta},~{}~{}~{}~{}k=1,2,\cdots;$
$\displaystyle\omega_{2k}$ $\displaystyle=\frac{2k\pi}{\beta},$ $\displaystyle
c_{2k}(\tau)$
$\displaystyle=\sqrt{\frac{2}{\beta}}\cos\frac{2k\pi\tau}{\beta},~{}~{}~{}~{}k=1,2,\cdots.$

Remarkably, the orthonormal basis formed by $\\{c_{k}(\tau)\\}_{k=0}^{\infty}$
spans the Hilbert space

$\mathbb{H}=L^{2}([0,\beta];\mathbb{R}^{d}).$ (2.8)

This foundation allows any continuous loop $x(\tau)$ to be uniquely expressed
in

$x(\tau)=\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau),~{}~{}~{}~{}\tau\in[0,\beta],$
(2.9)

where the coefficients $\\{\xi_{k}\\}_{k=0}^{\infty}\subset\mathbb{R}^{d}$ are
referred to as the normal mode coordinates. Consequently, the energy function
$\mathcal{E}(x)$ in (1.2) takes a precise form:

$\displaystyle\mathcal{E}(\xi)$
$\displaystyle=-\frac{1}{2}\int_{0}^{\beta}x(\tau)^{\mathrm{T}}\ddot{x}(\tau)\mathrm{d}\tau+\int_{0}^{\beta}V(x(\tau))\mathrm{d}\tau$
$\displaystyle=\frac{1}{2}\int_{0}^{\beta}\bigg{(}\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{)}^{\mathrm{T}}\bigg{(}\sum_{k=0}^{\infty}\omega_{k}^{2}\xi_{k}c_{k}(\tau)\bigg{)}\mathrm{d}\tau+\int_{0}^{\beta}V(x(\tau))\mathrm{d}\tau$
$\displaystyle=\frac{1}{2}\sum_{k=0}^{\infty}\omega_{k}^{2}|\xi_{k}|^{2}+\int_{0}^{\beta}V(x(\tau))\mathrm{d}\tau.$
(2.10)

For the convenience of analysis, we introduce the constant $a>0$ and rewrite
the energy function $\mathcal{E}(\xi)$ as

$\mathcal{E}(\xi)=\frac{1}{2}\sum_{k=0}^{\infty}(\omega_{k}^{2}+a^{2})|\xi_{k}|^{2}+\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau,$
(2.11)

where the potential function

$V^{a}(q)=V(q)-\frac{a^{2}}{2}|q|^{2}.$ (2.12)

Note that the constant $a>0$ ensures that the coefficient
$\omega_{k}^{2}+a^{2}$ in each normal mode is strictly positive.

To achieve a finite-dimensional approximation of the quantum thermal average
$\langle{O(\hat{q})}\rangle_{\beta}$, we introduce a finite integer parameter
$N\in\mathbb{N}$ to indicate the number of normal modes, and truncate the
continuous loop $x(\tau)$ to be

$x_{N}(\tau)=\sum_{k=0}^{N-1}\xi_{k}c_{k}(\tau).$ (2.13)

This truncation leads to the corresponding energy function

$\mathcal{E}_{N}(\xi)=\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k}^{2}+a^{2})|\xi_{k}|^{2}+\int_{0}^{\beta}V^{a}(x_{N}(\tau))\mathrm{d}\tau,~{}~{}~{}~{}\xi\in\mathbb{R}^{dN},$
(2.14)

and the Boltzmann distribution $\pi_{N}(\xi)$ in $\mathbb{R}^{dN}$ assumes the
form

$\pi_{N}(\xi)=\frac{1}{Z_{N}}\exp(-\mathcal{E}_{N}(\xi)),~{}~{}~{}~{}Z_{N}=\int_{\mathbb{R}^{dN}}\exp(-\mathcal{E}_{N}(\xi))\mathrm{d}\xi.$
(2.15)

Then the quantum thermal average $\langle{O(\hat{q})}\rangle_{\beta}$ finds
its approximation through

$\langle{O(\hat{q})}\rangle_{\beta}\approx\langle{O(\hat{q})}\rangle_{\beta,N}=\int_{\mathbb{R}^{dN}}\bigg{(}\frac{1}{\beta}\int_{0}^{\beta}O(x_{N}(\tau))\mathrm{d}\tau\bigg{)}\pi_{N}(\xi)\mathrm{d}\xi.$
(2.16)

The approximation (2.16) is referred to as the truncated CL-PIR.

While the calculation of $\langle{O(\hat{q})}\rangle_{\beta,N}$ can be
implemented with a finite-dimensional distribution, the CL-PIR presents
challenges in numerical computation because of the inconvenience to evaluate
the integrals

$\int_{0}^{\beta}V^{a}(x_{N}(\tau))\mathrm{d}\tau\mbox{~{}~{}and~{}~{}}\int_{0}^{\beta}O(x_{N}(\tau))\mathrm{d}\tau$

analytically. As a consequence, we seek for numerical integration techniques
to approximate these integrals. Let $D\in\mathbb{N}$ be the number of grid
points in the torus $[0,\beta]$, and define $\beta_{D}=\beta/D$, then these
integral terms are approximated by

$\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D}))\mbox{~{}~{}and~{}~{}}\beta_{D}\sum_{j=0}^{D-1}O(x_{N}(j\beta_{D})).$
(2.17)

The outcome is a discretized energy function

$\mathcal{E}_{N,D}(\xi)=\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k}^{2}+a^{2})|\xi_{k}|^{2}+\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D})),~{}~{}~{}~{}\xi\in\mathbb{R}^{dN}.$
(2.18)

By defining the corresponding Boltzmann distribution

$\pi_{N,D}(\xi)=\frac{1}{Z_{N,D}}\exp(-\mathcal{E}_{N,D}(\xi)),~{}~{}~{}~{}Z_{N,D}=\int_{\mathbb{R}^{dD}}\exp(-\mathcal{E}_{N,D}(\xi))\mathrm{d}\xi,$
(2.19)

the approximated quantum thermal average takes the form

$\langle{O(\hat{q})}\rangle_{\beta,N,D}=\int_{\mathbb{R}^{dN}}\bigg{(}\frac{1}{D}\sum_{j=0}^{D-1}O(x_{N}(j\beta_{D}))\bigg{)}\pi_{N,D}(\xi)\mathrm{d}\xi.$
(2.20)

The approximation (2.20) is referred to as the discretized truncated CL-PIR.

###  2.3 Relation between std-PIR and discretized truncated CL-PIR

We study the relation between the std-PIR and the discretized truncated PIR in
the case $N=D$, i.e., the number of normal modes $N$ is chosen to be the same
as the number of grid points $D$. In the discrerized truncated CL-PIR, the
potential function $\mathcal{E}_{D,D}(\xi)$ is given by

$\mathcal{E}_{D,D}(\xi)=\frac{1}{2}\sum_{k=0}^{D-1}(\omega_{k}^{2}+a^{2})|\xi_{k}|^{2}+\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{D}(j\beta_{D})),$
(2.21)

where $x_{D}(\tau)$ is the continuous loop

$x_{D}(\tau)=\sum_{k=0}^{D-1}\xi_{k}c_{k}(\tau),~{}~{}~{}~{}\tau\in[0,\beta].$
(2.22)

The $D$ grid values of the continuous loop $x_{D}(\tau)$ are given by

$x_{j}=x_{D}(j\beta_{D})=\sum_{k=0}^{D-1}\xi_{k}c_{k}(j\beta_{D}),~{}~{}~{}~{}j=0,1,\cdots,D-1.$
(2.23)

Then for given index $k\in\\{0,1,\cdots,N-1\\}$, we have

$\beta_{D}\sum_{j=0}^{D-1}x_{j}c_{k}(j\beta_{D})=\sum_{l=0}^{D-1}\xi_{l}\bigg{(}\beta_{D}\sum_{j=0}^{D-1}c_{l}(j\beta_{D})c_{k}(j\beta_{D})\bigg{)}=\sum_{l=0}^{D-1}\xi_{l}\delta_{lk}=\xi_{k},$
(2.24)

hence the coefficients $\\{\xi_{k}\\}_{k=0}^{D-1}$ can be represented by the
grid values $\\{x_{j}\\}_{j=0}^{D-1}$ via

$\xi_{k}=\beta_{D}\sum_{j=0}^{D-1}x_{j}c_{k}(j\beta_{D}),~{}~{}~{}~{}k=0,1,\cdots,D-1.$
(2.25)

Then we can represent the potential function
$\mathcal{E}_{D}^{\mathrm{std}}(x)$ in $\\{\xi_{k}\\}_{k=0}^{D-1}$:

$\displaystyle\mathcal{E}_{D}^{\mathrm{std}}(\xi)$
$\displaystyle=\frac{1}{2\beta_{D}}\sum_{j=0}^{D-1}|x_{j}-x_{j+1}|^{2}+\beta_{D}\sum_{j=0}^{D-1}V(x_{j})$
$\displaystyle=\frac{1}{2\beta_{D}}\sum_{j=0}^{D-1}|x_{j}-x_{j+1}|^{2}+\frac{a^{2}\beta_{D}}{2}\sum_{j=0}^{D-1}|x_{j}|^{2}+\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{j})$
$\displaystyle=\frac{1}{2}\sum_{k=0}^{D-1}(\omega_{k,D}^{2}+a^{2})|\xi_{k}|^{2}+\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{D}(j\beta_{D})),$
(2.26)

where the coefficients $\\{\omega_{k,D}\\}_{k=0}^{D-1}$ are given by

$\omega_{0,D}=0,~{}~{}~{}~{}\omega_{2k-1,D}=\omega_{2k,D}=2\sin\frac{k\pi}{\beta},~{}~{}~{}~{}k=1,2,\cdots.$
(2.27)

The energy function $\mathcal{E}_{D}^{\mathrm{std}}(\xi)$ is almost equivalent
to $\mathcal{E}_{D,D}(\xi)$ given in (2.21), except for the difference on the
coefficients $\omega_{k}$ and $\omega_{k,D}$. Furthermore, there holds

$\lim_{D\rightarrow\infty}\omega_{k,D}=\omega_{k},$ (2.28)

and thus the std-PIR and the discretized truncated CL-PIR are very close when
$D$ is large.

##  3 Convergence analysis of path integral representations

In this section we prove the convergence results of the std-PIR and the CL-
PIR. For the convenience, we list the the assumptions on the potential
function $V(q)$ and the observable function $O(q)$ as follows.

######  Assumption.

Given $a>0$, the potential function $V^{a}(q)=V(q)-a^{2}|q|^{2}/2$ satisfies

  1. (i)

$|V^{a}(0)|\leqslant M_{1}$, $V^{a}(q)\geqslant-M_{1}$, $|\nabla
V^{a}(q)|\leqslant M_{1}+M_{1}|q|$,

and the observable function $O(q)$ satisfies

  1. (ii)

$\max\\{|O(q)|,|\nabla O(q)|\\}\leqslant M_{2}$.

Using the fundamental theorem of calculus,

$V^{a}(q)=V^{a}(0)+q\cdot\int_{0}^{1}\nabla V^{a}(\theta q)\mathrm{d}\theta,$
(3.1)

then Assumption (i) implies

$\displaystyle V^{a}(q)$
$\displaystyle\leqslant|V^{a}(0)|+|q|\int_{0}^{1}|\nabla V^{a}(\theta
q)|\mathrm{d}q$ $\displaystyle\leqslant
M_{1}+|q|\int_{0}^{1}\big{(}M_{1}+\theta M_{1}|q|\big{)}\mathrm{d}\theta$
$\displaystyle=M_{1}+M_{1}|q|+\frac{M_{1}}{2}|q|^{2}\leqslant\frac{3}{2}M_{1}+M_{1}|q|^{2}.$
(3.2)

This section is organized as below. In Section 3.1, we study the Trotter
product formula, which is the key ingredient in the convergence analysis of
both the std-PIR and the CL-PIR. In Section 3.2, we prove the convergence of
the std-PIR. In Section 3.3, we validate the CL-PIR produces the accurate
quantum thermal average. In Section 3.4, we quantify the convergence rate of
the truncated CL-PIR. In Section 3.5, we quantify the convergence rate of the
discretized truncated CL-PIR.

###  3.1 Discussion on the Trotter product formula

Before delving into the details of the convergence analysis, we briefly
discuss on a key ingredient of the proof: the Trotter product formula. For the
convenience, we define the free particle Schrödinger operator $\hat{H}^{0}$
and the potential operator $\hat{V}$ by

$\hat{H}^{0}=\frac{\hat{p}^{2}}{2}=-\frac{1}{2}\Delta_{d},~{}~{}~{}~{}\hat{V}=V(\hat{q}),$
(3.3)

where $\Delta_{d}$ is the Laplace operator in $\mathbb{R}^{d}$. For the
constant $a>0$, define the quantum harmonic oscillator $\hat{H}^{a}$ and the
potential operator $\hat{V}^{a}$ by

$\hat{H}^{a}=\frac{\hat{p}^{2}}{2}+\frac{a^{2}}{2}\hat{q}^{2}=-\frac{1}{2}\Delta_{d}+\frac{a^{2}}{2}\sum_{i=1}^{d}\hat{q}_{i}^{2}.$
(3.4)

It is easy to observe the Hamiltonian $\hat{H}$ can be written as

$\hat{H}=\hat{H}^{0}+\hat{V}=\hat{H}^{a}+\hat{V}^{a}.$ (3.5)

The Trotter product formula is stated as

$\mathrm{Tr}[e^{-\beta\hat{H}}]=\lim_{D\rightarrow\infty}\mathrm{Tr}\big{[}\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}\big{]}=\lim_{D\rightarrow\infty}\mathrm{Tr}\big{[}\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}\big{]},$
(3.6)

where the LHS is the partition function appearing in the quantum thermal
average (1.3), while $(e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}})^{D}$
and $(e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}})^{D}$ in the RHS
correspond to the derivation of the std-PIR and the CL-PIR, respectively. If
we extract the kernel function from (3.6), the Trotter product formula is
given by

$\matrixelement{q}{e^{-\beta\hat{H}}}{q}=\lim_{D\rightarrow\infty}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}=\lim_{D\rightarrow\infty}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{q}.$
(3.7)

Here, we use the bra–ket notation to represent the kernel function, see
Chapter 3 of [23].

The Trotter product formula (3.6) is quite easy to understand, but the proof
is nontrivial. Most literature discuss the strong convergence of the operators
[23, 24], i.e.,

$\lim_{D\rightarrow\infty}\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}=e^{-\beta\hat{H}},~{}~{}~{}~{}\lim_{D\rightarrow\infty}\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}=e^{-\beta\hat{H}}$
(3.8)

in the strong sense, and the convergence in the kernel functions are rarely
mentioned. Therefore, we rely on the Feynman–Kac formula derived in Section
3.2 of [23], which represents the kernel functions in a special Wiener
measure. The details of the Feynman–Kac formula and the Wiener measure are
given in the proof of Lemma 3.1 and Lemma 3.3 in this paper.

###  3.2 Standard path integral representation

The convergence analysis of the std-PIR relies on the following Trotter
product formula.

######  Lemma 3.1.

Under Assumption (i), for any $q\in\mathbb{R}^{d}$, we have

$\lim_{D\rightarrow\infty}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}=\matrixelement{q}{e^{-\beta\hat{H}}}{q}.$
(3.9)

Also, there exist constants $A,\lambda>0$ independent of $D$ and $q$ such that

$\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}\leqslant
A\exp(-\lambda|q|^{2}).$ (3.10)

The proof is given in Appendix A. Note that (3.10) guarantees the exponential
decay of the kernel function, which allows the usage of the dominated
convergence theorem.

Now we state the main theorem.

######  Theorem 3.1.

Under Assumptions (i)(ii), we have

$\lim_{D\rightarrow\infty}\langle{O(\hat{q})}\rangle_{\beta,D}^{\mathrm{std}}=\langle{O(\hat{q})}\rangle_{\beta}.$
(3.11)

The proof of Theorem 3.1 mainly utilizes the equality

$\langle{O(\hat{q})}\rangle_{\beta,D}^{\mathrm{std}}=\frac{\mathrm{Tr}\big{[}\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}O(\hat{q})\big{]}}{\mathrm{Tr}\big{[}(e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}})^{D}\big{]}},$
(3.12)

which holds true for any integer $D\in\mathbb{N}$.

######  Proof.

With Lemma 3.1, we can apply the dominated convergence theorem to derive

$\lim_{D\rightarrow\infty}\int_{\mathbb{R}^{d}}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}O(q)\mathrm{d}q=\int_{\mathbb{R}^{d}}\matrixelement{q}{e^{-\beta\hat{H}}}{q}O(q)\mathrm{d}q,$
(3.13)

which is exactly

$\lim_{D\rightarrow\infty}\mathrm{Tr}\big{[}\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}O(\hat{q})\big{]}=\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})].$
(3.14)

Similarly, by choosing $O(q)\equiv 1$ in (3.14) we have

$\lim_{D\rightarrow\infty}\mathrm{Tr}\big{[}\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}\big{]}=\mathrm{Tr}[e^{-\beta\hat{H}}].$
(3.15)

By inserting the free positions $\\{x_{j}\\}_{j=0}^{D-1}$ in $\mathbb{R}^{d}$,
we can write the LHS of (3.14) as

$\displaystyle~{}~{}~{}~{}\mathrm{Tr}\big{[}\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}O(\hat{q})\big{]}$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathrm{d}x_{0}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{D-1}\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}}{x_{j+1}}O(x_{0})$
$\displaystyle=\frac{1}{(2\pi\beta_{D})^{\frac{dD}{2}}}\int_{\mathbb{R}^{d}}\mathrm{d}x_{0}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{D-1}\exp\bigg(-\frac{1}{2\beta_{D}}\sum_{j=0}^{D-1}|x_{j}-x_{j+1}|^{2}-\beta_{D}\sum_{j=0}^{D-1}V(x_{j})\bigg{missing})O(x_{0})$
$\displaystyle=\frac{1}{(2\pi\beta_{D})^{\frac{dD}{2}}}\int_{\mathbb{R}^{dD}}O(x_{0})\exp\big(-\mathcal{E}_{D}^{\mathrm{std}}(x)\big{missing})\mathrm{d}x.$
(3.16)

Using the symmetry of the expression (3.16), we can write

$\mathrm{Tr}\big{[}\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}O(\hat{q})\big{]}=\frac{1}{(2\pi\beta_{D})^{\frac{dD}{2}}}\int_{\mathbb{R}^{dD}}\bigg{(}\frac{1}{D}\sum_{j=0}^{D-1}O(x_{j})\bigg{)}\exp\big(-\mathcal{E}_{D}^{\mathrm{std}}(x)\big{missing})\mathrm{d}x.$
(3.17)

Similarly,

$\mathrm{Tr}\big{[}(e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}})^{D}\big{]}=\frac{1}{(2\pi\beta_{D})^{\frac{dD}{2}}}\int_{\mathbb{R}^{dD}}\exp\big(-\mathcal{E}_{D}^{\mathrm{std}}(x)\big{missing})\mathrm{d}x.$
(3.18)

Dividing (3.17) by (3.18), we obtain

$\frac{\mathrm{Tr}\big{[}\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}O(\hat{q})\big{]}}{\mathrm{Tr}\big{[}(e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}})^{D}\big{]}}=\int_{\mathbb{R}^{dD}}\bigg{(}\frac{1}{D}\sum_{j=0}^{D-1}O(x_{j})\bigg{)}\pi_{D}^{\mathrm{std}}(x)\mathrm{d}x=\langle{O(\hat{q})}\rangle_{\beta,D}^{\mathrm{std}}.$
(3.19)

Let the number of grid points $D$ tend to infinity, from (3.14) and (3.15) we
obtain

$\lim_{D\rightarrow\infty}\langle{O(\hat{q})}\rangle_{\beta,D}^{\mathrm{std}}=\lim_{D\rightarrow\infty}\frac{\mathrm{Tr}\big{[}\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}O(\hat{q})\big{]}}{\mathrm{Tr}\big{[}(e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}})^{D}\big{]}}=\frac{\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]}{\mathrm{Tr}[e^{-\beta\hat{H}}]}=\langle{O(\hat{q})}\rangle_{\beta}.$
(3.20)

$\square$

######  Remark.

Although Thoerem 3.1 guarantees the convergence of the std-PIR, the
quantification of
$\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,D}^{\mathrm{std}}\big{|}$
in terms of the number of grid points $D$ is still unknown.

###  3.3 Continuous loop path integral representation

To begin with, we prove the Hölder continuity of the continuous loop
$x(\tau)$. Recall that any continuous loop $x(\tau)\in\mathbb{H}$ can be
written in

$x(\tau)=\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau).$ (3.21)

Consider the Gaussian distribution $\nu_{0}$ of the normal mode coordinates
$\\{\xi_{k}\\}_{k=0}^{\infty}$ given as

$\xi_{k}\sim\mathcal{N}\bigg{(}0,\frac{I_{d}}{\omega_{k}^{2}+a^{2}}\bigg{)},~{}~{}~{}~{}k=0,1,\cdots,$
(3.22)

then we define $\nu$ to be the pushforward of the distribution $\nu_{0}$ in
the continuous loop mapping (3.21). It is clear that $\nu$ is a Gaussian
distribution in the Hilbert space
$\mathbb{H}=L^{2}([0,\beta];\mathbb{R}^{d})$. Note that the first eigenvalue
$\omega_{0}=0$, and the introduction of the constant $a>0$ ensures the well-
posedness of Gaussian distribution $\nu$. Now we study the properties of the
random continuous loop $x(\tau)$ in the distribution $\nu$.

######  Lemma 3.2.

The random continuous loop $x(\tau)$ with the distribution $\nu$ satisfies

$\mathbb{E}_{\nu}\bigg{[}\int_{0}^{\beta}|x(\tau)|^{2}\mathrm{d}\tau\bigg{]}=C_{0},$
(3.23)

where the constant $C_{0}=\frac{d\beta}{2a}\coth\frac{a\beta}{2}$. For any
$\tau_{1},\tau_{2}\in[0,\beta]$,

$\mathbb{E}_{\nu}\big{|}x(\tau_{1})-x(\tau_{2})\big{|}^{2}\leqslant
d(2\beta+1)|\tau_{1}-\tau_{2}|.$ (3.24)

For any constant $\gamma\in(0,\frac{1}{2})$, $x(\tau)$ is $\gamma$-Hölder
continuous in the torus $[0,\beta]$ almost surely.

The proof is given in Appendix A.

######  Remark.

The Hölder continuity of the continuous loop $x(\tau)$ implies the regularity
of $x(\tau)$ is the same as the standard Brownian process.

Using the Gaussian distribution $\nu$, we can interpret the formal Boltzmann
distribution $\pi(\xi)\propto\exp(-\mathcal{E}(\xi))$ as a probability
distribution in the Hilbert space $\mathbb{H}=L^{2}([0,\beta];\mathbb{R}^{d})$
defined by the Randon–Nikodym derivative

$\frac{\mathrm{d}\pi}{\mathrm{d}\nu}(\xi)=\frac{1}{Z}\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{missing}),$
(3.25)

where the constant $Z$ is the normalization constant defined by

$Z=\mathbb{E}_{\nu}\bigg{[}\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\bigg{missing})\mathrm{d}\tau\bigg{]}.$
(3.26)

Therefore, the PIR in (1.4) can be interpreted as the statiscal average

$\langle{O(\hat{q})}\rangle_{\beta}=\mathbb{E}_{\pi}\bigg{[}\frac{1}{\beta}\int_{0}^{\beta}O(x(\tau))\mathrm{d}\tau\bigg{]}.$
(3.27)

######  Remark.

Although the distribution $\nu$ depends on the parameter $a>0$, the
distribution $\pi$ does not depend on the parameter $a$. This is because
formally $\pi(\xi)\propto\exp(-\mathcal{E}(\xi))$, where the energy function
$\mathcal{E}(\xi)$ does not depend on the parameter $a$.

Using the Radon–Nikodym derivative in (3.25), we can equivalently rewrite the
CL-PIR (3.27) in the following result.

######  Theorem 3.2.

Under Assumption (i)(ii), we have

$\langle{O(\hat{q})}\rangle_{\beta}=\frac{\displaystyle\mathbb{E}_{\nu}\bigg{[}\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{missing})\times\bigg{(}\frac{1}{\beta}\int_{0}^{\beta}O(x(\tau))\mathrm{d}\tau\bigg{)}\bigg{]}}{\displaystyle\mathbb{E}_{\nu}\bigg{[}\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{missing})\bigg{]}}.$
(3.28)

Although the derivation of (3.28) in the paragraph above is natural, the
rigorous verification of Theorem 3.2 requires careful arguments on the Trotter
product formula. This is also the case for the proof of Theorem 3.1. To prove
Theorem 3.2, we state the following Trotter product formula, which is an
analogue of Lemma 3.1.

######  Lemma 3.3.

Under Assumption (i), for any $q,\tilde{q}\in\mathbb{R}^{d}$, we have

$\lim_{D\rightarrow\infty}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}}=\matrixelement{q}{e^{-\beta\hat{H}}}{\tilde{q}}.$
(3.29)

Also, there exist constants $A,\lambda>0$ independent of $D$ and $q$ such that

$\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{q}\leqslant
A\exp(-\lambda|q|^{2}).$ (3.30)

The proof is given in Appendix A. Now we present the proof of Theorem 3.2.

######  Proof.

The proof is accomplished in several steps.
1\. Simplification of the result
We claim that we only need to prove

$\frac{\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]}{\mathrm{Tr}[e^{-\beta\hat{H}^{a}}]}=\mathbb{E}_{\nu}\bigg{[}\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{missing})\times\bigg{(}\frac{1}{\beta}\int_{0}^{\beta}O(x(\tau))\mathrm{d}\tau\bigg{)}\bigg{]}.$
(3.31)

In particular, by choosing $O(q)\equiv 1$ in (3.31), we obtain

$\frac{\mathrm{Tr}[e^{-\beta\hat{H}}]}{\mathrm{Tr}[e^{-\beta\hat{H}^{a}}]}=\mathbb{E}_{\nu}\bigg{[}\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{missing})\bigg{]}$
(3.32)

Combining (3.31) and (3.32), we immediatelly obtain (3.28). Therefore, we
focus on the proof of (3.31). As a consequence of the uniform-in-$D$ bound in
Lemma 3.3, we can apply the dominated convergence theorem to deduce

$\lim_{D\rightarrow\infty}\int_{\mathbb{R}^{d}}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{q}O(q)\mathrm{d}q=\int_{\mathbb{R}^{d}}\matrixelement{q}{e^{-\beta\hat{H}}}{q}O(q)\mathrm{d}q,$
(3.33)

which is exactly equivalent to the Trotter product formula

$\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]=\lim_{D\rightarrow\infty}\mathrm{Tr}\big{[}\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}O(\hat{q})\big{]}.$
(3.34)

2\. Expansion of $\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]$ in the ring
polymer distribution
Using the Trotter product formula (3.34), we can conveniently approximate
$\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]$ in the ring polymer distribution.

$\displaystyle\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]$
$\displaystyle=\lim_{D\rightarrow\infty}\mathrm{Tr}\big{[}\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}O(\hat{q})\big{]}$
$\displaystyle=\lim_{D\rightarrow\infty}\int_{\mathbb{R}^{dD}}\mathrm{d}x\,\exp\bigg(-\beta_{D}\sum_{j=0}^{D-1}V(x_{j+1})\bigg{missing})O(x_{0})\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{a}}}{x_{j+1}}$

Using the symmetry of the expression in $O(x_{0})$, we obtain

$\displaystyle\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]=\lim_{D\rightarrow\infty}$
$\displaystyle\int_{\mathbb{R}^{dD}}\mathrm{d}x\,\exp\bigg(-\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{j})\bigg{missing})\,\times$
(3.35)
$\displaystyle\bigg{(}\frac{1}{D}\sum_{j=0}^{D-1}O(x_{j})\bigg{)}\times\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{a}}}{x_{j+1}}.$

Motivated by (3.35), we define the probability distribution of the ring
polymer by

$\Theta_{D}(x)=\frac{1}{Z_{D}}\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{a}}}{x_{j+1}},~{}~{}~{}~{}x\in\mathbb{R}^{dD}.$
(3.36)

where each $\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{a}}}{x_{j+1}}$ is
given by the Mehler kernel in Lemma A.1, and $Z_{D}$ is the normalization
constant given by

$Z_{D}=\int_{\mathbb{R}^{dD}}\mathrm{d}x\,\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{a}}}{x_{j+1}}.$
(3.37)

If we choose $V^{a}(q)\equiv 0$ and $O(q)\equiv 1$ in (3.35), then

$\mathrm{Tr}[e^{-\beta\hat{H}^{a}}]=\lim_{D\rightarrow\infty}\int_{\mathbb{R}^{dD}}\mathrm{d}x\,\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{a}}}{x_{j+1}}.$
(3.38)

Dividing (3.35) by (3.38), we obtain

$\frac{\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]}{\mathrm{Tr}[e^{-\beta\hat{H}^{a}}]}=\lim_{D\rightarrow\infty}\int_{\mathbb{R}^{dD}}\exp\bigg(-\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{j})\bigg{missing})\bigg{(}\frac{1}{D}\sum_{j=0}^{D-1}O(x_{j})\bigg{)}\Theta(x)\mathrm{d}x.$
(3.39)

In conclusion, we show that $\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]$ can be
accurately approximated by the statistical average of the ring polymer
distribution with $D$ grid points.
3\. Equivalence between the Gaussian distributions
We define the Gaussian distribution $\tilde{U}^{\beta}$ in the Hilbert space
$\mathbb{H}=L^{2}([0,\beta];\mathbb{R}^{d})$ by the following rule: for any
constants $0\leqslant\tau_{0}<\tau_{1}<\cdots<\tau_{P-1}<\beta$, the joint
distribution of the random variables
$x_{0}=x(\tau_{0}),x_{1}=x(\tau_{1}),\cdots,x_{P-1}=x(\tau_{P-1})$ is
proportional to

$\prod_{j=0}^{P-2}\matrixelement{x_{j}}{e^{-(\tau_{j+1}-\tau_{j})\hat{H}^{a}}}{x_{j+1}}\times\matrixelement{x_{P-1}}{e^{-(\tau_{0}-\tau_{P-1}+\beta)}}{x_{0}},$
(3.40)

which is the product of the Mehler kernels of the $n$ adjacent pairs in the
position coordinates $x(\tau_{0}),x(\tau_{1}),\cdots,x(\tau_{P-1})$. From the
Kolmogorov extension theorem, $\tilde{U}^{\beta}$ is indeed a well-defined
Gaussian distribution in $\mathbb{H}$. The difference between the Gaussian
distribution $\tilde{U}^{\beta}$ and the Wiener measure $U_{q,q}^{\beta}$
defined in (A.36) is that the endpoints of the continuous loop in
$\tilde{U}^{\beta}$ is flexible, while the endpoints in $U_{q,q}^{\beta}$ are
fixed at $q$.

Upon the definition of the Gaussian distribution $\tilde{U}^{\beta}$, the
distribution $\Theta_{P}(x_{1},\cdots,x_{P-1})$ can be viewed as the joint
distribution of the $D$ grid points $\\{x(j\beta_{D})\\}_{j=0}^{D-1}$. As a
consequence, we can rewrite (3.39) as

$\frac{\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]}{\mathrm{Tr}[e^{-\beta\hat{H}^{a}}]}=\lim_{D\rightarrow\infty}\int_{\mathbb{H}}\exp\bigg(-\beta_{D}\sum_{j=0}^{D-1}V^{a}(x(j\beta_{D}))\bigg{missing})\bigg{(}\frac{1}{D}\sum_{j=0}^{D-1}O(x(j\beta_{D}))\bigg{)}\mathrm{d}\tilde{U}^{\beta}.$
(3.41)

Then as the number of grid points $D\rightarrow\infty$, we can apply the
dominated convergence theorem on (3.41) to deduce

$\frac{\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]}{\mathrm{Tr}[e^{-\beta\hat{H}^{a}}]}=\int_{\mathbb{H}}\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{missing})\times\bigg{(}\frac{1}{\beta}\int_{0}^{\beta}O(x(\tau))\mathrm{d}\tau\bigg{)}\times\mathrm{d}\tilde{U}^{\beta}.$
(3.42)

The final step of the proof is to verify that the Gaussian distribution
$\tilde{U}^{\beta}$ defined in (3.40) and the distribution $\nu$ defined in
(3.22) are the same. Recall that the distribution $\nu$ is defined using the
normal mode coordinates,

$x(\tau)=\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau),~{}~{}~{}~{}\xi_{k}\sim\mathcal{N}\bigg{(}0,\frac{I_{d}}{\omega_{k}^{2}+a^{2}}\bigg{)},~{}~{}~{}~{}k=0,1,2,\cdots.$

Since both $\tilde{U}^{\beta}$ and $\nu$ are zero-mean Gaussian processes in
$[0,\beta]$, we only need to check their covariance functions are the same.

  1. 1.

In the distribution $\tilde{U}^{\beta}$, the joint distribution of $x=x(0)$
and $y=x(\tau)$ is given by the Mehler kernel introduced in Lemma A.1,

$\matrixelement{x}{e^{-\tau\hat{H}^{a}}}{y}\matrixelement{y}{e^{-(\beta-\tau)\hat{H}^{a}}}{x}=\bigg{(}\frac{a}{2\pi\sinh(a\beta)}\bigg{)}^{d}\exp\bigg(-\frac{A}{2}(|x^{2}|+|y|^{2})+Bx^{\mathrm{T}}y\bigg{missing}),$

where the constants $A$ and $B$ are given ny

$A=\frac{a}{\tanh(a\tau)}+\frac{a}{\tanh(a(\beta-\tau))},~{}~{}~{}~{}B=\frac{a}{\sinh(a\tau)}+\frac{a}{\sinh(a(\beta-\tau))}.$
(3.43)

From the Gaussian distribution of $(x,y)$ in
$\mathbb{R}^{d}\times\mathbb{R}^{d}$, the covariance function is

$\mathbb{E}_{\tilde{U}^{\beta}}\big{[}x(0)x(\tau)\big{]}=\frac{B^{2}}{A^{2}-B^{2}}.$
(3.44)

  2. 2.

In the distribution $\nu$, the covariance function can be calculated as

$\displaystyle\mathbb{E}_{\nu}\big{[}x(0)x(\tau)\big{]}$
$\displaystyle=\mathbb{E}_{\nu}\bigg{[}\sum_{k=0}^{\infty}\xi_{k}c_{k}(0)\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{]}=\sum_{k=0}^{\infty}\frac{c_{k}(0)c_{k}(\tau)}{\omega_{k}^{2}+a^{2}}$
$\displaystyle=\frac{1}{\beta}\cdot\frac{1}{a^{2}}+\frac{2}{\beta}\sum_{k=1}^{\infty}\frac{\cos\frac{2k\pi\tau}{\beta}}{(\frac{2k\pi}{\beta})^{2}+a^{2}}=\frac{1}{\beta}\sum_{k\in\mathbb{Z}}\frac{\cos\frac{2k\pi\tau}{\beta}}{\omega_{k}^{2}+a^{2}}.$
(3.45)

Note that the RHS of (3.45) can be explicitly calculated111See the answer in
https://math.stackexchange.com/a/4725694/402582., we can verify that
$\mathbb{E}_{\tilde{U}^{\beta}}\big{[}x(0)x(\tau)\big{]}$ and
$\mathbb{E}_{\nu}\big{[}x(0)x(\tau)\big{]}$ are exactly the same, which
implies the Gaussian processes $\tilde{U}^{\beta}$ and $\nu$ have the same
covariance function. Therefore, $\tilde{U}^{\beta}$ are $\nu$ are the same
distribution, and (3.42) directly yields (3.31). $\square$

###  3.4 Truncated continuous loop path integral representation

Similar to the expression in (3.28), we can write the statistical average
$\langle{O(\hat{q})}\rangle_{\beta,N}$ as

$\langle{O(\hat{q})}\rangle_{\beta,N}=\frac{\displaystyle\mathbb{E}_{\nu}\bigg{[}\exp\bigg(-\int_{0}^{\beta}V^{a}(x_{N}(\tau))\mathrm{d}\tau\bigg{missing})\times\bigg{(}\frac{1}{\beta}\int_{0}^{\beta}O(x_{N}(\tau))\mathrm{d}\tau\bigg{)}\bigg{]}}{\displaystyle\mathbb{E}_{\nu}\bigg{[}\exp\bigg(-\int_{0}^{\beta}V^{a}(x_{N}(\tau))\mathrm{d}\tau\bigg{missing})\bigg{]}},$
(3.46)

where $x_{N}(\tau)$ is the truncated continuous loop

$x_{N}(\tau)=\sum_{k=0}^{N-1}\xi_{k}c_{k}(\tau),~{}~{}~{}~{}\tau\in[0,\beta].$
(3.47)

From (3.28) and (2.16), we observe that the difference between
$\langle{O(\hat{q})}\rangle_{\beta}$ and
$\langle{O(\hat{q})}\rangle_{\beta,N}$ results from the difference between the
continuous loops $x(\tau)$ and $x_{N}(\tau)$. For the convenience of analysis,
we introduce the random variables

$\mathcal{A}=\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{missing}),~{}~{}~{}~{}\mathcal{B}=\frac{1}{\beta}\int_{0}^{\beta}O(x(\tau))\mathrm{d}\tau,$
(3.48)

and

$\mathcal{A}_{N}=\exp\bigg(-\int_{0}^{\beta}V^{a}(x_{N}(\tau))\mathrm{d}\tau\bigg{missing}),~{}~{}~{}~{}\mathcal{B}_{N}=\frac{1}{\beta}\int_{0}^{\beta}O(x_{N}(\tau))\mathrm{d}\tau,$
(3.49)

then $\langle{O(\hat{q})}\rangle_{\beta}$ and
$\langle{O(\hat{q})}\rangle_{\beta,N}$ can be expressed by

$\langle{O(\hat{q})}\rangle_{\beta}=\frac{\mathbb{E}_{\nu}[\mathcal{A}\mathcal{B}]}{\mathbb{E}_{\nu}[\mathcal{A}]},~{}~{}~{}~{}\langle{O(\hat{q})}\rangle_{\beta,N}=\frac{\mathbb{E}_{\nu}[\mathcal{A}_{N}\mathcal{B}_{N}]}{\mathbb{E}_{\nu}[\mathcal{A}_{N}]}.$
(3.50)

To estimate
$\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N}\big{|}$,
we only need to calculate $\mathbb{E}_{\nu}|\mathcal{A}-\mathcal{A}_{N}|$ and
$\mathbb{E}_{\nu}|\mathcal{B}-\mathcal{B}_{N}|$.

######  Lemma 3.4.

Under Assumptions (i)(ii), the random variables $\mathcal{A}_{N}$ and
$\mathcal{B}_{N}$ satisfy

$\displaystyle\mathcal{A}\leqslant\exp(\beta
M_{1}),~{}~{}~{}~{}\mathcal{A}_{N}\leqslant\exp(\beta
M_{1}),~{}~{}~{}~{}|\mathcal{B}|\leqslant
M_{2},~{}~{}~{}~{}|\mathcal{B}_{N}|\leqslant M_{2},$ (3.51)
$\displaystyle\mathbb{E}_{\nu}[\mathcal{A}]\geqslant\exp\Big(-\frac{3}{2}\beta
M_{1}-C_{0}M_{1}\Big{missing}),~{}~{}~{}~{}\mathbb{E}_{\nu}[\mathcal{A}_{N}]\geqslant\exp\Big(-\frac{3}{2}\beta
M_{1}-C_{0}M_{1}\Big{missing}),$ (3.52)

and

$\mathbb{E}_{\nu}|\mathcal{A}-\mathcal{A}_{N}|\leqslant\frac{K_{1}}{\sqrt{N}},~{}~{}~{}~{}\mathbb{E}_{\nu}|\mathcal{B}-\mathcal{B}_{N}|\leqslant\frac{K_{2}}{\sqrt{N}},$
(3.53)

where $C_{0}=\frac{d\beta}{2a}\coth\frac{a\beta}{2}$, and the constants
$K_{1}$ and $K_{2}$ are given by

$K_{1}=\beta\exp(\beta
M_{1})M_{1}\sqrt{\frac{d(\beta+2C_{0})}{2}},~{}~{}~{}~{}K_{2}=\frac{M_{2}}{2}\sqrt{d\beta}.$
(3.54)

The proof is given in Appendix A. Employing Lemma 3.4, it is direct to derive
the estimate of
$\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N}\big{|}$
in terms of the number of normal modes $N$.

######  Theorem 3.3.

Under Assumptions (i)(ii), the difference between
$\langle{O(\hat{q})}\rangle_{\beta}$ and
$\langle{O(\hat{q})}\rangle_{\beta,N}$ is estimated as

$\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N}\big{|}\leqslant\frac{K}{\sqrt{N}},$
(3.55)

where $C_{0}=\frac{d\beta}{2a}\coth\frac{a\beta}{2}$, and the constant $K$ is
given by

$K=\exp(6\beta M_{1}+2C_{0}M_{1})M_{2}\sqrt{2d(2\beta+3C_{0})}.$ (3.56)

The proof is given in Appendix A.

######  Remark.

In the numerical experiments of [22], it can be observed that the convergence
rate of
$\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N}\big{|}$
is actually 1 rather $1/2$ given in Theorem 3.3. In other words, the
convergence rate given in Theorem 3.3 is not optimal.

###  3.5 Discretized truncated continuous loop path integral representation

Similar to the expression in (3.46), we can write the statistical average
$\langle{O(\hat{q})}\rangle_{\beta,N,D}$ as

$\langle{O(\hat{q})}\rangle_{\beta,N,D}=\frac{\displaystyle\mathbb{E}_{\nu}\bigg{[}\exp\bigg(-\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D}))\bigg{missing})\times\bigg{(}\frac{1}{D}\sum_{j=0}^{D-1}O(x_{N}(j\beta_{D}))\bigg{)}\bigg{]}}{\displaystyle\mathbb{E}_{\nu}\bigg{[}\exp\bigg(-\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D}))\bigg{missing})\bigg{]}}.$
(3.57)

The difference between (3.46) and (3.57) is that the numerical integration in
(3.46) is replaced by the Riemann summation. For this reason, we define the
random variables

$\mathcal{A}_{N,D}=\exp\bigg(-\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D}))\bigg{missing}),~{}~{}~{}~{}\mathcal{B}_{N,D}=\frac{1}{D}\sum_{j=0}^{D-1}O(x_{N}(j\beta_{D})).$
(3.58)

Note that $\mathcal{A}_{N,D}$ and $\mathcal{B}_{N,D}$ are accurate
approximations to $\mathcal{A}_{N}$ and $\mathcal{B}_{N}$ as the number of
grid points $D\rightarrow\infty$. In the following we establish the estimates
the random variables $\mathcal{A}_{N,D}$ and $\mathcal{B}_{N,D}$ similar to
Lemma 3.4.

######  Lemma 3.5.

Under Assumptions (i)(ii), the random variables $\mathcal{A}_{N,D}$ and
$\mathcal{B}_{N,D}$ satisfy

$\displaystyle\mathcal{A}_{N,D}\leqslant\exp(\beta
M_{1}),~{}~{}~{}~{}|\mathcal{B}_{N,D}|\leqslant M_{2},$ (3.59)
$\displaystyle\mathbb{E}_{\nu}[\mathcal{A}_{N,D}]\geqslant\exp\Big(-\frac{3}{2}\beta
M_{1}-C_{0}M_{1}\Big{missing}),$ (3.60)

and

$\mathbb{E}_{\nu}|\mathcal{A}_{N}-\mathcal{A}_{N,D}|\leqslant\frac{L_{1}}{\sqrt{D}},~{}~{}~{}~{}\mathbb{E}_{\nu}|\mathcal{B}_{N}-\mathcal{B}_{N,D}|\leqslant\frac{L_{2}}{\sqrt{D}},$
(3.61)

where $C_{0}=\frac{d\beta}{2a}\coth\frac{a\beta}{2}$, and the constants
$L_{1}$ and $L_{2}$ are given by

$L_{1}=\beta\exp(\beta
M_{1})M_{1}\sqrt{2d(\beta+2C_{0})(2\beta+1)},~{}~{}~{}~{}L_{2}=M_{2}\sqrt{d\beta(2\beta+1)}.$
(3.62)

The proof is given in Appendix A. Employing Lemma 3.5, it is direct to derive
the estimate of
$\big{|}\langle{O(\hat{q})}\rangle_{\beta,N}-\langle{O(\hat{q})}\rangle_{\beta,N,D}\big{|}$
in terms of $N$ and $D$.

######  Theorem 3.4.

Under Assumptions (i)(ii), the difference between
$\langle{O(\hat{q})}\rangle_{\beta,N}$ and
$\langle{O(\hat{q})}\rangle_{\beta,N,D}$ is estimated as

$\big{|}\langle{O(\hat{q})}\rangle_{\beta,N}-\langle{O(\hat{q})}\rangle_{\beta,N,D}\big{|}\leqslant\frac{L}{\sqrt{D}},$
(3.63)

where $C_{0}=\frac{d\beta}{2a}\coth\frac{a\beta}{2}$, and the constant $L$ is
given by

$L=2\exp(6\beta M_{1}+2C_{0}M_{1})M_{2}\sqrt{2d(2\beta+1)(2\beta+3C_{0})}.$
(3.64)

The proof is given in Appendix A.

Combining the results in Theorem 3.3 and Theorem 3.4, we finally obtain

######  Corollary 3.5.

Under Assumptions (i)(ii), the difference between
$\langle{O(\hat{q})}\rangle_{\beta}$ and
$\langle{O(\hat{q})}\rangle_{\beta,N,D}$ is estimated as

$\displaystyle~{}~{}~{}~{}\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N,D}\big{|}$
(3.65) $\displaystyle\leqslant 2\exp(6\beta
M_{1}+2C_{0}M_{1})M_{2}\sqrt{2d(2\beta+3C_{0})}\bigg{(}\frac{1}{\sqrt{N}}+\frac{2\sqrt{2\beta+1}}{\sqrt{D}}\bigg{)},$

where the constant $C_{0}=\frac{d\beta}{2a}\coth\frac{a\beta}{2}$.

The result above shows that $\langle{O(\hat{q})}\rangle_{N,D}$ is indeed an
accurate approximation to the quantum thermal average
$\langle{O(\hat{q})}\rangle_{\beta}$ as the number of normal modes $N$ and the
number of grid points $D$ tend to infinity.

##  4 Conclusion

In this paper we study two kinds of path integral representations (PIR), the
std-PIR and the CL-PIR. We prove the convergence of the std-PIR, and quantify
the convergence of the truncated CL-PIR and the discretized truncated CL-PIR.
The proof is based on the Trotter product formula in the trace form. The
future studies focus on the intrinsic connection between the PIR and the
stochastic partial differential equations as well as other probabilistic
approaches.

## Acknowledgement

The work of Z. Zhou is partially supported by the National Key R&D Program of
China (Project No. 2020YFA0712000, 2021YFA1001200), and the National Natural
Science Foundation of China (Grant No. 12031013, 12171013).

X. Ye has used ChatGPT to improve the language in the introduction part. The
authors would like to thank Haitao Wang (SJTU) and Weijun Xu (PKU) for the
helpful discussions.

## Appendix

##  Appendix A Additional proofs for Section 3

######  Lemma A.1.

Given $a>0$, consider the quantum harmonic oscillator

$\hat{H}^{a}=\frac{\hat{p}^{2}}{2}+\frac{a^{2}}{2}\hat{q}^{2}=-\frac{1}{2}\Delta_{d}+\frac{a^{2}}{2}\sum_{i=1}^{d}\hat{q}_{i}^{2},$
(A.1)

then for any $q,\tilde{q}\in\mathbb{R}^{d}$, the kernel function
$\matrixelement{q}{e^{-\beta\hat{H}^{a}}}{\tilde{q}}$ is explicitly given by

$\matrixelement{q}{e^{-\beta\hat{H}^{a}}}{\tilde{q}}=\bigg{(}\frac{a}{2\pi\sinh(a\beta)}\bigg{)}^{\frac{d}{2}}\exp\bigg(-\frac{a}{\sinh(a\beta)}\bigg{(}\cosh(a\beta)\frac{|q|^{2}+|\tilde{q}|^{2}}{2}-q\cdot\tilde{q}\bigg{)}\bigg{missing}).$
(A.2)

The expression (A.2) is known as the Mehler kernel, and the derivation of the
result can be found in Problem 3-8 of [6].

######  Remark.

As the parameter $a\rightarrow 0$, $\hat{H}^{a}$ becomes the free particle
Schrödinger operator $\hat{H}^{0}$, and the kernel function
$\matrixelement{q}{e^{-\beta\hat{H}^{0}}}{\tilde{q}}$ is exactly the heat
kernel

$\matrixelement{q}{e^{-\beta\hat{H}^{0}}}{\tilde{q}}=\frac{1}{(2\pi\beta)^{\frac{d}{2}}}\exp\bigg(-\frac{|q-\tilde{q}|^{2}}{2\beta}\bigg{missing}).$
(A.3)

######  Proof (of Lemma 3.1).

The proof consists of two parts.
1\. Limit of
$\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}$
as $D\rightarrow\infty$
From Theorem 3.1.1 of [23], we define the Wiener measure $W_{q,q}^{\beta}$ of
the continuous loop $x(\tau)\in\mathbb{H}$ with the following rule: for
$0<t_{1}<\cdots<t_{P-1}<\beta$, the measure of the set

$\\{x(\tau)\in\mathbb{H}:x(0)=x(\beta)=q,~{}x(t_{j})\in
I_{j},~{}j=1,\cdots,P-1\\}$ (A.4)

in $\mathbb{R}^{d(P-1)}$ is defined by

$\int_{I_{1}}\mathrm{d}x_{1}\int_{I_{2}}\mathrm{d}x_{2}\cdots\int_{I_{P-1}}\mathrm{d}x_{P-1}\prod_{j=0}^{P-1}\matrixelement{x_{j}}{e^{-(\tau_{j+1}-\tau_{j})\hat{H}^{0}}}{x_{j+1}},$
(A.5)

where $I_{1},I_{2},\cdots,I_{P-1}$ are closed cuboids in $\mathbb{R}^{d}$, and
we presume

$t_{0}=0~{}~{}~{}~{}t_{P}=\beta,~{}~{}~{}~{}x_{0}=x_{P}=q\in\mathbb{R}^{d}.$
(A.6)

In other words, the continuous loop $x(\tau)$ in the Wiener measure
$W_{q,q}^{\beta}$ is a Brownian bridge with the endpoints fixed at
$q\in\mathbb{R}^{d}$. Using (3.1.10) of [23], we express
$\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}$
in the Feynman–Kac formula:

$\displaystyle~{}~{}~{}~{}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathrm{d}x_{1}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{D-1}\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}}{x_{j+1}}$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathrm{d}x_{1}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{D-1}\exp\bigg(-\beta_{D}\sum_{j=0}^{D-1}V(x_{j+1})\bigg{missing})\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{0}}}{x_{j+1}}$
$\displaystyle=\int_{\mathbb{H}}\exp\bigg(-\beta_{D}\sum_{j=1}^{D}V(x(j\beta_{D}))\bigg{missing})\mathrm{d}W_{q,q}^{\beta}.$
(A.7)

Here, for the continuous loop $x(\tau)$ in the Wiener measure
$W_{q,q}^{\beta}$, its marginal measure in the $D$ grid points
$x_{j}=x(j\beta_{D})$ is exactly

$\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{0}}}{x_{j+1}}.$
(A.8)

As the number of grid points $D\rightarrow\infty$, the dominated convergence
theorem implies

$\lim_{D\rightarrow\infty}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}}{q}=\int_{\mathbb{H}}\exp\bigg(-\int_{0}^{\beta}V(x(\tau))\mathrm{d}\tau\bigg{missing})\mathrm{d}W_{q,q}^{\beta}.$
(A.9)

Using the Feynman–Kac formula in Theorem 3.2.3 of [23], we have

$\matrixelement{q}{e^{-\beta\hat{H}}}{q}=\int_{\mathbb{H}}\exp\bigg(-\int_{0}^{\beta}V(x(\tau))\mathrm{d}\tau\bigg{missing})\mathrm{d}W_{q,q}^{\beta}.$
(A.10)

Combining (A.9) and (A.10), we obtain the desired result.
2\. Uniform-in-$D$ bound of
$\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}$
The potential function $V(q)$ satisfies

$V(q)=\frac{a^{2}}{2}|q|^{2}+V^{a}(q)\geqslant\frac{a^{2}}{2}|q|^{2}-M_{1},$
(A.11)

and thus we can write (A.7) as

$\displaystyle~{}~{}~{}~{}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathrm{d}x_{1}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{D-1}\exp\bigg(-\beta_{D}\sum_{j=0}^{D-1}V(x_{j+1})\bigg{missing})\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{0}}}{x_{j+1}}$
$\displaystyle\leqslant\exp(\beta
M_{1})\int_{\mathbb{R}^{d}}\mathrm{d}x_{1}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{D-1}\exp\bigg(-\frac{a^{2}\beta_{D}}{2}\sum_{j=1}^{D}|x_{j}|^{2}\bigg{missing})\prod_{j=0}^{D-1}\frac{1}{(2\pi\beta_{D})^{\frac{d}{2}}}e^{-\frac{|x_{j}-x_{j+1}|^{2}}{2\beta_{D}}}$
$\displaystyle=\frac{\exp(\beta
M_{1})}{(2\pi\beta_{D})^{\frac{dD}{2}}}\int_{\mathbb{R}^{d}}\mathrm{d}x_{1}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{D-1}F(x_{1},\cdots,x_{D-1}),$
(A.12)

where the function $F(x_{1},\cdots,x_{D-1})$ is given by

$F=\exp\bigg(-\frac{1}{2\beta_{D}}\sum_{j=0}^{D-1}|x_{j}-x_{j+1}|^{2}-\frac{a^{2}\beta_{D}}{2}\sum_{j=0}^{D-1}|x_{j}|^{2}\bigg{missing}).$
(A.13)

Now we consider the quantum harmonic oscillator

$\hat{H}^{a}=\frac{\hat{p}^{2}}{2}+\frac{a^{2}}{2}\hat{q}^{2}=-\frac{1}{2}\Delta_{d}+\frac{a^{2}}{2}\sum_{i=1}^{d}\hat{q}_{i}^{2},$
(A.14)

then using the Mehler kernel in Lemma A.1, we have

$\displaystyle~{}~{}~{}~{}\matrixelement{q}{e^{-\beta\hat{H}^{a}}}{q}=\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}\big{)}^{D}}{q}$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathrm{d}x_{1}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{D-1}\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{a}}}{x_{j+1}}$
$\displaystyle=\bigg{(}\frac{a}{2\pi\sinh(a\beta_{D})}\bigg{)}^{\frac{dD}{2}}\int_{\mathbb{R}^{d}}\mathrm{d}x_{1}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{D-1}G(x_{1},\cdots,x_{D-1}),$

where the function $G(x_{1},\cdots,x_{D-1})$ is given by

$G=\exp\bigg(-\frac{a}{\sinh(a\beta_{D})}\sum_{j=0}^{D-1}\bigg{(}\cosh(a\beta_{D})|x_{j}|^{2}-\sum_{j=0}^{D-1}x_{j}\cdot
x_{j+1}\bigg{)}\bigg{missing}).$ (A.15)

Comparing the expressions of $F$ and $G$, we observe that the coefficients
satisfy

$\frac{a}{\sinh(a\beta_{D})}\leqslant\frac{1}{\beta_{D}},~{}~{}~{}~{}a\cdot\frac{\cosh(a\beta_{D})-1}{\sinh(a\beta_{D})}\leqslant\frac{a^{2}\beta_{D}}{2},$
(A.16)

hence there is always $F\leqslant G$. Therefore, the inequality

$\int_{\mathbb{R}^{d}}\mathrm{d}x_{1}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{D-1}F(x_{1},\cdots,x_{D-1})\leqslant\int_{\mathbb{R}^{d}}\mathrm{d}x_{1}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{D-1}G(x_{1},\cdots,x_{D-1}).$
(A.17)

implies for any integer $D\in\mathbb{D}$,

$\frac{\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}}{\matrixelement{q}{e^{-\beta\hat{H}^{a}}}{q}}\leqslant\exp(\beta
M_{1})\bigg{(}\frac{\sinh(a\beta_{D})}{a\beta_{D}}\bigg{)}^{\frac{dD}{2}}.$
(A.18)

Note that as the integer $D\rightarrow\infty$,

$\frac{dD}{2}\log\bigg(\frac{\sinh(a\beta_{D})}{a\beta_{D}}\bigg{missing})\sim\frac{dD}{2}\bigg{(}\frac{\sinh(a\beta_{D})}{a\beta_{D}}-1\bigg{)}\sim\frac{dD}{2}\cdot\frac{(a\beta_{D})^{2}}{6}\sim
0,$ (A.19)

and thus there exists a constant $A$ such that

$A=\sup_{D\in\mathbb{N}}\bigg{(}\frac{\sinh(a\beta_{D})}{a\beta_{D}}\bigg{)}^{\frac{dD}{2}}<+\infty.$
(A.20)

Then we conclude

$\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}\leqslant
A\exp(\beta
M_{1})\matrixelement{q}{e^{-\beta\hat{H}^{a}}}{q},~{}~{}~{}~{}\forall
q\in\mathbb{R}^{d}.$ (A.21)

Again using the Mehler kernel, there exists constants $A,\lambda>0$ such that

$\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}\leqslant
A\exp(-\lambda|q|^{2}),~{}~{}~{}~{}\forall q\in\mathbb{R}^{d}.$ (A.22)

$\square$

######  Proof (of Lemma 3.2).

It is easy to calculate

$\mathbb{E}_{\nu}\bigg{[}\int_{0}^{\beta}|x(\tau)|^{2}\mathrm{d}\tau\bigg{]}=\mathbb{E}_{\nu}\bigg{[}\sum_{k=0}^{\infty}|\xi_{k}|^{2}\bigg{]}=\sum_{k=0}^{\infty}\frac{d}{\omega_{k}^{2}+a^{2}}=C_{0}.$
(A.23)

For any $\tau_{1},\tau_{2}\in[0,\beta]$, we assume
$\delta=|\tau_{1}-\tau_{2}|\leqslant\frac{\beta}{2}$ because $[0,\beta]$ is a
torus. Then

$\mathbb{E}_{\nu}\big{|}x(\tau_{1})-x(\tau_{2})\big{|}^{2}=\mathbb{E}_{\nu}\bigg{|}\sum_{k=0}^{\infty}\xi_{k}\big{(}c_{k}(\tau_{1})-c_{k}(\tau_{2})\big{)}\bigg{|}^{2}=\sum_{k=0}^{\infty}\mathbb{E}_{\nu}|\xi_{k}|^{2}\big{(}c_{k}(\tau_{1})-c_{k}(\tau_{2})\big{)}^{2}.$
(A.24)

Here, we have used the fact that $\\{\xi_{k}\\}_{k=0}^{\infty}$ are
independent random variables in the distribution $\nu$. Using
$\mathbb{E}|\xi_{k}|^{2}=\frac{d}{\omega_{k}^{2}+a^{2}}$, then (A.24)
immediately implies

$\displaystyle\mathbb{E}_{\nu}\big{|}x(\tau_{1})-x(\tau_{2})\big{|}^{2}$
$\displaystyle=\sum_{k=0}^{\infty}\frac{d}{\omega_{k}^{2}+a^{2}}\big{(}c_{k}(\tau_{1})-c_{k}(\tau_{2})\big{)}^{2}$
$\displaystyle=\frac{2d}{\beta}\sum_{k=1}^{\infty}\frac{(\sin\frac{2k\pi\tau_{1}}{\beta}-\sin\frac{2k\pi\tau_{2}}{\beta})^{2}}{(\frac{2k\pi}{\beta})^{2}+a^{2}}+\frac{2d}{\beta}\sum_{k=1}^{\infty}\frac{(\sin\frac{2k\pi\tau_{1}}{\beta}-\sin\frac{2k\pi\tau_{2}}{\beta})^{2}}{(\frac{2k\pi}{\beta})^{2}+a^{2}}$
$\displaystyle=\frac{8d}{\beta}\sum_{k=1}^{\infty}\frac{\sin^{2}\frac{k\pi(\tau_{1}-\tau_{2})}{\beta}}{(\frac{2k\pi}{\beta})^{2}+a^{2}}=\frac{8d}{\beta}\sum_{k\leqslant\frac{\beta}{\delta}}\frac{\sin^{2}\frac{k\pi(\tau_{1}-\tau_{2})}{\beta}}{(\frac{2k\pi}{\beta})^{2}+a^{2}}+\frac{8d}{\beta}\sum_{k>\frac{\beta}{\delta}}\frac{\sin^{2}\frac{k\pi(\tau_{1}-\tau_{2})}{\beta}}{(\frac{2k\pi}{\beta})^{2}+a^{2}}$
$\displaystyle\leqslant\frac{8d}{\beta}\sum_{k\leqslant\frac{\beta}{\delta}}\frac{(\frac{k\pi\delta}{\beta})^{2}}{(\frac{2k\pi}{\beta})^{2}}+\frac{8d}{\beta}\sum_{k>\frac{\beta}{\delta}}\frac{1}{(\frac{2k\pi}{\beta})^{2}}\leqslant\frac{8d}{\beta}\cdot\frac{\beta}{\delta}\cdot\frac{\delta^{2}}{4}+\frac{2d\beta}{\pi^{2}}\sum_{k>\frac{\beta}{\delta}}\frac{1}{k^{2}}$
$\displaystyle\leqslant
2d\beta\delta+\frac{2d\beta}{\pi^{2}}\cdot\frac{1}{\frac{\beta}{\delta}-1}\leqslant
2d\beta\delta+\frac{2d\beta}{\pi^{2}}\cdot\frac{2\delta}{\beta}\leqslant
d(2\beta+1)\delta.$ (A.25)

Next we prove the Hölder continuity of the continuous loop $x(\tau)$. For any
integer $m\in\mathbb{N}$,

$\displaystyle\mathbb{E}_{\nu}\big{|}x(\tau_{1})-x(\tau_{2})\big{|}^{2m}$
$\displaystyle=\mathbb{E}_{\nu}\bigg{|}\sum_{k=0}^{\infty}\xi_{k}\big{(}c_{k}(\tau_{1})-c_{k}(\tau_{2})\big{)}\bigg{|}^{2m}$
$\displaystyle=\mathbb{E}_{\nu}\bigg{(}\sum_{k=0}^{\infty}|\xi_{k}|^{2}\big{(}c_{k}(\tau_{1})-c_{k}(\tau_{2})\big{)}^{2}\bigg{)}^{m}.$
(A.26)

Here, we use the fact that the odd power of $\xi_{k}$ does not contribute to
the expectation. Expanding the RHS for the $m$ indices $k_{1},\cdots,k_{m}$,
we obtain

$\displaystyle\mathbb{E}_{\nu}\big{|}x(\tau_{1})-x(\tau_{2})\big{|}^{2m}=\sum_{k_{1}=0}^{\infty}\cdots\sum_{k_{m}=0}^{\infty}$
$\displaystyle\mathbb{E}_{\nu}\Big{(}|\xi_{k_{1}}|^{2}\cdots|\xi_{k_{m}}|^{2}\Big{)}$
(A.27)
$\displaystyle\big{(}c_{k_{1}}(\tau_{1})-c_{k_{1}}(\tau_{2})\big{)}^{2}\cdots\big{(}c_{k_{m}}(\tau_{1})-c_{k_{m}}(\tau_{2})\big{)}^{2}.$

For any indices $k_{1},\cdots,k_{m}$, the random variables
$\xi_{k_{1}},\cdots,\xi_{k_{m}}$ (possibly contain duplicate ones) are in the
Gaussian distribution, and thus there exists a constant $B_{m}$ such that

$\mathbb{E}_{\nu}\Big{(}|\xi_{k_{1}}|^{2}\cdots|\xi_{k_{m}}|^{2}\Big{)}\leqslant
B_{m}\,\mathbb{E}_{\nu}|\xi_{k_{1}}|^{2}\cdots\mathbb{E}_{\nu}|\xi_{k_{m}}|^{2}.$
(A.28)

Therefore we obtain

$\displaystyle\mathbb{E}_{\nu}\big{|}x(\tau_{1})-x(\tau_{2})\big{|}^{2m}$
$\displaystyle\leqslant
B_{m}\sum_{k=0}^{\infty}\bigg{(}\mathbb{E}_{\nu}|\xi_{k}|^{2}\big{(}c_{k}(\tau_{1})-c_{k}(\tau_{2})\big{)}^{2}\bigg{)}^{m}$
$\displaystyle\leqslant
B_{m}\big{(}d(2\beta+1)\big{)}^{m}|\tau_{1}-\tau_{2}|^{m}.$ (A.29)

Using the Kolmogorov continuity theorem, the continuous loop $x(\tau)$ is
$\gamma$-Hölder continuous for any $\gamma\in(0,\frac{m-1}{2m})$. Since the
integer $m$ can be sufficiently large, the constant $\gamma$ can be
arbitrarily close to $\frac{1}{2}$. Therefore $x(\tau)$ is $\gamma$-Hölder
continuous for any $\gamma\in(0,\frac{1}{2})$. $\square$

######  Proof (of Lemma 3.3).

The proof consists of three parts.
1\. Uniform-in-$D$ bound of
$\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}}$
Consider the quantum harmonic oscillator

$\hat{H}^{a}=\frac{\hat{p}^{2}}{2}+\frac{a^{2}}{2}\hat{q}^{2}=-\frac{1}{2}\Delta_{d}+\frac{a^{2}}{2}\sum_{i=1}^{d}\hat{q}_{i}^{2},$
(A.30)

introduced in Lemma A.1. From Theorem X.29 of [24], we deduce that both
$\hat{H}^{a}$ and $\hat{H}$ are essentially self-adjoint operators in
$C_{0}^{\infty}(\mathbb{R}^{d})$, which comprises all smooth functions in
$\mathbb{R}^{d}$ with compact support.

Now we aim to prove the Trotter product formula

$\lim_{D\rightarrow\infty}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}{e^{-\beta_{D}\hat{V}^{a}}}\big{)}^{D}}{\tilde{q}}=\matrixelement{q}{e^{-\beta\hat{H}}}{\tilde{q}}$
(A.31)

for any spatial coordinates $q,\tilde{q}\in\mathbb{R}^{d}$. Observe that
$e^{-\beta_{D}\hat{H}^{a}}$ and $e^{-\beta_{D}\hat{V}^{a}}$ are both
positivity-preserving operators, and from $V^{a}(q)\geqslant-M_{1}$ we have

$0\leqslant\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}}\leqslant\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{\beta_{D}M_{1}}\big{)}^{D}}{\tilde{q}}=e^{\beta
M_{1}}\matrixelement{q}{e^{-\beta\hat{H}^{a}}}{\tilde{q}}.$ (A.32)

Using the Mehler kernel in (A.2), we obtain the uniform-in-$D$ bound

$\displaystyle 0$
$\displaystyle\leqslant\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}}$
(A.33) $\displaystyle\leqslant e^{\beta M_{1}}\bigg{(}\frac{a}{2\pi
a\sinh(a\beta)}\bigg{)}^{\frac{d}{2}}\exp\bigg(-\frac{a}{\sinh(a\beta)}\bigg{(}\cosh(a\beta)\frac{|q|^{2}+|\tilde{q}|^{2}}{2}-q\cdot\tilde{q}\bigg{)}\bigg{missing}).$

As a consequence, there exists constants $A,\lambda>0$ such that

$\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{q}\leqslant
A\exp(-\lambda|q|^{2}),~{}~{}~{}~{}\forall q\in\mathbb{R}^{d}.$ (A.34)

2\. Limit of
$\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}}$
as $D\rightarrow\infty$
We show that
$\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}}$
can be represented in the Feynman–Kac formula. Similar to the Wiener measure
$W_{q,q}^{\beta}$ defined in (A.5), we define another kind of the Wiener
measure $U_{q,\tilde{q}}^{\beta}$ based on the Mehler kernel as follows. For
given $q,\tilde{q}\in\mathbb{R}^{d}$, let $\tilde{U}_{q,\tilde{q}}^{\beta}$ be
the Wiener measure of the continuous loop $x(\tau)\in\mathbb{H}$ defined the
following rule: for given $0<t_{1}<t_{2}<\cdots<t_{P-1}<\beta$, the measure of
the set

$\\{x(\tau)\in\mathbb{H}:x(0)=q,~{}x(\beta)=\tilde{q},~{}x(t_{j})\in
I_{j},~{}j=1,2,\cdots,P-1\\}$ (A.35)

is given by

$\int_{I_{1}}\mathrm{d}x_{1}\int_{I_{2}}\mathrm{d}x_{2}\cdots\int_{I_{P-1}}\mathrm{d}x_{P-1}\prod_{j=0}^{P-1}\matrixelement{x_{j}}{e^{-(\tau_{j+1}-\tau_{j})\hat{H}^{a}}}{x_{j+1}},$
(A.36)

where $I_{1},I_{2},\cdots,I_{P-1}$ are closed cuboids in $\mathbb{R}^{d}$, and
we presume

$t_{0}=0,~{}~{}~{}~{}x_{0}=q,~{}~{}~{}~{}t_{P}=\beta,~{}~{}~{}~{}x_{P}=\tilde{q}.$
(A.37)

Similar to (A.7), we can write
$\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}{e^{-\beta_{D}\hat{V}^{a}}}\big{)}^{D}}{\tilde{q}}$
as

$\displaystyle~{}~{}~{}~{}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}}$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathrm{d}x_{1}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{P-1}\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}}{x_{j+1}}$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathrm{d}x_{1}\cdots\int_{\mathbb{R}^{d}}\mathrm{d}x_{D-1}\exp\bigg(-\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{j+1})\bigg{missing})\prod_{j=0}^{D-1}\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{a}}}{x_{j+1}}$
$\displaystyle=\int_{\mathbb{H}}\exp\bigg(-\beta_{D}\sum_{j=1}^{D}V^{a}(x(j\beta_{D}))\bigg{missing})\mathrm{d}U_{q,\tilde{q}}^{\beta}.$
(A.38)

Here, the integration is taken over the continuous loop $x(\tau)$ in the
Wiener measure $U_{q,\tilde{q}}^{\beta}$. Let the number of grid points $D$
tend to infinity, we can apply the dominated convergence theorem on (A.38) to
deduce

$\lim_{D\rightarrow\infty}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}}=\int_{\mathbb{H}}\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{missing})\mathrm{d}U_{q,\tilde{q}}^{\beta}.$
(A.39)

3\. Feynman–Kac formula
There is one more step to obtain (A.31) from the limit (A.39). Multiply (A.39)
by the test function $\psi(\tilde{q})=\bra{q}\ket{\psi}\in
L^{2}(\mathbb{R}^{d})$, and integrate the expression over the variable
$q\in\mathbb{R}^{d}$,

$\displaystyle\lim_{D\rightarrow\infty}$
$\displaystyle\int_{\mathbb{R}^{d}}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}}\bra{\tilde{q}}\ket{\psi}\mathrm{d}\tilde{q}$
(A.40)
$\displaystyle=\bigg{\langle}\int_{\mathbb{H}}\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{missing})\mathrm{d}U_{q,\tilde{q}}^{\beta},\psi(\tilde{q})\bigg{\rangle}_{L^{2}(\mathbb{R}^{d})}.$

Here, $\langle{\cdot,\cdot}\rangle_{L^{2}(\mathbb{R}^{d})}$ is the inner
product in $L^{2}(\mathbb{R}^{d})$. Using the equality

$I=\int_{\mathbb{R}^{d}}\ket{\tilde{q}}\bra{\tilde{q}}\mathrm{d}\tilde{q},$
(A.41)

the limit (A.40) can be simplified as

$\displaystyle\lim_{D\rightarrow\infty}$
$\displaystyle\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\psi}$
(A.42)
$\displaystyle=\bigg{\langle}\int_{\mathbb{H}}\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{missing})\mathrm{d}U_{q,\tilde{q}}^{\beta},\psi(\tilde{q})\bigg{\rangle}_{L^{2}(\mathbb{R}^{d})}.$

Since $\hat{H}^{a}$ and $\hat{H}$ are both essentially self-adjoint operators
in $L^{2}(\mathbb{R}^{d})$, we can apply the Trotter product formula (Theorem
VIII.31 of [24]) to derive the strong limit

$\lim_{D\rightarrow\infty}\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}=e^{-\beta\hat{H}},~{}~{}~{}~{}\mbox{in
the strong $L^{2}(\mathbb{R}^{d})$ sense}.$ (A.43)

In particular, for the test function $\psi\in L^{2}(\mathbb{R}^{d})$, we have

$\lim_{D\rightarrow\infty}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\psi}=\matrixelement{q}{e^{-\beta\hat{H}}}{\psi},~{}~{}~{}~{}\mbox{in
the $L^{2}(\mathbb{R}^{d})$ sense}.$ (A.44)

Combining the limits (A.42) and (A.44), we obtain

$\matrixelement{q}{e^{-\beta\hat{H}}}{\psi}=\bigg{\langle}\int_{\mathbb{H}}\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{missing})\mathrm{d}U_{q,\tilde{q}}^{\beta},\psi(\tilde{q})\bigg{\rangle}_{L^{2}(\mathbb{R}^{d})}.$
(A.45)

Since $\psi(\tilde{q})$ can be any test function in $L^{2}(\mathbb{R}^{d})$,
we obtain the Feynman–Kac formula

$\matrixelement{q}{e^{-\beta\hat{H}}}{\tilde{q}}=\int_{\mathbb{H}}\exp\bigg(-\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{missing})\mathrm{d}U_{q,\tilde{q}}^{\beta}.$
(A.46)

Comparing (A.39) and (A.46), we finally obtain the Trotter product formula
(A.31), i.e.,

$\lim_{D\rightarrow\infty}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}}=\matrixelement{q}{e^{-\beta\hat{H}}}{\tilde{q}}.$

$\square$

######  Lemma A.2.

For any $\\{\xi_{k}\\}_{k=0}^{\infty}$ and $\\{\eta_{k}\\}_{k=0}^{\infty}$ in
$\mathbb{R}^{d}$, under Assumption (i), we have

$\displaystyle\bigg{|}\exp\bigg(-\int_{0}^{\beta}&V^{a}\bigg{(}\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{)}\mathrm{d}\tau\bigg{missing})-\exp\bigg(-\int_{0}^{\beta}V^{a}\bigg{(}\sum_{k=0}^{\infty}\eta_{k}c_{k}(\tau)\bigg{)}\mathrm{d}\tau\bigg{missing})\bigg{|}$
(A.47) $\displaystyle\leqslant\sqrt{2}\exp(\beta
M_{1})M_{1}\sqrt{\sum_{k=0}^{\infty}|\xi_{k}-\eta_{k}|^{2}\bigg{(}\beta+\sum_{k=0}^{\infty}|\xi_{k}|^{2}+\sum_{k=0}^{\infty}|\eta_{k}|^{2}\bigg{)}}.$

and under Assumption (ii), we have

$\bigg{|}\frac{1}{\beta}\int_{0}^{\beta}O\bigg{(}\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{)}-\frac{1}{\beta}\int_{0}^{\beta}O\bigg{(}\sum_{k=0}^{\infty}\eta_{k}c_{k}(\tau)\bigg{)}\bigg{|}\leqslant\frac{1}{\sqrt{\beta}}M_{2}\sqrt{\sum_{k=0}^{\infty}|\xi_{k}-\eta_{k}|^{2}}.$
(A.48)

######  Proof.

For any $q_{0},q_{1}\in\mathbb{R}^{d}$, the fundamental theorem of calculus
implies

$V(q_{1})-V(q_{0})=(q_{1}-q_{0})\cdot\int_{0}^{1}\nabla
V^{a}(q_{0}+\theta(q_{1}-q_{0}))\mathrm{d}\theta.$ (A.49)

By Assumption (i), we have

$\big{|}\nabla V^{a}(q_{0}+\theta(q_{1}-q_{0}))\big{|}\leqslant
M_{1}+M_{1}|q_{0}+\theta(q_{1}-q_{0})|\leqslant
M_{1}+M_{1}\max\\{|q_{0}|,|q_{1}|\\}.$ (A.50)

Then (A.49) and (A.50) imply the following estimate of
$|V^{a}(q_{1})-V^{a}(q_{0})|$:

$\displaystyle\big{|}V^{a}(q_{1})-V^{a}(q_{0})\big{|}$
$\displaystyle\leqslant|q_{1}-q_{0}|\int_{0}^{1}\big{|}\nabla
V^{a}(q_{0}+\theta(q_{1}-q_{0}))\big{|}\mathrm{d}\theta$
$\displaystyle\leqslant
M_{1}|q_{1}-q_{0}|\big{(}1+\max\\{|q_{0}|,|q_{1}|\\}\big{)}$
$\displaystyle\leqslant\sqrt{2}M_{1}|q_{1}-q_{0}|\sqrt{1+\max\\{|q_{0}|^{2},|q_{1}|^{2}\\}}$
$\displaystyle\leqslant\sqrt{2}M_{1}|q_{1}-q_{0}|\sqrt{1+|q_{0}|^{2}+|q_{1}|^{2}}.$
(A.51)

For given $\tau\in[0,\beta]$, we choose

$q_{0}=\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau),~{}~{}~{}~{}q_{1}=\sum_{k=0}^{\infty}\eta_{k}c_{k}(\tau),$
(A.52)

then (A.51) yields the estimate

$\displaystyle\bigg{|}$ $\displaystyle
V^{a}\bigg{(}\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{)}-V^{a}\bigg{(}\sum_{k=0}^{\infty}\eta_{k}c_{k}(\tau)\bigg{)}\bigg{|}$
(A.53)
$\displaystyle\leqslant\sqrt{2}M_{1}\bigg{|}\sum_{k=0}^{\infty}\big{(}\xi_{k}-\eta_{k}\big{)}c_{k}(\tau)\bigg{|}\sqrt{1+\bigg{|}\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{|}^{2}+\bigg{|}\sum_{k=0}^{\infty}\eta_{k}c_{k}(\tau)\bigg{|}^{2}}.$

Finally, using the Cauchy’s inequality,

$\displaystyle~{}~{}~{}~{}\int_{0}^{\beta}\bigg{|}V^{a}\bigg{(}\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{)}-V^{a}\bigg{(}\sum_{k=0}^{\infty}\eta_{k}c_{k}(\tau)\bigg{)}\bigg{|}\mathrm{d}\tau$
$\displaystyle\leqslant\sqrt{2}M_{1}\sqrt{\int_{0}^{\beta}\bigg{|}\sum_{k=0}^{\infty}\big{(}\xi_{k}-\eta_{k}\big{)}c_{k}(\tau)\bigg{|}^{2}\mathrm{d}\tau\int_{0}^{\beta}\bigg{(}1+\bigg{|}\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{|}^{2}+\bigg{|}\sum_{k=0}^{\infty}\eta_{k}c_{k}(\tau)\bigg{|}^{2}\bigg{)}\mathrm{d}\tau}$
$\displaystyle=\sqrt{2}M_{1}\sqrt{\sum_{k=0}^{\infty}|\xi_{k}-\eta_{k}|^{2}\bigg{(}\beta+\sum_{k=0}^{\infty}|\xi_{k}|^{2}+\sum_{k=0}^{\infty}|\eta_{k}|^{2}\bigg{)}}.$

Applying the variable substitution

$x=\int_{0}^{\beta}V^{a}\bigg{(}\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{)},~{}~{}~{}~{}y=\int_{0}^{\beta}V^{a}\bigg{(}\sum_{k=0}^{\infty}\eta_{k}c_{k}(\tau)\bigg{)},$
(A.54)

we have $x,y\geqslant-\beta M_{1}$ and thus

$|e^{-x}-e^{-y}|\leqslant|x-y|\max_{z\in[x,y]}|e^{-z}|\leqslant\exp(\beta
M_{1})|x-y|,$ (A.55)

which produces the inequality (A.47). From $|O(q_{0})-O(q_{1})|\leqslant
M_{2}|q_{0}-q_{1}|$ we derive

$\bigg{|}O\bigg{(}\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{)}-O\bigg{(}\sum_{k=0}^{\infty}\eta_{k}c_{k}(\tau)\bigg{)}\bigg{|}\leqslant
M_{2}\bigg{|}\sum_{k=0}^{\infty}\big{(}\xi_{k}-\eta_{k}\big{)}c_{k}(\tau)\bigg{|}.$
(A.56)

Using the Cauchy’s inequality,

$\displaystyle~{}~{}~{}~{}\int_{0}^{\beta}\bigg{|}O\bigg{(}\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{)}-O\bigg{(}\sum_{k=0}^{\infty}\eta_{k}c_{k}(\tau)\bigg{)}\bigg{|}\mathrm{d}\tau$
$\displaystyle\leqslant
M_{2}\int_{0}^{\beta}\bigg{|}\sum_{k=0}^{\infty}\big{(}\xi_{k}-\eta_{k}\big{)}c_{k}(\tau)\bigg{|}\mathrm{d}\tau$
$\displaystyle\leqslant
M_{2}\sqrt{\int_{0}^{\beta}\mathrm{d}\tau\int_{0}^{\beta}\bigg{|}\sum_{k=0}^{\infty}\big{(}\xi_{k}-\eta_{k}\big{)}c_{k}(\tau)\bigg{|}^{2}\mathrm{d}\tau}\leqslant\sqrt{\beta}M_{2}\sqrt{\sum_{k=0}^{\infty}|\xi_{k}-\eta_{k}|^{2}}.$

And thus we obtain the inequality (A.48). $\square$

######  Proof (of Lemma 3.4).

It is easy to see Assumption (i) implies $\mathcal{A}\leqslant\exp(\beta
M_{1})$ and $\mathcal{A}_{N}\leqslant\exp(\beta M_{1})$, while Assumption (ii)
implies $|\mathcal{B}|\leqslant M_{2}$ and $|\mathcal{B}_{N}|\leqslant M_{2}$.
Using the upper bound of $V^{a}(q)$

$V^{a}(q)\leqslant\frac{3}{2}M_{1}+M_{1}|q|^{2}$

derived in (3.2), we have

$\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\leqslant\int_{0}^{\beta}\bigg{(}\frac{3}{2}M_{1}+M_{1}|x(\tau)|^{2}\bigg{)}\mathrm{d}\tau=\frac{3}{2}\beta
M_{1}+M_{1}\sum_{k=0}^{\infty}|\xi_{k}|^{2}.$ (A.57)

Taking the expectation in both sides,

$\mathbb{E}_{\nu}\bigg{[}\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{]}\leqslant\frac{3}{2}\beta
M_{1}+M_{1}\sum_{k=0}^{\infty}\frac{d}{\omega_{k}^{2}+a^{2}}=\frac{3}{2}\beta
M_{1}+C_{0}M_{1}.$ (A.58)

Using the Jensen’s inequality, we obtain

$\mathbb{E}_{\nu}[\mathcal{A}]\geqslant\exp\bigg(-\mathbb{E}_{\nu}\bigg{[}\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{]}\bigg{missing})\geqslant\exp\Big(-\frac{3}{2}\beta
M_{1}-C_{0}M_{1}\Big{missing}).$ (A.59)

Note that for the continuous loop $x_{N}(\tau)$, we have the similar
inequality

$\mathbb{E}_{\nu}\bigg{[}\int_{0}^{\beta}V^{a}(x(\tau))\mathrm{d}\tau\bigg{]}\leqslant\frac{3}{2}\beta
M_{1}+M_{1}\sum_{k=0}^{N-1}\frac{d}{\omega_{k}^{2}+a^{2}}\leqslant\frac{3}{2}\beta
M_{1}+C_{0}M_{1},$ (A.60)

and thus $\mathcal{A}_{N}$ also satisfies

$\mathbb{E}_{\nu}[\mathcal{A}_{N}]\geqslant\exp\bigg(-\mathbb{E}_{\nu}\bigg{[}\int_{0}^{\beta}V^{a}(x_{N}(\tau))\mathrm{d}\tau\bigg{]}\bigg{missing})\geqslant\exp\Big(-\frac{3}{2}\beta
M_{1}-C_{0}M_{1}\Big{missing}).$ (A.61)

Now we calculate the difference between the continuous loops $x(\tau)$ and
$x_{N}(\tau)$, so that we can estimate
$\mathbb{E}_{\nu}|\mathcal{A}-\mathcal{A}_{N}|$ and
$\mathbb{E}_{\nu}|\mathcal{B}-\mathcal{B}_{N}|$. It is easy to calculate

$\int_{0}^{\beta}|x(\tau)-x_{N}(\tau)|^{2}\mathrm{d}\tau=\int_{0}^{\beta}\bigg{|}\sum_{k=N}^{\infty}\xi_{k}c_{k}(\tau)\bigg{|}^{2}\mathrm{d}\tau=\sum_{k=N}^{\infty}|\xi_{k}|^{2}.$
(A.62)

Taking the expectation in both sides,

$\mathbb{E}_{\nu}\bigg{[}\int_{0}^{\beta}|x(\tau)-x_{N}(\tau)|^{2}\mathrm{d}\tau\bigg{]}=\mathbb{E}_{\nu}\bigg{[}\sum_{k=N}^{\infty}|\xi_{k}|^{2}\bigg{]}=\sum_{k=N}^{\infty}\frac{d}{\omega_{k}^{2}+a^{2}}.$
(A.63)

Note the the eigenvalues $\\{\omega_{k}\\}_{k=0}^{\infty}$ satisfy

$\omega_{k}\geqslant\frac{k\pi}{\beta},~{}~{}~{}~{}k=0,1,2,\cdots,$ (A.64)

we have

$\sum_{k=N}^{\infty}\frac{d}{\omega_{k}^{2}+a^{2}}\leqslant
d\sum_{k=N}^{\infty}\bigg{(}\frac{\beta}{k\pi}\bigg{)}^{2}=\frac{d\beta^{2}}{\pi^{2}}\frac{1}{N-1}\leqslant\frac{d\beta^{2}}{\pi^{2}}\frac{2}{N}\leqslant\frac{d\beta^{2}}{4N},$
(A.65)

which implies

$\mathbb{E}_{\nu}\bigg{[}\sum_{k=N}^{\infty}|\xi_{k}|^{2}\bigg{]}\leqslant\frac{d\beta^{2}}{4N}.$
(A.66)

Applying Lemma A.2 on the two continuous loops $x(\tau)$ and $x_{N}(\tau)$, it
is easy to deduce

$|\mathcal{A}-\mathcal{A}_{N}|\leqslant\sqrt{2}\exp(\beta
M_{1})M_{1}\sqrt{\sum_{k=N}^{\infty}|\xi_{k}|^{2}\bigg{(}\beta+2\sum_{k=0}^{\infty}|\xi_{k}|^{2}\bigg{)}}.$
(A.67)

Then using (A.66) and the Cauchy’s inequality,

$\displaystyle\mathbb{E}_{\nu}|\mathcal{A}-\mathcal{A}_{N}|$
$\displaystyle\leqslant\sqrt{2}\exp(\beta
M_{1})M_{1}\sqrt{\mathbb{E}_{\nu}\bigg{[}\sum_{k=N}^{\infty}|\xi_{k}|^{2}\bigg{]}\mathbb{E}_{\nu}\bigg{[}\beta+2\sum_{k=0}^{\infty}|\xi_{k}|^{2}\bigg{]}}$
$\displaystyle\leqslant\sqrt{2}\exp(\beta
M_{1})M_{1}\sqrt{\frac{d\beta^{2}}{4N}\cdot(\beta+2C_{0})}$
$\displaystyle=\beta\exp(\beta M_{1})M_{1}\sqrt{\frac{d(\beta+2C_{0})}{2N}}.$
(A.68)

Also, the inequality

$|\mathcal{B}-\mathcal{B}_{N}|\leqslant\frac{1}{\sqrt{\beta}}M_{2}\sqrt{\sum_{k=N}^{\infty}|\xi_{k}|^{2}},$
(A.69)

implies

$\mathbb{E}_{\nu}|\mathcal{B}-\mathcal{B}_{N}|\leqslant\frac{1}{\sqrt{\beta}}M_{2}\sqrt{\mathbb{E}_{\nu}\bigg{[}\sum_{k=N}^{\infty}|\xi_{k}|^{2}\bigg{]}}\leqslant\frac{1}{\sqrt{\beta}}M_{2}\sqrt{\frac{d\beta^{2}}{4N}}=\frac{M_{2}}{2}\sqrt{\frac{d\beta}{N}}.$
(A.70)

$\square$

######  Proof (of Theorem 3.3).

Using the expressions of the quantum thermal average and the statistical
average,

$\langle{O(\hat{q})}\rangle_{\beta}=\frac{\mathbb{E}_{\nu}[\mathcal{A}\mathcal{B}]}{\mathbb{E}_{\nu}[\mathcal{A}]},~{}~{}~{}~{}\langle{O(\hat{q})}\rangle_{\beta,N}=\frac{\mathbb{E}_{\nu}[\mathcal{A}_{N}\mathcal{B}_{N}]}{\mathbb{E}_{\nu}[\mathcal{A}_{N}]},$
(A.71)

we can calculate

$\displaystyle~{}~{}~{}~{}\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N}\big{|}=\bigg{|}\frac{\mathbb{E}_{\nu}[\mathcal{A}\mathcal{B}]}{\mathbb{E}_{\nu}[\mathcal{A}]}-\frac{\mathbb{E}_{\nu}[\mathcal{A}_{N}\mathcal{B}_{N}]}{\mathbb{E}_{\nu}[\mathcal{A}_{N}]}\bigg{|}$
$\displaystyle\leqslant\frac{1}{\mathbb{E}_{\nu}[\mathcal{A}]\mathbb{E}_{\nu}[\mathcal{A}_{N}]}\Big{|}\mathbb{E}_{\nu}[\mathcal{A}\mathcal{B}]\mathbb{E}_{\nu}[A_{N}]-\mathbb{E}_{\nu}[\mathcal{A}_{N}\mathcal{B}_{N}]\mathbb{E}_{\nu}[\mathcal{A}]\Big{|}$
$\displaystyle\leqslant\exp\big(3\beta
M_{1}+2C_{0}M_{1}\big{missing})\Big{(}\big{|}\mathbb{E}_{\nu}[\mathcal{A}_{N}]\big{|}\big{|}\mathbb{E}_{\nu}[\mathcal{A}\mathcal{B}]-\mathbb{E}_{\nu}[\mathcal{A}_{N}\mathcal{B}_{N}]\big{|}+\big{|}\mathbb{E}_{\nu}[\mathcal{A}\mathcal{B}]\big{|}\big{|}\mathbb{E}_{\nu}[\mathcal{A}]-\mathbb{E}_{\nu}[\mathcal{A}_{N}]\big{|}\Big{)}$
$\displaystyle\leqslant\exp\big(3\beta
M_{1}+2C_{0}M_{1}\big{missing})\Big{(}\exp(\beta
M_{1})\mathbb{E}_{\nu}\big{|}\mathcal{A}\mathcal{B}-\mathcal{A}_{N}\mathcal{B}_{N}\big{|}+\exp(\beta
M_{1})M_{2}\mathbb{E}_{\nu}\big{|}\mathcal{A}-\mathcal{A}_{N}\big{|}\Big{)}.$

Furthermore, by Lemma 3.4,
$\mathbb{E}_{\nu}\big{|}\mathcal{A}\mathcal{B}-\mathcal{A}_{N}\mathcal{B}_{N}\big{|}$
is estimated by

$\mathbb{E}_{\nu}\big{|}\mathcal{A}\mathcal{B}-\mathcal{A}_{N}\mathcal{B}_{N}\big{|}\leqslant\exp(\beta
M_{1})\mathbb{E}_{\nu}\big{|}\mathcal{B}-\mathcal{B}_{N}\big{|}+M_{2}\mathbb{E}_{\nu}|\mathcal{A}-\mathcal{A}_{N}|\leqslant\frac{\exp(\beta
M_{1})K_{2}+M_{2}K_{1}}{\sqrt{N}},$

and thus we have

$\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N}\big{|}\leqslant\exp(3\beta
M_{1}+2C_{0}M_{1})\frac{\exp(2\beta M_{1})K_{2}+2\exp(\beta
M_{1})M_{2}K_{1}}{\sqrt{N}}.$ (A.72)

Therefore, the constant $K$ is explicitly given by

$\displaystyle K$ $\displaystyle=\exp(3\beta
M_{1}+2C_{0}M_{1})\Big{(}\exp(2\beta M_{1})K_{2}+2\exp(\beta
M_{1})M_{2}K_{1}\Big{)}$ $\displaystyle=\exp(3\beta
M_{1}+2C_{0}M_{1})\bigg{(}\exp(2\beta
M_{1})\cdot\frac{M_{2}}{2}\sqrt{d\beta}\,+$ $\displaystyle\hskip
142.26378pt2\exp(\beta M_{1})M_{2}\cdot\beta\exp(\beta
M_{1})M_{1}\sqrt{\frac{d(\beta+2C_{0})}{2}}\bigg{)}$
$\displaystyle=\exp(5\beta
M_{1}+2C_{0}M_{1})M_{2}\bigg{(}\frac{1}{2}\sqrt{d\beta}+\beta
M_{1}\sqrt{2d(\beta+2C_{0})}\bigg{)}$ $\displaystyle\leqslant\exp(6\beta
M_{1}+2C_{0}M_{1})M_{2}\bigg{(}\frac{1}{2}\sqrt{d\beta}+\sqrt{2d(\beta+2C_{0})}\bigg{)}$
$\displaystyle\leqslant\exp(6\beta
M_{1}+2C_{0}M_{1})M_{2}\sqrt{2d(2\beta+3C_{0})}.$

In the last inequality, we have used the algebraic inequality

$\frac{1}{2}\sqrt{x}+\sqrt{2x+4y}\leqslant\sqrt{4x+6y},~{}~{}~{}~{}\forall
x,y\geqslant 0.$ (A.73)

$\square$

######  Proof (of Lemma 3.5).

It is easy to see $\mathcal{A}_{N,D}\leqslant\exp(\beta M_{1})$ and
$|\mathcal{B}_{N,D}|\leqslant M_{2}$. Using

$V^{a}(q)\leqslant\frac{3}{2}M_{1}+M_{1}|q|^{2}$

derived in (3.2), we have

$\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D}))\leqslant\frac{3}{2}\beta
M_{1}+M_{1}\beta_{D}\sum_{j=0}^{D-1}|x_{N}(j\beta_{D})|^{2}.$ (A.74)

Taking the expectation in both sides, we obtain

$\mathbb{E}_{\nu}\bigg{[}\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D}))\bigg{]}\leqslant\frac{3}{2}\beta
M_{1}+M_{1}\mathbb{E}_{\nu}\bigg{[}\beta_{D}\sum_{j=0}^{D-1}|x_{N}(j\beta_{D})|^{2}\bigg{]}.$
(A.75)

For any $\tau\in[0,\beta]$, the value of the continuous loop at $\tau$ is

$x_{N}(\tau)=\sum_{k=0}^{N-1}\xi_{k}c_{k}(\tau),$ (A.76)

then using the independence of the random variables
$\\{\xi_{k}\\}_{k=0}^{\infty}$, we obtain

$\mathbb{E}_{\nu}\big{[}|x_{N}(\tau)|^{2}\big{]}=\sum_{k=0}^{N-1}\frac{d}{\omega_{k}^{2}+a^{2}}|c_{k}(\tau)|^{2}\leqslant\sum_{k=0}^{\infty}\frac{d}{\omega_{k}^{2}+a^{2}}|c_{k}(\tau)|^{2}.$
(A.77)

Note that the eigenfunctions $\\{c_{k}(\tau)\\}_{k=0}^{\infty}$ satisfy

$c_{0}(\tau)=\sqrt{\frac{1}{\beta}},~{}~{}~{}~{}|c_{2k-1}(\tau)|^{2}+|c_{2k}(\tau)|^{2}=\frac{2}{\beta},$
(A.78)

hence from (A.77) we have

$\mathbb{E}_{\nu}\big{[}|x_{N}(\tau)|^{2}\big{]}\leqslant\frac{1}{\beta}\sum_{k=0}^{\infty}\frac{d}{\omega_{k}^{2}+a^{2}}=\frac{C_{0}}{\beta}.$
(A.79)

As a consequence,

$\mathbb{E}_{\nu}\bigg{[}\beta_{D}\sum_{j=0}^{D-1}|x_{N}(j\beta_{D})|^{2}\bigg{]}\leqslant
C_{0}.$ (A.80)

Now from (A.75) and (A.80) we derive

$\mathbb{E}_{\nu}[\mathcal{A}_{N,D}]\geqslant\exp\bigg(-\mathbb{E}_{\nu}\bigg{[}\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D}))\bigg{]}\bigg{missing})\geqslant\exp\Big(-\frac{3}{2}\beta
M_{1}-C_{0}M_{1}\Big{missing}).$ (A.81)

Next we estimate the difference between $\mathcal{A}_{N}$ and
$\mathcal{A}_{N,D}$. By choosing

$x=\int_{0}^{\beta}V^{a}(x_{N}(\tau))\mathrm{d}\tau,~{}~{}~{}~{}y=\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D})),$
(A.82)

in the inequality (A.55), we have

$\displaystyle|\mathcal{A}_{N}-\mathcal{A}_{N,D}|$
$\displaystyle=\bigg{|}\exp\bigg(-\int_{0}^{\beta}V^{a}(x_{N}(\tau))\mathrm{d}\tau\bigg{missing})-\exp\bigg(-\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D}))\bigg{missing})\bigg{|}$
$\displaystyle\leqslant\exp(\beta
M_{1})\bigg{|}\int_{0}^{\beta}V^{a}(x_{N}(\tau))\mathrm{d}\tau-\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D}))\bigg{|}$
$\displaystyle\leqslant\exp(\beta
M_{1})\sum_{j=0}^{D-1}\bigg{|}\int_{j\beta_{D}}^{(j+1)\beta_{D}}\big{(}V^{a}(x_{N}(\tau))-V^{a}(x_{N}(j\beta_{D}))\big{)}\mathrm{d}\tau\bigg{|}$
$\displaystyle\leqslant\exp(\beta
M_{1})\sum_{j=0}^{N-1}\int_{j\beta_{D}}^{(j+1)\beta_{D}}\big{|}V^{a}(x_{N}(\tau))-V^{a}(x_{N}(j\beta_{D}))\big{|}\mathrm{d}\tau.$
(A.83)

Similar to the proof of Lemma A.2, we have

$\big{|}V^{a}(x_{N}(\tau_{1}))-V^{a}(x_{N}(\tau_{2}))\big{|}\leqslant\sqrt{2}M_{1}\big{|}x_{N}(\tau_{1})-x_{N}(\tau_{2})\big{|}\sqrt{1+|x_{N}(\tau_{1})|^{2}+|x_{N}(\tau_{2})|^{2}}.$
(A.84)

On the one hand, by Lemma 3.2 we have

$\mathbb{E}_{\nu}|x_{N}(\tau_{1})-x_{N}(\tau_{2})|^{2}\leqslant
d(2\beta+1)|\tau_{1}-\tau_{2}|.$ (A.85)

On the other hand, from (A.79) we have

$\mathbb{E}_{\nu}\big{[}1+|x_{N}(\tau_{1})|^{2}+|x_{N}(\tau_{2})|^{2}\big{]}\leqslant
1+\frac{2C_{0}}{\beta}.$ (A.86)

Therefore, taking the expectation in (A.84), we obtain

$\displaystyle~{}~{}~{}~{}\mathbb{E}_{\nu}\big{|}V^{a}(x_{N}(\tau_{1}))-V^{a}(x_{N}(\tau_{2}))\big{|}$
$\displaystyle\leqslant\sqrt{2}M_{1}\sqrt{\mathbb{E}_{\nu}|x_{N}(\tau_{1})-x_{N}(\tau_{2})|^{2}\,\mathbb{E}_{\nu}\big{[}1+|x_{N}(\tau_{1})|^{2}+|x_{N}(\tau_{2})|^{2}\big{]}}$
$\displaystyle\leqslant\sqrt{2}M_{1}\sqrt{d(2\beta+1)|\tau_{1}-\tau_{2}|\cdot\bigg{(}1+\frac{2C_{0}}{\beta}\bigg{)}}$
$\displaystyle=M_{1}\sqrt{\frac{2d}{\beta}(\beta+2C_{0})(2\beta+1)|\tau_{1}-\tau_{2}|}.$
(A.87)

Therefore from (A.83) we obtain the estimate of
$\mathbb{E}_{\nu}|\mathcal{A}_{N}-\mathcal{A}_{N,D}|$,

$\displaystyle~{}~{}~{}~{}\mathbb{E}_{\nu}|\mathcal{A}_{N}-\mathcal{A}_{N,D}|$
$\displaystyle\leqslant\exp(\beta
M_{1})M_{1}\sqrt{\frac{2d}{\beta}(\beta+2C_{0})(2\beta+1)}\sum_{j=0}^{D-1}\int_{j\beta_{D}}^{(j+1)\beta_{D}}\sqrt{|\tau-j\beta_{D}|}\mathrm{d}\tau$
$\displaystyle\leqslant\exp(\beta
M_{1})M_{1}\sqrt{\frac{2d}{\beta}(\beta+2C_{0})(2\beta+1)}\cdot\beta\sqrt{\beta_{D}}$
$\displaystyle\leqslant\beta\exp(\beta
M_{1})M_{1}\sqrt{\frac{2d(\beta+2C_{0})(2\beta+1)}{D}}.$ (A.88)

Also, $|\mathcal{B}_{N}-\mathcal{B}_{N,D}|$ is estimated by

$\displaystyle|\mathcal{B}_{N}-\mathcal{B}_{N,D}|$
$\displaystyle=\bigg{|}\frac{1}{\beta}\int_{0}^{\beta}O(x_{N}(\tau))\mathrm{d}\tau-\frac{1}{D}\sum_{j=0}^{D-1}O(x_{N}(j\beta_{D}))\bigg{|}$
$\displaystyle\leqslant\frac{1}{\beta}\sum_{j=0}^{D-1}\bigg{|}\int_{j\beta_{D}}^{(j+1)\beta_{D}}O(x_{N}(\tau))-O(x_{N}(j\beta_{D}))\bigg{|}\mathrm{d}\tau$
$\displaystyle\leqslant\frac{M_{2}}{\beta}\sum_{j=0}^{D-1}\int_{j\beta_{D}}^{(j+1)\beta_{D}}\big{|}x_{N}(\tau)-x_{N}(j\beta_{D})\big{|}.$
(A.89)

And thus from Lemma 3.2 we obtain

$\displaystyle\mathbb{E}_{\nu}|\mathcal{B}_{N}-\mathcal{B}_{N,D}|$
$\displaystyle\leqslant\frac{M_{2}}{\beta}\sum_{j=0}^{D-1}\int_{j\beta_{D}}^{(j+1)\beta_{D}}\sqrt{\mathbb{E}_{\nu}|x_{N}(\tau)-x_{N}(j\beta_{D})|^{2}}\mathrm{d}\tau$
$\displaystyle\leqslant\frac{M_{2}}{\beta}\sqrt{d(2\beta+1)}\sum_{j=0}^{D-1}\int_{j\beta_{D}}^{(j+1)\beta_{D}}\sqrt{|\tau-j\beta_{D}|}\mathrm{d}\tau$
$\displaystyle\leqslant\frac{M_{2}}{\beta}\sqrt{d(2\beta+1)}\cdot\beta\sqrt{\beta_{D}}$
$\displaystyle\leqslant M_{2}\sqrt{\frac{d\beta(2\beta+1)}{D}}.$ (A.90)

$\square$

######  Proof (of Theorem 3.4).

Similar to the proof of Theorem 3.3, we have

$\big{|}\langle{O(\hat{q})}\rangle_{\beta,N}-\langle{O(\hat{q})}\rangle_{\beta,N,D}\big{|}\leqslant\exp(3\beta
M_{1}+2C_{0}M_{1})\frac{\exp(2\beta M_{1})L_{2}+2\exp(\beta
M_{1})M_{2}L_{1}}{\sqrt{N}}.$ (A.91)

Therefore, the constant $L$ is explicitly given by

$\displaystyle L$ $\displaystyle=\exp(3\beta
M_{1}+2C_{0}M_{1})\Big{(}\exp(2\beta M_{1})L_{2}+2\exp(\beta
M_{1})M_{2}L_{1}\Big{)}$ $\displaystyle=\exp(3\beta
M_{1}+2C_{0}M_{1})\bigg{(}\exp(2\beta M_{1})\cdot
M_{2}\sqrt{d\beta(2\beta+1)}\,+$ $\displaystyle\hskip 85.35826pt2\exp(\beta
M_{1})M_{2}\cdot\beta\exp(\beta M_{1})M_{1}\sqrt{2d(\beta+2C_{0})}\bigg{)}$
$\displaystyle=\exp(5\beta
M_{1}+2C_{0}M_{1})M_{2}\Big{(}\sqrt{d\beta(2\beta+1)}+2\beta
M_{1}\sqrt{\beta+2C_{0}}\Big{)}$ $\displaystyle\leqslant 2\exp(6\beta
M_{1}+2C_{0}M_{1})\sqrt{2\beta+1}\bigg{(}\frac{1}{2}\sqrt{d\beta}+\sqrt{2d(\beta+2C_{0})}\bigg{)}$
$\displaystyle\leqslant 2\exp(6\beta
M_{1}+2C_{0}M_{1})M_{2}\sqrt{2d(2\beta+1)(2\beta+3C_{0})}.$

In particular, we have $L=2\sqrt{2\beta+1}K$. $\square$

## References

  * [1] Kerson Huang.  Statistical mechanics.  John Wiley & Sons, 2008.
  * [2] Donald A McQuarrie.  Statistical mechanics.  Sterling Publishing Company, 2000.
  * [3] Peter Atkins, Peter William Atkins, and Julio de Paula.  Atkins’ physical chemistry.  Oxford university press, 2014.
  * [4] Neil W Ashcroft and N David Mermin.  Solid state physics.  Cengage Learning, 2022.
  * [5] Subir Sachdev.  Quantum phase transitions.  Physics world, 12(4):33, 1999.
  * [6] Richard P Feynman, Albert R Hibbs, and Daniel F Styer.  Quantum mechanics and path integrals.  Courier Corporation, 2010.
  * [7] Mark Kac, Kenneth Baclawski, and Monroe David Donsker.  Mark kac: probability, number theory, and statistical physics: selected papers.  (No Title), 1979.
  * [8] MF Herman, EJ Bruskin, and BJ Berne.  On path integral monte carlo simulations.  The Journal of Chemical Physics, 76(10):5150–5155, 1982.
  * [9] Edwin L Pollock and David M Ceperley.  Simulation of quantum many-body systems by path-integral methods.  Physical Review B, 30(5):2555, 1984.
  * [10] Kenneth S Schweizer, Richard M Stratt, David Chandler, and Peter G Wolynes.  Convenient and accurate discretized path integral methods for equilibrium quantum mechanical calculations.  The Journal of Chemical Physics, 75(3):1347–1364, 1981.
  * [11] Michiel Sprik, Roger W Impey, and Michael L Klein.  Study of electron solvation in liquid ammonia using quantum path integral monte carlo calculations.  The Journal of chemical physics, 83(11):5802–5809, 1985.
  * [12] Stephen D Bond, Brian B Laird, and Benedict J Leimkuhler.  On the approximation of feynman–kac path integrals.  Journal of Computational Physics, 185(2):472–483, 2003.
  * [13] Jianfeng Lu, Yulong Lu, and Zhennan Zhou.  Continuum limit and preconditioned langevin sampling of the path integral molecular dynamics.  Journal of Computational Physics, 423:109788, 2020.
  * [14] Nawaf Bou-Rabee and Andreas Eberle.  Two-scale coupling for preconditioned hamiltonian monte carlo in infinite dimensions.  Stochastics and Partial Differential Equations: Analysis and Computations, 9:207–242, 2021.
  * [15] MF Herman, EJ Bruskin, and BJ Berne.  On path integral monte carlo simulations.  The Journal of Chemical Physics, 76(10):5150–5155, 1982.
  * [16] Philippe Sindzingre, Michael L Klein, and David M Ceperley.  Path-integral monte carlo study of low-temperature he 4 clusters.  Physical review letters, 63(15):1601, 1989.
  * [17] C Chakravarty, MC Gordillo, and DM Ceperley.  A comparison of the efficiency of fourier-and discrete time-path integral monte carlo.  The Journal of chemical physics, 109(6):2123–2134, 1998.
  * [18] Dominik Marx and Michele Parrinello.  Ab initio path integral molecular dynamics: Basic ideas.  The Journal of chemical physics, 104(11):4077–4082, 1996.
  * [19] Mark E Tuckerman, Dominik Marx, Michael L Klein, and Michele Parrinello.  Efficient and general algorithms for path integral car–parrinello molecular dynamics.  The Journal of chemical physics, 104(14):5579–5588, 1996.
  * [20] Michele Ceriotti, Michele Parrinello, Thomas E Markland, and David E Manolopoulos.  Efficient stochastic thermostatting of path integral molecular dynamics.  The Journal of chemical physics, 133(12), 2010.
  * [21] Jian Liu, Dezhang Li, and Xinzijian Liu.  A simple and accurate algorithm for path integral molecular dynamics with the langevin thermostat.  The Journal of chemical physics, 145(2), 2016.
  * [22] Xuda Ye and Zhennan Zhou.  Exact calculation of quantum thermal average from continuous loop path integral molecular dynamics.  arXiv preprint arXiv:2307.06510, 2023.
  * [23] James Glimm and Arthur Jaffe.  Quantum physics: a functional integral point of view.  Springer Science & Business Media, 2012.
  * [24] Michael Reed and Barry Simon.  I: Functional analysis, volume 1.  Academic press, 1981.