This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

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: abneryepku@pku.edu.cn    Zhennan Zhou Beijing International Center for Mathematical Research, Peking University, Beijing, 100871, P. R. China.  Email: zhennan@bicmr.pku.edu.cn
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 ℝd\mathbb{R}^{d}:

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

Here, q^\hat{q} and p^\hat{p} denote position and momentum operators in ℝd\mathbb{R}^{d}, and V​(q)V(q) is a real-valued potential function. When the system exists at a constant temperature T=1/βT=1/\beta, its state is described by the canonical ensemble with the density operator e−β​H^e^{-\beta\hat{H}}, and thus the partition function is given by 𝒵=Tr​[e−β​H^]\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,β][0,\beta], denoted as x​(τ)∈C​([0,β];ℝd)x(\tau)\in C([0,\beta];\mathbb{R}^{d}). The central to this approach is the energy function

ℰ​(x)=∫0β[12​|x′​(τ)|2+V​(x​(τ))]​dτ,\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−β​H^e^{-\beta\hat{H}} is the formal Boltzmann distribution π​(x)∝exp⁡(−ℰ​(x))\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​(q^)O(\hat{q}) in the canonical ensemble e−β​H^e^{-\beta\hat{H}}, where O​(q)O(q) is real-valued function. The expression takes shape as follows:

⟨O​(q^)⟩β=Tr​[e−β​H^​O​(q^)]Tr​[e−β​H^].\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:

⟨O​(q^)⟩β=∫(1β​∫0βO​(x​(τ))​dτ)​π​(x)​𝒟​[x],\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 𝒟​[x]\mathcal{D}[x] embodies the formal Lebesgue measure in the function space C​([0,β];ℝd)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 π​(x)\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​(τ)x(\tau) undergoes the approximation by utilizing its grid values, denoted as xj=x​(j​βD)x_{j}=x(j\beta_{D}), where the index j=0,1,⋯,D−1j=0,1,\cdots,D-1. Here, D∈ℕD\in\mathbb{N} signifies the number of grid points, and βD=β/D\beta_{D}=\beta/D. Consequently, the energy function ℰ​(x)\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 ⟨O​(q^)⟩β,Dstd\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∈ℕN\in\mathbb{N}, we arrive at the truncated CL-PIR, and the resulting approximated quantum thermal average is denoted as ⟨O​(q^)⟩β,N\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 ⟨O​(q^)⟩β,N,D\langle{O(\hat{q})}\rangle_{\beta,N,D}, where D∈ℕD\in\mathbb{N} denotes the number of grid points.

Note that we have employed the symbol D∈ℕD\in\mathbb{N} to denote the number of grid points in the interval [0,β][0,\beta] for both the std-PIR and the CL-PIR. In the special case N=DN=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.

    limD→∞⟨O​(q^)⟩β,Dstd=⟨O​(q^)⟩β\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.

    |⟨O​(q^)⟩β−⟨O​(q^)⟩β,N|≲1N\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.

    |⟨O​(q^)⟩β−⟨O​(q^)⟩β,N,D|≲1N+1D\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—NN and DD—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​(τ)x(\tau) as a Gaussian stochastic process in the torus [0,β][0,\beta]. Notably, we establish a significant link between the quantum thermal average ⟨O​(q^)⟩β\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 ℰ​(x)\mathcal{E}(x) from (1.2) comes to the forefront. Let D∈ℕD\in\mathbb{N} be the number of grid points in [0,β][0,\beta], and we employ the grid values xj=x​(j​βD)x_{j}=x(j\beta_{D}) of the continuous loop to represent the energy ℰ​(x)\mathcal{E}(x), where we presume βD=β/D\beta_{D}=\beta/D. Utilizing the finite-difference approximation

x′​(j​βD)≈x​((j+1)​βD)−x​(j​βD)βD=xj+1−xjβD,j=0,1,⋯,D−1,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 ℰD​(x)\mathcal{E}_{D}(x) can be approximated as

ℰDstd​(x)=12​βD​∑j=0D−1|xj−xj+1|2+βD​∑j=0D−1V​(xj),x∈ℝd​D.\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, {xj}j=0D−1\{x_{j}\}_{j=0}^{D-1} represent the grid values expected to align with the continuous loop x​(τ)x(\tau), characterized by the approximation:

xj≈x​(j​βD),j=0,1,⋯,D−1.x_{j}\approx x(j\beta_{D}),~{}~{}~{}~{}j=0,1,\cdots,D-1. (2.3)

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

πDstd​(x)=1ZDstd​exp⁡(−ℰDstd​(x)​missing),ZDstd=∫ℝd​Dexp⁡(−ℰDstd​(x)​missing)​dx.\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 ⟨O​(q^)⟩β\langle{O(\hat{q})}\rangle_{\beta} finds approximation through the expression:

⟨O​(q^)⟩β≈⟨O​(q^)⟩β,Dstd=∫ℝd​D(1D​∑j=0D−1O​(x​(j​βD)))​πDstd​(x)​dx.\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→∞D\rightarrow\infty, the formal continuum limit of the std-PIR (2.5) is given by (1.4).

To compute the statistical average ⟨O​(q^)⟩β,Dstd\langle{O(\hat{q})}\rangle_{\beta,D}^{\mathrm{std}} numerically, one usually simulates an ergodic Langevin process whose invariant distribution is exactly πDstd​(x)\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​(τ)x(\tau) finds representation through normal mode coordinates. Beginning with the eigenvalue problem associated with the second-order differential operator,

−c¨k​(τ)=ωk2​ck​(τ),τ∈[0,β],k=0,1,⋯,-\ddot{c}_{k}(\tau)=\omega_{k}^{2}c_{k}(\tau),~{}~{}~{}~{}\tau\in[0,\beta],~{}~{}~{}~{}k=0,1,\cdots, (2.6)

the eigenvalues {ωk}k=0∞\{\omega_{k}\}_{k=0}^{\infty} and the eigenfunctions {ck​(τ)}k=0∞\{c_{k}(\tau)\}_{k=0}^{\infty} manifest as

ω0\displaystyle\omega_{0} =0,\displaystyle=0, c0​(τ)\displaystyle c_{0}(\tau) =1β;\displaystyle=\sqrt{\frac{1}{\beta}}; (2.7)
ω2​k−1\displaystyle\omega_{2k-1} =2​k​πβ,\displaystyle=\frac{2k\pi}{\beta}, c2​k−1​(τ)\displaystyle c_{2k-1}(\tau) =2β​sin⁡2​k​π​τβ,k=1,2,⋯;\displaystyle=\sqrt{\frac{2}{\beta}}\sin\frac{2k\pi\tau}{\beta},~{}~{}~{}~{}k=1,2,\cdots;
ω2​k\displaystyle\omega_{2k} =2​k​πβ,\displaystyle=\frac{2k\pi}{\beta}, c2​k​(τ)\displaystyle c_{2k}(\tau) =2β​cos⁡2​k​π​τβ,k=1,2,⋯.\displaystyle=\sqrt{\frac{2}{\beta}}\cos\frac{2k\pi\tau}{\beta},~{}~{}~{}~{}k=1,2,\cdots.

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

ℍ=L2​([0,β];ℝd).\mathbb{H}=L^{2}([0,\beta];\mathbb{R}^{d}). (2.8)

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

x​(τ)=∑k=0∞ξk​ck​(τ),τ∈[0,β],x(\tau)=\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau),~{}~{}~{}~{}\tau\in[0,\beta], (2.9)

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

ℰ​(ξ)\displaystyle\mathcal{E}(\xi) =−12​∫0βx​(τ)T​x¨​(τ)​dτ+∫0βV​(x​(τ))​dτ\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
=12​∫0β(∑k=0∞ξk​ck​(τ))T​(∑k=0∞ωk2​ξk​ck​(τ))​dτ+∫0βV​(x​(τ))​dτ\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
=12​∑k=0∞ωk2​|ξk|2+∫0βV​(x​(τ))​dτ.\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>0a>0 and rewrite the energy function ℰ​(ξ)\mathcal{E}(\xi) as

ℰ​(ξ)=12​∑k=0∞(ωk2+a2)​|ξk|2+∫0βVa​(x​(τ))​dτ,\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

Va​(q)=V​(q)−a22​|q|2.V^{a}(q)=V(q)-\frac{a^{2}}{2}|q|^{2}. (2.12)

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

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

xN​(τ)=∑k=0N−1ξk​ck​(τ).x_{N}(\tau)=\sum_{k=0}^{N-1}\xi_{k}c_{k}(\tau). (2.13)

This truncation leads to the corresponding energy function

ℰN​(ξ)=12​∑k=0N−1(ωk2+a2)​|ξk|2+∫0βVa​(xN​(τ))​dτ,ξ∈ℝd​N,\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 πN​(ξ)\pi_{N}(\xi) in ℝd​N\mathbb{R}^{dN} assumes the form

πN​(ξ)=1ZN​exp⁡(−ℰN​(ξ)),ZN=∫ℝd​Nexp⁡(−ℰN​(ξ))​dξ.\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 ⟨O​(q^)⟩β\langle{O(\hat{q})}\rangle_{\beta} finds its approximation through

⟨O​(q^)⟩β≈⟨O​(q^)⟩β,N=∫ℝd​N(1β​∫0βO​(xN​(τ))​dτ)​πN​(ξ)​dξ.\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 ⟨O​(q^)⟩β,N\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

∫0βVa​(xN​(τ))​dτ​ and ​∫0βO​(xN​(τ))​dτ\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∈ℕD\in\mathbb{N} be the number of grid points in the torus [0,β][0,\beta], and define βD=β/D\beta_{D}=\beta/D, then these integral terms are approximated by

βD​∑j=0D−1Va​(xN​(j​βD))​ and ​βD​∑j=0D−1O​(xN​(j​βD)).\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

ℰN,D​(ξ)=12​∑k=0N−1(ωk2+a2)​|ξk|2+βD​∑j=0D−1Va​(xN​(j​βD)),ξ∈ℝd​N.\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

