Geometric ergodicity of Langevin dynamics with Coulomb interactions

TL;DR

This paper proves that Langevin dynamics for Coulomb gases in multi-dimensional space converge exponentially to a unique equilibrium, using a novel Lyapunov function to analyze long-term behavior.

Contribution

It introduces a new Lyapunov function approach to establish exponential ergodicity for Coulomb gas Langevin dynamics.

Findings

Exponential convergence to the Boltzmann-Gibbs measure.

Construction of a novel Lyapunov function for Coulomb systems.

Proof of ergodicity in multi-dimensional Coulomb gases.

Abstract

This paper is concerned with the long time behavior of Langevin dynamics of {\em Coulomb gases} in $\mathbf{R}^d$ with $d\geq 2$, that is a second order system of Brownian particles driven by an external force and a pairwise repulsive Coulomb force. We prove that the system converges exponentially to the unique Boltzmann-Gibbs invariant measure under a weighted total variation distance. The proof relies on a novel construction of Lyapunov function for the Coulomb system.

Equations

$$\frac{d\mathbf{q}}{dt}$$

$$\frac{d\mathbf{p}}{dt}$$

$$H(\mathbf{q},\mathbf{p}):=\frac{1}{2}\sum _{i=1}^N|p_i{|}^2+\sum _{i=1}^{N}V(q_i)+\frac{1}{2N}\sum _{1\le i\ne j\le N}K(q_i-q_j),$$

$$U(\mathbf{q}):=\sum _{i=1}^{N}V(q_i)+\frac{1}{2N}\sum _{1\le i\ne j\le N}K(q_i-q_j).$$

