<EMAIL_ADDRESS>(X. Ye)
<EMAIL_ADDRESS>(Z. Zhou)

# Dimension-free Ergodicity of Path Integral Molecular Dynamics

Xuda Ye  1  Zhennan Zhou  2  11affiliationmark:  Beijing International Center
for Mathematical Research, Peking University, Beijing 100871, China.
22affiliationmark:  Institute for Theoretical Sciences, Westlake University,
Hangzhou, 310030, China.

###### Abstract

The quantum thermal average plays a central role in describing the
thermodynamic properties of a quantum system. Path integral molecular dynamics
(PIMD) is a prevailing approach for computing quantum thermal averages by
approximating the quantum partition function as a classical isomorphism on an
augmented space, enabling efficient classical sampling, but the theoretical
knowledge of the ergodicity of the sampling is lacking. Parallel to the
standard PIMD with $N$ ring polymer beads, we also study the Matsubara mode
PIMD, where the ring polymer is replaced by a continuous loop composed of $N$
Matsubara modes. Utilizing the generalized $\Gamma$ calculus, we prove that
both the Matsubara mode PIMD and the standard PIMD have uniform-in-$N$
ergodicity, i.e., the convergence rate towards the invariant distribution does
not depend on the number of modes or beads $N$.

###### keywords:

quantum thermal average, path integral molecular dynamics, Matsubara modes,
ergodicity, generalized $\Gamma$ calculus.

37A30, 82B31, 81S40

##  1 Introduction

Calculation of the quantum thermal average plays an important role in quantum
physics and quantum chemistry, not only because it fully characterizes the
canonical ensemble of the quantum system, but also because of its wide
applications in describing the thermal properties of complex quantum systems,
including the idea quantum gases [1], chemical reaction rates [2, 3], density
of states of crystals [4], quantum phase transitions [5], etc. However, since
the computational cost of direct discretization methods (finite difference,
pseudospectral methods, etc.) grows exponentially with the spatial dimension
[6, 7], the exact calculation of the quantum thermal average is hardly
affordable in high dimensions. In the past decades, there have been numerous
methods committed to compute the quantum thermal average approximately, and
the path integral molecular dynamics (PIMD) is among the most prevailing ones.

The PIMD is a computational framework to obtain accurate quantum thermal
averages. Using the imaginary time-slicing approach in Feynman path integral
[8], the PIMD maps the quantum system in $\mathbb{R}^{d}$ to a ring polymer of
$N$ beads in $\mathbb{R}^{dN}$, where each bead represents a classical
duplicate of the original quantum system, and adjacent beads are connected by
a harmonic spring potential. When the number of beads $N$ is large enough, the
classical Boltzmann distribution of the ring polymer in $\mathbb{R}^{dN}$ is
expected to yield an accurate approximation of the quantum thermal average.
Since the development of the PIMD in 1970s, this framework has been widely
used in the calculations in the chemical reaction rates [9, 10, 11],
transition state theory [12] and tunneling splittings [13, 14]. Several
variants of the PIMD, the ring polymer molecular dynamics (RPMD) [15, 16], the
centroid molecular dynamics (CMD) [17, 18] and the Matsubara dynamics [19,
20], have been employed to compute the quantum correlation function. Recently,
there have been fruitful studies on designing efficient and accurate
algorithms to enhance the numerical performance of the PIMD [21, 22, 23, 24].
The PIMD can also be applied in multi-electronic-state quantum systems, see
[25, 26] for instance. The general procedure to compute the quantum thermal
average in the PIMD is demonstrated as follows:

  1. 1.

Ring polymer approximation. Pick a positive integer $N$ and map the quantum
system to a ring polymer system of $N$ beads, where each bead is a classical
duplicate of the original system.

  2. 2.

Construction of sampling. Equip the ring polymer system with a stochastic
thermostat (e.g., Langevin, Andersen and Nosé–Hoover) so that the resulting
stochastic process samples the Boltzmann distribution of the ring polymer.

  3. 3.

Stochastic simulation. Evolve the stochastic process with a proper time
discretization method, and approximate the quantum thermal average by the time
average in the long-time simulation.

The PIMD framework based on the ring polymer beads is referred to as the
standard PIMD in this paper.

Although the standard PIMD has become a mature framework to compute the
quantum thermal average, a rigorous mathematical understanding of its
convergence is still sorely lacking. A fundamental question is the dimension-
free ergodicity, i.e., does the stochastic process sampling the Boltzmann
distribution has a convergence rate which does not depend on the number of
beads $N$? The property of uniform-in-$N$ ergodicity has been partially
validated when the potential function is quadratic [27, 28], but a rigorous
justification in the general case is still to be investigated.

An important goal of this paper is to justify the dimension-free ergodicity of
the standard PIMD. Nevertheless, our proof mainly relies on an alternative
PIMD framework, the Matsubara mode PIMD, which is more convenient for
analysis. As its name indicates, this framework is based on the Matsubara
modes [29, 19] rather than the discrete beads of the ring polymer. In the
Matsubara mode PIMD, the ring polymer is expressed as the superpositions of
$N$ Matsubara modes, resulting in a continuous imaginary-time ring polymer
loop rather than a ragged one (see Figure 3.2 of [19] for example). The
Matsubara mode PIMD has a close relation to the standard PIMD, because the
Matsubara mode PIMD with suitable discretization is equivalent to the standard
PIMD [29, 30], except for the minor difference in the mode frequencies. In
fact, the standard PIMD utilizes the normal mode frequencies [19], whose
continuum limit give rise to the Matsubara frequencies. It is also pointed in
[29, 30] that the standard PIMD has better performance than the Matsubara mode
PIMD in computing the quantum thermal average.

The main contribution of this paper is that we prove the uniform-in-$N$
ergodicity of the Matsubara mode PIMD and the standard PIMD, namely, the
convergence rate towards the invariant distribution does not depend on $N$,
which stands for the number of modes in the Matsubara mode PIMD and the number
of beads in the standard PIMD. The authors believe that this is the first
dimension-free egordicity result of the PIMD in the form of the underdamped
Langevin dynamics. Our proof strategies include the Bakry–Émery calculus [31]
and the recently developed generalized $\Gamma$ calculus [32, 33]. It should
be stressed that the convergence rate can still depend on the potential
function $V(x)$ and the temperature, and in particular the convergence rate is
exponentially small in the low temperature limit. This a a common feature of
the Bakry–Émery calculus with a non-convex potential function, referring to
Proposition 5.1.6 of [31].

The paper is organized as follows. In Section 2 we derive the Matsubara mode
PIMD from the Feynman path integral representation of the quantum thermal
average. In Section 3 we prove the uniform-in-$N$ ergodicity the Matsubara
mode PIMD and the standard PIMD. In Section 4 we implement the numerical tests
to show the performance of the Matsubara mode PIMD and the standard PIMD in
computing the quantum thermal average and validate the dimension-free
ergodicity.

##  2 Derivation of the Matsubara mode PIMD

In this section we derive the Matsubara mode PIMD from the Feynman path
integral representation of the quantum thermal average. First, we show that
the quantum equilibrium can be equivalently interpreted as the classical
Boltzmann distribution of a continuous loop representing the ring polymer in
the position space. Second, we represent the energy function
$\mathcal{E}(\xi)$ of the ring polymer loop in terms of the Matsubara
coordinates $\xi$. Finally, we construct an infinite-dimensional underdamped
Langevin dynamics to sample the distribution $\exp{-\mathcal{E}(\xi)}$,
formulating the Matsubara mode PIMD.

Note that the our derivation does not rely on the normal mode transform, and
is thus different from the traditional ways to generate the Matsubara modes
[19].

###  2.1 Quantum equilibrium as the Boltzmann distribution of a continuous
loop

We consider the quantum system in $\mathbb{R}^{d}$ given by the Hamiltonian

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

where $\hat{q}$ and $\hat{p}$ are the position and momentum operators in
$\mathbb{R}^{d}$, and $V(\cdot)$ is a real-valued potential function in
$\mathbb{R}^{d}$. When the quantum system is at a constant temperature
$T=1/\beta$, the state of the system can be described by the canonical
ensemble with the density operator $e^{-\beta\hat{H}}$, and thus the partition
function $\mathcal{Z}=\mathrm{Tr}[e^{-\beta\hat{H}}]$. The quantum thermal
average of the system means the canonical average of an observable operator
$O(\hat{q})$, namely

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

Here, we assume the observable operator $O(\hat{q})$ depends only on the
position operator $\hat{q}$, where $O(\cdot)$ is a real-valued function in
$\mathbb{R}^{d}$.

Utilizing the Feynman path integral [34, 8], the partition function can be
interpreted as the integral with respect to a continuous loop $x(\tau)$ in
$\mathbb{R}^{d}$ parameterized by $\tau\in[0,\beta]$:

$\mathcal{Z}=\mathrm{Tr}[e^{-\beta\hat{H}}]=\int\exp\\{-\mathcal{E}(x)\\}\mathcal{D}[x],$

where $\mathcal{D}[x]$ is the formal Lebesgue measure of the continuous loop
in $\mathbb{R}^{d}$, and the energy functional $\mathcal{E}(x)$ is defined by

$\mathcal{E}(x)=\frac{1}{2}\int_{0}^{\beta}|x^{\prime}(\tau)|^{2}\mathrm{d}\tau+\int_{0}^{\beta}V(x(\tau))\mathrm{d}\tau.$
(3)

As a consequence, the quantum canonical ensemble $e^{-\beta\hat{H}}$ is
exactly mapped to a classical Boltzmann distribution of the continuous loop
$x(\tau)$ with the energy $\mathcal{E}(x)$. Finally, the quantum thermal
average $\langle{O(\hat{q})}\rangle_{\beta}$ can be formally written as

$\langle{O(\hat{q})}\rangle_{\beta}=\frac{1}{\mathcal{Z}}\int\bigg{[}\frac{1}{\beta}\int_{0}^{\beta}O(x(\tau))\mathrm{d}\tau\bigg{]}\exp\\{-\mathcal{E}(x)\\}\mathcal{D}[x].$
(4)

###  2.2 Represent the energy function with Matsubara modes

Consider the following eigenvalue problem with periodic boundary conditions:

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

where the eigenvalues and eigenfunctions are explicitly given by

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

The eigenfunctions $\\{c_{k}(\tau)\\}_{k=0}^{\infty}$ form the orthonormal
Fourier basis of $L^{2}([0,\beta];\mathbb{R})$, and thus any continuous loop
$x(\cdot)$ can be uniquely represented as

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

where $\\{\xi_{k}\\}_{k=0}^{\infty}$ in $\mathbb{R}^{d}$ are the coordinates
of $x(\cdot)$ in different Matsubara modes. Conversely, for a given continuous
loop $x(\cdot)$, the Matsubara coordinates are calculated from

$\xi_{k}=\int_{0}^{\beta}x(\tau)c_{k}(\tau)\mathrm{d}\tau,~{}~{}~{}~{}k=0,1,2,\cdots.$

Utilizing the Matsubara coordinates, the energy function $\mathcal{E}(x)$ in
(3) is written as

$\mathcal{E}(\xi)=\frac{1}{2}\sum_{k=0}^{\infty}\omega_{k}^{2}|\xi_{k}|^{2}+\mathcal{V}(\xi),~{}~{}~{}\text{with}~{}~{}\mathcal{V}(\xi):=\int_{0}^{\beta}V\bigg{(}\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{)}\mathrm{d}\tau,$
(5)

and the target Boltzmann distribution is formally given by
$\exp\\{-\mathcal{E}(\xi)\\}$.

###  2.3 Formulation of the Matsubara mode PIMD

In the following, we construct an infinite-dimensional Langevin dynamics to
sample the Boltzmann distribution $\exp\\{-\mathcal{E}(\xi)\\}$. Without any
preconditioning, the vanilla underdamped Langevin dynamics for sampling
$\exp\\{-\mathcal{E}(\xi)\\}$ reads