πN,D​(ξ)=1ZN,D​exp⁡(−ℰN,D​(ξ)),ZN,D=∫ℝd​Dexp⁡(−ℰN,D​(ξ))​dξ,\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

⟨O​(q^)⟩β,N,D=∫ℝd​N(1D​∑j=0D−1O​(xN​(j​βD)))​πN,D​(ξ)​dξ.\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=DN=D, i.e., the number of normal modes NN is chosen to be the same as the number of grid points DD. In the discrerized truncated CL-PIR, the potential function ℰD,D​(ξ)\mathcal{E}_{D,D}(\xi) is given by

ℰD,D​(ξ)=12​∑k=0D−1(ωk2+a2)​|ξk|2+βD​∑j=0D−1Va​(xD​(j​βD)),\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 xD​(τ)x_{D}(\tau) is the continuous loop

xD​(τ)=∑k=0D−1ξk​ck​(τ),τ∈[0,β].x_{D}(\tau)=\sum_{k=0}^{D-1}\xi_{k}c_{k}(\tau),~{}~{}~{}~{}\tau\in[0,\beta]. (2.22)

The DD grid values of the continuous loop xD​(τ)x_{D}(\tau) are given by

xj=xD​(j​βD)=∑k=0D−1ξk​ck​(j​βD),j=0,1,⋯,D−1.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∈{0,1,⋯,N−1}k\in\{0,1,\cdots,N-1\}, we have

βD​∑j=0D−1xj​ck​(j​βD)=∑l=0D−1ξl​(βD​∑j=0D−1cl​(j​βD)​ck​(j​βD))=∑l=0D−1ξl​δl​k=ξk,\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 {ξk}k=0D−1\{\xi_{k}\}_{k=0}^{D-1} can be represented by the grid values {xj}j=0D−1\{x_{j}\}_{j=0}^{D-1} via

ξk=βD​∑j=0D−1xj​ck​(j​βD),k=0,1,⋯,D−1.\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 ℰDstd​(x)\mathcal{E}_{D}^{\mathrm{std}}(x) in {ξk}k=0D−1\{\xi_{k}\}_{k=0}^{D-1}:

ℰDstd​(ξ)\displaystyle\mathcal{E}_{D}^{\mathrm{std}}(\xi) =12​βD​∑j=0D−1|xj−xj+1|2+βD​∑j=0D−1V​(xj)\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})
=12​βD​∑j=0D−1|xj−xj+1|2+a2​βD2​∑j=0D−1|xj|2+βD​∑j=0D−1Va​(xj)\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})
=12​∑k=0D−1(ωk,D2+a2)​|ξk|2+βD​∑j=0D−1Va​(xD​(j​βD)),\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 {ωk,D}k=0D−1\{\omega_{k,D}\}_{k=0}^{D-1} are given by

ω0,D=0,ω2​k−1,D=ω2​k,D=2​sin⁡k​πβ,k=1,2,⋯.\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 ℰDstd​(ξ)\mathcal{E}_{D}^{\mathrm{std}}(\xi) is almost equivalent to ℰD,D​(ξ)\mathcal{E}_{D,D}(\xi) given in (2.21), except for the difference on the coefficients ωk\omega_{k} and ωk,D\omega_{k,D}. Furthermore, there holds

limD→∞ωk,D=ωk,\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 DD 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)V(q) and the observable function O​(q)O(q) as follows.

Assumption.

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

  1. (i)

    |Va​(0)|⩽M1|V^{a}(0)|\leqslant M_{1}, Va​(q)⩾−M1V^{a}(q)\geqslant-M_{1}, |∇Va​(q)|⩽M1+M1​|q||\nabla V^{a}(q)|\leqslant M_{1}+M_{1}|q|,

and the observable function O​(q)O(q) satisfies

  1. (ii)

    max⁡{|O​(q)|,|∇O​(q)|}⩽M2\max\{|O(q)|,|\nabla O(q)|\}\leqslant M_{2}.

Using the fundamental theorem of calculus,

Va​(q)=Va​(0)+q⋅∫01∇Va​(θ​q)​dθ,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

Va​(q)\displaystyle V^{a}(q) ⩽|Va​(0)|+|q|​∫01|∇Va​(θ​q)|​dq\displaystyle\leqslant|V^{a}(0)|+|q|\int_{0}^{1}|\nabla V^{a}(\theta q)|\mathrm{d}q
⩽M1+|q|​∫01(M1+θ​M1​|q|)​dθ\displaystyle\leqslant M_{1}+|q|\int_{0}^{1}\big{(}M_{1}+\theta M_{1}|q|\big{)}\mathrm{d}\theta
=M1+M1​|q|+M12​|q|2⩽32​M1+M1​|q|2.\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 H^0\hat{H}^{0} and the potential operator V^\hat{V} by

H^0=p^22=−12​Δd,V^=V​(q^),\hat{H}^{0}=\frac{\hat{p}^{2}}{2}=-\frac{1}{2}\Delta_{d},~{}~{}~{}~{}\hat{V}=V(\hat{q}), (3.3)

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

H^a=p^22+a22​q^2=−12​Δd+a22​∑i=1dq^i2.\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 H^\hat{H} can be written as

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

The Trotter product formula is stated as

Tr​[e−β​H^]=limD→∞Tr​[(e−βD​H^0​e−βD​V^)D]=limD→∞Tr​[(e−βD​H^a​e−βD​V^a)D],\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−βD​H^0​e−βD​V^)D(e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}})^{D} and (e−βD​H^a​e−βD​V^a)D(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

⟨q|e−β​H^|q⟩=limD→∞⟨q|(e−βD​H^0​e−βD​V^)D|q⟩=limD→∞⟨q|(e−βD​H^a​e−βD​V^a)D|q⟩.\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.,

limD→∞(e−βD​H^0​e−βD​V^)D=e−β​H^,limD→∞(e−βD​H^a​e−βD​V^a)D=e−β​H^\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∈ℝdq\in\mathbb{R}^{d}, we have

limD→∞⟨q|(e−βD​H^0​e−βD​V^)D|q⟩=⟨q|e−β​H^|q⟩.\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,λ>0A,\lambda>0 independent of DD and qq such that

⟨q|(e−βD​H^0​e−βD​V^)D|q⟩⩽A​exp⁡(−λ​|q|2).\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

limD→∞⟨O​(q^)⟩β,Dstd=⟨O​(q^)⟩β.\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

⟨O​(q^)⟩β,Dstd=Tr​[(e−βD​H^0​e−βD​V^)D​O​(q^)]Tr​[(e−βD​H^0​e−βD​V^)D],\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∈ℕD\in\mathbb{N}.

Proof.

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

limD→∞∫ℝd⟨q|(e−βD​H^0​e−βD​V^)D|q⟩​O​(q)​dq=∫ℝd⟨q|e−β​H^|q⟩​O​(q)​dq,\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

limD→∞Tr​[(e−βD​H^0​e−βD​V^)D​O​(q^)]=Tr​[e−β​H^​O​(q^)].\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)≡1O(q)\equiv 1 in (3.14) we have

limD→∞Tr​[(e−βD​H^0​e−βD​V^)D]=Tr​[e−β​H^].\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 {xj}j=0D−1\{x_{j}\}_{j=0}^{D-1} in ℝd\mathbb{R}^{d}, we can write the LHS of (3.14) as

Tr​[(e−βD​H^0​e−βD​V^)D​O​(q^)]\displaystyle~{}~{}~{}~{}\mathrm{Tr}\big{[}\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}O(\hat{q})\big{]}
=∫ℝddx0​⋯​∫ℝddxD−1​∏j=0D−1⟨xj|e−βD​H^0​e−βD​V^|xj+1⟩​O​(x0)\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})
=1(2​π​βD)d​D2​∫ℝddx0​⋯​∫ℝddxD−1​exp⁡(−12​βD​∑j=0D−1|xj−xj+1|2−βD​∑j=0D−1V​(xj)​missing)​O​(x0)\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})
=1(2​π​βD)d​D2​∫ℝd​DO​(x0)​exp⁡(−ℰDstd​(x)​missing)​dx.\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

Tr​[(e−βD​H^0​e−βD​V^)D​O​(q^)]=1(2​π​βD)d​D2​∫ℝd​D(1D​∑j=0D−1O​(xj))​exp⁡(−ℰDstd​(x)​missing)​dx.\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,

Tr​[(e−βD​H^0​e−βD​V^)D]=1(2​π​βD)d​D2​∫ℝd​Dexp⁡(−ℰDstd​(x)​missing)​dx.\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

Tr​[(e−βD​H^0​e−βD​V^)D​O​(q^)]Tr​[(e−βD​H^0​e−βD​V^)D]=∫ℝd​D(1D​∑j=0D−1O​(xj))​πDstd​(x)​dx=⟨O​(q^)⟩β,Dstd.\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 DD tend to infinity, from (3.14) and (3.15) we obtain

limD→∞⟨O​(q^)⟩β,Dstd=limD→∞Tr​[(e−βD​H^0​e−βD​V^)D​O​(q^)]Tr​[(e−βD​H^0​e−βD​V^)D]=Tr​[e−β​H^​O​(q^)]Tr​[e−β​H^]=⟨O​(q^)⟩β.\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 |⟨O​(q^)⟩β−⟨O​(q^)⟩β,Dstd|\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 DD is still unknown.

3.3 Continuous loop path integral representation

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

x​(τ)=∑k=0∞ξk​ck​(τ).x(\tau)=\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau). (3.21)

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

ξk∼𝒩​(0,Idωk2+a2),k=0,1,⋯,\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 ν0\nu_{0} in the continuous loop mapping (3.21). It is clear that ν\nu is a Gaussian distribution in the Hilbert space ℍ=L2​([0,β];ℝd)\mathbb{H}=L^{2}([0,\beta];\mathbb{R}^{d}). Note that the first eigenvalue ω0=0\omega_{0}=0, and the introduction of the constant a>0a>0 ensures the well-posedness of Gaussian distribution ν\nu. Now we study the properties of the random continuous loop x​(τ)x(\tau) in the distribution ν\nu.