$$K(x)=\{\begin{array}{ll}-\log |x|& \text{ if }d=2,\\ |x{|}^{2-d}& \text{ if }d\ge 3.\end{array}$$

$$-△K=c_d{\delta }_0,\text{ with }{c}_d=\{\begin{array}{ll}2\pi & \text{ if }d=2\\ (d-2)|{\mathbb{S}}^{d-1}|=\frac{d(d-2){\pi }^{d/2}}{\Gamma (1+d/2)}& \text{ if }d\ge 3.\end{array}$$

$$\pi (d\mathbf{q}d\mathbf{p})=P(d\mathbf{p})\otimes Q(d\mathbf{q})=\frac{1}{Z}{e}^{-\beta H(\mathbf{q},\mathbf{p})}d\mathbf{q}d\mathbf{p},$$

$$\displaystyle P(d\mathbf{p})=\frac{1}{Z_{p}}\exp\Big{(}-\frac{\beta}{2}|\mathbf{p}|^{2}\Big{)}\text{ and }Q(d\mathbf{q})=\frac{1}{Z_{q}}\exp\Big{(}-\beta U(\mathbf{q})\Big{)}.$$

$$\frac{d\mathbf{q}}{dt}=-\nabla U(\mathbf{q})+\sqrt{2{\beta }^{-1}}\frac{d\mathbf{B}}{dt}$$

$$\frac{d{q}_i}{dt}=-\nabla V(q_i)-\frac{1}{N}\sum _{j\ne i}^N\nabla K(q_i-q_j)+\sqrt{2{\beta }^{-1}}\frac{d{B}_i}{dt},i=1,2,\cdots ,N.$$

$$Q(d\mathbf{q})\propto\exp\Big{(}-\beta\Big{(}\sum_{i=1}^{N}V(q_{i})+\frac{1}{2N}\sum_{1\leq i\neq j\leq N}K(q_{i}-q_{j})\Big{)}\Big{)}d\mathbf{q}.$$

$$\frac{d\mathbf{q}}{dt}$$

$$\frac{d\mathbf{p}}{dt}$$

$$U(\mathbf{q})=\sum _{i=1}^{N}V(q_i)+\frac{1}{2N}\sum _{1\le i\ne j\le N}K(q_i-q_j).$$

$$H(\mathbf{q},\mathbf{p})=\frac{1}{2}|\mathbf{p}{|}^2+\sum _{i=1}^{N}V(q_i)+\frac{1}{2N}\sum _{1\le i\ne j\le N}K(q_i-q_j).$$

$$\pi (d\mathbf{q}d\mathbf{p})=\frac{1}{Z}{e}^{-\beta H(\mathbf{q},\mathbf{p})}d\mathbf{q}d\mathbf{p}.$$

$$V(q)\ge c_1|q{|}^2-M\forall q\in {\mathbf{R}}^d,$$

$$c_{2}V(q)-M\le \nabla V(q)\cdot q\le c_{3}V(q)+M\forall q\in {\mathbf{R}}^d.$$

$$|\nabla V(q)|\le \epsilon V(q)+M_{\epsilon }\forall q\in {\mathbf{R}}^d.$$

$$Z={\iint }_{{\mathbf{R}}^d\times {\mathbf{R}}^d}{e}^{-\beta H(\mathbf{q},\mathbf{p})}d\mathbf{q}d\mathbf{p}<\infty .$$

$$\mathcal{D}:=\{\mathbf{q}\in ({\mathbf{R}}^d{)}^N:U(\mathbf{q})<\infty \}.$$

$$\mathcal{D}=\{\mathbf{q}=(q_1,q_2,\cdots ,q_N)\in ({\mathbf{R}}^d{)}^N:q_i\ne q_j\text{ if }i\ne j\text{ and }{q}_i<\infty \forall i=1,2,\cdots ,N\}.$$

$$\partial \mathcal{D}:=\{\infty \}\cup \bigcup _{1\le i\ne j\le N}\{\mathbf{q}=(q_1,q_2,\cdots ,q_N)\in ({\mathbf{R}}^d{)}^N:q_i=q_j\}.$$

$$\partial \mathcal{X}:=\{\partial \mathcal{D}\times ({\mathbf{R}}^d{)}^N\}\cup \{\mathcal{D}\times \{\infty \}\}.$$

$$(P_t\phi )(\mathbf{q},\mathbf{p}):={\mathbb{E}}_{(\mathbf{q},\mathbf{p})}[\phi (\mathbf{q}(t),\mathbf{p}(t))]$$

$$(\nu P_t)(A)={\int }_{\mathcal{X}}(P_t{1}_A)(\mathbf{q},\mathbf{p})\nu (d\mathbf{q}d\mathbf{p}),$$

$$\mathcal{L}f=\mathbf{p}\cdot {\nabla }_{\mathbf{q}}f-\gamma \mathbf{p}\cdot {\nabla }_{\mathbf{p}}f-{\nabla }_{\mathbf{q}}U(\mathbf{q})\cdot {\nabla }_{\mathbf{p}}f+\gamma {\beta }^{-1}{△}_{\mathbf{p}}f.$$

$$\mathcal{L}W\le -\lambda W+C.$$

$${\mathcal{P}}_W:=\{\nu \in {\mathcal{P}}_{\mathcal{X}}:{\int }_{\mathcal{X}}W(\mathbf{x})\nu (d\mathbf{x})<\infty \}.$$

$${\rho }_W({\nu }_1,{\nu }_2)=\underset{|\phi {|}_W\le 1}{\sup }{\int }_{\mathcal{X}}\phi (x)({\nu }_1(d\mathbf{x})-{\nu }_2(d\mathbf{x})),$$

$$|\phi {|}_W=\underset{\mathbf{x}\in \mathcal{X}}{\sup }\frac{|\phi (\mathbf{x})|}{1+W(\mathbf{x})}.$$

Full text

1. Introduction

1.1. Model description

In this paper, we study a classical Langevin dynamics for a system of $N$ particles in $\mathbf{R}^{d},d\geq 2$ subject to a confining external force and interacting through the Coulomb force. The evolution of particles is described by the following stochastic differential equations (SDEs):

[TABLE]

where $\mathbf{p}=(p_{1},\cdots,p_{N})\in(\mathbf{R}^{d})^{N}$ and $\mathbf{q}=(q_{1},\cdots,q_{N})\in(\mathbf{R}^{d})^{N}$ are vectors of positions and momentum of the $N-$particles in the Euclidean space $\mathbf{R}^{d}$ with $d\geq 2$ and $\mathbf{B}(t)$ is a vector of $N$ independent $d$-dimensional Brownian motions. Here the constant $\gamma>0$ is a friction parameter, $\beta$ is an inverse temperature parameter.

The Hamiltonian (or energy) $H(\mathbf{q},\mathbf{p})$ is defined by

[TABLE]

which is the sum of the kinetic energy $\frac{1}{2}\sum_{i=1}^{N}|p_{i}|^{2}$ and the potential energy $U(\mathbf{q})$ defined by

[TABLE]

The potential energy consists of two parts: an external energy and an interacting energy. Observe that the $N$-dependent prefactors appearing in the interaction term in the definition of the Hamiltonian $H$ keep them balanced in the same order $\mathcal{O}(N)$. The function $V:\mathbf{R}^{d}\rightarrow\mathbf{R}$ is an external confining potential, growing fast enough at infinity, and satisfying some extra conditions to be specified later. The pairwise interacting kernel $K:\mathbf{R}^{d}\rightarrow\mathbf{R}$ is the Coulomb kernel in $\mathbf{R}^{d}$:

[TABLE]

The Coulomb kernel $K$ is the fundamental solution of the Poisson equation in the free physical space $\mathbf{R}^{d}$, i.e.

[TABLE]

It is easy to check that the Boltzmann-Gibbs measure $\pi(d\mathbf{q}\,d\mathbf{p})$ is an invariant measure of the particle system (1.1):

[TABLE]

where

[TABLE]

The primary aim of this paper is to prove that the particle system (1.1) converges exponentially fast to the (unique) invariant measure $\pi$ as $t\rightarrow\infty$. This result will be a consequence of a classical Harris ergodic theorem which assumes the existence for a minorizing measure and a Lyapunov function in the spirit of [24, 25, 14]. The minorizing measure is the source of the probabilistic mixing while the Lyapunov function ensures the dynamics does not escape to infinite. The system is hypoellitpic which quickly implies the minorization condition [23]. The main contribution of this paper is to construct a Lyapunov function in the Coulomb setting. The construction used in [23], does not seem to work for the Coulomb system. The approach used here is a further refinement of the ideas developed in [15]. This is further discussed in Section 1.3.

1.2. Motivating and connected problems

The dynamics (1.1) is usually referred to as the inertial Langevin equation in the physical-chemistry literature. It models for example the time evolution of a system of confining charged particles (e.g. electrons) interacting each other through the repulsive Coulomb force in a thermal environment. When $\gamma\rightarrow\infty$, the Langevin dynamics (1.1), after proper accelerating in time, reduces to the following first order equation

[TABLE]

or more concretely

[TABLE]

The dynamics (1.5) is usually called the overdamped limit of the Langevin equation (1.1). A rigorous derivation of the overdamped limit can be found in [31]. Note that the equilibrium measure of the overdamped Langevin equation (1.5) is given by

[TABLE]

When $\beta=\beta_{N}=N\tilde{\beta}$, the equilibrium measure $Q(d\mathbf{q})$ is usually refered to as the Coulomb gases model. To clarify the terminology, we will refer (1.1) and (1.5) as the dynamical Coulomb gases, and the measure (1.7) as the steady Coulomb gases. In recent years, there has been a growing interest in the probabilistic properties of steady Coulomb gases in meso/microscopic regimes; see e.g. [28, 27, 26, 29, 21] and the references therein. However, the long time behavior of dynamical Coulomb gases was less understood. For the overdamped model (1.5) with a logarithmic $K$ in one dimension, the exponential convergence to the equilibrium can be established by using the Bakry-Émery argument thanks to the convexity of the negative logarithmic function in 1D; see e.g. [8]. The recent paper [1] studied the convergence to the equilibrium in two dimensions via the Poincaré inequality.

The Coulomb gas model (1.7) has many important connections with other models in mathematical physics. In 1D or 2D, the Coulomb gas model $Q$ is closely related to some famous random matrix models. In fact, if the potential $V$ is quadratic, e.g. $V(q)=\frac{|q|^{2}}{2\beta}$, then under the scaling $\beta=N\tilde{\beta}$, the measure $Q$ corresponds to the famous $\beta$-Hermite and $\beta$-Ginibre ensemble in 1D and 2D respectively. In particular, when $\tilde{\beta}=2$, the gas describes the eigenvalue distribution of a random $N\times N$ matrix of i.i.d. real and complex normal entries, whose limiting distribution are the Wigner’s semicircle law and the circular law respectively. Because of the link between the equilibrium measure $Q$ and the eigenvalue distribution of random matrix ensembles, the overdamped dynamics (1.6) and the under-damped dynamics (1.6) can be used for Monte Carlo simulation of random matrices. The recent paper [3] proposed an efficient Hybrid Monte Carlo (HMC) method for sampling Coulomb gases, based on a Metropolis-Hastings algorithm with a proposal produced by discretizing the second order dynamics (1.1). However, the convergences of the sampling algorithm proposed there and the associated continuous dynamics (1.1) have not been analyzed yet. This paper provides a theoretical underpinning for the HMC algorithm of [3] by showing that the (1.1) converges to the equilibrium Gibbs measure $\pi$ exponentially fast. Finally, we would like to point out a link between the dynamical Coulomb gas and partial differential equations (PDEs). At the dynamical level, one expects that as $N\rightarrow\infty$ the empirical measure of the SDEs (1.1) converges to its mean field limit which is the so-called Vlasov-Possion-Fokker-Planck equation. In the case where the kernel is Lipschitz, this picture has been well-understood. However, rigourous justifications of the mean field limits of dynamics of interacting particles with singular kernels, especially the Coulomb kernel is highly non-trivial and far from understood. We will not discuss this problem in the present paper, but we refer the interested readers to [10, 19, 30, 22, 18, 2, 17] for some recent progress.

1.3. Discussion on ergodicity for the Langevin equation

The Langevin equation is an SDE with degenerate noise. In the last few decades, both probabilistic and analytical methods have been developed for studying the long time behavior of the Langevin dynamics and the corresponding kinetic Fokker Planck equation. At the PDE level, an important analytical tool is based on the so-called hypocoercivity initiated by Villani [35] and then later advanced by Dolbeault, Mouhot and Schmeiser [6, 7]. The framework of hypercoercivity was also employed recently by [4, 11] to prove ergodicity for Langevin particle systems with some mild singular potentials excluding the Coulomb potential.

At the probabilistic level, the geometric ergodicity of the Langevin dynamics with regular drifts has been proved by [34, 23, 36] based on the techniques of Meyn and Tweedie [24, 25] and Khasminskii [20]. These techniques have been formalized now a standard route of establishing ergodicity of Markov processes: the validity of a local minorization condition and the existence of a Lyapunov function; see [14]. The local minorization condition leads to contractivity of any two coupled stochastic trajectories in some average sense if they both start from a local compact set $\mathcal{C}$. The existence of a Lyapunov function guarantees that the process return to the compact set $\mathcal{C}$ infinitely many often and usually provides quantitative estimates on the excursion time outside the set $\mathcal{C}$. In the case where the underlying diffusion is degenerate but hypoelliptic, such as the Langevin dynamics considered here, the minorization condition can be justified by a standard control-type argument [23]. However, constructing Lyapunov functions for degenerate diffusions is more technical and is problem-dependent.

For Langevin dynamics, a natural candidate of Lyapunov function would be the Hamiltonian. Yet, due to the absence of diffusion in the position-coordinate, the Hamiltonian fails to be a real Lyapunov function since it gives no dissipation in the position-coordinate. In the case where the drift of Langevin dynamics is regular, adding a $cp\cdot q$ term into the Hamiltonian with a suitable small $c>0$ enables transferring the dissipation from the momentum variable to the position and hence leads to a Lyapunov function, see [34, 23]. However, this $pq$ trick does not work when the drift of Langevin dynamics is singular. Recently, Herzog and Mattingly [15] (see also [5]) constructed a Lyapunov function for a Langevin particle system with a class of admissible singular interaction kernels, including Lennard-Jones type kernels. The key idea underlying their construction is to understand the high energy limit of the dynamics via apppriate space-time scaling. Lyapunov functions can then be built by identifying the effective/dominant generator in the high energy limit. However, the Coulomb kernel (at least in 1D/2D) can not be treated in the framework of [15, 5] since the admissible condition of [15] (see also below) is violated.

In this work, we revisit the idea of [15] and tailor an explicit Lyapunov function for the dynamical Coulomb gas (1.1) in $\mathbf{R}^{d}$ with any $d\geq 2$. In fact, the same Lyapunov function also works for dynamical Log gases [9] in one dimensions, but we restrict our attension here on $d\geq 2$ only for ease of exposition.

1.4. Plan of the paper

The rest of the paper is organized as follows. In the next section, we set up the model with some assumptions on the potential function and state the main results of the paper. In Section 3, we develop some heuristics by building a Lyapunov function for a simple two-body Coulomb system. In Section 4, we extend this heuristic idea to general many body Coulomb systems. The proofs of main results are given in Section 5. Finally, we conclude in Section 6 with several extensions and outlook.

2. Set-up and Main results

Let us recall the SDE (1.1) and rewrite it as

[TABLE]

where

[TABLE]

Recall also the Hamiltonian $H(\mathbf{q},\mathbf{p})$ defined by (1.2):

[TABLE]

We aim to show that the dynamics (2.1) converges to the Boltzmann-Gibbs measure

[TABLE]

To this end, we need to make some technical assumptions on the external potential $V$.

Assumptions 2.1**.**

The function $V:\mathbf{R}^{d}\rightarrow\mathbf{R}$ is at least twice differentiable and satisfies

(i) $V(q)\geq 0$ for all $q\in\mathbf{R}^{d}$ and there exist constants $c_{1},c_{2},c_{3},M>0$ such that

[TABLE]

and that

[TABLE]

(ii) For any $\varepsilon>0$, there exists $M_{\varepsilon}>0$ such that

[TABLE]

Remark 2.2*.*

The quadratic growth condition (2.2) implies that $V$ is strongly confining and it is necessary for the $pq$ trick to work in constructing Lyapunov functions for Langevin dynamics, see e.g. [23]. It also implies that

[TABLE]

The inequalities in (2.3) are standard dissipative conditions on $V$. The condition (2.4) controlling $\nabla V$ by $V$ itself will be useful for our explicit construction of Lyapunov function for the Coulomb dynamics. Assumption 2.1 is by no means sharp, but proves to be sufficient for the validity of our main results. In particular, the standard double well function $V(q)=\frac{1}{4}(1-|q|^{2})^{2}$ fullfils Assumption 2.1.

We denote by $\mathcal{D}$ the domain of $U$ in $(\mathbf{R}^{d})^{N}$, i.e.

[TABLE]

We also define the state space $\mathcal{X}:=\mathcal{D}\times(\mathbf{R}^{d})^{N}$. Due to the singularity of $K$ and the growth condition on $V$ at infinity, it is clear that

[TABLE]

The boundary of $\mathcal{D}$ is given by

[TABLE]

Similarly, we have the boundary of $\mathcal{X}$ given by

[TABLE]

In addition, we have the following lemma whose proof is trivial and hence omitted.

Lemma 2.3**.**

(i) The set $\mathcal{D}$ is open and path connected;

(ii) The set $\mathcal{D}_{R}:=\{\mathbf{q}\in(\mathbf{R}^{d})^{N}:U(\mathbf{q})<R\}$ is precompact in $(\mathbf{R}^{d})^{N}$ for every $R>0$.

Due to the singularity of the Coulomb kernel, the existence and uniqueness of a global pathwise solution to the SDE (2.1) does not automatically follow from a standard solution theory of SDE based on local Lipschitz continuity. In particular, it is not guaranteed a priori that two particles could not collide. The next proposition shows the global existence and uniqueness of the process $\mathbf{x}=(\mathbf{q},\mathbf{p})\in\mathcal{X}$ satisfying (2.1) and thus demonstrates that collisions between particles never occur in finite time.

Proposition 2.4**.**

Given any initial condition $\mathbf{x}(0)=(\mathbf{q}(0),\mathbf{p}(0))\in\mathcal{X}$, there exists a unique pathwise solution $\mathbf{x}(t)=(\mathbf{q}(t),\mathbf{p}(t))$ satisfying the equation (2.1). Moreover, $\mathbf{x}(t)\in\mathcal{X}$ for any finite time $t>0$ almost surely.

Proposition 2.4 can be proved by just using the Hamiltonian $H$ as the Lyapunov function, but we will give a proof in Section 5 based on a new Lyapunov function tailored for the dynamics (2.1); see Proposition 5.1.

Thanks to Proposition 2.4, one can define the Markov semigroup $P_{t}$ associated to the SDE (2.1) by

[TABLE]

where $(\mathbf{q}(t),\mathbf{p}(t))$ is the solution process of (2.1) starting from $(\mathbf{q},\mathbf{p})\in\mathcal{X}$ and $\varphi:\mathcal{X}\rightarrow\mathbf{R}$ is a bounded measurable function. Let $\mathscr{P}_{\mathcal{X}}$ be the space of probability measures on $\mathcal{X}$. We also define the dual operation of $P_{t}$ on probability measures by setting for a given probability measure $\nu\in\mathscr{P}_{X}$,

[TABLE]

where $A$ is a Boreal measurable subset of $\mathcal{X}$ and $\mathbf{1}_{A}$ is the indicator function of $A$.

The generator $\mathcal{L}$ associated to the semigroup $P_{t}$ is given for every smooth test function $f:\mathcal{X}\rightarrow\mathbf{R}$ by

[TABLE]

With the definition of the generator $\mathcal{L}$, we introduce an important concept of Lyapunov function. We say a function $W:\mathcal{X}\rightarrow(0,\infty)$ is a Lyapunov function of the equation (2.1) if $W(\mathbf{q},\mathbf{p})\rightarrow\infty$ as $H(\mathbf{q},\mathbf{p})\rightarrow\infty$ and there exist constants $\lambda,C>0$ such that

[TABLE]

The main result of the paper addresses the convergence to equilibrium for the dynamical system of Coulomb gas (2.1). Motivated by [14, 15], we define a weighted total variation distance of probability measures. To be more specific, for any measurable function $W:\mathcal{X}\rightarrow(0,\infty)$, we define the set $\mathscr{P}_{W}$ by

[TABLE]

The weighted total variation metric $\rho_{W}$ on $\mathscr{P}_{W}$ is defined for any $\nu_{1},\nu_{2}\in\mathscr{P}_{W}$ by

[TABLE]

where the semi-norm $|\varphi|_{W}$ is defined by

[TABLE]

The convergence to the equilibrium result is stated in the theorem below.

Theorem 2.5**.**

The Boltzmann-Gibbs measure $\pi(d\mathbf{q}\,d\mathbf{p})$ is the unique invariant measure of the dynamics (2.1). Moreover, for every $a\in(0,\beta)$, there exists $W\in C^{2}(\mathcal{X},(0,\infty))$ and constants $\lambda,C>0$ such that

[TABLE]

with $\Psi(\mathbf{x})=o(H(\mathbf{x}))$ as $H(\mathbf{x})\rightarrow\infty$ and such that for every $\nu\in\mathscr{P}_{W}$ and $t\geq 0$

[TABLE]

Remark 2.6*.*

The function $W$ in the above statement is in fact the Lyapunov function to be constructed below. Theorem 2.5 implies an exponential convergence to the equilibrium of the Markov semigroup $P_{t}$ for almost optimal observables with sharp integrability. Indeed, first Theorem 2.5 states that (2.8) holds with any initial condition $\nu$ satisfying

[TABLE]

Furthermore, consider $\nu=\delta_{\mathbf{x}}$ and a test function $\varphi(x):\mathcal{X}\rightarrow\mathbf{R}$ with $|\varphi|_{W}<\infty$. Then (2.8) translates into

[TABLE]

One can not expect this hold when $a\geq\beta$ since in this case there exist test functions such that $\pi\varphi=\infty$ while $\sup_{\mathbf{x}}|\varphi(\mathbf{x})|/(e^{aH(\mathbf{x})})<\infty$.

The proof of Theorem 2.5 is based on two ingredients: a local minorization condition on the probability density $P_{t}(x,\cdot)$ and the existence of a suitable Lyapunov function. The proof of the former is quite standard and follows from the support theorem [33, 32]. This relies on the fact that diffusion of consideration here is hypoelliptic so that the probability density $P_{t}(x,\cdot)$ is smooth by the Hörmander’s theorem [16].

The second ingredient, i.e. the existence of a suitable Lyapunov is established in Proposition 2.7. Its proof requires a more delicate analysis and is the major contribution of the present paper; see the next two sections. The complete proof of the Theorem 2.5 can be found in Section 5.

Proposition 2.7**.**

For any $a\in(0,\beta)$, there exists $W,\Psi\in C^{2}(\mathcal{X},(0,\infty))$ satisfying the following:

(i) $W(\mathbf{x})=e^{aH(\mathbf{x})+\Psi(\mathbf{x})}$ and $\Psi(\mathbf{x})=o(H(\mathbf{x}))$ as $H(\mathbf{x})\rightarrow\infty$. As a consequence, $W(\mathbf{x})\rightarrow\infty$ as $H(\mathbf{x})\rightarrow\infty$.

(ii) There exists $\lambda,C>0$ depending on $N$ such that

[TABLE]

3. Construction of Lyapunov function: Building intuition from a simple example

In this section, we would like to build some intuition of constructing a Lyapunov function for (2.1) through examining carefully the structure of a simple two-body Coulomb system. To be more concrete, we consider a simple system of Coulomb gases with only two particles ($N=2$ in (2.1)). To simplify the dynamics even further, we artificially freeze one of the particle at the origin by setting $q_{1}=p_{1}=0$ and only study the dynamics of $(q,p)=(q_{2},p_{2})\in\mathcal{X}=(\mathbf{R}^{d})^{2}\setminus(0,0)$ with $d\geq 2$:

[TABLE]

with $U(q)=V(q)+K(q)$. Recall the Coulomb kernel $K$ in (1.4). In this case, the generator $\mathcal{L}$ associated to above is

[TABLE]

where $\mathcal{H}$ is the conservative part (corresponding to the conservative Hamiltonian transport dynamics) and $\mathcal{O}$ is the dissipative part (corresponding to the fluctuation/dissipation of the Orstein-Uhlenbeck diffusion process on the momentum):

[TABLE]

3.1. Revisting old ideas

An initial natural guess of the Lyapunov function $W$ would simply be the Hamiltonian $H(q,p)=\frac{1}{2}|p|^{2}+V(q)+K(q)$. However, a simply calculation yields that

[TABLE]

According to the assumptions on $V$ and $K$, the functions $p\rightarrow|p|^{2}$ and $(q,p)\rightarrow H(q,p)$ are not comparable when either $|p|\rightarrow\infty$ or $|p|\rightarrow 0$ while $p$ remains bounded. Consequently, the function $|p|^{2}$ can not be bounded from below by $H$ and hence the inequality (2.7) fails to hold. This is mainly because the dissipative operator $\mathcal{O}$ only acts in the $p$-direction. To transfer the dissipation from $p$ to $q$ one needs to modify the Lyapunov function from $W=H$ to something else.

In the case where $U(q)$ has no singularity and grows at least quadratically at infinity, a Lyapunov function of the form

[TABLE]

was constructed previously (see e.g [34, 23]). By choosing $c>0$ small enough, one can justify that the function $H(q,p)$ and $W(q,p)$ defined above are comparable when $|p|^{2}+|q|^{2}\rightarrow\infty$. The trick of adding a $cp\cdot q$ term is essential to the proof of hypercoercivity of kinetic equations [35, 7] which shows that the generator is coercive with respect to some modified norms that are equivalent to $H^{1}$ or $L^{2}$ (weighted against the invariant measure $\pi$) although it is not coercive with respect to the corresponding standard norms.

Why adding a correction term $cp\cdot q$ term with $H$ leads to a Lyapunov function was explained in detail in [5, 13]. The key observation there is that at high energy level ($H\gg 1$), the system (3.1) evolves approximately around the deterministic orbit of the Hamiltonian dynamics ($\gamma=0$ in (3.1)). Therefore the average dissipative effect is approximately given by

[TABLE]

where $\langle|p|^{2}\rangle(h)$ is the average of $|p|^{2}$ along the deterministic Hamiltonian orbits at energy level $h$. To correct $-|p|^{2}+d\gamma\beta^{-1}$ on the right side of (3.2) to the average dissipation as above, a Lyapunov function of the form $W=H+\Psi$ was constructed where the corrector $\Psi$ solves the following Poisson equation:

[TABLE]

When $U(q)$ has no singularity and grows at least quadratically at infinity, it was shown that $\Psi(q,p)\sim cp\cdot q$ as $|q|\rightarrow\infty$. The same averaging idea was further exploited by [5] to construct Lyapunov function for the Langevin dynamics with a Lennard-Jones type singularity

[TABLE]

for some constants $a,b>0$. In this case, the operator $\mathcal{H}$ used in the Poisson equation (3.3) to define the corrector has to be replaced by another effective operator to take account of the blow-up of $H$ as $q\rightarrow 0$; see [5, Section 4] for more details. However, the idea of reducing stochastic dynamics to the deterministic Hamiltonian dynamics seems to be difficult to be generalized to more complicated singular potentials $U(q)$, such as (1.3).

Recently, the authors of [15] proposed a different idea to construct a Lyapunov function and used it to prove the geometric ergodicity of Langevin dynamics with a larger class of admissible singular potentials $U$. More specifically, their admissible condition on $U$ states that for any sequence $\{q_{k}\}\subset\mathcal{D}$ with $U(q_{k})\rightarrow\infty$,

[TABLE]

Note that this condition is fullffiled by both the Lennard-Jones type potential defined in (3.4) and the Coulomb potential in $\mathbf{R}^{d}$ with $d\geq 3$. However, the second limit fails to hold for two dimensional Coulomb potential and the Log potential in one dimension. Like in our the present paper, the choice of Lyapunov function in [15] was inspired by the dominating scaling at the singularity of the kernel (see the discussion below) and the scaling analysis which inspired the hypocoercivity constructions in [7]. However, the choice made in this paper is able to treat all Coulomb potentials in $\mathbf{R}^{d}$ for any dimension $d\geq 2$ uniformly. Moreover, the new Lyapunov function constructed here satisfies a Lyapunov condition with constants that have explicit dependence in $N$ (see Remark 4.1). The latter may be helpful for studying scaling limit problems as $N\rightarrow\infty$. A more geometric understanding of why this new Lyapunov function is better still eludes us.

Since our construction of Lyapunov function for the Coulomb system follows closely that of [15], it is worth to briefly review the techniques therein. Their construction relies on the following key observation (which was described earlier in [7]): under the following parabolic scaling

[TABLE]

as $\lambda\rightarrow 0$, the dynamics (3.1) formally reduces to its overdamped limit

[TABLE]

when the potential $U(q)$ is also properly scaled up to infinity as $|q|\rightarrow\infty$ or $|q|\rightarrow 0$. This observation is nice because the overdamped Langevin dynamics above has a natural Lyapunov function, namely the potential $U(q)$. As a result, effectively it suffices to consider building a Lyapunov function by focusing on the regime where $U$ (instead of $H$) is large while $p$ remains bounded. Based on such observation, the authors of [15] proposed a Lyapunov function of the form

[TABLE]

where $b>0$ and $\Psi(q,p)=o(H(q,p))$ as $H(q,p)\rightarrow\infty$. With this definition, assuming $\Psi\in C^{2}(\mathcal{X})$, one obtains that

[TABLE]

In order to let $W$ satisfy the Lyapunov condition (2.7), the corrector $\Psi$ has to be chosen appropriately in order to gain dissipation in $q$ when $U(q)$ is large. In particular, it is necessary to have for some $C>0$

[TABLE]

For this end, one need to consider scaling large $U$ in two different routes: $|q|\rightarrow\infty$ and $|q|\rightarrow 0$. For the Lennard-Jones potential defined in (3.4), it was shown that in both scaling regimes the dominating term in the generator is given by the following operator

[TABLE]

Hence effectively the condition (3.8) becomes

[TABLE]

A typical choice of $\Psi$ suggested by [15] is a function $\Psi\in C^{2}(\mathcal{X})$ and

[TABLE]

for some large constant $\kappa>0$ which leads to $\mathcal{A}\Psi=-\kappa$ on $\mathcal{X}\cap\{U\geq R\}$. Therefore on $\mathcal{X}\cap\{U\geq R\}$, it holds that

[TABLE]

where the inequality above follows from the Young’s inequality. Now we fix $\kappa=3d\gamma\beta^{-1}$ and $b<\beta^{-1}$. Then by condition (3.5), there exists $R>0$ large enough such that $\delta>0$ and that

[TABLE]

which implies that

[TABLE]

On the other hand, on $\mathcal{X}\cap\{U\leq R\}$, the right side of (3.9) is bounded. By combining the above estimates and (3.7), one sees that the function $W$ defined in (3.6) satisfies the Lyapunov condition (2.7). Again this relies on the crucial assumption (3.5) on $U(q)$. However, this assumption fails for two dimensional Coulomb potential or one dimensional Log potential.

3.2. Lyapunov function for the two-body system

We now build a new Lyapunov function for the two-body system (3.1) with a Coulomb potential $K$ in $\mathbf{R}^{d}$ for any $d\geq 2$. Similar to the previous discussion, we would like to consider scalings for $q$ at large $U$. Since the main issue comes from the singularity of the kernel, we first restrict our attention on the regime $|q|\rightarrow 0$. In particular, we consider the scaling

[TABLE]

with $\lambda\rightarrow\infty$. Then the generator $\mathcal{L}$, when rewritten in the $(Q,P)$-coordinates, takes the form

[TABLE]

where in the last line we have used the assumption that $V(q)\in C^{2}(\mathbf{R}^{d})$ and hence finite when $|q|$ is bounded and that the gradient of the Coulomb kernel is $\nabla K(q)=\frac{q}{|q|^{d}}$ for any $d\geq 2$. Therefore the dominating term of the generator in the regime $|q|\rightarrow 0$ is

[TABLE]

Motivated by this, we define a Lyapunov function similar to (3.6):

[TABLE]

with $a>0$. We would like to pick $\Psi$ so that $\mathcal{L}\Psi\approx\mathcal{A}\Psi$ becomes very negative when $U(q)\rightarrow\infty$ in the route $|q|\rightarrow 0$. In view of the prefactors in front of the derivation operators in (3.11), it is natural to consider $\Psi$ with the property that

[TABLE]

for some constants $b,\alpha>0$. In fact, with this choice, we have from (3.11) that

[TABLE]

In addition, by setting $0<\alpha\leq 1$, one sees that the first term on the right side of above in the case $d=2$ is always non-positive and hence obtains that

[TABLE]

In the subsequent discussion, we adopt the choice $\alpha=1$ for the purpose of simplifying calculations.

To obtain similar behavior about $\mathcal{L}\Psi$ in the other regime of large $U(q)$, namely $|q|\rightarrow\infty$, we simply add a $cp\cdot q$ term in our $\Psi$ since by the discussion in the last section this term gives extra dissipation in $q$-direction when $|q|$ is large. To sum up, the considerations above lead to the following new choice of corrector

[TABLE]

Notice that we have set $\alpha=1$ in (3.13). With this $\Psi$, we verify that the function $W$ defined (3.11) indeed satisfies the Lyapunov condition (2.7) with appropriate parameters $a,b,c$. Indeed, recall the Hamiltonian $H(q,p)=\frac{1}{2}|p|^{2}+V(q)+K(q)$. By definition,

[TABLE]

In view of the growth condition Assumption 2.1 (ii) on $V(q)$, it is clear that when $c>0$ is sufficiently small, $\Psi(q,p)=o(H(q,p))$ as $H\rightarrow\infty$ and thus $W(q,p)=\exp(H(q,p)(1+o(1)))$. This in particular implies that $W\rightarrow\infty$ as $H\rightarrow\infty$.

Next we show that $W$ satisfies the inequality (2.7). A straightforward calculation yields

[TABLE]

To evaluate the terms appearing in the above, we need the following ansatze:

[TABLE]

Consequently,

[TABLE]

Let us fix $a\in(0,\beta)$ and $b>0$. Then $I_{1}\leq 0$. Moreover, if $c<-a\gamma(1-a\beta^{-1})/2$, then

[TABLE]

In addition, thanks to Assumption (2.1), we have

[TABLE]

By Young’s inequality, it holds for any $\varepsilon>0$ that

[TABLE]

Lastly, one has

[TABLE]

From above, one sees that if one first sets $\varepsilon<\frac{a\gamma(1-a\beta^{-1})}{2}$ and then chooses $c>0$ sufficiently small so that the terms involving $|q|^{2}$ in $I_{4}$ can be absorbed by those in $I_{3}$, then one obtains that

[TABLE]

for some $c_{5},c_{6},c_{7}>0$. It is clear that $W(q,p)\rightarrow\infty$ leads to one of the following options: $|p|\rightarrow\infty,|q|\rightarrow\infty$ or $|q|\rightarrow 0$. The latter implies by (3.15) that there exists $\lambda>0$ such that

[TABLE]

when $R$ is sufficiently large. On the other hand, it follows from (3.15) that $\mathcal{L}W$ is bounded when $W(q,p)\leq R$. Hence we have shown that there exist constants $\lambda,K>0$ such that

[TABLE]

4. Construction of Lyapunov function: General case

This section devotes to the proof of Proposition 2.7, i.e. building a Lyapunov function for the dynamics of Coulomb gases with any finite number of particles ($N\geq 2$). Informed by the analysis of the two-body case in the previous section, we define with constants $b,c>0$ the correction term

[TABLE]

Note that the function above reduces to (3.14) when $(\mathbf{q},\mathbf{p})=((q,0),(p,0))$ and $N=2$. Hence, we can understand this choice as the many-particle analog of model problem presented in Section 3. We also define with $a>0$ the Lyapunov function

[TABLE]

By the same reasoning as in the previous section, one sees that when $c>0$ is sufficiently small, $\Psi(\mathbf{q},\mathbf{p})=o(H(\mathbf{q},\mathbf{p}))$ and hence $W(\mathbf{q},\mathbf{p})=\exp(aH(\mathbf{q},\mathbf{p})(1+o(1)))$. This implies that $W\rightarrow\infty$ as $H\rightarrow\infty$.

Next we continue to show that $W$ satisfies (2.7). In fact,

[TABLE]

We need the following derivatives:

[TABLE]

It follows that

[TABLE]

Moreover,

[TABLE]

By Assumption 2.1, the first term on the second line of above can be bounded by

[TABLE]

where we have used (2.4) with $\varepsilon=\frac{cc_{2}}{4b}$ in the penultimate line and (2.2) in the last line above. Thanks to Lemma 4.2, the third line of (4.3) can be bounded by

[TABLE]

It then follow from the last two estimates that

[TABLE]

In addition, it holds that

[TABLE]

Combining the estimates above yields

[TABLE]

Now we choose the parameters $a,b,c,\varepsilon_{1},\varepsilon_{2}>0$ such that

[TABLE]

Note that the last two inequality hold in particular when $\varepsilon_{1}+\varepsilon_{2}\leq\gamma a(1-a\beta^{-1})/2$ and $c$ is sufficiently small. With above choice of parameters, we have obtained that there exist $\alpha,C=C(N)>0$ such that

[TABLE]

Observe from the derivarion above that the constant $C(N)=\mathcal{O}(N^{d-1})$. By the same reasoning as in the previous section, the above inequality implies that there exists $\lambda>0$ such that when $R\gg 1$,

[TABLE]

On the other hand, it follows from (4.5) that $\mathcal{L}W$ is bounded on the set $\mathcal{X}\cap\{W(\mathbf{q},\mathbf{p})\geq R\}$. To conclude, we have

[TABLE]

for some constant $C$ depending on $N$.

Remark 4.1*.*

The dependence on $N$ of the constants appearing in the Lyapunov condition above can be made explicit. In fact, let us define for $\kappa>0$ the set

[TABLE]

From the estimate (4.5), one sees that if $(\mathbf{q},\mathbf{p})\in\mathcal{X}\setminus\mathcal{C}_{\kappa}$ with some $\kappa\gg 1$, then

[TABLE]

holds with $\lambda(N)=C_{\kappa}N^{d-1}$ for some $C_{\kappa}>0$. However, due to the exponential form of $W$, it follows from (4.5) that

[TABLE]

with $\alpha>0$ depending on the growth of $V$. Combining above leads to (4.6) with $\lambda=\lambda(N)=\mathcal{O}(N^{d-1})$ and $C=C(N)=\mathcal{O}(e^{N^{\alpha+d+2}})$.

Lemma 4.2**.**

For any $\mathbf{q}\in\mathcal{D}=\{\mathbf{q}=(q_{1},\cdots,q_{n})^{T}\in(\mathbf{R}^{d})^{N}:q_{i}\neq q_{j}\ \forall i,j=1,2\cdots,n\text{ with }i\neq j\}$, the following inequality holds:

[TABLE]

Proof.

It is straightforward to see that

[TABLE]

where the second equality follows from swapping the indices $k$ and $j$ in the summation in the second line of above.

We show that $J_{2}$ is non-negative. In fact, $J_{2}$ can be expanded as

[TABLE]

By relabelling the indices $i\rightarrow j,j\rightarrow i,k\rightarrow k$ in the second term on the right side and $i\rightarrow k,j\rightarrow i,k\rightarrow j$ in the third term, one gets

[TABLE]

We claim that the individual terms inside the summation above are non-negative. Indeed, let us define for fixed indices $i<j<k$ the angles $\{\theta_{m}\}_{m=1}^{3}$ as follows:

[TABLE]

It suffices to show that $S\geq 0$ with $S$ defined by

[TABLE]

We prove this by considering the following two cases.

Case 1: The three points $\{q_{i},q_{j},q_{k}\}$ lie on a common line. Without loss of generality we may assume that $q_{i}$ lies in the middle of the segment connecting $q_{j}$ and $q_{k}$. In this case, we have $\cos\theta_{1}=-1,\cos\theta_{2}=\cos\theta_{3}=1$. This leads to

[TABLE]

Case 2: The three points $\{q_{i},q_{j},q_{k}\}$ do not lie on a common line. In this case, the there points form a triangle with edges of lengths $|q_{i}-q_{j}|,|q_{i}-q_{k}|$ and $|q_{j}-q_{k}|$ and with angles $\{\theta_{m}\}_{m=1}^{3}$. Hence it holds that $\theta_{1}+\theta_{2}+\theta_{3}=\pi$ and $0<\theta_{i}<\pi,i=1,2,3$. If the triangle is an acute or right triangle, then it is obvious that $S\geq 0$ since $\cos(\theta_{i})\geq 0$ for every $i=1,2,3$. On the other hand, if the triangle is an obtuse triangle, then one of the angles, say $\theta_{1}$, is larger than $\pi/2$. As a result, we have $\theta_{2}+\theta_{3}<\pi/2$ and thus $\cos\theta_{2}\geq\cos(\theta_{2}+\theta_{3})=-\cos(\theta_{1})>0$ and $\cos\theta_{2}\geq\cos(\theta_{2}+\theta_{3})-\cos(\theta_{1})>0$. Consequently,

[TABLE]

Summarizing the discussion above yields $J_{2}\geq 0$. It follows that

[TABLE]

∎

A similar result to Lemma 4.2 was proved in [1, Lemma 5.2].

5. Proof of Main Results

5.1. Proof of Proposition 2.4

Proposition 2.4 follows immediately from the following stronger result.

Proposition 5.1**.**

For any initial point $\mathbf{x}=(\mathbf{q},\mathbf{p})\in\mathcal{X}$, there exits a unique strong pathwise solution $\mathbf{x}(t)$ to the equation (2.1). Furthermore, for $a\in(0,\beta)$, let $W$ be the Lyapunov function given in Proposition 2.7 with constants $\lambda,C$ defined in (2.9). Then for every $\mathbf{x}=(\mathbf{q},\mathbf{p})\in\mathcal{X}$ and $t\geq 0$,

[TABLE]

Proof.

For $R>0$, let us define stopping times

[TABLE]

Thanks to the coercivity of the function $W(\mathbf{x})$ from Lemma 5.3, we have the inclusion $\{\tau=\infty\}\subset\{\tau_{\partial\mathcal{X}}=\infty\}$. As a result, equation (2.1) has a strong pathwise solution $\mathbf{x}_{t}$ up to the stopping time $\tau_{R}$. We want to show that $\mathbf{P}_{\mathbf{x}}(\tau=\infty)=1$.

In fact, an application of Ito’s formula to the function $f(\mathbf{x},t)=e^{\lambda t}(W(\mathbf{x})-\frac{C}{\lambda})$ yields

[TABLE]

Now consider a fixed $t>0$ and letting $R\rightarrow\infty$. One can obtain from above and the fact that $W\geq 0$ that

[TABLE]

which implies that

[TABLE]

Thus we have obtained that $\mathbf{P}_{\mathbf{x}}(\tau_{R}>t)=1$ for any finite $t>0$. This in turn implies that the particles system $\mathbf{x}_{t}$ do not explode or collide in finite time.

Finally, using again the estimate (5.2), one has

[TABLE]

which leads to (5.1) by letting $R\rightarrow\infty$ since $\lim_{R\rightarrow\infty}\mathbf{P}_{\mathbf{x}}({\tau_{R}\geq t})=1$.

∎

Remark 5.2*.*

A similar stopping time argument using the Hamiltonian $H$ as a Lyapunov function instead of $W$ would lead to

[TABLE]

and hence also prevents explosion as well as collision leading to global well-posedness of solution. However, an estimate of the form (5.1) would not be available.

Lemma 5.3** (Coercivity of $W$ and $H$).**

Let $H(\mathbf{q},\mathbf{p})$ and $W(\mathbf{q},\mathbf{p})$ be given by (1.2) and (4.2) with parameters chosen as in (4.4). Then

[TABLE]

Proof.

Since we have seen in the previous section that $W(\mathbf{q},\mathbf{p})=e^{aH(\mathbf{q},\mathbf{p})(1+o(1))}$, we only need to show the above is valid for $H$. For this, it suffices to show that for any $R>0$, there exists $R_{1}>0,R_{2}>0,r>0$ such that $H(\mathbf{q},\mathbf{p})\geq R$ provided that one of the following holds: $|\mathbf{p}|\geq R_{1}$, or $|\mathbf{q}|\geq R_{2}$ or $\min_{i\neq j}|q_{i}-q_{j}|\leq r$.

This is obvious for $d\geq 3$. In fact, since both $K$ and $V$ are non-negative and $V$ grows at least quadratically at infinity by Assumption 2.1, in view of the expression

[TABLE]

the claim above is valid if one picks for example $R_{1}=\sqrt{2R},R_{2}=\sqrt{\frac{M+R}{c_{1}}},r=(\frac{1}{2NR})^{\frac{1}{d-2}}$.

For $d=2$, we follow [1, Lemma 3.2]. In fact, it follows from the simple inequality $\log x\leq x$ for any $x>0$ that

[TABLE]

where we used the Young’s inequality in the last line. Consequently, $H\geq R$ if either $|\mathbf{p}|\geq R_{1}=\sqrt{2(R+M+\frac{N}{2c_{1}})}$ or $|\mathbf{q}|\geq R_{2}=\sqrt{\frac{2}{c_{1}}(R+M+\frac{N}{2c_{1}})}$. Finally, we only need to show that $H\geq R$ when $|\mathbf{p}|\leq R_{1},|\mathbf{q}|\leq R_{2}$ and $|q_{i}-q_{j}|\leq r$ for some small $r>0$ and for some $1\leq i\neq j\leq N$. To see this, observe that for any $i\neq j$,

[TABLE]

where the last line follows from

[TABLE]

Hence if $|q_{i}-q_{j}|\leq r=e^{-N(R+\sqrt{N}R_{2})}$, then $H\geq R$. This completes the proof of the lemma. ∎

5.2. Proof of Theorem 2.5

The proof of Theorem 2.5 is based on Proposition 2.7 and a local minorization condition on the transition density $P_{t}(\mathbf{x},\cdot)$ as shown in the next proposition.

Proposition 5.4**.**

For any $a\in(0,\beta)$, let $W$ be the Lyapunov function defined in Proposition 2.7. Define for $R>0$ the compact set $\mathcal{C}_{R}=\{\mathbf{x}\in\mathcal{X}:W(\mathbf{x})\leq R\}$. Then for any sufficiently large $R>0$ and any $t_{0}>0$, there exists a probability measure $\nu$ on $\mathcal{X}$ and a constant $\alpha>0$ such that for every Borel set $A\subset\mathcal{X}$ and every $\mathbf{x}\in\mathcal{C}_{R}$,

[TABLE]

Proof.

The proof of the proposition is the same as the proof of [15, Corollary 5.12] (see also [23, Lemma 2.3]) and hence omitted. ∎

Proposition 5.5**.**

For every $\mathbf{x}\in\mathcal{X}$ and $t>0$, we have $\text{supp }P_{t}(\mathbf{x},\cdot)=\mathcal{X}$, where $\text{supp }\nu$ denotes the support of the measure $\nu$. In addition, the measure $P_{t}(\mathbf{x},\cdot)$ is absolutely continuous with respect to the Lebesgue measure on $\mathcal{X}$. Let $p_{t}(\mathbf{x},\mathbf{y})$ be the associated probability density. Then the mapping $(t,\mathbf{x},\mathbf{y})\mapsto p_{t}(\mathbf{x},\mathbf{y}):(0,\infty)\times\mathcal{X}\times\mathcal{X}\rightarrow[0,\infty)$ is continuous.

Proof.

The proof of the proposition is the same as the proof of [15, Proposition] and hence omitted. ∎

Proof of Theorem 2.5.

First it is straightforward to check that the measure $\pi$ is an invariant measure. The uniqueness follows from [12, Theorem 3.16] and the fact that $P_{t}$ is strong Feller and that $\text{supp }\nu=\mathcal{X}$ for any invariant measure $\nu$ as shown in Proposition 5.5.

To show the exponential convergence to the equilibrium we apply [14, Theorem 1.2] to the sampled Markov chain on $\mathcal{X}$ with transition density $P_{n}(\mathbf{x},\cdot):=P_{n\triangle t}(\mathbf{x},\cdot)$ with $\triangle t>0$ and $n\in\mathbb{N}$. In fact, Proposition 5.1 and Proposition 5.4 imply Assumption 1 and Assumption 2 in [14]. As a result of [14, Theorem 1.2], one obtains that there exists $C>0$ and $\delta\in(0,1)$ such that

[TABLE]

for all $\nu_{1},\nu_{2}\in\mathcal{P}_{W}$. To see the above bound holds for any $t$, let us write $t=n\triangle t+r$ with $r\in(0,\triangle t)$. For any measurable function $\varphi:\mathcal{X}\rightarrow\mathbf{R}$ with $|\varphi|_{W}\leq 1$, it follows from Proposition 5.1 that

[TABLE]

for some $C^{\prime}>0$ independent of $r>0$. As a consequence of Fubini theorem and Chapman-Kolmogorov equation, it holds that

[TABLE]

By defining $\eta=-\frac{\log\delta}{\triangle t}$ and $\tilde{C}=CC^{\prime}/\delta^{1/\triangle t}$, one obtains the desired bound

[TABLE]

The inequality (2.8) follows by setting $\nu_{2}=\pi$ and $\pi P_{t}=\pi$.

∎

6. Extension and outlook

Although the present paper focuses on the Langevin particle system with Coulomb kernel, the main theorem 2.5 is still valid for Riesz-type repulsive kernels

[TABLE]

The restriction $s<d$ is natural because otherwise the invariant measure is not normalizable and hence not a probability measure. In the Riesz case, Proposition 2.7 still holds with the same Lyapunov function (4.2) and the same corrector $\Psi$ given by (4.1). The proof proceeds in the same manner as in Section 4 as soon as one observes that the key inequality (4.7) in Lemma 4.2 is still valid with $d$ replaced by $s+2$.

For the same reason, our results also extend trivially to Langevin particles system (1.1) in 1D in which the kernel $K$ has Riesz or Log singularity. The only differece of the one dimensional system from systems in high dimensions is that the domain $\mathcal{D}$ of $U$ in this case is not pathwise connected. Therefore, the dynamics of the particles need to be restricted to one of the connected components (depending on the initial configuration). For example, one may set the state space $\mathcal{X}=\mathcal{D}\times\mathbf{R}^{N}$ with

[TABLE]

The convergence to the equilibrium result holds with the equilibrium measure $\pi$ restricted on $\mathcal{X}$ as well.

The convergence of Langevin dynamics with Coulomb kernel is studied via a probabilistic approach in this paper based on a novel construction of Lyapunov function. The Lyapunov function satisfies the Lyapunov condition (2.9) with the constants $\lambda,C$ depending on $N$ explicitly. However, the $N$-dependence of the constant $\alpha$ in the minorization condition (5.3) is very implicit. It remains an interesting question to obtain an explicit estimate for the minorization condition so that one can obtain an $N$-explicit convergence rate for the Langevin dynamics.

Acknowledgement

Both authors acknowledge the support of the National Science Foundation through the grant DMS-1613337. YL also acknowledges the partial support from the Duke Trinity College of Arts & Sciences Office of the Dean.