$\left\\{\begin{aligned} \dot{\xi}_{k}&=\eta_{k},\\\
\dot{\eta}_{k}&=-\omega_{k}^{2}\xi_{k}-\int_{0}^{\beta}\nabla
V(x(\tau))c_{k}(\tau)\mathrm{d}\tau-\gamma\eta_{k}+\sqrt{2\gamma}\dot{B}_{k},\end{aligned}\right.$
(6)

where $\gamma>0$ is the damping rate on each mode,
$\\{\eta_{k}\\}_{k=0}^{\infty}$ are the auxiliary velocity variables in
$\mathbb{R}^{d}$, and $\\{B_{k}\\}_{k=0}^{\infty}$ are independent Brownian
motions in $\mathbb{R}^{d}$. However, the high-frequency modes in the energy
function $\mathcal{E}(\xi)$ poses the infamous stiffness problem [27, 35] in
the time discretization of (6): the time step needs to be extremely small to
stabilize of the high-frequency part of the dynamics.

In this paper, we employ the preconditioning [27, 36] to scale the internal
frequencies of the Matsubara modes. Introduce $a>0$ and rewrite the energy
function $\mathcal{E}(\xi)$ as

$\mathcal{E}(\xi)=\frac{1}{2}\sum_{k=0}^{\infty}(\omega_{k}^{2}+a)|\xi_{k}|^{2}+\mathcal{V}^{a}(\xi),~{}~{}~{}\text{with}~{}~{}\mathcal{V}^{a}(\xi):=\int_{0}^{\beta}V^{a}\bigg{(}\sum_{k=0}^{\infty}\xi_{k}c_{k}(\tau)\bigg{)}\mathrm{d}\tau,$
(7)

where the potential function $V^{a}(q):=V(q)-a|q|^{2}/2$. The positive
coefficient $\omega_{k}^{2}+a$ on each mode indicates that it is reasonable to
apply the preconditioning

$\left\\{\begin{aligned} \dot{\xi}_{k}&=\eta_{k},\\\
\dot{\eta}_{k}&=-\xi_{k}-\frac{1}{\omega_{k}^{2}+a}\int_{0}^{\beta}\nabla
V^{a}(x(\tau))c_{k}(\tau)\mathrm{d}\tau-\gamma\eta_{k}+\sqrt{\frac{2\gamma}{\omega_{k}^{2}+a}}\dot{B}_{k},\end{aligned}\right.$
(8)

which is an infinite-dimensional Langevin dynamics of the mode coordinates
$\\{\xi_{k}\\}_{k=0}^{\infty}$.

######  Remark 2.1.

If we neglect the gradient term in (8), the dynamics of $(\xi_{k},\eta_{k})$
reads

$\left\\{\begin{aligned} \dot{\xi}_{k}&=\eta_{k},\\\
\dot{\eta}_{k}&=-\xi_{k}-\gamma\eta_{k}+\sqrt{\frac{2\gamma}{\omega_{k}^{2}+a}}\dot{B}_{k},\end{aligned}\right.$

which is a underdamped Langevin dynamics in
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ with the invariant distribution

$\exp\\{-\frac{\omega_{k}^{2}+a}{2}\big{(}|\xi_{k}|^{2}+|\eta_{k}|^{2}\big{)}\bigg{\\}}.$

Clearly, the preconditioning removes the stiffness of all Matsubara modes, and
$\mathcal{O}(1)$ time step is applicable for the time discretization.

######  Remark 2.2.

The constant $a>0$ is to ensure the frequency of the centroid mode $(k=0)$ is
nonzero, so that the preconditioning can be applied on all Matsubara modes. An
alternative preconditioning scheme for (6) without introducing $a$ is defined
as:

$\left\\{\begin{aligned} \xi_{0}&=\eta_{0},\\\
\eta_{0}&=-\frac{1}{\sqrt{\beta}}\int_{0}^{\beta}\nabla
V(x(\tau))\mathrm{d}\tau-\gamma\eta_{0}+\sqrt{2\gamma}\dot{B}_{0},\end{aligned}\right.$
(9) $\left\\{\begin{aligned} \xi_{k}&=\eta_{k},\\\
\eta_{k}&=-\xi_{k}-\frac{1}{\omega_{k}^{2}}\int_{0}^{\beta}\nabla
V(x(\tau))c_{k}(\tau)\mathrm{d}\tau-\gamma\eta_{k}+\frac{\sqrt{2\gamma}}{\omega_{k}}\dot{B}_{k},\end{aligned}\right.~{}~{}~{}~{}k=1,2,\cdots.$
(10)

Moreover, the preconditioned Matsubara mode PIMD (9) and (10) can be employed
in quantum systems defined on a periodic torus space $\mathbb{T}^{d}$.

###  2.4 Implementation of the Matsubara mode PIMD in finite dimensions

At this stage, the infinite-dimensional Matsubara mode PIMD (8) is completely
formal, and neither its existence or well-posedness is not justified. For the
purpose of numerical simulation, we need to truncate the number of Matsubara
modes in (8) to a finite integer $N$. In this case the continuous loop
$x(\tau)$ is truncated as

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

and the potential energy of the loop in (7) is replaced by

$\mathcal{V}_{N}^{a}(\xi):=\int_{0}^{\beta}V^{a}(x_{N}(\tau))\mathrm{d}\tau=\int_{0}^{\beta}V^{a}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(\tau)\bigg{)}\mathrm{d}\tau,~{}~{}~{}~{}\xi\in\mathbb{R}^{dN}.$
(11)

However, even for a finite $N$, the integral in (11) does not have an explicit
expression, and the numerical integration is required to approximate
$\mathcal{V}_{N}^{a}(\xi)$. Let $D\in\mathbb{N}$ be the discretization size in
the interval $[0,\beta]$, then we obtain the approximation

$\mathcal{V}_{N}^{a}(\xi)\approx\mathcal{V}_{N,D}^{a}(\xi):=\beta_{D}\sum_{j=0}^{D-1}V^{a}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{D})\bigg{)},~{}~{}~{}~{}\xi\in\mathbb{R}^{dN},$
(12)

where $\beta_{D}=\beta/D$. The gradients of $\mathcal{V}^{a}_{N,D}(\xi)$ with
respect to each $\xi_{k}$ are correspondingly

$\nabla_{\xi_{k}}\mathcal{V}_{N,D}^{a}(\xi)=\beta_{D}\sum_{j=0}^{N-1}\nabla
V^{a}(x_{N}(j\beta_{D}))c_{k}(j\beta_{D}),~{}~{}~{}~{}k=0,1,\cdots,N-1.$

######  Remark 2.3.

We do not select the high-order numerical integration schemes mainly because
the continuous loop $x_{N}(\tau)$ becomes ragged when the number of modes $N$
is large. The low-regularity of $x_{N}(\tau)$ makes it ineffectual to apply
high-order integration.

With the potential function $\mathcal{V}_{N,D}^{a}(\xi)$, the Matsubara mode
PIMD (8) is approximated as

$\left\\{\begin{aligned} \dot{\xi}_{k}&=\eta_{k},\\\
\dot{\eta}_{k}&=-\xi_{k}-\frac{\beta_{D}}{\omega_{k}^{2}+a}\sum_{j=0}^{D-1}\nabla
V^{a}(x_{N}(j\beta_{D}))c_{k}(j\beta_{D})-\gamma\eta_{k}+\sqrt{\frac{2\gamma}{\omega_{k}^{2}+a}}\dot{B}_{k},\end{aligned}\right.~{}~{}~{}k=0,1,\cdots,N-1,$
(13)

whose invariant distribution is the classical Boltzmann distribution:

$\pi_{N,D}(\xi)\propto\exp\\{-\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k}^{2}+a)|\xi_{k}|^{2}-\beta_{D}\sum_{j=0}^{D-1}V^{a}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{D})\bigg{)}\bigg{\\}}.$
(14)

A typical choice of the discretization size $D$ is the number of modes $N$, so
that the computational cost per time step is equal to $\mathcal{O}(N\log N)$
in the Fast Fourier Transform. It is natural from (4) that the quantum thermal
average $\langle{O(\hat{q})}\rangle_{\beta}$ is approximated as

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

The long-time simulation of the Matsubara mode PIMD (13) at finite $N$ and $D$
then provides an accurate estimate of the statistical average in the
distribution $\pi_{N,D}(\xi)$, which is denoted by
$\langle{O(\hat{q})}\rangle_{\beta,N,D}$.

###  2.5 Relation to the standard PIMD

We demonstrate the connection between the Matsubara mode PIMD and the standard
PIMD. Suppose the number of modes $N$ is an odd integer and we view
$\\{\xi_{k}\\}_{k=0}^{N-1}$ as the normal mode coordinates [37] in the
standard PIMD, then the Boltzmann distribution of the ring polymer beads in
the standard PIMD is expressed as

$\exp\\{-\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k,N}^{2}+a)|\xi_{k}|^{2}-\beta_{N}\sum_{j=0}^{N-1}V^{a}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{N})\bigg{)}\bigg{\\}},$
(15)

where $\\{\omega_{k,N}\\}_{k=0}^{N-1}$ are the normal mode frequencies

$\omega_{0}=0,~{}~{}~{}~{}\omega_{2k-1,N}=\omega_{2k,N}=\frac{2}{\beta_{N}}\sin\bigg(\frac{k\pi}{N}\bigg{missing}),~{}~{}~{}~{}k=1,\cdots,\frac{N-1}{2}.$

Furthermore, the preconditioned Langevin dynamics for sampling (15) is written
as

$\left\\{\begin{aligned} \dot{\xi}_{k}&=\eta_{k},\\\
\dot{\eta}_{k}&=-\xi_{k}-\frac{\beta_{D}}{\omega_{k,N}^{2}+a}\sum_{j=0}^{D-1}\nabla
V^{a}(x_{N}(j\beta_{D}))c_{k}(j\beta_{D})-\gamma\eta_{k}+\sqrt{\frac{2\gamma}{\omega_{k,N}^{2}+a}}\dot{B}_{k},\end{aligned}\right.~{}~{}~{}k=0,1,\cdots,N-1.$
(16)

Comparing the dynamics (13) and (16), we observe the standard PIMD and the
Matsubara mode PIMD (with odd $D=N$) are equivalent except for the mode
frequencies. Since the normal mode frequencies satisfy

$\lim_{N\rightarrow\infty}\omega_{k,N}=\omega_{k},~{}~{}~{}~{}k=0,1,2,\cdots,$

we conclude that the Matsubara mode PIMD (13) and the standard PIMD (16) have
the identical continuum limit, namely the infinite-dimensional Langevin
dynamics (8). See Appendix C for the detailed derivation of the normal mode
coordinates.

##  3 Uniform-in-$N$ ergodicity of the Matsubara mode PIMD

In this section, we prove the uniform-in-$N$ ergodicity of the Matsubara mode
PIMD. First, we present the assumptions and ergodicity results in Section 3.1.
Then, we prove the uniform-in-$N$ ergodicity in the overdamped and underdamped
cases in Sections 3.2 and 3.3, respectively. Finally, in Section 3.4, we
demonstrate that our proof of the uniform-in-$N$ ergodicity in the Matsubara
mode PIMD applies to the standard PIMD.

###  3.1 Assumptions and ergodicity results

For the convenience of analysis, in the Matsubara mode PIMD (13) we choose the
damping rate $\gamma=1$, yielding the following underdamped Langevin dynamics:

$\left\\{\begin{aligned} \dot{\xi}_{k}&=\eta_{k},\\\
\dot{\eta}_{k}&=-\xi_{k}-\frac{\beta_{D}}{\omega_{k}^{2}+a}\sum_{j=0}^{D-1}\nabla
V^{a}(x_{N}(j\beta_{D}))c_{k}(j\beta_{D})-\eta_{k}+\sqrt{\frac{2}{\omega_{k}^{2}+a}}\dot{B}_{k}.\end{aligned}\right.$
(17)

A closely related dynamics is the overdamped version of (17):

$\dot{\xi}_{k}=-\xi_{k}-\frac{\beta_{D}}{\omega_{k}^{2}+a}\sum_{j=0}^{D-1}\nabla
V^{a}(x_{N}(j\beta_{D}))c_{k}(j\beta_{D})+\sqrt{\frac{2}{\omega_{k}^{2}+a}}\dot{B}_{k},$
(18)

which can be viewed as the overdamped limit of (17) as
$\gamma\rightarrow\infty$ [38]. The invariant distribution of (18) is
$\pi_{N,D}(\xi)$ defined in (14), and the invariant distribution of (17) is

$\mu_{N,D}(\xi,\eta)\propto\exp\bigg(-\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k}^{2}+a)(|\xi_{k}|^{2}+|\eta_{k}|^{2})-\beta_{D}\sum_{j=0}^{D-1}V^{a}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{D})\bigg{)}\bigg{missing}).$
(19)

Clearly, $\pi_{N,D}(\xi)$ is the marginal distribution of $\mu_{N,D}(\xi)$ in
the variable $\xi\in\mathbb{R}^{dN}$.

Under appropriate conditions on the potential $V(q)$, we prove that both the
overdamped and underdamped Matsubara mode PIMD in (18) and (17) possess the
uniform-in-$N$ ergodicity, namely, the convergence rate to the invariant
distribution does not depend on the number of modes $N$.

Before we conduct a detailed discussion on these results, we enumerate all the
assumptions required in the proof, mainly on the potential function $V(q)$ in
$\mathbb{R}^{d}$.

######  Assumption.

Given $a>0$, the potential function

$V^{a}(q)=V(q)-\frac{a}{2}|q|^{2},~{}~{}~{}~{}q\in\mathbb{R}^{d}$

is twice differentiable in $\mathbb{R}^{d}$, and for some constants
$M_{1},M_{2}\geqslant 0$:

  1. (i)

$V^{a}(q)$ can be decomposed as $V^{c}(q)+V^{b}(q)$, where
$\nabla^{2}V^{c}(q)\succcurlyeq O_{d}$ and $|V^{b}(q)|\leqslant M_{1}$ for any
$q\in\mathbb{R}^{d}$.

  2. (ii)