Lemma 3.2.

The random continuous loop x​(τ)x(\tau) with the distribution ν\nu satisfies

𝔼ν​[∫0β|x​(τ)|2​dτ]=C0,\mathbb{E}_{\nu}\bigg{[}\int_{0}^{\beta}|x(\tau)|^{2}\mathrm{d}\tau\bigg{]}=C_{0}, (3.23)

where the constant C0=d​β2​a​coth⁡a​β2C_{0}=\frac{d\beta}{2a}\coth\frac{a\beta}{2}. For any τ1,τ2∈[0,β]\tau_{1},\tau_{2}\in[0,\beta],

𝔼ν​|x​(τ1)−x​(τ2)|2⩽d​(2​β+1)​|τ1−τ2|.\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 γ∈(0,12)\gamma\in(0,\frac{1}{2}), x​(τ)x(\tau) is γ\gamma-Hölder continuous in the torus [0,β][0,\beta] almost surely.

The proof is given in Appendix A.

Remark.

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

Using the Gaussian distribution ν\nu, we can interpret the formal Boltzmann distribution π​(ξ)∝exp⁡(−ℰ​(ξ))\pi(\xi)\propto\exp(-\mathcal{E}(\xi)) as a probability distribution in the Hilbert space ℍ=L2​([0,β];ℝd)\mathbb{H}=L^{2}([0,\beta];\mathbb{R}^{d}) defined by the Randon–Nikodym derivative

d​πd​ν​(ξ)=1Z​exp⁡(−∫0βVa​(x​(τ))​dτ​missing),\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 ZZ is the normalization constant defined by

Z=𝔼ν​[exp⁡(−∫0βVa​(x​(τ))​missing)​d​τ].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

⟨O​(q^)⟩β=𝔼π​[1β​∫0βO​(x​(τ))​dτ].\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>0a>0, the distribution π\pi does not depend on the parameter aa. This is because formally π​(ξ)∝exp⁡(−ℰ​(ξ))\pi(\xi)\propto\exp(-\mathcal{E}(\xi)), where the energy function ℰ​(ξ)\mathcal{E}(\xi) does not depend on the parameter aa.

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

⟨O​(q^)⟩β=𝔼ν​[exp⁡(−∫0βVa​(x​(τ))​dτ​missing)×(1β​∫0βO​(x​(τ))​dτ)]𝔼ν​[exp⁡(−∫0βVa​(x​(τ))​dτ​missing)].\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,q~∈ℝdq,\tilde{q}\in\mathbb{R}^{d}, we have

limD→∞⟨q|(e−βD​H^a​e−βD​V^a)D|q~⟩=⟨q|e−β​H^|q~⟩.\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,λ>0A,\lambda>0 independent of DD and qq such that

⟨q|(e−βD​H^a​e−βD​V^a)D|q⟩⩽A​exp⁡(−λ​|q|2).\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

Tr​[e−β​H^​O​(q^)]Tr​[e−β​H^a]=𝔼ν​[exp⁡(−∫0βVa​(x​(τ))​dτ​missing)×(1β​∫0βO​(x​(τ))​dτ)].\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)≡1O(q)\equiv 1 in (3.31), we obtain

Tr​[e−β​H^]Tr​[e−β​H^a]=𝔼ν​[exp⁡(−∫0βVa​(x​(τ))​dτ​missing)]\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-DD bound in Lemma 3.3, we can apply the dominated convergence theorem to deduce

limD→∞∫ℝd⟨q|(e−βD​H^a​e−βD​V^a)D|q⟩​O​(q)​dq=∫ℝd⟨q|e−β​H^|q⟩​O​(q)​dq,\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

Tr​[e−β​H^​O​(q^)]=limD→∞Tr​[(e−βD​H^a​e−βD​V^a)D​O​(q^)].\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 Tr​[e−β​H^​O​(q^)]\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 Tr​[e−β​H^​O​(q^)]\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})] in the ring polymer distribution.

Tr​[e−β​H^​O​(q^)]\displaystyle\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})] =limD→∞Tr​[(e−βD​H^a​e−βD​V^a)D​O​(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{]}
=limD→∞∫ℝd​Ddx​exp⁡(−βD​∑j=0D−1V​(xj+1)​missing)​O​(x0)​∏j=0D−1⟨xj|e−βD​H^a|xj+1⟩\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​(x0)O(x_{0}), we obtain

Tr​[e−β​H^​O​(q^)]=limD→∞\displaystyle\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})]=\lim_{D\rightarrow\infty} ∫ℝd​Ddxexp⁡(−βD​∑j=0D−1Va​(xj)​missing)×\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)
(1D​∑j=0D−1O​(xj))×∏j=0D−1⟨xj|e−βD​H^a|xj+1⟩.\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

ΘD​(x)=1ZD​∏j=0D−1⟨xj|e−βD​H^a|xj+1⟩,x∈ℝd​D.\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 ⟨xj|e−βD​H^a|xj+1⟩\matrixelement{x_{j}}{e^{-\beta_{D}\hat{H}^{a}}}{x_{j+1}} is given by the Mehler kernel in Lemma A.1, and ZDZ_{D} is the normalization constant given by

ZD=∫ℝd​Ddx​∏j=0D−1⟨xj|e−βD​H^a|xj+1⟩.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 Va​(q)≡0V^{a}(q)\equiv 0 and O​(q)≡1O(q)\equiv 1 in (3.35), then

Tr​[e−β​H^a]=limD→∞∫ℝd​Ddx​∏j=0D−1⟨xj|e−βD​H^a|xj+1⟩.\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

Tr​[e−β​H^​O​(q^)]Tr​[e−β​H^a]=limD→∞∫ℝd​Dexp⁡(−βD​∑j=0D−1Va​(xj)​missing)​(1D​∑j=0D−1O​(xj))​Θ​(x)​dx.\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 Tr​[e−β​H^​O​(q^)]\mathrm{Tr}[e^{-\beta\hat{H}}O(\hat{q})] can be accurately approximated by the statistical average of the ring polymer distribution with DD grid points.
3. Equivalence between the Gaussian distributions
We define the Gaussian distribution U~β\tilde{U}^{\beta} in the Hilbert space ℍ=L2​([0,β];ℝd)\mathbb{H}=L^{2}([0,\beta];\mathbb{R}^{d}) by the following rule: for any constants 0⩽τ0<τ1<⋯<τP−1<β0\leqslant\tau_{0}<\tau_{1}<\cdots<\tau_{P-1}<\beta, the joint distribution of the random variables x0=x​(τ0),x1=x​(τ1),⋯,xP−1=x​(τP−1)x_{0}=x(\tau_{0}),x_{1}=x(\tau_{1}),\cdots,x_{P-1}=x(\tau_{P-1}) is proportional to

∏j=0P−2⟨xj|e−(τj+1−τj)​H^a|xj+1⟩×⟨xP−1|e−(τ0−τP−1+β)|x0⟩,\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 nn adjacent pairs in the position coordinates x​(τ0),x​(τ1),⋯,x​(τP−1)x(\tau_{0}),x(\tau_{1}),\cdots,x(\tau_{P-1}). From the Kolmogorov extension theorem, U~β\tilde{U}^{\beta} is indeed a well-defined Gaussian distribution in ℍ\mathbb{H}. The difference between the Gaussian distribution U~β\tilde{U}^{\beta} and the Wiener measure Uq,qβU_{q,q}^{\beta} defined in (A.36) is that the endpoints of the continuous loop in U~β\tilde{U}^{\beta} is flexible, while the endpoints in Uq,qβU_{q,q}^{\beta} are fixed at qq.

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

Tr​[e−β​H^​O​(q^)]Tr​[e−β​H^a]=limD→∞∫ℍexp⁡(−βD​∑j=0D−1Va​(x​(j​βD))​missing)​(1D​∑j=0D−1O​(x​(j​βD)))​dU~β.\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→∞D\rightarrow\infty, we can apply the dominated convergence theorem on (3.41) to deduce

Tr​[e−β​H^​O​(q^)]Tr​[e−β​H^a]=∫ℍexp⁡(−∫0βVa​(x​(τ))​dτ​missing)×(1β​∫0βO​(x​(τ))​dτ)×dU~β.\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 U~β\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​(τ)=∑k=0∞ξk​ck​(τ),ξk∼𝒩​(0,Idωk2+a2),k=0,1,2,⋯.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 U~β\tilde{U}^{\beta} and ν\nu are zero-mean Gaussian processes in [0,β][0,\beta], we only need to check their covariance functions are the same.

  1. 1.

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

    ⟨x|e−τ​H^a|y⟩​⟨y|e−(β−τ)​H^a|x⟩=(a2​π​sinh⁡(a​β))d​exp⁡(−A2​(|x2|+|y|2)+B​xT​y​missing),\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 AA and BB are given ny

    A=atanh⁡(a​τ)+atanh⁡(a​(β−τ)),B=asinh⁡(a​τ)+asinh⁡(a​(β−τ)).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)(x,y) in ℝd×ℝd\mathbb{R}^{d}\times\mathbb{R}^{d}, the covariance function is

    𝔼U~β​[x​(0)​x​(τ)]=B2A2−B2.\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

    𝔼ν​[x​(0)​x​(τ)]\displaystyle\mathbb{E}_{\nu}\big{[}x(0)x(\tau)\big{]} =𝔼ν​[∑k=0∞ξk​ck​(0)​∑k=0∞ξk​ck​(τ)]=∑k=0∞ck​(0)​ck​(τ)ωk2+a2\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}}
    =1β⋅1a2+2β​∑k=1∞cos⁡2​k​π​τβ(2​k​πβ)2+a2=1β​∑k∈ℤcos⁡2​k​π​τβωk2+a2.\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 𝔼U~β​[x​(0)​x​(τ)]\mathbb{E}_{\tilde{U}^{\beta}}\big{[}x(0)x(\tau)\big{]} and 𝔼ν​[x​(0)​x​(τ)]\mathbb{E}_{\nu}\big{[}x(0)x(\tau)\big{]} are exactly the same, which implies the Gaussian processes U~β\tilde{U}^{\beta} and ν\nu have the same covariance function. Therefore, U~β\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 ⟨O​(q^)⟩β,N\langle{O(\hat{q})}\rangle_{\beta,N} as

⟨O​(q^)⟩β,N=𝔼ν​[exp⁡(−∫0βVa​(xN​(τ))​dτ​missing)×(1β​∫0βO​(xN​(τ))​dτ)]𝔼ν​[exp⁡(−∫0βVa​(xN​(τ))​dτ​missing)],\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 xN​(τ)x_{N}(\tau) is the truncated continuous loop

xN​(τ)=∑k=0N−1ξk​ck​(τ),τ∈[0,β].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 ⟨O​(q^)⟩β\langle{O(\hat{q})}\rangle_{\beta} and ⟨O​(q^)⟩β,N\langle{O(\hat{q})}\rangle_{\beta,N} results from the difference between the continuous loops x​(τ)x(\tau) and xN​(τ)x_{N}(\tau). For the convenience of analysis, we introduce the random variables

𝒜=exp⁡(−∫0βVa​(x​(τ))​dτ​missing),ℬ=1β​∫0βO​(x​(τ))​dτ,\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

𝒜N=exp⁡(−∫0βVa​(xN​(τ))​dτ​missing),ℬN=1β​∫0βO​(xN​(τ))​dτ,\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 ⟨O​(q^)⟩β\langle{O(\hat{q})}\rangle_{\beta} and ⟨O​(q^)⟩β,N\langle{O(\hat{q})}\rangle_{\beta,N} can be expressed by

⟨O​(q^)⟩β=𝔼ν​[𝒜​ℬ]𝔼ν​[𝒜],⟨O​(q^)⟩β,N=𝔼ν​[𝒜N​ℬN]𝔼ν​[𝒜N].\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 |⟨O​(q^)⟩β−⟨O​(q^)⟩β,N|\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N}\big{|}, we only need to calculate 𝔼ν​|𝒜−𝒜N|\mathbb{E}_{\nu}|\mathcal{A}-\mathcal{A}_{N}| and 𝔼ν​|ℬ−ℬN|\mathbb{E}_{\nu}|\mathcal{B}-\mathcal{B}_{N}|.

Lemma 3.4.

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

𝒜⩽exp⁡(β​M1),𝒜N⩽exp⁡(β​M1),|ℬ|⩽M2,|ℬN|⩽M2,\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)
𝔼ν​[𝒜]⩾exp⁡(−32​β​M1−C0​M1​missing),𝔼ν​[𝒜N]⩾exp⁡(−32​β​M1−C0​M1​missing),\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

𝔼ν​|𝒜−𝒜N|⩽K1N,𝔼ν​|ℬ−ℬN|⩽K2N,\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 C0=d​β2​a​coth⁡a​β2C_{0}=\frac{d\beta}{2a}\coth\frac{a\beta}{2}, and the constants K1K_{1} and K2K_{2} are given by

K1=β​exp⁡(β​M1)​M1​d​(β+2​C0)2,K2=M22​d​β.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 |⟨O​(q^)⟩β−⟨O​(q^)⟩β,N|\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N}\big{|} in terms of the number of normal modes NN.

Theorem 3.3.

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

|⟨O​(q^)⟩β−⟨O​(q^)⟩β,N|⩽KN,\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N}\big{|}\leqslant\frac{K}{\sqrt{N}}, (3.55)

where C0=d​β2​a​coth⁡a​β2C_{0}=\frac{d\beta}{2a}\coth\frac{a\beta}{2}, and the constant KK is given by

K=exp⁡(6​β​M1+2​C0​M1)​M2​2​d​(2​β+3​C0).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 |⟨O​(q^)⟩β−⟨O​(q^)⟩β,N|\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N}\big{|} is actually 1 rather 1/21/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 ⟨O​(q^)⟩β,N,D\langle{O(\hat{q})}\rangle_{\beta,N,D} as

⟨O​(q^)⟩β,N,D=𝔼ν​[exp⁡(−βD​∑j=0D−1Va​(xN​(j​βD))​missing)×(1D​∑j=0D−1O​(xN​(j​βD)))]𝔼ν​[exp⁡(−βD​∑j=0D−1Va​(xN​(j​βD))​missing)].\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

𝒜N,D=exp⁡(−βD​∑j=0D−1Va​(xN​(j​βD))​missing),ℬN,D=1D​∑j=0D−1O​(xN​(j​βD)).\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 𝒜N,D\mathcal{A}_{N,D} and ℬN,D\mathcal{B}_{N,D} are accurate approximations to 𝒜N\mathcal{A}_{N} and ℬN\mathcal{B}_{N} as the number of grid points D→∞D\rightarrow\infty. In the following we establish the estimates the random variables 𝒜N,D\mathcal{A}_{N,D} and ℬN,D\mathcal{B}_{N,D} similar to Lemma 3.4.

Lemma 3.5.

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

𝒜N,D⩽exp⁡(β​M1),|ℬN,D|⩽M2,\displaystyle\mathcal{A}_{N,D}\leqslant\exp(\beta M_{1}),~{}~{}~{}~{}|\mathcal{B}_{N,D}|\leqslant M_{2}, (3.59)
𝔼ν​[𝒜N,D]⩾exp⁡(−32​β​M1−C0​M1​missing),\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

𝔼ν​|𝒜N−𝒜N,D|⩽L1D,𝔼ν​|ℬN−ℬN,D|⩽L2D,\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 C0=d​β2​a​coth⁡a​β2C_{0}=\frac{d\beta}{2a}\coth\frac{a\beta}{2}, and the constants L1L_{1} and L2L_{2} are given by

L1=β​exp⁡(β​M1)​M1​2​d​(β+2​C0)​(2​β+1),L2=M2​d​β​(2​β+1).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 |⟨O​(q^)⟩β,N−⟨O​(q^)⟩β,N,D|\big{|}\langle{O(\hat{q})}\rangle_{\beta,N}-\langle{O(\hat{q})}\rangle_{\beta,N,D}\big{|} in terms of NN and DD.

Theorem 3.4.

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

|⟨O​(q^)⟩β,N−⟨O​(q^)⟩β,N,D|⩽LD,\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 C0=d​β2​a​coth⁡a​β2C_{0}=\frac{d\beta}{2a}\coth\frac{a\beta}{2}, and the constant LL is given by

L=2​exp⁡(6​β​M1+2​C0​M1)​M2​2​d​(2​β+1)​(2​β+3​C0).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 ⟨O​(q^)⟩β\langle{O(\hat{q})}\rangle_{\beta} and ⟨O​(q^)⟩β,N,D\langle{O(\hat{q})}\rangle_{\beta,N,D} is estimated as

|⟨O​(q^)⟩β−⟨O​(q^)⟩β,N,D|\displaystyle~{}~{}~{}~{}\big{|}\langle{O(\hat{q})}\rangle_{\beta}-\langle{O(\hat{q})}\rangle_{\beta,N,D}\big{|} (3.65)
⩽2​exp⁡(6​β​M1+2​C0​M1)​M2​2​d​(2​β+3​C0)​(1N+2​2​β+1D),\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 C0=d​β2​a​coth⁡a​β2C_{0}=\frac{d\beta}{2a}\coth\frac{a\beta}{2}.

The result above shows that ⟨O​(q^)⟩N,D\langle{O(\hat{q})}\rangle_{N,D} is indeed an accurate approximation to the quantum thermal average ⟨O​(q^)⟩β\langle{O(\hat{q})}\rangle_{\beta} as the number of normal modes NN and the number of grid points DD 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>0a>0, consider the quantum harmonic oscillator

H^a=p^22+a22​q^2=−12​Δd+a22​∑i=1dq^i2,\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,q~∈ℝdq,\tilde{q}\in\mathbb{R}^{d}, the kernel function ⟨q|e−β​H^a|q~⟩\matrixelement{q}{e^{-\beta\hat{H}^{a}}}{\tilde{q}} is explicitly given by

⟨q|e−β​H^a|q~⟩=(a2​π​sinh⁡(a​β))d2​exp⁡(−asinh⁡(a​β)​(cosh⁡(a​β)​|q|2+|q~|22−q⋅q~)​missing).\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→0a\rightarrow 0, H^a\hat{H}^{a} becomes the free particle Schrödinger operator H^0\hat{H}^{0}, and the kernel function ⟨q|e−β​H^0|q~⟩\matrixelement{q}{e^{-\beta\hat{H}^{0}}}{\tilde{q}} is exactly the heat kernel

⟨q|e−β​H^0|q~⟩=1(2​π​β)d2​exp⁡(−|q−q~|22​β​missing).\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 ⟨q|(e−βD​H^0​e−βD​V^)D|q⟩\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q} as D→∞D\rightarrow\infty
From Theorem 3.1.1 of [23], we define the Wiener measure Wq,qβW_{q,q}^{\beta} of the continuous loop x​(τ)∈ℍx(\tau)\in\mathbb{H} with the following rule: for 0<t1<⋯<tP−1<β0<t_{1}<\cdots<t_{P-1}<\beta, the measure of the set

{x​(τ)∈ℍ:x​(0)=x​(β)=q,x​(tj)∈Ij,j=1,⋯,P−1}\{x(\tau)\in\mathbb{H}:x(0)=x(\beta)=q,~{}x(t_{j})\in I_{j},~{}j=1,\cdots,P-1\} (A.4)

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

∫I1dx1​∫I2dx2​⋯​∫IP−1dxP−1​∏j=0P−1⟨xj|e−(τj+1−τj)​H^0|xj+1⟩,\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 I1,I2,⋯,IP−1I_{1},I_{2},\cdots,I_{P-1} are closed cuboids in ℝd\mathbb{R}^{d}, and we presume