$-M_{2}I_{d}\preccurlyeq\nabla^{2}V^{a}(q)\preccurlyeq M_{2}I_{d}$ for any
$q\in\mathbb{R}^{d}$.

A special assumption required in the underdamped case is:

  1. (iii)

The number of modes $N$ is no larger than the discretization size $D$, namely,
$N\leqslant D$.

Here, $I_{d}$ and $O_{d}$ are the identity and zero matrix in
$\mathbb{R}^{d\times d}$. Assumption (i) can be shortly interpreted as:
$V^{c}(q)$ is globally convex and $V^{b}(q)$ is globally bounded.

Next, we display the ergodicity results of the Matsubara mode PIMD in Table 1.

Uniform-in-$N$ ergodicity of the Matsubara mode PIMD
---
dynamics | overdamped (18) | underdamped (17)
assumption | (i) $\Rightarrow$ Theorem 3.3 | (i)(ii)(iii) $\Rightarrow$ Theorem 3.5
distribution |  $\pi_{N,D}(\xi)$ in (14) |  $\mu_{N,D}(\xi,\eta)$ in (19)
ergodicity | $\mathrm{\,Ent}_{\pi_{N,D}}(P_{t}f)\leqslant e^{-2\lambda_{1}t}\mathrm{\,Ent}_{\pi_{N,D}}(f)$ | $W_{\mu_{N,D}}(P_{t}f)\leqslant e^{-2\lambda_{2}t}W_{\mu_{N,D}}(f)$

Table 1: The uniform-in-$N$ ergodicity of the Matsubara mode PIMD. The
relative entropy $\mathrm{\,Ent}_{\pi_{N,D}}(f)$ and the entropy-like quantity
$W_{\mu_{N,D}}(f)$ are defined in (20) and (21).

In Table 1, $\mathrm{\,Ent}_{\pi_{N,D}}(f)$ is the relative entropy of the
density function $f(\xi)$ in $\mathbb{R}^{dN}$:

$\mathrm{\,Ent}_{\pi_{N,D}}(f):=\int_{\mathbb{R}^{dN}}f\log
f\mathrm{d}\pi_{N,D}-\int_{\mathbb{R}^{dN}}f\mathrm{d}\pi_{N,D}\log\int_{\mathbb{R}^{dN}}f\mathrm{d}\pi_{N,D},$
(20)

and $W_{\mu_{N,D}}(f)$ is the entropy-like quantity of the density function
$f(\xi,\eta)$ in $\mathbb{R}^{2dN}$:

$W_{\mu_{N,D}}(f):=\bigg{(}\frac{M_{2}^{2}}{a^{2}}+1\bigg{)}\mathrm{\,Ent}_{\mu_{N,D}}(f)+\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}+a}\int_{\mathbb{R}^{2dN}}\frac{|\nabla_{\eta_{k}}f-\nabla_{\xi_{k}}f|^{2}+|\nabla_{\eta_{k}}f|^{2}}{f}\mathrm{d}\mu_{N,D}.$
(21)

Furthermore, the convergence rates $\lambda_{1}$ and $\lambda_{2}$ are given
by

$\lambda_{1}=\exp(-4\beta
M_{1}),~{}~{}~{}~{}\lambda_{2}=\frac{a^{2}}{3M_{2}^{2}+5a^{2}}\exp(-4\beta
M_{1}),$

######  Remark 3.1.

The convergence rates $\lambda_{2}<\lambda_{1}$, and they are exponentially
small in the low temperature limit $\beta\rightarrow\infty$ if the constant
$M_{1}>0$, which characterizes the non-convexity of the potential $V(q)$. This
is because we utilize the bounded perturbation of the log-Sobolev inequalities
(Proposition 5.1.6 of [31]). However, $\lambda_{2}<\lambda_{1}$ does not imply
the overdamped Matsubara mode PIMD converges faster than the underdamped one
in practical simulation. With the hypocoercivity methods, it is proved that
introducing auxiliary velocity variables accelerates the convergence of the
overdamped Langevin dynamics [38].

######  Remark 3.2.

Assumption (iii) comes from the discrete orthogonal condition (36):

$\sum_{j=0}^{D-1}c_{k}(j\beta_{D})c_{l}(j\beta_{D})=0,~{}~{}~{}~{}0\leqslant
k<l\leqslant N-1,$

which in general is incorrect when $N>D$. The authors conjecture that the
uniform-in-$N$ ergodicity of the underdamped Matsubara mode PIMD (17) also
holds true when $N>D$, nevertheless no viable strategy is available.

###  3.2 Uniform ergodicity of overdamped Matsubara mode PIMD

We prove the uniform-in-$N$ ergodicity of (18) in the relative entropy using
the Bakry–Émery calculus [31]. The log-Sobolev inequality for the distribution
$\pi_{N,D}(\xi)$ produces an explicit convergence rate in high-dimensions.

######  Theorem 3.3.

Under Assumption (i), let $(P_{t})_{t\geqslant 0}$ be the Markov semigroup of
the overdamped Matsubara mode PIMD (18), then for any positive smooth function
$f(\xi)$ in $\mathbb{R}^{dN}$,

$\mathrm{\,Ent}_{\pi_{N,D}}(P_{t}f)\leqslant\exp(-2\lambda_{1}t)\mathrm{\,Ent}_{\pi_{N,D}}(f),~{}~{}~{}~{}\forall
t\geqslant 0,$

where the convergence rate $\lambda_{1}=\exp(-4\beta M_{1})$.

######  Proof 3.4.

We study the overdamped Matsubara mode PIMD driven by the convex potential
$V^{c}(q)$, and prove the corresponding uniform-in-$N$ log-Sobolev inequality.
Utilizing the bounded perturbation, we obtain the uniform-in-$N$ log-Sobolev
inequality for (18).
1\. Ergodicity of overdamped Matsubara mode PIMD driven by $V^{c}(q)$
For notational convenience, introduce the potential function of the convex
part $V^{c}(q)$:

$\mathcal{V}_{N,D}^{c}(\xi):=\beta_{D}\sum_{j=0}^{D-1}V^{c}(x_{N}(j\beta_{D}))=\beta_{D}\sum_{j=0}^{D-1}V^{c}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{D})\bigg{)},$

Since $V^{c}(q)$ is globally convex in $\mathbb{R}^{d}$, it is easy to deduce
$\mathcal{V}_{N,D}^{c}(\xi)$ is globally convex in the mode coordinates
$\xi=\\{\xi_{k}\\}_{k=0}^{N-1}$ (see Lemma A.1 in Appendix A). Next consider
the overdamped Matsubara mode PIMD driven by $\mathcal{V}_{N,D}^{c}(\xi)$:

$\dot{\xi}_{k}=-\xi_{k}-\frac{1}{\omega_{k}^{2}+a}\nabla_{\xi_{k}}\mathcal{V}_{N,D}^{c}(\xi)+\sqrt{\frac{2}{\omega_{k}^{2}+a}}\dot{B}_{k},~{}~{}~{}~{}k=0,1,\cdots,N-1,$
(22)

It is easy to see the generator of (22) is given by

$L^{c}=-\sum_{k=0}^{N-1}\bigg{(}\xi_{k}+\frac{1}{\omega_{k}^{2}+a}\nabla_{\xi_{k}}\mathcal{V}_{N,D}^{c}(\xi)\bigg{)}\cdot\nabla_{\xi_{k}}+\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}+a}\Delta_{\xi_{k}},$

and the invariant distribution of (22) is

$\pi_{N,D}^{c}(\xi)\propto\exp\\{-\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k}^{2}+a)|\xi_{k}|^{2}-\beta_{D}\sum_{j=0}^{D-1}V^{c}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{D})\bigg{)}\bigg{\\}}.$
(23)

To establish the log-Sobolev inequality for the distribution
$\pi_{N,D}^{c}(\xi)$, we compute the carré du champ operator $\Gamma_{1}(f,f)$
and the iterated operator $\Gamma_{2}(f,f)$ corresponding to the generator
$L^{c}$ (see Definition 1.4.2 and Equation (1.16.1) in [31]). By direct
calculation,

$\displaystyle\Gamma_{1}(f,g)$
$\displaystyle=\sum_{k=0}^{N-1}\frac{\nabla_{\xi_{k}}f\cdot\nabla_{\xi_{k}}g}{\omega_{k}^{2}+a},$
$\displaystyle\Gamma_{2}(f,g)$
$\displaystyle=\sum_{k,l=0}^{N-1}\frac{\nabla_{\xi_{k}\xi_{l}}^{2}f:\nabla_{\xi_{k}\xi_{l}}^{2}g}{(\omega_{k}^{2}+a)(\omega_{l}^{2}+a)}+\sum_{k=0}^{N-1}\frac{\nabla_{\xi_{k}}f\cdot\nabla_{\xi_{k}}g}{\omega_{k}^{2}+a}+\sum_{k,l=0}^{N-1}\frac{\nabla_{\xi_{k}}f\cdot\nabla^{2}_{\xi_{k}\xi_{l}}\mathcal{V}_{N,D}^{c}(\xi)\cdot\nabla_{\xi_{l}}g}{(\omega_{k}^{2}+a)(\omega_{l}^{2}+a)}.$

Here, the dot product $u\cdot B\cdot v$ for $u\in\mathbb{R}^{d}$,
$B\in\mathbb{R}^{d\times d}$ and $v\in\mathbb{R}^{d}$ means

$u\cdot B\cdot v=u^{\mathrm{T}}Bv=\sum_{p,q=1}^{d}u_{p}B_{pq}v_{q},$

and the double dot product $A:B$ for $A\in\mathbb{R}^{d\times d}$ and
$B\in\mathbb{R}^{d\times d}$ means

$A:B=\mathrm{Tr}[A^{\mathrm{T}}B]=\sum_{p,q=1}^{d}A_{pq}B_{pq}.$

Utilizing the convexity of the potential function
$\mathcal{V}_{N,D}^{c}(\xi)$, we obtain

$\Gamma_{2}(f,f)\geqslant\sum_{k=0}^{N-1}\frac{|\nabla_{\xi_{k}}f|^{2}}{\omega_{k}^{2}+a}+\sum_{k,j=0}^{N-1}\frac{\nabla_{\xi_{k}}f\cdot\nabla^{2}_{\xi_{k}\xi_{j}}\mathcal{V}_{N,D}^{c}(\xi)\cdot\nabla_{\xi_{j}}f}{(\omega_{k}^{2}+a)(\omega_{j}^{2}+a)}\geqslant\sum_{k=0}^{N-1}\frac{|\nabla_{\xi_{k}}f|^{2}}{\omega_{k}^{2}+a}=\Gamma_{1}(f,f),$

hence the log-Sobolev constant for $\pi_{N,D}^{c}(\xi)$ is 1 according to
Proposition 5.7.1 of [31]. The log-Sobolev inequality for $\pi_{N,D}^{c}(\xi)$
then reads

$\mathrm{\,Ent}_{\pi_{N,D}^{c}}(f)\leqslant\frac{1}{2}\int_{\mathbb{R}^{dN}}\frac{\Gamma_{1}(f,f)}{f}\mathrm{d}\pi_{N,D}^{c}=\frac{1}{2}\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}+a}\int_{\mathbb{R}^{dN}}\frac{|\nabla_{\xi_{k}}f|^{2}}{f}\mathrm{d}\pi_{N,D}^{c}.$
(24)

Note that the relative entropy here is defined with respect to the
distribution $\pi_{N,D}^{c}(\xi)$ with $V^{c}(q)$ rather than the distribution
$\pi_{N,D}(\xi)$ with $V^{a}(q)$.
3\. Ergodicity of overdamped Matsubara mode PIMD driven by $V^{a}(q)$
Let $Z_{N,D}$ and $Z_{N,D}^{c}$ be the normalization constants of the
distributions $\pi_{N,D}(\xi)$ and $\pi_{N,D}^{c}(\xi)$,

$\displaystyle Z_{N,D}$
$\displaystyle=\int_{\mathbb{R}^{dN}}\exp\bigg(-\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k}^{2}+a)|\xi_{k}|^{2}-\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D}))\bigg{missing})\mathrm{d}\xi,$
$\displaystyle Z_{N,D}^{c}$
$\displaystyle=\int_{\mathbb{R}^{dN}}\exp\bigg(-\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k}^{2}+a)|\xi_{k}|^{2}-\beta_{D}\sum_{j=0}^{D-1}V^{c}(x_{N}(\beta_{D}))\bigg{missing})\mathrm{d}\xi.$

Using the inequality $|V^{a}(q)-V^{c}(q)|=|V^{b}(q)|\leqslant M_{1}$, we have

$|\mathcal{V}_{N,D}^{a}(\xi)-\mathcal{V}_{N,D}^{c}(\xi)|\leqslant\beta_{D}\sum_{j=0}^{D-1}\big{|}V^{b}(x_{N}(j\beta_{D}))\big{|}\leqslant\beta
M_{1},$

and thus the constants $Z_{N,D}$ and $Z_{N,D}^{c}$ satisfy