t0=0tP=β,x0=xP=q∈ℝd.t_{0}=0~{}~{}~{}~{}t_{P}=\beta,~{}~{}~{}~{}x_{0}=x_{P}=q\in\mathbb{R}^{d}. (A.6)

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

⟨q|(e−βD​H^0​e−βD​V^)D|q⟩\displaystyle~{}~{}~{}~{}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}
=∫ℝddx1​⋯​∫ℝddxD−1​∏j=0D−1⟨xj|e−βD​H^0​e−βD​V^|xj+1⟩\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}}
=∫ℝddx1​⋯​∫ℝddxD−1​exp⁡(−βD​∑j=0D−1V​(xj+1)​missing)​∏j=0D−1⟨xj|e−βD​H^0|xj+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}}
=∫ℍexp⁡(−βD​∑j=1DV​(x​(j​βD))​missing)​dWq,qβ.\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​(τ)x(\tau) in the Wiener measure Wq,qβW_{q,q}^{\beta}, its marginal measure in the DD grid points xj=x​(j​βD)x_{j}=x(j\beta_{D}) is exactly

∏j=0D−1⟨xj|e−βD​H^0|xj+1⟩.\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→∞D\rightarrow\infty, the dominated convergence theorem implies

limD→∞⟨q|(e−βD​H^0​e−βD​V^)|q⟩=∫ℍexp⁡(−∫0βV​(x​(τ))​dτ​missing)​dWq,qβ.\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

⟨q|e−β​H^|q⟩=∫ℍexp⁡(−∫0βV​(x​(τ))​dτ​missing)​dWq,qβ.\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-DD bound of ⟨q|(e−βD​H^0​e−βD​V^)D|q⟩\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}
The potential function V​(q)V(q) satisfies

V​(q)=a22​|q|2+Va​(q)⩾a22​|q|2−M1,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

⟨q|(e−βD​H^0​e−βD​V^)D|q⟩\displaystyle~{}~{}~{}~{}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{0}}e^{-\beta_{D}\hat{V}}\big{)}^{D}}{q}
=∫ℝddx1​⋯​∫ℝddxD−1​exp⁡(−βD​∑j=0D−1V​(xj+1)​missing)​∏j=0D−1⟨xj|e−βD​H^0|xj+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}}
⩽exp⁡(β​M1)​∫ℝddx1​⋯​∫ℝddxD−1​exp⁡(−a2​βD2​∑j=1D|xj|2​missing)​∏j=0D−11(2​π​βD)d2​e−|xj−xj+1|22​βD\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}}}
=exp⁡(β​M1)(2​π​βD)d​D2​∫ℝddx1​⋯​∫ℝddxD−1​F​(x1,⋯,xD−1),\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​(x1,⋯,xD−1)F(x_{1},\cdots,x_{D-1}) is given by

F=exp⁡(−12​βD​∑j=0D−1|xj−xj+1|2−a2​βD2​∑j=0D−1|xj|2​missing).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

H^a=p^22+a22​q^2=−12​Δd+a22​∑i=1dq^i2,\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

⟨q|e−β​H^a|q⟩=⟨q|(e−βD​H^a)D|q⟩\displaystyle~{}~{}~{}~{}\matrixelement{q}{e^{-\beta\hat{H}^{a}}}{q}=\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}\big{)}^{D}}{q}
=∫ℝddx1​⋯​∫ℝddxD−1​∏j=0D−1⟨xj|e−βD​H^a|xj+1⟩\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}}
=(a2​π​sinh⁡(a​βD))d​D2​∫ℝddx1​⋯​∫ℝddxD−1​G​(x1,⋯,xD−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​(x1,⋯,xD−1)G(x_{1},\cdots,x_{D-1}) is given by

G=exp⁡(−asinh⁡(a​βD)​∑j=0D−1(cosh⁡(a​βD)​|xj|2−∑j=0D−1xj⋅xj+1)​missing).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 FF and GG, we observe that the coefficients satisfy

asinh⁡(a​βD)⩽1βD,a⋅cosh⁡(a​βD)−1sinh⁡(a​βD)⩽a2​βD2,\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⩽GF\leqslant G. Therefore, the inequality

∫ℝddx1​⋯​∫ℝddxD−1​F​(x1,⋯,xD−1)⩽∫ℝddx1​⋯​∫ℝddxD−1​G​(x1,⋯,xD−1).\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∈𝔻D\in\mathbb{D},

⟨q|(e−βD​H^0​e−βD​V^)D|q⟩⟨q|e−β​H^a|q⟩⩽exp⁡(β​M1)​(sinh⁡(a​βD)a​βD)d​D2.\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→∞D\rightarrow\infty,

d​D2​log⁡(sinh⁡(a​βD)a​βD​missing)∼d​D2​(sinh⁡(a​βD)a​βD−1)∼d​D2⋅(a​βD)26∼0,\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 AA such that

A=supD∈ℕ(sinh⁡(a​βD)a​βD)d​D2<+∞.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

⟨q|(e−βD​H^0​e−βD​V^)D|q⟩⩽A​exp⁡(β​M1)​⟨q|e−β​H^a|q⟩,∀q∈ℝd.\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,λ>0A,\lambda>0 such that

⟨q|(e−βD​H^0​e−βD​V^)D|q⟩⩽A​exp⁡(−λ​|q|2),∀q∈ℝd.\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

𝔼ν​[∫0β|x​(τ)|2​dτ]=𝔼ν​[∑k=0∞|ξk|2]=∑k=0∞dωk2+a2=C0.\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 τ1,τ2∈[0,β]\tau_{1},\tau_{2}\in[0,\beta], we assume δ=|τ1−τ2|⩽β2\delta=|\tau_{1}-\tau_{2}|\leqslant\frac{\beta}{2} because [0,β][0,\beta] is a torus. Then

𝔼ν​|x​(τ1)−x​(τ2)|2=𝔼ν​|∑k=0∞ξk​(ck​(τ1)−ck​(τ2))|2=∑k=0∞𝔼ν​|ξk|2​(ck​(τ1)−ck​(τ2))2.\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 {ξk}k=0∞\{\xi_{k}\}_{k=0}^{\infty} are independent random variables in the distribution ν\nu. Using 𝔼​|ξk|2=dωk2+a2\mathbb{E}|\xi_{k}|^{2}=\frac{d}{\omega_{k}^{2}+a^{2}}, then (A.24) immediately implies

𝔼ν​|x​(τ1)−x​(τ2)|2\displaystyle\mathbb{E}_{\nu}\big{|}x(\tau_{1})-x(\tau_{2})\big{|}^{2} =∑k=0∞dωk2+a2​(ck​(τ1)−ck​(τ2))2\displaystyle=\sum_{k=0}^{\infty}\frac{d}{\omega_{k}^{2}+a^{2}}\big{(}c_{k}(\tau_{1})-c_{k}(\tau_{2})\big{)}^{2}
=2​dβ​∑k=1∞(sin⁡2​k​π​τ1β−sin⁡2​k​π​τ2β)2(2​k​πβ)2+a2+2​dβ​∑k=1∞(sin⁡2​k​π​τ1β−sin⁡2​k​π​τ2β)2(2​k​πβ)2+a2\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}}
=8​dβ​∑k=1∞sin2⁡k​π​(τ1−τ2)β(2​k​πβ)2+a2=8​dβ​∑k⩽βδsin2⁡k​π​(τ1−τ2)β(2​k​πβ)2+a2+8​dβ​∑k>βδsin2⁡k​π​(τ1−τ2)β(2​k​πβ)2+a2\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}}
⩽8​dβ​∑k⩽βδ(k​π​δβ)2(2​k​πβ)2+8​dβ​∑k>βδ1(2​k​πβ)2⩽8​dβ⋅βδ⋅δ24+2​d​βπ2​∑k>βδ1k2\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}}
⩽2​d​β​δ+2​d​βπ2⋅1βδ−1⩽2​d​β​δ+2​d​βπ2⋅2​δβ⩽d​(2​β+1)​δ.\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​(τ)x(\tau). For any integer m∈ℕm\in\mathbb{N},

𝔼ν​|x​(τ1)−x​(τ2)|2​m\displaystyle\mathbb{E}_{\nu}\big{|}x(\tau_{1})-x(\tau_{2})\big{|}^{2m} =𝔼ν​|∑k=0∞ξk​(ck​(τ1)−ck​(τ2))|2​m\displaystyle=\mathbb{E}_{\nu}\bigg{|}\sum_{k=0}^{\infty}\xi_{k}\big{(}c_{k}(\tau_{1})-c_{k}(\tau_{2})\big{)}\bigg{|}^{2m}
=𝔼ν​(∑k=0∞|ξk|2​(ck​(τ1)−ck​(τ2))2)m.\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 ξk\xi_{k} does not contribute to the expectation. Expanding the RHS for the mm indices k1,⋯,kmk_{1},\cdots,k_{m}, we obtain

𝔼ν​|x​(τ1)−x​(τ2)|2​m=∑k1=0∞⋯​∑km=0∞\displaystyle\mathbb{E}_{\nu}\big{|}x(\tau_{1})-x(\tau_{2})\big{|}^{2m}=\sum_{k_{1}=0}^{\infty}\cdots\sum_{k_{m}=0}^{\infty} 𝔼ν​(|ξk1|2​⋯​|ξkm|2)\displaystyle\mathbb{E}_{\nu}\Big{(}|\xi_{k_{1}}|^{2}\cdots|\xi_{k_{m}}|^{2}\Big{)} (A.27)
(ck1​(τ1)−ck1​(τ2))2​⋯​(ckm​(τ1)−ckm​(τ2))2.\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 k1,⋯,kmk_{1},\cdots,k_{m}, the random variables ξk1,⋯,ξkm\xi_{k_{1}},\cdots,\xi_{k_{m}} (possibly contain duplicate ones) are in the Gaussian distribution, and thus there exists a constant BmB_{m} such that