$\frac{Z_{N,D}}{Z_{N,D}^{c}}\in\big{[}\exp(-\beta M_{1}),\exp(\beta
M_{1})\big{]}.$

As a result, the density functions $\pi_{N,D}(\xi)$ and $\pi_{N,D}^{c}(\xi)$
satisfy

$\frac{\pi_{N,D}^{c}(\xi)}{\pi_{N,D}(\xi)}=\frac{Z_{N,D}}{Z_{N,D}^{c}}\exp\bigg(\beta_{D}\sum_{j=0}^{D-1}V^{b}(x_{N}(j\beta_{D}))\bigg{missing})\in\big{[}\exp(-2\beta
M_{1}),\exp(2\beta M_{1})\big{]}.$

Using the bounded perturbation (Proposition 5.1.6 of [31]), we obtain from
(24) that

$\exp(-4\beta
M_{1})\mathrm{\,Ent}_{\pi_{N,D}}(f)\leqslant\frac{1}{2}\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}+a}\int_{\mathbb{R}^{dN}}\frac{|\nabla_{\xi_{k}}f|^{2}}{f}\mathrm{d}\pi_{N,D}.$

Hence for the rate $\lambda_{1}=\exp(-4\beta M_{1})$, the relative entropy has
exponential decay,

$\mathrm{\,Ent}_{\pi_{N,D}}(P_{t}f)\leqslant\exp(-2\lambda_{1}t)\mathrm{\,Ent}_{\pi_{N,D}}(f),~{}~{}~{}~{}\forall
t\geqslant 0,$

for any positive smooth function $f(\xi)$.

In particular, the convergence rate $\lambda_{1}$ does not depend on the
number of modes $N$ or the discretization size $D$, hence we conclude the
uniform-in-$N$ ergodicity of the overdamped Matsubara mode PIMD (18).

###  3.3 Uniform ergodicity of underdamped Matsubara mode PIMD

We prove the uniform-in-$N$ ergodicity of (17) in the entropy-like quantity
(21), and the main technique is the generalized $\Gamma$ calculus developed in
[32, 33]. The generalized $\Gamma$ calculus is an extension of the Bakry–Émery
calculus and can be applied on stochastic processes with degenerate
diffusions. The generalized $\Gamma$ calculus is largely inspired from the
hypocoercivity theory [39] of Villani, and is able to produce an explicit
convergence rate in the relative entropy rather than $H^{1}$ or $L^{2}$. For
convenience, we present a brief review of the theory of the generalized
$\Gamma$ calculus in Appendix C.

Recall that the invariant distribution of the underdamped Matsubara mode PIMD
(17) is $\mu_{N,D}(\xi,\eta)$ defined in (19).

######  Theorem 3.5.

Under Assumptions (i)(ii)(iii), let $(P_{t})_{t\geqslant 0}$ be the Markov
semigroup of the underdamped Matsubara mode PIMD (17), then for any positive
smooth function $f(\xi,\eta)$ in $\mathbb{R}^{2dN}$,

$W_{\mu_{N,D}}(P_{t}f)\leqslant\exp(-2\lambda_{2}t)W_{\mu_{N,D}}(f),~{}~{}~{}~{}\forall
t\geqslant 0,$

where the convergence rate
$\lambda_{2}=\frac{a^{2}}{3M_{2}^{2}+5a^{2}}\exp(-4\beta M_{1})$.

######  Proof 3.6.

The proof is accomplished in several steps. First, we establish the uniform-
in-$N$ log-Sobolev inequality for the distribution $\mu_{N,D}(\xi,\eta)$.
Second, introduce the functions

$\Phi_{1}(f)=f\log
f,~{}~{}~{}~{}\Phi_{2}(f)=\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}+a}\frac{|\nabla_{\eta_{k}}f-\nabla_{\xi_{k}}f|^{2}+|\nabla_{\eta_{k}}f|^{2}}{f}$

and compute the generalized $\Gamma$ operators $\Gamma_{\Phi_{1}}(f)$ and
$\Gamma_{\Phi_{2}}(f)$. Finally, by validating the generalized curvature-
dimension condition

$\Gamma_{\Phi_{2}}(f)-\frac{1}{2}\Phi_{2}(f)+\bigg{(}\frac{M_{2}^{2}}{a^{2}}+1\bigg{)}\Gamma_{\Phi_{1}}(f)\geqslant
0,$

we can apply Theorem C.9 to prove the exponential decay of the quantity
$W(P_{t}f)$.
1\. Uniform-in-$N$ log-Sobolev inequality for $\mu_{N}(\xi,\eta)$
In Theorem 3.3, the log-Sobolev inequality for the distribution
$\pi_{N,D}(\xi)$ holds true:

$\lambda_{1}\mathrm{\,Ent}_{\pi_{N,D}}(f)\leqslant\frac{1}{2}\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}+a}\int_{\mathbb{R}^{dN}}\frac{|\nabla_{\xi_{k}}f|^{2}}{f}\mathrm{d}\pi_{N,D},$

where the convergence rate $\lambda_{1}=\exp(-4\beta M_{1})$. For the velocity
variables $\\{\eta_{k}\\}_{k=0}^{N-1}$ in $\mathbb{R}^{dN}$, define the
Gaussian distribution $\nu_{N}(\eta)$ by its density function:

$\nu_{N}(\eta)\propto\exp\bigg(-\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k}^{2}+a)|\eta_{k}|^{2}\bigg{missing}),~{}~{}~{}~{}\eta\in\mathbb{R}^{dN},$

then log-Sobolev inequality for $\nu_{N}$ holds true,

$\mathrm{\,Ent}_{\nu_{N}}(f)\leqslant\frac{1}{2}\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}+a}\int_{\mathbb{R}^{dN}}\frac{|\nabla_{\eta_{k}}f|^{2}}{f}\mathrm{d}\nu_{N}.$

Since the distribution
$\mu_{N,D}(\xi,\eta)=\pi_{N,D}(\xi)\otimes\nu_{N}(\eta)$ is the tensor
product, Proposition 5.2.7 of [31] yields the log-Sobolev inequality for the
product distribution:

$\lambda_{1}\mathrm{Ent}_{\mu_{N,D}}(f)\leqslant\frac{1}{2}\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}+a}\int_{\mathbb{R}^{2dN}}\frac{|\nabla_{\xi_{k}}f|^{2}+|\nabla_{\eta_{k}}f|^{2}}{f}\mathrm{d}\mu_{N,D},$
(25)

where the convergence rate is determined by the smaller one of the rates
$\lambda_{1}=\exp(-4\beta M_{1})$ and $1$, which is $\lambda_{1}$ itself. Note
that (25) does not imply the ergodicity of (17) directly, because (17) has
degenerate diffusion in the $\eta$ variable.
2\. Calculation of the generalized $\Gamma$ operators for $\Phi_{1}(f)$ and
$\Phi_{2}(f)$
We still use the notation

$\mathcal{V}_{N,D}^{a}(\xi)=\beta_{D}\sum_{j=0}^{D-1}V^{a}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{D})\bigg{)},~{}~{}~{}~{}\xi\in\mathbb{R}^{dN},$

and the generator of (17) is given by

$L=\sum_{k=0}^{N-1}\eta_{k}\cdot\nabla_{\xi_{k}}-\sum_{k=0}^{N-1}\bigg{(}\xi_{k}+\eta_{k}+\frac{1}{\omega_{k}^{2}+a}\nabla_{\xi_{k}}\mathcal{V}_{N,D}^{a}(\xi)\bigg{)}\cdot\nabla_{\eta_{k}}+\sum_{k=0}^{N-1}\frac{\Delta_{\xi_{k}}}{\omega_{k}^{2}+a}.$

By direct calculation, the commutators $[L,\nabla_{\xi_{k}}]$ and
$[L,\nabla_{\eta_{k}}]$ are given by

$[L,\nabla_{\xi_{k}}]=\nabla_{\eta_{k}}+\sum_{l=0}^{N-1}\frac{1}{\omega_{l}^{2}+a}\nabla_{\xi_{k}\xi_{l}}^{2}\mathcal{V}_{N,D}^{a}(\xi)\cdot\nabla_{\eta_{l}},~{}~{}~{}~{}[L,\nabla_{\eta_{k}}]=\nabla_{\eta_{k}}-\nabla_{\xi_{k}}.$

Inspired from Example 3 of [33], define the functions

$\Phi_{1}(f)=f\log
f,~{}~{}~{}\Phi_{2}(f)=\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}+a}\frac{|\nabla_{\eta_{k}}f-\nabla_{\xi_{k}}f|^{2}+|\nabla_{\eta_{k}}f|^{2}}{f},$
(26)

where $\Phi_{2}(f)$ can be viewed as a twisted form of

$\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}+a}\frac{|\nabla_{\xi_{k}}f|^{2}+|\nabla_{\eta_{k}}f|^{2}}{f},$

which appears in the RHS of the log-Sobolev inequality (25). From Example C.4
in Appendix C, the generalized $\Gamma$ operator $\Gamma_{\Phi_{1}}(f)$ is
given by

$\Gamma_{\Phi_{1}}(f)=\frac{1}{2}\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}+a}\frac{|\nabla_{\eta_{k}}f|^{2}}{f}.$
(27)

To compute $\Gamma_{\Phi_{2}}(f)$, we write
$\Phi_{2}(f)=\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}+a}\Phi_{2,k}(f)$, where

$\Phi_{2,k}(f)=\frac{|\nabla_{\eta_{k}}f-\nabla_{\xi_{k}}f|^{2}+|\nabla_{\eta_{k}}f|^{2}}{f},~{}~{}~{}~{}k=0,1,\cdots,N-1.$

According to Example C.5 in Appendix C, we have

$\displaystyle f\cdot\Gamma_{\Phi_{2,k}}(f)$
$\displaystyle\geqslant(\nabla_{\eta_{k}}f-\nabla_{\xi_{k}}f)\cdot[L,\nabla_{\eta_{k}}-\nabla_{\xi_{k}}]f+\nabla_{\eta_{k}}f\cdot[L,\nabla_{\eta_{k}}]f$
$\displaystyle=|\nabla_{\eta_{k}}f-\nabla_{\xi_{k}}f|^{2}-(\nabla_{\eta_{k}}f-\nabla_{\xi_{k}}f)\cdot\sum_{l=0}^{N-1}\frac{1}{\omega_{l}^{2}+a}\nabla_{\xi_{k}\xi_{l}}^{2}\mathcal{V}_{N,D}^{a}(\xi)\cdot\nabla_{\eta_{l}}f.$

Taking the summation over $k=0,1,\cdots,N-1$, we obtain

$f\cdot\Gamma_{\Phi_{2}}(f)\geqslant\sum_{k=0}^{N-1}\frac{|\nabla_{\eta_{k}}f-\nabla_{\xi_{k}}f|^{2}}{\omega_{k}^{2}+a}-\sum_{k,l=0}^{N-1}\frac{\nabla_{\eta_{k}}f-\nabla_{\xi_{k}}f}{\omega_{k}^{2}+a}\cdot\nabla_{\xi_{k}\xi_{l}}^{2}\mathcal{V}_{N,D}^{a}(\xi)\cdot\frac{\nabla_{\eta_{l}}f}{\omega_{l}^{2}+a}.$
(28)

To further simplify the expression of $\Gamma_{\Phi_{2}}(f)$, define the
vectors $X,Y\in\mathbb{R}^{dN}$ by

$X=\bigg{\\{}\frac{\nabla_{\eta_{k}}f-\nabla_{\xi_{k}}f}{\sqrt{\omega_{k}^{2}+a}}\bigg{\\}}_{k=0}^{N-1}\in\mathbb{R}^{dN},~{}~{}~{}~{}Y=\bigg{\\{}\frac{\nabla_{\eta_{k}}f}{\sqrt{\omega_{k}^{2}+a}}\bigg{\\}}_{k=0}^{N-1}\in\mathbb{R}^{dN},$
(29)

then the inequality (28) can be equivalently written as

$\Gamma_{\Phi_{2}}(f)\geqslant\frac{|X|^{2}-X^{\mathrm{T}}\Sigma Y}{f},$ (30)

where the symmetric matrix $\Sigma\in\mathbb{R}^{dN\times dN}$ is given by

$\Sigma_{kl}=\frac{1}{\sqrt{(\omega_{k}^{2}+a)(\omega_{l}^{2}+a)}}\nabla_{\xi_{k}\xi_{l}}^{2}\mathcal{V}_{N,D}^{a}(\xi)\in\mathbb{R}^{d\times
d},~{}~{}~{}~{}k,l=0,1,\cdots,N-1.$

By Lemma A.3 we have
$-\frac{M_{2}}{a}I_{dN}\preccurlyeq\Sigma\preccurlyeq\frac{M_{2}}{a}I_{dN}$,
hence (30) directly produces

$\Gamma_{\Phi_{2}}(f)\geqslant\frac{|X|^{2}-\frac{M_{2}}{a}|X||Y|}{f}.$ (31)

In conclusion, the functions $\Phi_{1}(f)$, $\Phi_{2}(f)$ and their
generalized $\Gamma$ operators satisfy

$\Phi_{1}(f)=f\log
f,~{}~{}~{}~{}\Phi_{2}(f)=\frac{|X|^{2}+|Y|^{2}}{f},~{}~{}~{}~{}\Gamma_{\Phi_{1}}(f)=\frac{|Y|^{2}}{2f},~{}~{}~{}~{}\Gamma_{\Phi_{2}}(f)\geqslant\frac{|X|^{2}-\frac{M_{2}}{a}|X||Y|}{f},$

where the vectors $X,Y\in\mathbb{R}^{dN}$ are given in (29).
3\. Generalized curvature-dimension condition produces ergodicity
Let us summarize the functional inequalities. The log-Sobolev inequality (25)
implies

$\lambda_{1}\mathrm{\,Ent}_{\mu_{N,D}}(f)\leqslant\frac{1}{2}\int_{\mathbb{R}^{2dN}}\frac{|X+Y|^{2}+|Y|^{2}}{f}\mathrm{d}\mu_{N,D}\leqslant\frac{3}{2}\int_{\mathbb{R}^{2dN}}\frac{|X|^{2}+|Y|^{2}}{f}\mathrm{d}\mu_{N,D},$

hence with the expressions of $\Phi_{1}(f)$ and $\Phi_{2}(f)$, we can
equivalently write

$\frac{2\lambda_{1}}{3}\bigg{(}\int_{\mathbb{R}^{2dN}}\Phi_{1}(f)\mathrm{d}\mu_{N,D}-\Phi_{1}\bigg{(}\int_{\mathbb{R}^{2dN}}f\mathrm{d}\mu_{N,D}\bigg{)}\bigg{)}\leqslant\int_{\mathbb{R}^{2dN}}\Phi_{2}(f)\mathrm{d}\mu_{N,D}.$
(32)

The inequality (31) implies

$\Gamma_{\Phi_{2}}(f)-\frac{1}{2}\Phi_{2}(f)\geqslant\frac{|X|^{2}-2\frac{M_{2}}{a}|X||Y|-|Y|^{2}}{2f},$

and thus using the expression of $\Gamma_{\Phi_{1}}(f)$ we have

$\Gamma_{\Phi_{2}}(f)-\frac{1}{2}\Phi_{2}(f)+\bigg{(}\frac{M_{2}^{2}}{a^{2}}+1\bigg{)}\Gamma_{\Phi_{1}}(f)\geqslant\frac{(|X|-\frac{M_{2}}{a}|Y|)^{2}}{2f}\geqslant
0.$ (33)

Utilizing Theorem C.9 in Appendix C, define the entropy-like quantity by

$W_{\mu_{N,D}}(f)=\bigg{(}\frac{M_{2}^{2}}{a^{2}}+1\bigg{)}\bigg{(}\int_{\mathbb{R}^{2dN}}\Phi_{1}(f)\mathrm{d}\mu_{N,D}-\Phi_{1}\bigg{(}\int_{\mathbb{R}^{2dN}}f\mathrm{d}\mu_{N,D}\bigg{)}\bigg{)}+\int_{\mathbb{R}^{2dN}}\Phi_{2}(f)\mathrm{d}\mu_{N,D}.$

From the functional inequalities (32) and (33), we derive the exponential
decay

$W_{\mu_{N,D}}(P_{t}f)\leqslant\exp\bigg(-\frac{t}{1+\frac{3(M_{2}^{2}/a^{2}+1)}{2\lambda_{1}}}\bigg{missing})W_{\mu_{N,D}}(f)\leqslant\exp(-2\lambda_{2}t)W_{\mu_{N,D}}(f),~{}~{}~{}~{}\forall
t\geqslant 0,$

which completes the proof.

The convergence rate $\lambda_{2}$ does not depend on the number of modes $N$
or the discretization size $D$. Nevertheless, Theorem 3.5 requires the special
Assumption (iii), i.e., $N\leqslant D$.

###  3.4 Extension of the the standard PIMD

In Theorem 3.5, we have proved the uniform-in-$N$ ergodicity of the
underdamped Matsubara mode PIMD (17) under Assumptions (i)(ii)(iii). In
particular, Assumption (iii) is satisfied when the discretization size $D=N$
or $D=\infty$. In these two cases, the Boltzmann distribution of the
continuous loop is given by

$\displaystyle(D=N)~{}:~{}\mu_{N,N}(\xi)$
$\displaystyle\propto\exp\\{-\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k}^{2}+a)(|\xi_{k}|^{2}+|\eta_{k}|^{2})-\beta_{N}\sum_{j=0}^{N-1}V^{a}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{D})\bigg{)}\bigg{\\}},$
$\displaystyle(D=\infty)~{}:~{}~{}~{}~{}\mu_{N}(\xi)$
$\displaystyle\propto\exp\\{-\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k}^{2}+a)(|\xi_{k}|^{2}+|\eta_{k}|^{2})-\int_{0}^{\beta}V^{a}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(\tau)\bigg{)}\mathrm{d}\tau\bigg{\\}}.$

Theorem 3.5 directly shows the Matsubara mode PIMD for sampling
$\pi_{N,N}(\xi)$ and $\pi_{N}(\xi)$ has the uniform-in-$N$ convergence rate
$\lambda_{2}=\frac{a^{2}}{3M_{2}^{2}+5a^{2}}\exp(-4\beta M_{1})$.

Given the ergodicity results of the Matsubara mode PIMD, we can conveniently
extend the uniform-in-$N$ ergodicity to the standard PIMD. Let $N$ be an odd
integer, we employ the $N$ normal mode coordinates $\\{\xi_{k}\\}_{k=0}^{N-1}$
to represent the $N$ beads of the ring polymer (see Appendix B). The Boltzmann
distribution of the $N$ beads is given in (15):

$\mu_{N,N}^{\mathrm{std}}(\xi,\eta)\propto\exp\\{-\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k,N}^{2}+a)(|\xi_{k}|^{2}+|\eta_{k}|^{2})-\beta_{N}\sum_{j=0}^{N-1}V^{a}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{N})\bigg{)}\bigg{\\}},$

and the standard PIMD for sampling $\mu_{N,N}^{\mathrm{std}}(\xi,\eta)$ is
given by

$\left\\{\begin{aligned} \dot{\xi}_{k}&=\eta_{k},\\\
\dot{\eta}_{k}&=-\xi_{k}-\frac{\beta_{N}}{\omega_{k,N}^{2}+a}\sum_{j=0}^{N-1}\nabla
V^{a}(x_{N}(j\beta_{N}))c_{k}(j\beta_{N})-\eta_{k}+\sqrt{\frac{2}{\omega_{k,N}^{2}+a}}\dot{B}_{k},\end{aligned}\right.$
(34)

where we choose the damping rate $\gamma=1$. Similar to the quantity (21) in
the Matsubara mode PIMD, we define the entropy-like quantity of $f(\xi,\eta)$
in $\mathbb{R}^{2dN}$:

$W_{\mu_{N,N}^{\mathrm{std}}}(f):=\bigg{(}\frac{M_{2}^{2}}{a^{2}}+1\bigg{)}\mathrm{\,Ent}_{\mu_{N,N}^{\mathrm{std}}}(f)+\sum_{k=0}^{N-1}\frac{1}{\omega_{k,N}^{2}+a}\int_{\mathbb{R}^{2dN}}\frac{|\nabla_{\eta_{k}}f-\nabla_{\xi_{k}}f|^{2}+|\nabla_{\eta_{k}}f|^{2}}{f}\mathrm{d}\mu_{N,N}^{\mathrm{std}}.$

Utilizing the same approaches with Theorem 3.5, we can prove the uniform-
in-$N$ ergodicity of the standard PIMD (34).

######  Theorem 3.7.

Let $N$ be an odd integer. Under Assumptions (i)(ii), let $(P_{t})_{t\geqslant
0}$ be the Markov semigroup of the standard PIMD (34), then for any positive
smooth function $f(\xi,\eta)$ in $\mathbb{R}^{2dN}$,

$W_{\mu_{N,N}^{\mathrm{std}}}(P_{t}f)\leqslant\exp(-2\lambda_{2}t)W_{\mu_{N,N}^{\mathrm{std}}}(f),~{}~{}~{}~{}\forall
t\geqslant 0,$

where the convergence rate
$\lambda_{2}=\frac{a^{2}}{3M_{2}^{2}+5a^{2}}\exp(-4\beta M_{1})$.

From Theorem 3.7, we conclude that the standard PIMD has uniform-in-$N$
ergodicity.

######  Remark 3.8.

The uniform-in-$N$ ergodicity of the standard PIMD (34) should also hold true
when $N$ is an even integer, and it requires a slight modification on the
normal mode transform (38).

We note that Theorem 3.5 of [36] also provides a result of the dimension-free
ergodicity of the PIMD, and the major difference is that [36] studies the
preconditioned Hamiltonian Monte Carlo (pHMC) rather than the preconditioned
Langevin dynamics, which is a more popular thermostat used in the PIMD [37].
Also, it is required in (3.10) that the duration parameter in the pHMC needs
to be small enough.

##  4 Numerical tests

###  4.1 Setup of parameters and the time average estimator

For several examples of quantum systems, we compute the quantum thermal
average employing the Matsubara mode PIMD (17) and the standard PIMD (34). For
simplicity, we choose the preconditioning parameter $a=1$ and the damping rate
$\gamma=1$. Also, the discretization size $D$ equal to the number of modes $N$
is an odd integer.

In the Matsubara mode PIMD, the quantum thermal average
$\langle{O(\hat{q})}\rangle_{\beta}$ is approximated by th statistical average
$\langle{O(\hat{q})}\rangle_{\beta,N,N}$:

$\langle{O(\hat{q})}\rangle_{\beta}\approx\langle{O(\hat{q})}\rangle_{\beta,N,N}:=\int_{\mathbb{R}^{dN}}\bigg{[}\frac{1}{N}\sum_{j=0}^{N-1}O\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{N})\bigg{)}\bigg{]}\pi_{N,N}(\xi)\mathrm{d}\xi,$

where $\pi_{N,N}(\xi)$ is the Boltzmann distribution of the continuous loop
defined in (14). Similarly, in the standard PIMD we have the approximation

$\langle{O(\hat{q})}\rangle_{\beta}\approx\langle{O(\hat{q})}\rangle_{\beta,N,N}^{\mathrm{std}}:=\int_{\mathbb{R}^{dN}}\bigg{[}\frac{1}{N}\sum_{j=0}^{N-1}O\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{N})\bigg{)}\bigg{]}\pi_{N,N}^{\mathrm{std}}(\xi)\mathrm{d}\xi,$

Finally, let $\xi_{k}(t)$ and $\eta_{k}(t)$ be the solution to the Matsubara
mode PIMD (17) or the standard PIMD (34), then the quantum thermal average is
computed from the time average

$A(\beta,N,T):=\frac{1}{T}\int_{0}^{T}\bigg{[}\frac{1}{N}\sum_{j=0}^{N-1}O\bigg{(}\sum_{k=0}^{N-1}\xi_{k}(t)c_{k}(j\beta_{N})\bigg{)}\bigg{]},$

where $T>0$ is the simulation time. The accuracy of the simulation can thus be
characterized by the time average error
$A(\beta,N,T)-\langle{O(\hat{q})}\rangle_{\beta}$, which depends on the number
of modes $N$ and the simulation time $T$.

For the discretization of the Langevin dynamics, we employ the standard BAOAB
integrator [40], which is widely applied in the molecular dynamics [40, 41].

###  4.2 1D model potential

Let the potential function $V(q)$ and the observable function $O(q)$ be given
by

$V(q)=\frac{1}{2}q^{2}+q\cos
q,~{}~{}~{}~{}O(q)=\sin\bigg(\frac{\pi}{2}q\bigg{missing}),~{}~{}~{}~{}q\in\mathbb{R}^{1}.$
(35)

The exact value of the quantum thermal average
$\langle{O(\hat{q})}\rangle_{\beta}$ is computed from the spectral method with
the Gauss–Hermite quadrature. In the numerical tests, fix the time step
$\Delta t=1/16$ and the simulation time $T=5\times 10^{6}$.

####  4.2.1 Convergence of the time average

We plot the time average error in computing the quantum thermal average
$\langle{O(\hat{q})}\rangle_{\beta}$ in Figure 1. The inverse temperature
$\beta=1,2,4,8$, and the number of modes $N=9,17,33,65,129$. The left and
right columns show the result of the Matsubara mode PIMD and the standard
PIMD, respectively.

Figure 1: Time average error in computing the quantum thermal average for the
1D model potential (35). Left: Matsubara mode PIMD. Right: standard PIMD. Top
to bottom: the inverse temperature $\beta=1,2,4,8$.

Figure 1 shows that the standard PIMD has a better accuracy than the Matsubara
mode PIMD in all temperatures, and requires a smaller number of modes $N$ for
convergence.

####  4.2.2 Correlation function of mode coordinates