𝔼ν​(|ξk1|2​⋯​|ξkm|2)⩽Bm​𝔼ν​|ξk1|2​⋯​𝔼ν​|ξkm|2.\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

𝔼ν​|x​(τ1)−x​(τ2)|2​m\displaystyle\mathbb{E}_{\nu}\big{|}x(\tau_{1})-x(\tau_{2})\big{|}^{2m} ⩽Bm​∑k=0∞(𝔼ν​|ξk|2​(ck​(τ1)−ck​(τ2))2)m\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}
⩽Bm​(d​(2​β+1))m​|τ1−τ2|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​(τ)x(\tau) is γ\gamma-Hölder continuous for any γ∈(0,m−12​m)\gamma\in(0,\frac{m-1}{2m}). Since the integer mm can be sufficiently large, the constant γ\gamma can be arbitrarily close to 12\frac{1}{2}. Therefore x​(τ)x(\tau) is γ\gamma-Hölder continuous for any γ∈(0,12)\gamma\in(0,\frac{1}{2}). □\square

Proof (of Lemma 3.3).

The proof consists of three parts.
1. Uniform-in-DD bound of ⟨q|(e−βD​H^a​e−βD​V^a)D|q~⟩\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}}
Consider the quantum harmonic oscillator

H^a=p^22+a22​q^2=−12​Δd+a22​∑i=1dq^i2,\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 H^a\hat{H}^{a} and H^\hat{H} are essentially self-adjoint operators in C0∞​(ℝd)C_{0}^{\infty}(\mathbb{R}^{d}), which comprises all smooth functions in ℝd\mathbb{R}^{d} with compact support.

Now we aim to prove the Trotter product formula

limD→∞⟨q|(e−βD​H^a​e−βD​V^a)D|q~⟩=⟨q|e−β​H^|q~⟩\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,q~∈ℝdq,\tilde{q}\in\mathbb{R}^{d}. Observe that e−βD​H^ae^{-\beta_{D}\hat{H}^{a}} and e−βD​V^ae^{-\beta_{D}\hat{V}^{a}} are both positivity-preserving operators, and from Va​(q)⩾−M1V^{a}(q)\geqslant-M_{1} we have

0⩽⟨q|(e−βD​H^a​e−βD​V^a)D|q~⟩⩽⟨q|(e−βD​H^a​eβD​M1)D|q~⟩=eβ​M1​⟨q|e−β​H^a|q~⟩.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-DD bound

0\displaystyle 0 ⩽⟨q|(e−βD​H^a​e−βD​V^a)D|q~⟩\displaystyle\leqslant\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}} (A.33)
⩽eβ​M1​(a2​π​a​sinh⁡(a​β))d2​exp⁡(−asinh⁡(a​β)​(cosh⁡(a​β)​|q|2+|q~|22−q⋅q~)​missing).\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,λ>0A,\lambda>0 such that

⟨q|(e−βD​H^a​e−βD​V^a)D|q⟩⩽A​exp⁡(−λ​|q|2),∀q∈ℝd.\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 ⟨q|(e−βD​H^a​e−βD​V^a)D|q~⟩\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}} as D→∞D\rightarrow\infty
We show that ⟨q|(e−βD​H^a​e−βD​V^a)D|q~⟩\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 Wq,qβW_{q,q}^{\beta} defined in (A.5), we define another kind of the Wiener measure Uq,q~βU_{q,\tilde{q}}^{\beta} based on the Mehler kernel as follows. For given q,q~∈ℝdq,\tilde{q}\in\mathbb{R}^{d}, let U~q,q~β\tilde{U}_{q,\tilde{q}}^{\beta} be the Wiener measure of the continuous loop x​(τ)∈ℍx(\tau)\in\mathbb{H} defined the following rule: for given 0<t1<t2<⋯<tP−1<β0<t_{1}<t_{2}<\cdots<t_{P-1}<\beta, the measure of the set

{x​(τ)∈ℍ:x​(0)=q,x​(β)=q~,x​(tj)∈Ij,j=1,2,⋯,P−1}\{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

∫I1dx1​∫I2dx2​⋯​∫IP−1dxP−1​∏j=0P−1⟨xj|e−(τj+1−τj)​H^a|xj+1⟩,\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 I1,I2,⋯,IP−1I_{1},I_{2},\cdots,I_{P-1} are closed cuboids in ℝd\mathbb{R}^{d}, and we presume

t0=0,x0=q,tP=β,xP=q~.t_{0}=0,~{}~{}~{}~{}x_{0}=q,~{}~{}~{}~{}t_{P}=\beta,~{}~{}~{}~{}x_{P}=\tilde{q}. (A.37)

Similar to (A.7), we can write ⟨q|(e−βD​H^a​e−βD​V^a)D|q~⟩\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}{e^{-\beta_{D}\hat{V}^{a}}}\big{)}^{D}}{\tilde{q}} as

⟨q|(e−βD​H^a​e−βD​V^a)D|q~⟩\displaystyle~{}~{}~{}~{}\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\tilde{q}}
=∫ℝddx1​⋯​∫ℝddxP−1​∏j=0D−1⟨xj|e−βD​H^a​e−βD​V^a|xj+1⟩\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}}
=∫ℝddx1​⋯​∫ℝddxD−1​exp⁡(−βD​∑j=0D−1Va​(xj+1)​missing)​∏j=0D−1⟨xj|e−βD​H^a|xj+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}}
=∫ℍexp⁡(−βD​∑j=1DVa​(x​(j​βD))​missing)​dUq,q~β.\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​(τ)x(\tau) in the Wiener measure Uq,q~βU_{q,\tilde{q}}^{\beta}. Let the number of grid points DD tend to infinity, we can apply the dominated convergence theorem on (A.38) to deduce

limD→∞⟨q|(e−βD​H^a​e−βD​V^a)D|q~⟩=∫ℍexp⁡(−∫0βVa​(x​(τ))​dτ​missing)​dUq,q~β.\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 ψ​(q~)=⟨q|ψ⟩∈L2​(ℝd)\psi(\tilde{q})=\bra{q}\ket{\psi}\in L^{2}(\mathbb{R}^{d}), and integrate the expression over the variable q∈ℝdq\in\mathbb{R}^{d},

limD→∞\displaystyle\lim_{D\rightarrow\infty} ∫ℝd⟨q|(e−βD​H^a​e−βD​V^a)D|q~⟩​⟨q~|ψ⟩​dq~\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)
=⟨∫ℍexp⁡(−∫0βVa​(x​(τ))​dτ​missing)​dUq,q~β,ψ​(q~)⟩L2​(ℝd).\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, ⟨⋅,⋅⟩L2​(ℝd)\langle{\cdot,\cdot}\rangle_{L^{2}(\mathbb{R}^{d})} is the inner product in L2​(ℝd)L^{2}(\mathbb{R}^{d}). Using the equality

I=∫ℝd|q~⟩​⟨q~|​dq~,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

limD→∞\displaystyle\lim_{D\rightarrow\infty} ⟨q|(e−βD​H^a​e−βD​V^a)D|ψ⟩\displaystyle\matrixelement{q}{\big{(}e^{-\beta_{D}\hat{H}^{a}}e^{-\beta_{D}\hat{V}^{a}}\big{)}^{D}}{\psi} (A.42)
=⟨∫ℍexp⁡(−∫0βVa​(x​(τ))​dτ​missing)​dUq,q~β,ψ​(q~)⟩L2​(ℝd).\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 H^a\hat{H}^{a} and H^\hat{H} are both essentially self-adjoint operators in L2​(ℝd)L^{2}(\mathbb{R}^{d}), we can apply the Trotter product formula (Theorem VIII.31 of [24]) to derive the strong limit

limD→∞(e−βD​H^a​e−βD​V^a)D=e−β​H^,in the strong L2​(ℝd) sense.\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 ψ∈L2​(ℝd)\psi\in L^{2}(\mathbb{R}^{d}), we have

limD→∞⟨q|(e−βD​H^a​e−βD​V^a)D|ψ⟩=⟨q|e−β​H^|ψ⟩,in the L2​(ℝd) sense.\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

⟨q|e−β​H^|ψ⟩=⟨∫ℍexp⁡(−∫0βVa​(x​(τ))​dτ​missing)​dUq,q~β,ψ​(q~)⟩L2​(ℝd).\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 ψ​(q~)\psi(\tilde{q}) can be any test function in L2​(ℝd)L^{2}(\mathbb{R}^{d}), we obtain the Feynman–Kac formula

⟨q|e−β​H^|q~⟩=∫ℍexp⁡(−∫0βVa​(x​(τ))​dτ​missing)​dUq,q~β.\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.,

limD→∞⟨q|(e−βD​H^a​e−βD​V^a)D|q~⟩=⟨q|e−β​H^|q~⟩.\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 {ξk}k=0∞\{\xi_{k}\}_{k=0}^{\infty} and {ηk}k=0∞\{\eta_{k}\}_{k=0}^{\infty} in ℝd\mathbb{R}^{d}, under Assumption (i), we have

|exp⁡(−∫0β&​Va​(∑k=0∞ξk​ck​(τ))​dτ​missing)−exp⁡(−∫0βVa​(∑k=0∞ηk​ck​(τ))​dτ​missing)|\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)
⩽2​exp⁡(β​M1)​M1​∑k=0∞|ξk−ηk|2​(β+∑k=0∞|ξk|2+∑k=0∞|ηk|2).\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

|1β​∫0βO​(∑k=0∞ξk​ck​(τ))−1β​∫0βO​(∑k=0∞ηk​ck​(τ))|⩽1β​M2​∑k=0∞|ξk−ηk|2.\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 q0,q1∈ℝdq_{0},q_{1}\in\mathbb{R}^{d}, the fundamental theorem of calculus implies

V​(q1)−V​(q0)=(q1−q0)⋅∫01∇Va​(q0+θ​(q1−q0))​dθ.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