The correlation function of the mode coordinates $\\{\xi_{k}\\}_{k=0}^{N-1}$
is computed from

$C_{k}(\beta,N,\Delta T):=\frac{\langle{\xi_{k}(t)\xi_{k}(t+\Delta
T)}\rangle}{\langle{\xi_{k}(t)\xi_{k}(t)}\rangle},~{}~{}~{}~{}k=0,1,\cdots,N-1,$

where
$\langle{f(t)}\rangle:=\lim_{T\rightarrow\infty}T^{-1}\int_{0}^{T}f(t)\mathrm{d}t$
denotes the time average of the function $f$, and $\Delta t$ is the time
interval of two instants recording the coordinate $\xi_{k}$. The exponential
decay of $C_{k}(\beta,N,\Delta)$ is an important criterion of the geometric
ergodicity.

We compute the correlation functions of the first five mode coordinates
$\\{\xi_{k}\\}_{k=0}^{4}$, where the inverse temperature $\beta=1,2,4,8$, and
the number of modes $N=9,17,33,65,129$. The correlation functions
$C_{k}(\beta,N,\Delta T)$ are plotted as the function of $\Delta T$ in Figure
2.

Figure 2: Correlation functions in the numerical simulation of the 1D
potential 35. Left: Matsubara mode PIMD. Right: standard PIMD. Top to bottom:
$\beta=1,2,4,8$. The centroid mode $\xi_{0}$ is colored in blue,
$\xi_{1},\xi_{2}$ are colored in red, and $\xi_{3},\xi_{4}$ are colored in
yellow.

The coincidence of the correlation functions for various $N$ shows that both
the Matsubara mode PIMD and the standard PIMD have uniform-in-$N$ ergodicity.
Meanwhile, the separation of the correlation functions for various $k$
verifies the convergence rates on the different modes can be divergent, and in
low temperatures ($\beta$ is large), high frequency modes tend to have a
longer correlation time.

###  4.3 3D spherical potential

Consider the 3D spherical potential

$V(q)=\frac{1}{2}|q|^{2}+\frac{1}{\sqrt{|q|^{2}+0.2^{2}}},~{}~{}~{}~{}|q|=\sqrt{q_{1}^{2}+q_{2}^{2}+q_{3}^{2}},$

and we aim to capture the probability distribution of $|q|$, namely, the
Euclidean distance from the origin in $\mathbb{R}^{3}$. Utilizing the density
operator $e^{-\beta\hat{H}}$, the distribution of $|q|$ can be expressed via
the density function

$\rho(r)=\frac{1}{\mathcal{Z}}\int_{\mathbb{R}^{3}}\matrixelement{q}{e^{-\beta\hat{H}}}{q}\delta(|q|-r)\mathrm{d}q,~{}~{}~{}~{}r\geqslant
0,~{}~{}~{}~{}\mathcal{Z}=\Tr[e^{-\beta\hat{H}}\big{]}.$

The distribution $\rho(r)$ is able to characterize the radial observable
functions. Choose inverse temperature $\beta=4$, the time step $\Delta
t=1/32$, and the simualtion time $T=5\times 10^{6}$. In Figure 3, we plot the
density of $|q|$ while the number of modes $N$ varies in $3,5,9,17,33$.

Figure 3: Probability density of $|q|$ in the simulation of the PIMD. Left:
Matsubara mode PIMD. Right: standard PIMD. The top and bottom graphs use
different scales.

Figure 3 shows that as the number of modes $N$ increases, the density function
shifts from the classical distribution (black dashed curve) to the quantum
limit. Both the Matsubara mode PIMD and the standard PIMD correctly computes
the correct density function $\rho(r)$, however the standard PIMD has better
accuracy than the Matsubara mode PIMD when the number of modes $N$ is small.

##  5 Conclusion

We prove the uniform-in-$N$ ergodicity of the Matsubara mode PIMD and the
standard PIMD for a general potential function $V(q)$, i.e., the convergence
towards the invariant distribution does not depend on the number of modes (for
the Matsubara mode PIMD) or the number of beads (for the standard PIMD).

## Acknowledgments

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).

The authors would like to thank Haitao Wang (SJTU), Xin Chen (SJTU), Wei Liu
(WHU) for the helpful discussions. The numerical tests are supported by High-
performance Computing Platform of Peking University.

##  Appendix A Additional proofs for Section 3

######  Lemma A.1.

Suppose $V^{c}(q)$ is convex in $\mathbb{R}^{d}$, then for any coordinates
$\xi=\\{\xi_{k}\\}_{k=0}^{N-1}$ in $\mathbb{R}^{dN}$,

$\mathcal{V}_{N,D}^{c}(\xi)=\beta_{D}\sum_{j=0}^{D-1}V^{c}(x_{N}(j\beta_{D}))$

is convex in $\mathbb{R}^{dN}$, where
$x_{N}(\tau)=\sum_{k=0}^{N-1}\xi_{k}c_{k}(\tau)$ is the continuous loop with
coordinates $\\{\xi_{k}\\}_{k=0}^{N-1}$.

######  Proof A.2.

The Hessian matrix of $V_{N,D}^{c}(\xi)$ is given by

$\nabla_{\xi_{k}\xi_{l}}^{2}\mathcal{V}_{N,D}^{c}(\xi)=\beta_{D}\sum_{j=0}^{D-1}\nabla^{2}V^{c}(x_{N}(j\beta_{D}))c_{k}(j\beta_{D})c_{l}(j\beta_{D})\in\mathbb{R}^{d\times
d},~{}~{}~{}~{}k,l=0,1,\cdots,N-1.$

To prove that $V_{N,D}^{c}(\xi)$ is convex, we only need to show for any
constants $\\{q_{k}\\}_{k=0}^{N-1}$,

$\mathcal{S}(q,q):=\sum_{k,j=0}^{N-1}q_{k}\cdot\nabla_{\xi_{k}\xi_{l}}^{2}\mathcal{V}_{N,D}^{c}(\xi)\cdot
q_{l}\geqslant 0.$

By direct calculation, we have

$\mathcal{S}(q,q)=\beta_{D}\sum_{j=0}^{D-1}\bigg{[}\bigg{(}\sum_{k=0}^{N-1}q_{k}c_{k}(j\beta_{D})\bigg{)}\cdot\nabla^{2}V^{c}(x_{N}(j\beta_{D}))\cdot\bigg{(}\sum_{k=0}^{N-1}q_{l}c_{l}(j\beta_{D})\bigg{)}\bigg{]}.$

Introduce the constants $v_{j}=\sum_{k=0}^{N-1}q_{k}c_{k}(j\beta_{D})$ in
$\mathbb{R}^{d}$ for $j=0,1,\cdots,D-1$, then we can write

$\mathcal{S}(q,q)=\beta_{D}\sum_{j=0}^{D-1}v_{j}\cdot\nabla^{2}V^{c}(x_{N}(j\beta_{D}))\cdot
v_{j}.$

Since each $\nabla^{2}V^{c}(x_{N}(j\beta_{D}))\in\mathbb{R}^{d\times d}$ is
positive semidefinite, we conclude $\mathcal{S}(q,q)\geqslant 0$.

######  Lemma A.3.

Let positive integers $N,D$ satisfy $N\leqslant D$, and the potential
$V^{a}(q)$ satisfies $-M_{2}I_{d}\preccurlyeq\nabla^{2}V^{a}(q)\preccurlyeq
M_{2}I_{d}$ in $\mathbb{R}^{d}$. For any mode coordinates
$\\{\xi_{k}\\}_{k=0}^{N-1}$, define the potential function

$\mathcal{V}_{N,D}^{a}(\xi)=\beta_{D}\sum_{j=0}^{D-1}V^{a}(x_{N}(j\beta_{D})),~{}~{}~{}~{}\xi\in\mathbb{R}^{dN}$

and the matrix $\Sigma\in\mathbb{R}^{dN}$ by

$\Sigma_{kl}=\frac{1}{\sqrt{(\omega_{k}^{2}+a)(\omega_{l}^{2}+a)}}\nabla_{\xi_{k}\xi_{l}}^{2}\mathcal{V}_{N,D}^{a}(\xi)\in\mathbb{R}^{d\times
d},~{}~{}~{}~{}k,l=0,1,\cdots,N-1,$

then $\Sigma$ satisfies
$-\frac{M_{2}}{a}I_{dN}\preccurlyeq\Sigma\preccurlyeq\frac{M_{2}}{a}I_{dN}$.

######  Proof A.4.

By direct calculation, the Hessian of $\mathcal{V}_{N,D}^{a}(\xi)$ is given by

$\nabla_{\xi_{k}\xi_{l}}^{2}\mathcal{V}_{N,D}^{a}(\xi)=\beta_{D}\sum_{j=0}^{D-1}\nabla^{2}V^{a}(x_{N}(j\beta_{D}))c_{k}(j\beta_{D})c_{l}(j\beta_{D})\in\mathbb{R}^{d\times
d},$

then for any constants $\\{q_{k}\\}_{k=0}^{N-1}$ in $\mathbb{R}^{d}$, we have

$\sum_{k,l=0}^{N-1}q_{k}\cdot\Sigma_{kl}\cdot
q_{l}=\beta_{D}\sum_{j=0}^{D-1}\bigg{(}\sum_{k=0}^{N-1}\frac{q_{k}}{\sqrt{\omega_{k}^{2}+a}}c_{k}(j\beta_{D})\bigg{)}\cdot\nabla^{2}V^{a}(x_{N}(j\beta_{D}))\cdot\bigg{(}\sum_{l=0}^{N-1}\frac{q_{l}}{\sqrt{\omega_{l}^{2}+a}}c_{l}(j\beta_{D})\bigg{)}.$

Next, it is convenient to introduce the vectors

$v_{j}=\sum_{k=0}^{N-1}\frac{q_{k}}{\sqrt{\omega_{k}^{2}+a}}c_{k}(j\beta_{D})\in\mathbb{R}^{d},~{}~{}~{}~{}j=0,1,\cdots,D-1,$

then we obtain the relation

$\sum_{k,l=0}^{N-1}q_{k}\cdot\Sigma_{kl}\cdot
q_{l}=\beta_{D}\sum_{j=0}^{D-1}v_{j}\cdot\nabla^{2}V^{a}(x_{N}(j\beta_{D}))\cdot
v_{j}\Longrightarrow\bigg{|}\sum_{k,l=0}^{N-1}q_{k}\cdot\Sigma_{kl}\cdot
q_{l}\bigg{|}\leqslant\beta_{D}M_{2}\sum_{j=0}^{D-1}|v_{j}|^{2}.$

To estimate the summation $\sum_{j=0}^{D-1}|v_{j}|^{2}$, we write

$\sum_{j=0}^{D-1}|v_{j}|^{2}=\sum_{k,l=0}^{N-1}\frac{q_{k}\cdot
q_{l}}{\sqrt{(\omega_{k}^{2}+a)(\omega_{l}^{2}+a)}}\sum_{j=0}^{D-1}c_{k}(j\beta_{D})c_{l}(j\beta_{D}).$

Since the integers $k,l\leqslant D$, we have the discrete orthogonal condition

$\beta_{D}\sum_{j=0}^{D-1}c_{k}(j\beta_{D})c_{l}(j\beta_{D})=\left\\{\begin{aligned}
&0~{}\mathrm{or}~{}1,&&\mbox{if $k=j$,}\\\ &0,&&\mbox{if $k\neq
j$.}\end{aligned}\right.$ (36)

Note that LHS of (36) can be 0 when $k=j=D$ is an even integer. Finally, we
obtain

$\beta_{D}\sum_{j=0}^{D-1}|v_{j}|^{2}\leqslant\frac{1}{a}\sum_{k=0}^{N-1}|q_{k}|^{2}\Longrightarrow\bigg{|}\sum_{k,l=0}^{N-1}q_{k}\cdot\Sigma_{kl}\cdot
q_{l}\bigg{|}\leqslant\frac{1}{a}\sum_{k=0}^{N-1}|q_{k}|^{2}.$ (37)

Since (37) holds true for any constants $\\{q_{k}\\}_{k=0}^{N-1}$ in
$\mathbb{R}^{d}$, we have
$-\frac{M_{2}}{a}I_{dN}\preccurlyeq\Sigma\preccurlyeq\frac{M_{2}}{a}I_{dN}$.

##  Appendix B Normal mode coordinates in path integral representation

We employ a slightly different approach from [37] to derive the normal mode
coordinates. For simplicity, we assume the number of modes $N$ be an odd
integer. Let $\\{\xi_{k}\\}_{k=0}^{N-1}$ the mode coordinates, and define the
continuous loop by

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

where $\\{c_{k}(\cdot)\\}_{k=0}^{N-1}$ are the Fourier basis functions defined
in Section 2.2. Suppose the discretization size in $[0,\beta]$ is also $N$,
then the $N$ bead positions of the loop are

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

From (36), the grid values of $c_{k}(\cdot)$ satisfy the normalization
condition

$\beta_{N}\sum_{j=0}^{N-1}c_{k}(j\beta_{N})c_{l}(j\beta_{N})=\delta_{k,l},~{}~{}~{}~{}k,l=0,1,\cdots,N-1,$

where $\delta_{k,l}$ is the Kronecker delta. Then we have the equality

$\sum_{j=0}^{N-1}|x_{j}|^{2}=\sum_{k=0}^{N-1}|\xi_{k}|^{2}\sum_{j=0}^{N-1}c_{k}^{2}(j\beta_{N})=\frac{1}{\beta_{N}}\sum_{k=0}^{N-1}|\xi_{k}|^{2}.$
(39)

The same procedure can be used to compute the spring potential

$\sum_{j=0}^{N-1}|x_{j}-x_{j+1}|^{2}=\sum_{j=0}^{N-1}x_{j}\cdot(2x_{j}-x_{j-1}-x_{j+1})={4}\sum_{k=0}^{N-1}\sin^{2}\bigg{(}\frac{\lceil\frac{k}{2}\rceil\pi}{N}\bigg{)}|\xi_{k}|^{2}.$

By defining the normal mode frequencies

$\omega_{k,N}=\frac{2}{\beta_{N}}\sin\bigg(\frac{\lceil\frac{k}{2}\rceil\pi}{N}\bigg{missing}),~{}~{}~{}~{}k=0,1,\cdots,N-1,$

we can conveniently write the spring potential as

$\frac{1}{\beta_{N}^{2}}\sum_{j=0}^{N-1}|x_{j}-x_{j+1}|^{2}=\sum_{k=0}^{N-1}\omega_{k,N}^{2}|\xi_{k}|^{2}.$
(40)

Recall that in the standard PIMD with $N$ beads, the energy of the ring
polymer is

$\mathcal{E}^{\mathrm{std}}(x)=\frac{1}{2\beta_{N}^{2}}\sum_{j=0}^{N-1}|x_{j}-x_{j+1}|^{2}+\beta_{N}\sum_{j=0}^{N-1}V(x_{j}),$

then we employ the equalities (39) and (40) to rewrite
$\mathcal{E}^{\mathrm{std}}(x)$ as

$\displaystyle\mathcal{E}^{\mathrm{std}}(x)$
$\displaystyle=\frac{1}{2\beta_{N}^{2}}\sum_{j=0}^{N-1}|x_{j}-x_{j+1}|^{2}+\frac{a\beta_{N}}{2}\sum_{j=0}^{N-1}|x_{j}|^{2}+\beta_{N}\sum_{j=0}^{N-1}V^{a}(x_{j})$
$\displaystyle=\frac{1}{2}\sum_{j=0}^{N-1}(\omega_{k,N}^{2}+a)|\xi_{k}|^{2}+\beta_{N}\sum_{j=0}^{N-1}V^{a}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{N})\bigg{)}.$

Therefore, the classical Boltzmann distribution in the standard PIMD is given
by

$\exp\\{-\frac{1}{2}\sum_{k=0}^{N-1}(\omega_{k,N}^{2}+a)|\xi_{k}|^{2}-\beta_{N}\sum_{j=0}^{N-1}V^{a}\bigg{(}\sum_{k=0}^{N-1}\xi_{k}c_{k}(j\beta_{N})\bigg{)}\bigg{\\}}.$

##  Appendix C Generalized $\Gamma$ calculus

We review the generalized $\Gamma$ calculus developed in [32, 33], which deals
with the ergodicity of the Markov processes with degenerate diffusions, for
example, the underdamped Langevin dynamics. The generalized $\Gamma$ calculus
is based on the Bakry–Émery theory [31].

Let $\\{X_{t}\\}_{t\geqslant 0}$ be a reversible Markov process in
$\mathbb{R}^{d}$, and $(P_{t})_{t\geqslant 0}$ be the Markov semigroup. Let
$L$ be the infinitesimal generator of $\\{X_{t}\\}_{t\geqslant 0}$, and $\pi$
be the invariant distribution. Then $L$ is self-adjoint in $L^{2}(\pi)$. In
the classical Bakry–Émery theory, the carré du champ operator
$\Gamma_{1}(f,g)$ and the iterated operator $\Gamma_{2}(f,g)$ are defined by

$\displaystyle\Gamma_{1}(f,g)$ $\displaystyle=\frac{1}{2}(L(fg)-gLf-fLg),$
$\displaystyle\Gamma_{2}(f,g)$
$\displaystyle=\frac{1}{2}(L\Gamma_{1}(f,g)-\Gamma_{1}(f,Lg)-\Gamma_{1}(g,Lf)).$

The curvature-dimension condition $CD(\rho,\infty)$ is known as the function
inequality

$\Gamma_{2}(f,f)\geqslant\rho\Gamma_{1}(f,f),~{}~{}~{}~{}\mbox{for any smooth
function $f$}.$

For a given positive smooth function $f$ in $\mathbb{R}^{d}$, define the
relative entropy of $f$ with respect to the invariant distribution $\pi$ by

$\mathrm{\,Ent}_{\pi}(f)=\int_{\mathbb{R}^{d}}f\log
f\mathrm{d}\pi-\int_{\mathbb{R}^{d}}f\mathrm{d}\pi\log\int_{\mathbb{R}^{d}}f\mathrm{d}\pi.$

If $\rho>0$, then $CD(\rho,\infty)$ implies the log-Sobolev inequality (see
Equation (5.7.1) of [31])

$\mathrm{\,Ent}_{\pi}(f)\leqslant\frac{1}{2\rho}\int_{\mathbb{R}^{d}}\frac{\Gamma_{1}(f,f)}{f}\mathrm{d}\pi,~{}~{}~{}~{}\mbox{for
any positive smooth function $f$}.$ (41)

and thus the exponential decay of the relative entropy (see Theorem 5.2.1 of
[31])

$\mathrm{\,Ent}_{\pi}(P_{t}f)\leqslant e^{-2\rho
t}\mathrm{\,Ent}_{\pi}(f),~{}~{}~{}~{}\forall t\geqslant 0.$ (42)

Inspired from the operators $\Gamma_{1}$ and $\Gamma_{2}$, we define the
generalized $\Gamma$ operator as follows.

######  Definition C.1.

Suppose $f$ is a smooth function, and $\Phi(f)$ is a function of $f$ and the
derivatives of $f$. For a stochastic process $\\{X_{t}\\}_{t\geqslant 0}$ with
generator $L$, define the generalized $\Gamma$ operator by

$\Gamma_{\Phi}(f)=\frac{1}{2}(L\Phi(f)-\mathrm{d}\Phi(f)\cdot Lf),$

where $\mathrm{d}\Phi(f)\cdot g$ for two smooth functions $f,g$ are given by

$\mathrm{d}\Phi(f)\cdot g=\lim_{s\rightarrow 0}\frac{\Phi(f+sg)-\Phi(f)}{s}.$

The expression of $\Gamma_{\Phi}(f)$ can be obtained via the following result
(Lemma 5 of [33]).

######  Lemma C.2.

Suppose $C_{1},C_{2}$ are two linear operators and $\Phi(f)=C_{1}f\cdot
C_{2}f$, then

$\Gamma_{\Phi}(f)=\Gamma_{1}(C_{1}f,C_{2}f)+\frac{1}{2}C_{1}f\cdot[L,C_{2}]f+\frac{1}{2}[L,C_{1}]f\cdot
C_{2}f,$

where $\Gamma_{1}(\cdot,\cdot)$ is the classical carré du champ operator.

In the following we assume the stochastic process $\\{X_{t}\\}_{t\geqslant 0}$
in $\mathbb{R}^{d}$ is solved by

$\mathrm{d}X_{t}=b(X_{t})\mathrm{d}t+\sigma\mathrm{d}B_{t},~{}~{}~{}~{}t\geqslant
0,$

where $b:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is the drift force,
$\sigma\in\mathbb{R}^{d\times m}$ is a constant matrix, and
$\\{B_{t}\\}_{t\geqslant 0}$ is the standard Brownian motion in
$\mathbb{R}^{m}$. Then the generator of $\\{X_{t}\\}_{t\geqslant 0}$ is given
by

$Lf(x)=b(x)\cdot f(x)+\nabla\cdot(D\nabla f),$ (43)

where $D=\frac{1}{2}\sigma\sigma^{\mathrm{T}}\in\mathbb{R}^{d\times d}$ is the
constant diffusion matrix. For the generator $L$ given in (43), we calculate
the $\Gamma$ operator $\Gamma_{\Phi}(f)$ with some classical choices of
$\Phi$.

######  Example C.3.

If $\Phi(f)=|f|^{2}$, then we can take $C_{1}=C_{2}=1$ in Lemma C.2 and obtain

$\Gamma_{\Phi}(f)=\Gamma_{1}(f,f)=\frac{1}{2}\big{[}\nabla(D\nabla(f^{2}))-2f\nabla(D\nabla
f)\big{]}=(\nabla f)^{\mathrm{T}}D\nabla f.$

In particular, since $D\in\mathbb{R}^{d\times d}$ is positive semidefinite, we
always have $\Gamma_{1}(f)\geqslant 0$.

######  Example C.4.

If $\Phi(f)=f\log f$, then by direct calculation we have

$\Gamma_{\Phi}(f)=\frac{1}{2}\Big{[}\nabla\cdot\big{(}D(\log f+1)\nabla
f\big{)}-(\log f+1)Lf\Big{]}=\frac{(\nabla f)^{\mathrm{T}}D\nabla f}{2f}.$

######  Example C.5.

If $\Phi(f)=|Cf|^{2}/f$ for some linear operator $C$, then

$\Gamma_{\Phi}(f)\geqslant\frac{Cf\cdot[L,Cf]}{f}.$ (44)

The proof below is given in Lemma 7 of [33].

######  Proof C.6.

It is easy to verify for any smooth functions $f,g$, there holds

$L(fg)=gLf+fLg+2\Gamma_{1}(f,g).$

By replacing $f\rightarrow|Cf|^{2}$ and $g\rightarrow 1/f$, we have

$\displaystyle L\bigg{(}\frac{|Cf|^{2}}{f}\bigg{)}$
$\displaystyle=\frac{1}{f}L(|Cf|^{2})+|Cf|^{2}L\bigg{(}\frac{1}{f}\bigg{)}+2\Gamma_{1}\bigg{(}|Cf|^{2},\frac{1}{f}\bigg{)}$
$\displaystyle=\frac{1}{f}L(|Cf|^{2})+|Cf|^{2}\bigg{(}-\frac{Lf}{f^{2}}+\frac{2}{f^{3}}\Gamma_{1}(f)\bigg{)}+\frac{4Cf\cdot\Gamma_{1}(Cf,f)}{f^{2}}.$
(45)

Note that

$\mathrm{d}\bigg{(}\frac{|Cf|^{2}}{f}\bigg{)}\cdot|Cf|^{2}=\frac{\mathrm{d}(|Cf|^{2})\cdot
Lf}{f^{2}}-|Cf|^{2}\frac{Lf}{f^{2}}.$ (46)

Hence from (45)(46) and the definition of the generalized $\Gamma$ operator in
(C.1), we have

$\Gamma_{\Phi}(f)\geqslant\frac{\Gamma_{|C\cdot|^{2}}(f)}{f}+\frac{1}{f^{3}}|Cf|^{2}\Gamma_{1}(f)+\frac{2Cf\cdot\Gamma_{1}(Cf,f)}{f^{2}},$
(47)

where $\Gamma_{|C\cdot|^{2}}$ is the generalized $\Gamma$ operator induced by
the function $|Cf|^{2}$. Since the matrix $D\in\mathbb{R}^{d\times d}$ is
positive semidefinite, we have the Cauchy-Swarchz inequality

$|\Gamma_{1}(f,Cf)|^{2}\leqslant\Gamma_{1}(f,f)\Gamma_{1}(Cf,Cf),~{}~{}~{}~{}\mbox{for
any smooth function $f$}.$ (48)

From (47) and (48) we obtain

$\displaystyle\Gamma_{\Phi}(f)$
$\displaystyle\geqslant\frac{\Gamma_{1}(Cf,Cf)+Cf\cdot[L,C]f}{f}+\frac{|Cf|^{2}\Gamma_{1}(f,f)}{f^{3}}$
$\displaystyle~{}~{}~{}~{}-2\sqrt{\frac{|Cf|^{2}\Gamma_{1}(f,f)}{f^{3}}}\sqrt{\frac{\Gamma_{1}(Cf,Cf)}{f}}\geqslant\frac{Cf\cdot[L,C]f}{f}.$

Hence we obtain the desired result.

Now we establish the curvature-dimension conditions for the generalized
$\Gamma$ operators. The following result relates the time derivative of
$\Phi(f)$ with $\Gamma_{\Phi}(f)$.

######  Lemma C.7.

Given the constant $t>0$, for any $s\in[0,t]$, we have the equality

$\frac{\mathrm{d}}{\mathrm{d}s}\big{[}P_{s}\Phi(P_{t-s}f)(x)\big{]}=2P_{s}\Gamma_{\Phi}(P_{t-s}f)(x).$

As a consequence, for any $t\geqslant 0$,