|∇Va​(q0+θ​(q1−q0))|⩽M1+M1​|q0+θ​(q1−q0)|⩽M1+M1​max⁡{|q0|,|q1|}.\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 |Va​(q1)−Va​(q0)||V^{a}(q_{1})-V^{a}(q_{0})|:

|Va​(q1)−Va​(q0)|\displaystyle\big{|}V^{a}(q_{1})-V^{a}(q_{0})\big{|} ⩽|q1−q0|​∫01|∇Va​(q0+θ​(q1−q0))|​dθ\displaystyle\leqslant|q_{1}-q_{0}|\int_{0}^{1}\big{|}\nabla V^{a}(q_{0}+\theta(q_{1}-q_{0}))\big{|}\mathrm{d}\theta
⩽M1​|q1−q0|​(1+max⁡{|q0|,|q1|})\displaystyle\leqslant M_{1}|q_{1}-q_{0}|\big{(}1+\max\{|q_{0}|,|q_{1}|\}\big{)}
⩽2​M1​|q1−q0|​1+max⁡{|q0|2,|q1|2}\displaystyle\leqslant\sqrt{2}M_{1}|q_{1}-q_{0}|\sqrt{1+\max\{|q_{0}|^{2},|q_{1}|^{2}\}}
⩽2​M1​|q1−q0|​1+|q0|2+|q1|2.\displaystyle\leqslant\sqrt{2}M_{1}|q_{1}-q_{0}|\sqrt{1+|q_{0}|^{2}+|q_{1}|^{2}}. (A.51)

For given τ∈[0,β]\tau\in[0,\beta], we choose

q0=∑k=0∞ξk​ck​(τ),q1=∑k=0∞ηk​ck​(τ),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{|} Va(∑k=0∞ξkck(τ))−Va(∑k=0∞ηkck(τ))|\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)
⩽2​M1​|∑k=0∞(ξk−ηk)​ck​(τ)|​1+|∑k=0∞ξk​ck​(τ)|2+|∑k=0∞ηk​ck​(τ)|2.\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,

∫0β|Va​(∑k=0∞ξk​ck​(τ))−Va​(∑k=0∞ηk​ck​(τ))|​dτ\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
⩽2​M1​∫0β|∑k=0∞(ξk−ηk)​ck​(τ)|2​dτ​∫0β(1+|∑k=0∞ξk​ck​(τ)|2+|∑k=0∞ηk​ck​(τ)|2)​dτ\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}
=2​M1​∑k=0∞|ξk−ηk|2​(β+∑k=0∞|ξk|2+∑k=0∞|ηk|2).\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=∫0βVa​(∑k=0∞ξk​ck​(τ)),y=∫0βVa​(∑k=0∞ηk​ck​(τ)),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⩾−β​M1x,y\geqslant-\beta M_{1} and thus

|e−x−e−y|⩽|x−y|​maxz∈[x,y]⁡|e−z|⩽exp⁡(β​M1)​|x−y|,|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​(q0)−O​(q1)|⩽M2​|q0−q1||O(q_{0})-O(q_{1})|\leqslant M_{2}|q_{0}-q_{1}| we derive

|O​(∑k=0∞ξk​ck​(τ))−O​(∑k=0∞ηk​ck​(τ))|⩽M2​|∑k=0∞(ξk−ηk)​ck​(τ)|.\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,

∫0β|O​(∑k=0∞ξk​ck​(τ))−O​(∑k=0∞ηk​ck​(τ))|​dτ\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
⩽M2​∫0β|∑k=0∞(ξk−ηk)​ck​(τ)|​dτ\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
⩽M2​∫0βdτ​∫0β|∑k=0∞(ξk−ηk)​ck​(τ)|2​dτ⩽β​M2​∑k=0∞|ξk−ηk|2.\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 𝒜⩽exp⁡(β​M1)\mathcal{A}\leqslant\exp(\beta M_{1}) and 𝒜N⩽exp⁡(β​M1)\mathcal{A}_{N}\leqslant\exp(\beta M_{1}), while Assumption (ii) implies |ℬ|⩽M2|\mathcal{B}|\leqslant M_{2} and |ℬN|⩽M2|\mathcal{B}_{N}|\leqslant M_{2}. Using the upper bound of Va​(q)V^{a}(q)

Va​(q)⩽32​M1+M1​|q|2V^{a}(q)\leqslant\frac{3}{2}M_{1}+M_{1}|q|^{2}

derived in (3.2), we have

∫0βVa​(x​(τ))​dτ⩽∫0β(32​M1+M1​|x​(τ)|2)​dτ=32​β​M1+M1​∑k=0∞|ξk|2.\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,

𝔼ν​[∫0βVa​(x​(τ))​dτ]⩽32​β​M1+M1​∑k=0∞dωk2+a2=32​β​M1+C0​M1.\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

𝔼ν​[𝒜]⩾exp⁡(−𝔼ν​[∫0βVa​(x​(τ))​dτ]​missing)⩾exp⁡(−32​β​M1−C0​M1​missing).\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 xN​(τ)x_{N}(\tau), we have the similar inequality

𝔼ν​[∫0βVa​(x​(τ))​dτ]⩽32​β​M1+M1​∑k=0N−1dωk2+a2⩽32​β​M1+C0​M1,\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 𝒜N\mathcal{A}_{N} also satisfies

𝔼ν​[𝒜N]⩾exp⁡(−𝔼ν​[∫0βVa​(xN​(τ))​dτ]​missing)⩾exp⁡(−32​β​M1−C0​M1​missing).\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​(τ)x(\tau) and xN​(τ)x_{N}(\tau), so that we can estimate 𝔼ν​|𝒜−𝒜N|\mathbb{E}_{\nu}|\mathcal{A}-\mathcal{A}_{N}| and 𝔼ν​|ℬ−ℬN|\mathbb{E}_{\nu}|\mathcal{B}-\mathcal{B}_{N}|. It is easy to calculate

∫0β|x​(τ)−xN​(τ)|2​dτ=∫0β|∑k=N∞ξk​ck​(τ)|2​dτ=∑k=N∞|ξk|2.\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,

𝔼ν​[∫0β|x​(τ)−xN​(τ)|2​dτ]=𝔼ν​[∑k=N∞|ξk|2]=∑k=N∞dωk2+a2.\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 {ωk}k=0∞\{\omega_{k}\}_{k=0}^{\infty} satisfy

ωk⩾k​πβ,k=0,1,2,⋯,\omega_{k}\geqslant\frac{k\pi}{\beta},~{}~{}~{}~{}k=0,1,2,\cdots, (A.64)

we have

∑k=N∞dωk2+a2⩽d​∑k=N∞(βk​π)2=d​β2π2​1N−1⩽d​β2π2​2N⩽d​β24​N,\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

𝔼ν​[∑k=N∞|ξk|2]⩽d​β24​N.\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​(τ)x(\tau) and xN​(τ)x_{N}(\tau), it is easy to deduce

|𝒜−𝒜N|⩽2​exp⁡(β​M1)​M1​∑k=N∞|ξk|2​(β+2​∑k=0∞|ξk|2).|\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,

𝔼ν​|𝒜−𝒜N|\displaystyle\mathbb{E}_{\nu}|\mathcal{A}-\mathcal{A}_{N}| ⩽2​exp⁡(β​M1)​M1​𝔼ν​[∑k=N∞|ξk|2]​𝔼ν​[β+2​∑k=0∞|ξk|2]\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{]}}
⩽2​exp⁡(β​M1)​M1​d​β24​N⋅(β+2​C0)\displaystyle\leqslant\sqrt{2}\exp(\beta M_{1})M_{1}\sqrt{\frac{d\beta^{2}}{4N}\cdot(\beta+2C_{0})}
=β​exp⁡(β​M1)​M1​d​(β+2​C0)2​N.\displaystyle=\beta\exp(\beta M_{1})M_{1}\sqrt{\frac{d(\beta+2C_{0})}{2N}}. (A.68)

Also, the inequality

|ℬ−ℬN|⩽1β​M2​∑k=N∞|ξk|2,|\mathcal{B}-\mathcal{B}_{N}|\leqslant\frac{1}{\sqrt{\beta}}M_{2}\sqrt{\sum_{k=N}^{\infty}|\xi_{k}|^{2}}, (A.69)

implies

𝔼ν​|ℬ−ℬN|⩽1β​M2​𝔼ν​[∑k=N∞|ξk|2]⩽1β​M2​d​β24​N=M22​d​βN.\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,

⟨O​(q^)⟩β=𝔼ν​[𝒜​ℬ]𝔼ν​[𝒜],⟨O​(q^)⟩β,N=𝔼ν​[𝒜N​ℬN]𝔼ν​[𝒜N],\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

|⟨O​(q^)⟩β−⟨O​(q^)⟩β,N|=|𝔼ν​[𝒜​ℬ]𝔼ν​[𝒜]−𝔼ν​[𝒜N​ℬN]𝔼ν​[𝒜N]|\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{|}
⩽1𝔼ν​[𝒜]​𝔼ν​[𝒜N]​|𝔼ν​[𝒜​ℬ]​𝔼ν​[AN]−𝔼ν​[𝒜N​ℬN]​𝔼ν​[𝒜]|\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{|}
⩽exp⁡(3​β​M1+2​C0​M1​missing)​(|𝔼ν​[𝒜N]|​|𝔼ν​[𝒜​ℬ]−𝔼ν​[𝒜N​ℬN]|+|𝔼ν​[𝒜​ℬ]|​|𝔼ν​[𝒜]−𝔼ν​[𝒜N]|)\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{)}
⩽exp⁡(3​β​M1+2​C0​M1​missing)​(exp⁡(β​M1)​𝔼ν​|𝒜​ℬ−𝒜N​ℬN|+exp⁡(β​M1)​M2​𝔼ν​|𝒜−𝒜N|).\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, 𝔼ν​|𝒜​ℬ−𝒜N​ℬN|\mathbb{E}_{\nu}\big{|}\mathcal{A}\mathcal{B}-\mathcal{A}_{N}\mathcal{B}_{N}\big{|} is estimated by