$\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathbb{R}^{d}}\Phi(P_{t}f)\mathrm{d}\pi=-2\int_{\mathbb{R}^{d}}\Gamma_{\Phi}(P_{t}f)\mathrm{d}\pi.$

######  Proof C.8.

Using the chain rule, we have

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}s}\big{[}P_{s}\Phi(P_{t-s}f)\big{]}$
$\displaystyle=LP_{s}\Phi(P_{t-s}f)+P_{s}\frac{\mathrm{d}}{\mathrm{d}s}\big{[}\Phi(P_{t-s}f)\big{]}$
$\displaystyle=LP_{s}\Phi(P_{t-s}f)+P_{s}\lim_{r\rightarrow
0}\frac{\Phi(P_{t-s-r}f)-\Phi(P_{t-s}f)}{r}$
$\displaystyle=LP_{s}\Phi(P_{t-s}f)-P_{s}\mathrm{d}\Phi(P_{t-s}f)\cdot
LP_{t-s}f$
$\displaystyle=P_{s}\big{(}L\Phi_{t-s}f-\mathrm{d}\Phi(P_{t-s}f)\cdot
LP_{t-s}f\big{)}$ $\displaystyle=2P_{s}\Gamma_{\Phi}(P_{t-s}f).$

Integrating the equality over the distribution $\pi$, we obtain

$\frac{\mathrm{d}}{\mathrm{d}s}\int_{\mathbb{R}^{d}}P_{s}\Phi(P_{t-s}f)\mathrm{d}\pi=2\int_{\mathbb{R}^{d}}\Gamma_{\Phi}(P_{t-s}f)\mathrm{d}\pi.$
(49)

Since $\pi$ is the invariant distribution, replacing $t-s$ by $s$ in (49), we
obtain

$\frac{\mathrm{d}}{\mathrm{d}s}\int_{\mathbb{R}^{d}}\Phi(P_{t}f)\mathrm{d}\pi=-2\int_{\mathbb{R}^{d}}\Gamma_{\Phi}(P_{t}f)\mathrm{d}\pi,$

which completes the proof.

By choosing $\Phi(f)=f\log f$, Lemma C.7 implies

$\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{\,Ent}_{\pi}(P_{t}f)=\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathbb{R}^{d}}\Phi(P_{t}f)\mathrm{d}t=-\int_{\mathbb{R}^{d}}\frac{(\nabla
f)^{\mathrm{T}}D\nabla
f}{f}\mathrm{d}\pi=-\int_{\mathbb{R}^{d}}\frac{\Gamma_{1}(f,f)}{f}\mathrm{d}\pi.$
(50)

If we have the log-Sobolev inequality (41), from (50) we have

$\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{\,Ent}_{\pi}(P_{t}f)\leqslant-2\rho\mathrm{\,Ent}_{\pi}(f),$

which implies $\mathrm{\,Ent}_{\pi}(P_{t}f)\leqslant e^{-2\rho
t}\mathrm{\,Ent}_{\pi}(f)$, and we recover the result in (42).

Now we state the main theorem, which provides the generalized curvature-
dimension condition for degenerate diffusion processes.

######  Theorem C.9.

Let $\\{X_{t}\\}_{t\geqslant 0}$ be an ergodic stochastic process with the
invariant distribution $\pi$. If for two functions $\Phi_{1}(f)$ and
$\Phi_{2}(f)$, there hold the functional inequalities

$0\leqslant\int_{\mathbb{R}^{d}}\Phi_{1}(f)\mathrm{d}\pi-\Phi_{1}\bigg{(}\int_{\mathbb{R}^{d}}f\mathrm{d}\pi\bigg{)}\leqslant
c\int_{\mathbb{R}^{d}}\Phi_{2}(f)\mathrm{d}\pi,$ (51)
$\Gamma_{\Phi_{2}}(f)\geqslant\rho\Phi_{2}(f)-m\Gamma_{\Phi_{1}}(f),$ (52)

for some constants $c,\rho,m>0$, then by defining the entropy-like quantity

$W_{\pi}(f)=m\bigg{(}\int_{\mathbb{R}^{d}}\Phi_{1}(f)\mathrm{d}\pi-\Phi_{1}\bigg{(}\int_{\mathbb{R}^{d}}f\mathrm{d}\pi\bigg{)}\bigg{)}+\int_{\mathbb{R}^{d}}\Phi_{2}(f)\mathrm{d}\pi,$

we have the exponential decay

$W_{\pi}(P_{t}f)\leqslant\exp\bigg(-\frac{2\rho
t}{1+mc}\bigg{missing})W_{\pi}(f),~{}~{}~{}~{}\forall t\geqslant 0.$

The proof below is given in Lemma 3 of [33].

######  Proof C.10.

Using Lemma C.7 and (52), we have

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}W_{\pi}(P_{t}f)$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}t}\bigg{[}m\int_{\mathbb{R}^{d}}\Phi_{1}(P_{t}f)\mathrm{d}\pi+\int_{\mathbb{R}^{d}}\Phi_{2}(P_{t}f)\mathrm{d}\pi\bigg{]}$
$\displaystyle=-2\int_{\mathbb{R}^{d}}\big{(}m\Gamma_{\Phi_{1}}+\Gamma_{\Phi_{2}}\big{)}(P_{t}f)\mathrm{d}\pi\leqslant-2\rho\int_{\mathbb{R}^{d}}\Phi_{2}(P_{t}f)\mathrm{d}\pi.$
(53)

Using (51) and the definition of $W(f)$, we have

$W_{\pi}(f)=m\bigg{(}\int_{\mathbb{R}^{d}}\Phi_{1}(f)\mathrm{d}\pi-\Phi_{1}\bigg{(}\int_{\mathbb{R}^{d}}f\mathrm{d}\pi\bigg{)}\bigg{)}+\int_{\mathbb{R}^{d}}\Phi_{2}(f)\mathrm{d}\pi\leqslant(1+mc)\int_{\mathbb{R}^{d}}\Phi_{2}(f)\mathrm{d}\pi.$

Hence (53) implies

$\frac{\mathrm{d}}{\mathrm{d}t}W_{\pi}(P_{t}f)\leqslant-\frac{2\rho}{1+mc}W_{\pi}(P_{t}f),~{}~{}~{}~{}\forall
t\geqslant 0,$

yielding the desired result.

## References

  * [1] Kerson Huang.  Statistical mechanics.  John Wiley & Sons, 2008.
  * [2] Donald A. McQuarrie.  Statistical mechanics.  Sterling Publishing Company, 2000.
  * [3] Peter W. Atkins, Julio De Paula, and James Keeler.  Atkins’ physical chemistry.  Oxford university press, 2023.
  * [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] Steven A. Orszag.  Comparison of pseudospectral and spectral approximation.  Studies in Applied Mathematics, 51(3):253–259, 1972.
  * [7] Jie Shen, Tao Tang, and Li-Lian Wang.  Spectral methods: algorithms, analysis and applications.  Springer, 41, 2011.
  * [8] Richard P. Feynman, Albert R. Hibbs, and Daniel F. Styer.  Quantum mechanics and path integrals.  Courier Corporation, 2010.
  * [9] William H. Miller.  Path integral representation of the reaction rate constant in quantum mechanical transition state theory.  The Journal of Chemical Physics, 63(3):1166–1172, 1975.
  * [10] Ian R. Craig and David E. Manolopoulos.  Chemical reaction rates from ring polymer molecular dynamics.  The Journal of chemical physics, 122(8):084106, 2005.
  * [11] Xuecheng Tao, Philip Shushkov, and Thomas F. Miller.  Microcanonical rates from ring-polymer molecular dynamics: Direct-shooting, stationary-phase, and maximum-entropy approaches.  The Journal of Chemical Physics, 152(12), 2020.
  * [12] Gregory A. Voth.  Feynman path integral formulation of quantum mechanical transition-state theory.  The Journal of Physical Chemistry, 97(32):8365–8377, 1993.
  * [13] Edit Mátyus, David J. Wales, and Stuart C. Althorpe.  Quantum tunneling splittings from path-integral molecular dynamics.  The Journal of chemical physics, 144(11):114108, 2016.
  * [14] Christophe L. Vaillant, David J. Wales, and Stuart C. Althorpe.  Tunneling splittings from path-integral molecular dynamics using a Langevin thermostat.  The Journal of chemical physics, 148(23):234102, 2018.
  * [15] Scott Habershon, David E. Manolopoulos, Thomas E. Markland, and Thomas F. Miller III.  Ring-polymer molecular dynamics: Quantum effects in chemical dynamics from classical trajectories in an extended phase space.  Annual review of physical chemistry, 64:387–413, 2013.
  * [16] Ian R. Craig and David E. Manolopoulos.  Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics.  The Journal of chemical physics, 121(8):3368–3373, 2004.
  * [17] Seogjoo Jang and Gregory A. Voth.  A derivation of centroid molecular dynamics and other approximate time evolution methods for path integral centroid variables.  The Journal of chemical physics, 111(6):2371–2384, 1999.
  * [18] Jianshu Cao and Gregory A. Voth.  The formulation of quantum statistical mechanics based on the Feynman path centroid density. III. Phase space formalism and analysis of centroid molecular dynamics.  The Journal of chemical physics, 101(7):6157–6167, 1994.
  * [19] Michael J. Willatt.  Matsubara dynamics and its practical implementation.  PhD thesis, Cambridge University, 2017.
  * [20] Timothy J. H. Hele, Michael J. Willatt, Andrea Muolo, and Stuart C. Althorpe.  Boltzmann-conserving classical dynamics in quantum time-correlation functions: “Matsubara dynamics”.  The Journal of Chemical Physics, 142(13):134103, 2015.
  * [21] Jianshu Cao and Glenn J. Martyna.  Adiabatic path integral molecular dynamics methods. II. Algorithms.  The Journal of chemical physics, 104(5):2028–2035, 1996.
  * [22] 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.
  * [23] Thomas E. Markland and David E. Manolopoulos.  An efficient ring polymer contraction scheme for imaginary time path integral simulations.  The Journal of chemical physics, 129(2):024105, 2008.
  * [24] Bruce J. Berne and D. Thirumalai.  On the simulation of quantum systems: path integral methods.  Annual Review of Physical Chemistry, 37(1):401–424, 1986.
  * [25] Jianfeng Lu and Zhennan Zhou.  Path integral molecular dynamics with surface hopping for thermal equilibrium sampling of nonadiabatic systems.  The Journal of Chemical Physics, 146(15), 2017.
  * [26] Xinzijian Liu and Jian Liu.  Path integral molecular dynamics for exact quantum statistics of multi-electronic-state systems.  The Journal of Chemical Physics, 148(10), 2018.
  * [27] 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.
  * [28] Roman Korol, Nawaf Bou-Rabee, and Thomas F. Miller.  Cayley modification for strongly stable path-integral and ring-polymer molecular dynamics.  The Journal of chemical physics, 151(12), 2019.
  * [29] Rob D. Coalson.  On the connection between Fourier coefficient and Discretized Cartesian path integration.  The Journal of chemical physics, 85(2):926–936, 1986.
  * [30] Charusita Chakravarty, Maria C. Gordillo, and David M. 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.
  * [31] Dominique Bakry, Ivan Gentil, Michel Ledoux, et al.  Analysis and geometry of Markov diffusion operators.  Springer, 103, 2014.
  * [32] Pierre Monmarché.  Hypocoercivity in metastable settings and kinetic simulated annealing.  Probability Theory and Related Fields, 172(3-4):1215–1248, 2018\.
  * [33] Pierre Monmarché.  Generalized $\Gamma$ calculus and application to interacting particles on a graph.  Potential Analysis, 50(3):439–466, 2019.
  * [34] James Glimm and Arthur Jaffe.  Quantum physics: a functional integral point of view.  Springer Science & Business Media, 2012.
  * [35] Xuda Ye and Zhennan Zhou.  Efficient sampling of thermal averages of interacting quantum particle systems with random batches.  The Journal of Chemical Physics, 154(20), 2021.
  * [36] 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.
  * [37] 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):024103, 2016.
  * [38] Yu Cao, Jianfeng Lu, and Lihan Wang.  On explicit $L^{2}$-convergence rate estimate for underdamped Langevin dynamics.  arXiv preprint arXiv:1908.04746, 2019.
  * [39] Cédric Villani.  Hypocoercivity.  Number 949-951. American Mathematical Soc., 2009.
  * [40] Ben Leimkuhler and Charles Matthews.  Molecular dynamics.  Interdisciplinary applied mathematics, 39:443, 2015.
  * [41] Marco Lauricella, Letizia Chiodo, Fabio Bonaccorso, Mihir Durve, Andrea Montessori, Adriano Tiribocchi, Alessandro Loppini, Simonetta Filippi, and Sauro Succi.  Multiscale Hybrid Modeling of Proteins in Solvent: SARS-CoV2 Spike Protein as test case for Lattice Boltzmann–All Atom Molecular Dynamics Coupling.  arXiv preprint arXiv:2305.05025, 2023.