𝔼ν​|𝒜​ℬ−𝒜N​ℬN|⩽exp⁡(β​M1)​𝔼ν​|ℬ−ℬN|+M2​𝔼ν​|𝒜−𝒜N|⩽exp⁡(β​M1)​K2+M2​K1N,\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

|⟨O​(q^)⟩β−⟨O​(q^)⟩β,N|⩽exp⁡(3​β​M1+2​C0​M1)​exp⁡(2​β​M1)​K2+2​exp⁡(β​M1)​M2​K1N.\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 KK is explicitly given by

K\displaystyle K =exp⁡(3​β​M1+2​C0​M1)​(exp⁡(2​β​M1)​K2+2​exp⁡(β​M1)​M2​K1)\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{)}
=exp⁡(3​β​M1+2​C0​M1)(exp⁡(2​β​M1)⋅M22d​β+\displaystyle=\exp(3\beta M_{1}+2C_{0}M_{1})\bigg{(}\exp(2\beta M_{1})\cdot\frac{M_{2}}{2}\sqrt{d\beta}\,+
2exp⁡(β​M1)M2⋅βexp⁡(β​M1)M1d​(β+2​C0)2)\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{)}
=exp⁡(5​β​M1+2​C0​M1)​M2​(12​d​β+β​M1​2​d​(β+2​C0))\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{)}
⩽exp⁡(6​β​M1+2​C0​M1)​M2​(12​d​β+2​d​(β+2​C0))\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{)}
⩽exp⁡(6​β​M1+2​C0​M1)​M2​2​d​(2​β+3​C0).\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

12​x+2​x+4​y⩽4​x+6​y,∀x,y⩾0.\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 𝒜N,D⩽exp⁡(β​M1)\mathcal{A}_{N,D}\leqslant\exp(\beta M_{1}) and |ℬN,D|⩽M2|\mathcal{B}_{N,D}|\leqslant M_{2}. Using

Va​(q)⩽32​M1+M1​|q|2V^{a}(q)\leqslant\frac{3}{2}M_{1}+M_{1}|q|^{2}

derived in (3.2), we have

βD​∑j=0D−1Va​(xN​(j​βD))⩽32​β​M1+M1​βD​∑j=0D−1|xN​(j​βD)|2.\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

𝔼ν​[βD​∑j=0D−1Va​(xN​(j​βD))]⩽32​β​M1+M1​𝔼ν​[βD​∑j=0D−1|xN​(j​βD)|2].\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 τ∈[0,β]\tau\in[0,\beta], the value of the continuous loop at τ\tau is

xN​(τ)=∑k=0N−1ξk​ck​(τ),x_{N}(\tau)=\sum_{k=0}^{N-1}\xi_{k}c_{k}(\tau), (A.76)

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

𝔼ν​[|xN​(τ)|2]=∑k=0N−1dωk2+a2​|ck​(τ)|2⩽∑k=0∞dωk2+a2​|ck​(τ)|2.\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 {ck​(τ)}k=0∞\{c_{k}(\tau)\}_{k=0}^{\infty} satisfy

c0​(τ)=1β,|c2​k−1​(τ)|2+|c2​k​(τ)|2=2β,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

𝔼ν​[|xN​(τ)|2]⩽1β​∑k=0∞dωk2+a2=C0β.\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,

𝔼ν​[βD​∑j=0D−1|xN​(j​βD)|2]⩽C0.\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

𝔼ν​[𝒜N,D]⩾exp⁡(−𝔼ν​[βD​∑j=0D−1Va​(xN​(j​βD))]​missing)⩾exp⁡(−32​β​M1−C0​M1​missing).\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 𝒜N\mathcal{A}_{N} and 𝒜N,D\mathcal{A}_{N,D}. By choosing

x=∫0βVa​(xN​(τ))​dτ,y=βD​∑j=0D−1Va​(xN​(j​βD)),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

|𝒜N−𝒜N,D|\displaystyle|\mathcal{A}_{N}-\mathcal{A}_{N,D}| =|exp⁡(−∫0βVa​(xN​(τ))​dτ​missing)−exp⁡(−βD​∑j=0D−1Va​(xN​(j​βD))​missing)|\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{|}
⩽exp⁡(β​M1)​|∫0βVa​(xN​(τ))​dτ−βD​∑j=0D−1Va​(xN​(j​βD))|\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{|}
⩽exp⁡(β​M1)​∑j=0D−1|∫j​βD(j+1)​βD(Va​(xN​(τ))−Va​(xN​(j​βD)))​dτ|\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{|}
⩽exp⁡(β​M1)​∑j=0N−1∫j​βD(j+1)​βD|Va​(xN​(τ))−Va​(xN​(j​βD))|​dτ.\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

|Va​(xN​(τ1))−Va​(xN​(τ2))|⩽2​M1​|xN​(τ1)−xN​(τ2)|​1+|xN​(τ1)|2+|xN​(τ2)|2.\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

𝔼ν​|xN​(τ1)−xN​(τ2)|2⩽d​(2​β+1)​|τ1−τ2|.\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

𝔼ν​[1+|xN​(τ1)|2+|xN​(τ2)|2]⩽1+2​C0β.\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

𝔼ν​|Va​(xN​(τ1))−Va​(xN​(τ2))|\displaystyle~{}~{}~{}~{}\mathbb{E}_{\nu}\big{|}V^{a}(x_{N}(\tau_{1}))-V^{a}(x_{N}(\tau_{2}))\big{|}
⩽2​M1​𝔼ν​|xN​(τ1)−xN​(τ2)|2​𝔼ν​[1+|xN​(τ1)|2+|xN​(τ2)|2]\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{]}}
⩽2​M1​d​(2​β+1)​|τ1−τ2|⋅(1+2​C0β)\displaystyle\leqslant\sqrt{2}M_{1}\sqrt{d(2\beta+1)|\tau_{1}-\tau_{2}|\cdot\bigg{(}1+\frac{2C_{0}}{\beta}\bigg{)}}
=M1​2​dβ​(β+2​C0)​(2​β+1)​|τ1−τ2|.\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 𝔼ν​|𝒜N−𝒜N,D|\mathbb{E}_{\nu}|\mathcal{A}_{N}-\mathcal{A}_{N,D}|,

𝔼ν​|𝒜N−𝒜N,D|\displaystyle~{}~{}~{}~{}\mathbb{E}_{\nu}|\mathcal{A}_{N}-\mathcal{A}_{N,D}|
⩽exp⁡(β​M1)​M1​2​dβ​(β+2​C0)​(2​β+1)​∑j=0D−1∫j​βD(j+1)​βD|τ−j​βD|​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
⩽exp⁡(β​M1)​M1​2​dβ​(β+2​C0)​(2​β+1)⋅β​βD\displaystyle\leqslant\exp(\beta M_{1})M_{1}\sqrt{\frac{2d}{\beta}(\beta+2C_{0})(2\beta+1)}\cdot\beta\sqrt{\beta_{D}}
⩽β​exp⁡(β​M1)​M1​2​d​(β+2​C0)​(2​β+1)D.\displaystyle\leqslant\beta\exp(\beta M_{1})M_{1}\sqrt{\frac{2d(\beta+2C_{0})(2\beta+1)}{D}}. (A.88)

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

|ℬN−ℬN,D|\displaystyle|\mathcal{B}_{N}-\mathcal{B}_{N,D}| =|1β​∫0βO​(xN​(τ))​dτ−1D​∑j=0D−1O​(xN​(j​β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{|}
⩽1β​∑j=0D−1|∫j​βD(j+1)​βDO​(xN​(τ))−O​(xN​(j​βD))|​d​τ\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
⩽M2β​∑j=0D−1∫j​βD(j+1)​βD|xN​(τ)−xN​(j​βD)|.\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

𝔼ν​|ℬN−ℬN,D|\displaystyle\mathbb{E}_{\nu}|\mathcal{B}_{N}-\mathcal{B}_{N,D}| ⩽M2β​∑j=0D−1∫j​βD(j+1)​βD𝔼ν​|xN​(τ)−xN​(j​βD)|2​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
⩽M2β​d​(2​β+1)​∑j=0D−1∫j​βD(j+1)​βD|τ−j​βD|​dτ\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
⩽M2β​d​(2​β+1)⋅β​βD\displaystyle\leqslant\frac{M_{2}}{\beta}\sqrt{d(2\beta+1)}\cdot\beta\sqrt{\beta_{D}}
⩽M2​d​β​(2​β+1)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

|⟨O​(q^)⟩β,N−⟨O​(q^)⟩β,N,D|⩽exp⁡(3​β​M1+2​C0​M1)​exp⁡(2​β​M1)​L2+2​exp⁡(β​M1)​M2​L1N.\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 LL is explicitly given by

L\displaystyle L =exp⁡(3​β​M1+2​C0​M1)​(exp⁡(2​β​M1)​L2+2​exp⁡(β​M1)​M2​L1)\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{)}
=exp⁡(3​β​M1+2​C0​M1)(exp⁡(2​β​M1)⋅M2d​β​(2​β+1)+\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)}\,+
2exp⁡(β​M1)M2⋅βexp⁡(β​M1)M12​d​(β+2​C0))\displaystyle\hskip 85.35826pt2\exp(\beta M_{1})M_{2}\cdot\beta\exp(\beta M_{1})M_{1}\sqrt{2d(\beta+2C_{0})}\bigg{)}
=exp⁡(5​β​M1+2​C0​M1)​M2​(d​β​(2​β+1)+2​β​M1​β+2​C0)\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{)}
⩽2​exp⁡(6​β​M1+2​C0​M1)​2​β+1​(12​d​β+2​d​(β+2​C0))\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{)}
⩽2​exp⁡(6​β​M1+2​C0​M1)​M2​2​d​(2​β+1)​(2​β+3​C0).\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​2​β+1​KL=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.