Langevin Dynamics with Variable Coefficients and Nonconservative Forces: From Stationary States to Numerical Methods

Abstract

Langevin dynamics is a versatile stochastic model used in biology, chemistry, engineering, physics and computer science. Traditionally, in thermal equilibrium, one assumes (i) the forces are given as the gradient of a potential and (ii) a fluctuation-dissipation relation holds between stochastic and dissipative forces; these assumptions ensure that the system samples a prescribed invariant Gibbs-Boltzmann distribution for a specified target temperature. In this article, we relax these assumptions, incorporating variable friction and temperature parameters and allowing nonconservative force fields, for which the form of the stationary state is typically not known a priori. We examine theoretical issues such as stability of the steady state and ergodic properties, as well as practical aspects such as the design of numerical methods for stochastic particle models. Applications to nonequilibrium systems with thermal gradients and active particles are discussed.

1. Introduction

Langevin dynamics is a system of stochastic differential equations of the form

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}\mathit{q}}{\mathrm{d}t}& =& {\mathit{M}}^{-1}\mathit{p},\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}\mathit{p}}{\mathrm{d}t}& =& \mathit{F}\left(\mathit{q}\right)-\mathbf{\Gamma }\left(\mathit{q}\right)\mathit{p}+\mathbf{\Sigma }\left(\mathit{q}\right)\dot{W}.\hfill \end{array}$$

(2)

Here $\mathit{q},\mathit{p}$ represent vectors of the positions and momenta of the particles comprising together $2n\in \mathbb{N}$ degrees of freedom. The mass matrix $\mathit{M}$ is symmetric positive definite. The force $\mathit{F}$ is usually taken to be conservative, i.e., $\mathit{F}\left(\mathit{q}\right)=-{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)$ for some potential energy function U. Langevin dynamics describes a physical system of particles moving under prescribed interaction forces and subject to collisions with particles of a “heat bath.” The friction matrix $\mathbf{\Gamma }$, which is here allowed to vary with position, typically models drag on a set of distinguished particles due to the interactions with the surrounding environment, whereas the matrix $\mathbf{\Sigma }$, also potentially varying with position, characterizes the stochastic effects of collisions. $\mathit{W}$ represents a vector of independent and uncorrelated Wiener processes in ${\mathbb{R}}^m$, $\mathit{W}={[W_1,\dots W_m]}^T$, where the components of its formal time derivative $\dot{W}$ are independent white noise processes such that

$$\mathbb{E}\left[{\dot{W}}_i\left(t\right){\dot{W}}_j\left(t^{\prime }\right)\right]={\delta }_{i,j}\delta \left(t-t^{\prime }\right).$$

(3)

and

$$\mathbb{E}\left[{\dot{W}}_i\left(t\right)\right]=0.$$

(4)

In case $\mathbf{\Gamma }\left(\mathit{q}\right)\equiv \gamma \mathbf{I}$ where $\gamma$ is a nonnegative scalar and constant, one may take $\mathbf{\Sigma }\left(\mathit{q}\right)\equiv \sigma \mathbf{I}$, where $\sigma =\sqrt{2\gamma k_{\mathrm{B}}T}$, T is the temperature and $k_{\mathrm{B}}$ is Boltzmann’s constant. The resulting system can then be shown, under conditions reviewed in Section 2, to have a unique invariant distribution with density ${\rho }_{\beta }(\mathit{q},\mathit{p})\propto \exp (-\beta H(\mathit{q},\mathit{p}))$, where $H(\mathit{q},\mathit{p})={\mathit{p}}^T{\mathit{M}}^{-1}\mathit{p}/2+U\left(\mathit{q}\right)$ represents the total energy of the isolated deterministic system and ${\beta }^{-1}=k_{\mathrm{B}}T$.

The emphasis in this article is on the general form (1) and (2), for which we will relax several typical assumptions. First, the system is not assumed to admit a universal fluctuation-dissipation relation, instead we assume only certain nondegeneracy conditions. Second, we do not assume a conservative force field. Such generalized forms of Langevin dynamics can be used to model diffusion (or thermal) gradients by particle simulation [1,2,3], as well as a variety of models for flocking [4,5,6], protein folding [7] and bacterial suspensions [8]. Problems with nonconservative forces are considered often in the physics literature, but the precise characterization of the convergence to stationary conditions is rarely discussed. In the general case, the form of the stationary distribution is non-obvious and the uniqueness of the steady state and its asymptotic stability are not well understood. An important contribution of this article is to show the existence of a steady state for the model (1) and (2) under relaxed conditions and to explicitly construct a Lyapunov function, thus allowing us to conclude the geometric convergence of observable averages.

Our approach to proving ergodicity is dependent on the fact that the friction and noise terms in the equations depend only on positions. On the other hand, within this class, and subject to the nondegeneracy condition on the friction and noise tensors, our treatment is general and has a considerable range of applications. Formally (1) and (2) includes dissipative particle dynamics (DPD), which is a momentum-conserving, 2nd order, gradient-type system [9,10], but we specifically exclude in this article cases for which $\mathbf{\Gamma }$ is not positive definite, which occur as a direct consequence of momentum conservation in DPD; ergodic properties for DPD systems in one dimension were discussed in previous work of Shardlow and Yan [11]. See also [12,13,14] for some special cases of Langevin-type systems with velocity dependent coefficients to which our formulation is not applicable.

To perform numerical discretization of the system (1) and (2), we proceed by splitting. This involves performing an additive decomposition of the stochastic vector field into several component parts, each of which can be exactly integrated (in the sense of distributions). The resulting stochastic maps can be composed to approximate the solution in a single timestep. Even for a particular choice of decomposition, there are many ways to combine the sequence of steps. The properties of various choices of the splitting method have been examined in detail for constant $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$ in [15]. In that work it has been shown that the numerical methods constructed using certain components inherit the ergodicity properties of the SDE system, whereas the invariant measure of the numerical method differs from that of the SDE. The error in the invariant measure can be characterized in terms of the order of accuracy of averages of test functions. In general, it is found that symmetric compositions are preferred as they exhibit, under mild conditions, even order of accuracy for the invariant measure (thus one obtains second order accuracy for the same computational work as a first order scheme).

In the case of the more general systems considered in this article, the calculation of the Ornstein-Uhlenbeck solution, which is required at each step of our splitting-based numerical methods, becomes potentially demanding from a computational standpoint. In the methods of [15,16] this is done exactly, however in the current setting the solution of a Lyapunov equation would need to be obtained at each step in time, a calculation that would dominate the computational load in a large scale simulation. Therefore, we offer an efficient numerical procedure based on multiple timestepping [17], relying on a further splitting of the OU equation at each interior timestep.

In the last section of this article, we illustrate the theory in three numerical examples. The first is a particle system with an imposed diffusion gradient, as in [7,18,19,20,21]. Here the issue is to capture the correct density fluctuations as a function of position. In the second model, we incorporate a stirring force which alters the equilibrium state of a physical model. In the last example, we discuss the use of a blended Langevin model which draws on ideas from the literature of flocking, in particular the Cucker-Smale system. We use our theory to infer attractive steady states for the system, and characterize flocking tendencies by the use of two order parameters: one modelling the formation of consensus and the other characterizing the peculiarity which can be viewed as the average internal energy of isolated clumps of matter.

2. Stationary States of SDEs and Their Stability

In this section, we briefly outline the general theory on which our analysis of the ergodic properties of the system (1) and (2) in the next section is based. We are following in large parts the presentation in the review articles [22,23] and we refer the reader to these articles for a more detailed and comprehensive presentation.

Let ${\left(\mathit{x}\left(t\right)\right)}_{t\ge 0}$ denote the solution of an Itô diffusion process of the form

$$\dot{\mathit{x}}=\mathit{a}\left(\mathit{x}\right)+\mathit{B}\left(\mathit{x}\right)\dot{W},\mathit{x}\left(0\right)\sim {\mu }_0,$$

(5)

taking values in a suitable domain. In this general presentation, $\mathit{x}\left(t\right)\in {\Omega }_x$ may represent either a position vector or else the combined vector of positions and momenta, i.e., $\mathit{x}\left(t\right)=(\mathit{q}\left(t\right),\mathit{p}\left(t\right))\in {\Omega }_q\times {\mathbb{R}}^n$. In the latter case, we would typically either assume a compact configurational domain, as for example when periodic boundary conditions are used in the position space, i.e., ${\Omega }_q=L{\mathbb{T}}^n$, where $L>0$ and $\mathbb{T}=\mathbb{R}/\mathbb{Z}$ denotes the 1-torus, or else an unbounded domain, e.g., ${\Omega }_q={\mathbb{R}}^n$. ${\left(\mathit{W}\left(t\right)\right)}_{t\ge 0}$ represents a standard Wiener process in ${\mathbb{R}}^m$ and the coefficients $\mathit{a}$, $\mathit{B}$ are assumed to be smooth, i.e., $\mathit{a}\in {\mathcal{C}}^{\infty }({\Omega }_x,{\mathbb{R}}^{2n})$, and $\mathit{B}\in {\mathcal{C}}^{\infty }({\Omega }_x,{\mathbb{R}}^{2n\times m})$. The initial state of $\mathit{x}$ is specified by the probability measure ${\mu }_0$, which throughout this article we assume to be such that $\mathit{x}\left(0\right)$ has finite mean and variance, i.e.,

$${\int }_{{\Omega }_x}{\parallel \mathit{x}\parallel }_2^2{\mu }_0\left(\mathrm{d}\mathit{x}\right)<\infty .$$

2.1. The Associated Semigroup of Evolution Operators and Their Adjoints

The time evolution of the expectations of observables under the dynamics of the SDE is described by the semigroup of evolution operators ${\left({\mathcal{P}}_t\right)}_{t\ge 0}$

$$\left({\mathcal{P}}_t\phi \right)\left(\mathit{x}\right):=\mathbb{E}\left[\phi \left(\mathit{x}\left(t\right)\right)\right|\mathit{x}\left(0\right)=\mathit{x}],$$

for $\phi \in \mathcal{S},\mathit{x}\in {\Omega }_x$, where the expectation is taken with respect to the Wiener measure associated with the driving noise process ${\left(\mathit{W}\left(t\right)\right)}_{t\ge 0}$ and $\mathcal{S}\subseteq \mathcal{M}({\Omega }_x,\mathbb{R})$ denotes a set of test functions or observables, which is contained in $\mathcal{M}({\Omega }_x,\mathbb{R})$, the set of all real valued measurable functions defined on the domain ${\Omega }_x$. If the test function set $\mathcal{S}$ is chosen appropiately, (e.g., $\mathcal{S}={\mathcal{C}}_b^{\infty }({\Omega }_x,\mathbb{R})$, where ${\mathcal{C}}_b^{\infty }({\Omega }_x,\mathbb{R})$ denotes the set of all smooth bounded real valued functions defined on ${\Omega }_x$) the action of ${\mathcal{P}}_t$ corresponds to the solution of an initial value problem, i.e., $\left({\mathcal{P}}_t\phi \right)\left(\mathit{x}\right)=u(\mathit{x},t)$, where u solves

$$\frac{\partial }{\partial t}u(\mathit{x},t)=\mathcal{L}u(\mathit{x},t),u(\mathit{x},0)=\phi \left(\mathit{x}\right).$$

(6)

The operator $\mathcal{L}$, defined such that

$$\left(\mathcal{L}\phi \right)\left(\mathit{x}\right)=\underset{\tau \to 0}{\lim }\frac{\mathbb{E}\left[\phi \right(\mathit{x}\left(\tau \right)\left)\right|\mathit{x}\left(0\right)=\mathit{x}]-\phi (\mathit{x})}{\tau },$$

for all $\phi \in \mathcal{S}$, is referred to as the infinitesimal generator. As the action of the operators ${\left({\mathcal{P}}_t\right)}_{t\ge 0}$ is given as the solution of the differential Equation (6), it is common to use the notation ${\mathcal{P}}_t:=e^{t\mathcal{L}}$. In the case of an Itô diffusion process (5), and with the regularity assumptions on the coefficients $\mathit{a},\mathit{B}$ stated above, it can be shown that the infinitesimal generator takes the form

$$\mathcal{L}=\mathit{a}·\nabla +\frac{1}{2}{\mathit{B}}^T\mathit{B}:{\nabla }^2,$$

where “$:$” denotes the Frobenius product, i.e.,

$${\mathit{B}}^T\mathit{B}:{\nabla }^2=\sum _{i,j=1}^n{\left[{\mathit{B}}^T\mathit{B}\right]}_{i,j}{\partial }_{x_i}{\partial }_{x_j},$$

See, e.g., [24] for the case $\mathcal{S}={\mathcal{C}}_b^{\infty }({\Omega }_x,\mathbb{R})$, and [25] for extensions from that core to more general test function sets such as the ones considered below.

It is often convenient to adopt a dual perspective by considering the evolution of the density of the law ${\mu }_t$ of $\mathit{x}\left(t\right)$ in time. Let ${\mu }_t\left(\mathrm{d}\mathit{x}\right)=\rho (\mathit{x},t)\mathrm{d}\mathit{x}$, for $t\ge 0$, in the sense of distributions. The corresponding semigroup ${({\mathcal{P}}_t^{†})}_{t\ge 0}$ of transfer operators is then naturally defined as the formal adjoint operators of ${\mathcal{P}}_t$, i.e.,

$${\int }_{{\Omega }_x}\left({\mathcal{P}}_t\phi \right)\left(\mathit{x}\right){\rho }_0\left(\mathit{x}\right)\mathrm{d}\mathit{x}={\int }_{{\Omega }_x}\phi \left(\mathit{x}\right){\mathcal{P}}_t^{†}{\rho }_0\left(\mathit{x}\right)\mathrm{d}\mathit{x},$$

assuming $\mathit{x}\left(0\right)\sim {\mu }_0$.

Similarly as for its adjoint, the action of ${\mathcal{P}}_t^{†}$ corresponds to the solution of an initial value problem, which in this case is known as the Fokker-Planck equation. More specifically,

$$\frac{\partial }{\partial t}\rho (\mathit{x},t)={\mathcal{L}}^{†}\rho (\mathit{x},t),\rho (\mathit{x},0)={\rho }_0\left(\mathit{x}\right),$$

(7)

where $\rho (·,t)$ denotes the probability density of ${\mu }_t$ and ${\mathcal{L}}^{†}$—the Fokker-Planck operator—can be shown to correspond to the $L^2$-adjoint of the infinitesimal generator $\mathcal{L}$, i.e.,

$${\mathcal{L}}^{†}\rho =-\nabla ·\left(\mathit{a}\rho \right)+{\nabla }^2:\left(\frac{1}{2}{\mathit{B}}^T\mathit{B}\rho \right),$$

thus

$$e^{t{\mathcal{L}}^{†}}={\mathcal{P}}_t^{†}.$$

2.2. Hypoellipticity and Existence of a Smooth Transition Kernel

Note that (7) is in general to be interpreted in a weak sense as, up to this point, we have not made any assumptions on the regularity of $\rho$. However, within the scope of this article it is sufficient to consider the case where the differential operator ${\partial }_t-{\mathcal{L}}^{†}$ is hypoelliptic. A differential operator A is said to be hypoelliptic, if for any g solving the differential equation $Ag=f$, it follows that g is of higher regularity than f in the sense that

$$f\in H_s^{\mathrm{loc}}\Rightarrow g\in H_{s+ϵ}^{\mathrm{loc}},$$

with $ϵ>0$, where $H_s^{\mathrm{loc}}$ denotes the local Sobolev space of order $s\in \mathbb{N}$. This means that if ${\partial }_t-{\mathcal{L}}^{†}$ is hypoelliptic, then the solution of (7) is smooth in the sense that $\rho \in {\mathcal{C}}^{\infty }({\Omega }_x,(0,\infty ))$ irrespective of the regularity of ${\rho }_0$. A common way to establish hypoellipticity of a differential operator is via Hörmander’s theorem ([26], Theorem 22.2.1, on p. 353):

Theorem 1.

Let A be a differential operator of the form

$$A={\mathit{a}}_0·\nabla +\sum _{i=1}^M{({\mathit{a}}_i·\nabla )}^{†}({\mathit{a}}_i·\nabla ),$$

where ${\mathit{a}}_i,0\le i\le M$ are ${\mathcal{C}}^{\infty }$ vector fields in ${\mathbb{R}}^n$ and ^† indicates the formal $L^2$ adjoint. Iteratively define a collection of vector fields by

$${\mathcal{V}}_0=\{{\mathit{a}}_i:i\ge 0\},{\mathcal{V}}_{k+1}={\mathcal{V}}_k\cup \{[\mathit{v},{\mathit{a}}_i]:\mathit{v}\in {\mathcal{V}}_k,0\le i\le M\},$$

(8)

where

$$[\mathit{X},\mathit{Y}]=(\nabla \mathit{Y})\mathit{X}-(\nabla \mathit{X})\mathit{Y},$$

denotes the commutator of vector fields $\mathit{X},\mathit{Y}\in {\mathcal{C}}^{\infty }({\Omega }_x,{\mathbb{R}}^n)$ and $(\nabla \mathit{X}),(\nabla \mathit{Y})$ their Jacobian matrices. If

$$\forall \mathit{x}\in {\mathbb{R}}^n,\mathrm{lin}\left\{\mathit{v}\left(\mathit{x}\right):\mathit{v}\in \bigcup _{k\in \mathbb{N}}{\mathcal{V}}_k\right\}={\mathbb{R}}^n,$$

(9)

then the operator A is hypoelliptic.

The condition (9) is commonly referred to as Hörmander’s condition. In the particular case of $A={\partial }_t-{\mathcal{L}}^{†}$ one can easily verify that (9) is exactly satisfied if

$${\mathcal{V}}_0=\{{\mathit{B}}_i:i\ge 1\},{\mathcal{V}}_{k+1}={\mathcal{V}}_k\cup \{[\mathit{v},{\mathit{B}}_i]:\mathit{v}\in {\mathcal{V}}_k,0\le i\le M\}.$$

with ${\mathit{B}}_0=\mathit{a}$ and ${\mathit{B}}_i$ refers to the i-th column of the diffusion tensor $\mathit{B}$ in (5). This particular version of Hörmander’s condition adapted to the parabolic PDE of the form (7) is referred to as the parabolic Hörmander condition.

The direct consequence of the parabolic Hörmander condition is the smoothness of the underlying transition kernel describing the evolution of the probability measure associated to the SDE.

2.3. Ergodicity and Convergence in Law

Let now $\mu \left(\mathrm{d}\mathit{x}\right)=\rho \left(\mathit{x}\right)\mathrm{d}\mathit{x}$ be a probability measure with density $\rho \in {\mathcal{C}}^2({\Omega }_x,[0,\infty ))$. When we say that $\mu$ is an invariant measure of the SDE (5) or that the latter preserves the probability measure $\mu$, we mean that the probability density $\rho$ is a stationary solution of the Fokker-Planck equation associated to the SDE, i.e.,

$${\mathcal{L}}^{†}\rho =0.$$

(10)

Note that for (10) to be well posed in a strong sense it is sufficient that $\mathcal{L}$ is hypoelliptic, i.e., the Fokker-Planck operator ${\mathcal{L}}^{†}$ satisfies the Hörmander condition. It is also straightforward to see that every convex combination of two invariant measures is again an invariant measure of the respective SDE which by definition means that the set of invariant measures of a particular SDE is convex.

Let for $t>0$ and some suitable observable $\phi$,

$${\overline{\phi }}_t:=\frac{1}{t}{\int }_0^t\phi \left(\mathit{x}\left(s\right)\right)\mathrm{d}s,$$

denote the finite time average of $\phi$ evaluated along the path of the solution $\mathit{x}$ of (5). Similarly, we write

$$\mu \left[\phi \right]:={\int }_{{\Omega }_x}\phi \left(\mathit{x}\right)\mu \left(\mathrm{d}\mathit{x}\right).$$

for the $\mu$-weighted spatial average of $\phi$. The process $\mathit{x}$ is said to be ergodic with respect to the invariant probability measure $\mu$ if for all $\phi \in L^1\left(\mu \right)$ and for almost all realizations of the Wiener process $\mathit{W}$, and ${\mu }_0$-almost all initial values $\mathit{x}\left(0\right)$, trajectory averages coincide with expectations with respect to the measure $\mu$ in the asymptotic limit $t\to \infty$, i.e.,

$$\underset{t\to \infty }{\lim }{\overline{\phi }}_t=\mu \left[\phi \right],\mathrm{a}.\mathrm{s}.$$

(11)

It can be shown (see [27]) that ergodicity follows if (i) there exists an invariant measure with positive smooth density and (ii) ${\partial }_t-{\mathcal{L}}^{†}$ is hypoelliptic.

Assume that $\mathit{x}\left(t\right)$ converges in law towards a unique invariant measure $\mu$, i.e.,

$$\underset{t\to \infty }{\lim }\mathbb{E}\left[\phi \left(\mathit{x}\left(t\right)\right)\right|\mathit{x}\left(0\right)=\mathit{x}]=\mu \left[\phi \right],$$

(12)

for all $\phi \in {\mathcal{C}}_b^{\infty }({\Omega }_x,\mathbb{R})$, ${\mu }_0$-almost all $\mathit{x}\in {\Omega }_x$. A common way to characterize the convergence of the expectation $\mathbb{E}\left[\phi \right(\mathit{x}\left(t\right)\left)\right|\mathit{x}\left(0\right)=\mathit{x}]$, or, in semi group notation the convergence of $e^{t\mathcal{L}}\phi \left(\mathit{x}\right)$ to the $\mu$-weighted average $\mu \left[\phi \right]$, is via functional decay estimates of the semi-group operators $e^{t\mathcal{L}}$. For this purpose a set of test functions $\mathcal{S}$ is fixed and equipped with a norm ${\parallel ·\parallel }_{\mathcal{S}}$ such that $E:=(\mathcal{S},{\parallel ·\parallel }_{\mathcal{S}})$ forms a Banach-space. Of particular interest in this context is exponential convergence of $e^{t\mathcal{L}}\phi$ towards $\mu \left[\phi \right]$ in the respective norm, i.e.,

$${\parallel e^{t\mathcal{L}}\phi -\mu \left[\phi \right]\parallel }_{\mathcal{S}}\le C{e}^{-\kappa t}{\parallel \phi -\mu \left[\phi \right]\parallel }_{\mathcal{S}},$$

(13)

where $C,\kappa$ are positive constants, the latter corresponding to the spectral gap of the generator $\mathcal{L}$ in the functional space $E_0=({\mathcal{S}}_0,{\parallel ·\parallel }_{\mathcal{S}})$, where ${\mathcal{S}}_0\subseteq \mathcal{S}$ denotes the subset of test functions with vanishing mean, i.e.,

$${\mathcal{S}}_0=\left\{\phi \in \mathcal{S}:\mu \left[\phi \right]=0\right\}.$$

Let the operator $\Pi$ denote the orthogonal projection from $\mathcal{S}$ onto ${\mathcal{S}}_0$, i.e.,

$$\Pi \phi =\phi -\mu \left[\phi \right],\phi \in \mathcal{S}.$$

Denote further by ${\parallel A\parallel }_{\mathcal{B}\left(E\right)}$ the operator norm

$${\parallel A\parallel }_{\mathcal{B}\left(E\right)}:=\underset{\begin{array}{c}\phi \in E\\ \phi \ne 0\end{array}}{\sup }\frac{{\parallel A\phi \parallel }_{\mathcal{S}}}{{\parallel \phi \parallel }_{\mathcal{S}}},$$

of an operator $A:E\to E$. Equation (13) implies that $e^{t\mathcal{L}}\Pi$ when considered as an operator on E is bounded in the operator norm $\mathcal{B}\left(E\right)$ as

$${\parallel e^{t\mathcal{L}}\Pi \parallel }_{\mathcal{B}\left(E\right)}\le C{e}^{-\kappa t}.$$

(14)

2.4. Finite-Time Averages and the Central Limit Theorem

Ergodicity of an SDE ensures that time averages of infinitely long trajectories almost surely coincide with spatial averages with respect to the target measure. However, ergodicity per se does not address the statistical properties of finite-time averages apart from the convergence of the time average ${\overline{\phi }}_t$ as $t\to \infty$. For practical applications, where for a given unique invariant measure, the aim is to approximate the $\mu$-weighted average $\mu \left[\phi \right]$ of a test function $\phi \in \mathcal{S}$, it is important that fluctuations of finite-time averages ${\overline{\phi }}_t$ (i.e., the Monte Carlo error in the finite-time approximations) around the infinite-time value ${\lim }_{t\to \infty }{\overline{\phi }}_t=\mu \left[\phi \right]$ can be quantified. For this purpose it typically regarded as necessary that a central limit theorem holds, i.e.,

$$\sqrt{t}\left({\overline{\phi }}_t-\mu \left[\phi \right]\right)\sim \mathcal{N}(0,{\sigma }_{\phi }^2),\mathrm{as}t\to \infty ,$$

(15)

where ${\sigma }_{\phi }^2$ is commonly referred to as the asymptotic variance of the observable $\phi$. A sufficient condition for a central limit theorem of the form (15) to hold for the solution $\mathit{x}$ of (5) and the observable $\phi$ is that the Poisson equation

$$-\mathcal{L}\Phi =\Pi \phi ,$$

(16)

has a solution which belongs to $L^2\left(\mu \right)$ (see [28]). The asymptotic variance then takes the form

$${\sigma }_{\phi }^2=-2{\int }_{{\Omega }_x}\left({\mathcal{L}}^{-1}\Pi \phi \right)\Pi \phi \mathrm{d}\mu .$$

Let $L_0^2\left(\mu \right)$ denote the subspace of $L^2\left(\mu \right)$ consisting of functions with vanishing mean. Note that for $\Phi \in L^2\left(\mu \right)$ it is a priori not clear how to interpret (16) since only under additional regularity assumptions (16) can be interpreted in a weak sense. A common way to makes sense of (16) is by deriving bounds for the operator ${\mathcal{L}}^{-1}$ in $\mathcal{B}\left(E_0\right)$ where $E_0$ is some subspace of $L_0^2\left(\mu \right)$. If ${\mathcal{L}}^{-1}$ is bounded in $\mathcal{B}\left(E_0\right)$, this then directly implies $\Phi \in L^2\left(\mu \right)$ in (16) for $\phi \in E_0$. The relationship between the spectral properties of the operator ${\mathcal{L}}^{-1}\Pi$ and the convergence properties of the solution process $\mathit{x}$, or more specifically the decay properties of the semi-group operator $e^{t\mathcal{L}}\Pi$ as $t\to \infty$ can be made more precise via the formal identity

$${\mathcal{L}}^{-1}\Pi =-{\int }_0^{\infty }{e}^{t\mathcal{L}}\Pi \mathrm{d}t,$$

(17)

which follows from

$$-{\int }_0^{\infty }\mathcal{L}{e}^{t\mathcal{L}}\phi \mathrm{d}t=-{\int }_0^{\infty }\left(\frac{\mathrm{d}}{\mathrm{d}t}{e}^{t\mathcal{L}}\phi \right)\mathrm{d}t=\phi ,$$

(18)

for $\phi \in \{\varphi \in {\mathcal{S}}_0:\mathcal{L}\varphi \in {\mathcal{S}}_0\}$. Using the identity (17), one directly finds that (14) is a sufficient condition for ${\mathcal{L}}^{-1}\Pi$ to be bounded since

$$\begin{array}{cc}\hfill \parallel {{\mathcal{L}}^{-1}\Pi \phi \parallel }_{\mathcal{S}}& ={\parallel {\int }_0^{\infty }{e}^{t\mathcal{L}}\Pi \phi \mathrm{d}t\parallel }_{\mathcal{S}}\hfill \\ & \le {\int }_0^{\infty }{\parallel e^{t\mathcal{L}}\Pi \phi \parallel }_{\mathcal{S}}\mathrm{d}t\hfill \\ & \le C{\int }_0^{\infty }{e}^{-\kappa t}\parallel {\Pi \phi \parallel }_{\mathcal{S}}\mathrm{d}t\hfill \\ & \le \frac{C}{\kappa }\parallel {\Pi \phi \parallel }_{\mathcal{S}}.\hfill \end{array}$$

(19)

We conclude that the exponential decay (13) implies the central limit theorem (15) for $\phi$ in the corresponding function space $E_0$. Estimates for $E_0=H^1\left(\mu \right)\cap L_0^2\left(\mu \right)$ can be obtained using the framework of hypocoercivity as presented in [29]. In [30] techniques are introduced to show the decay estimate (13) for $E_0=L_0^2\left(\mu \right)$. In this article we will use Lyapunov function-based techniques which allow to show exponential convergence as in (13) in some weighted $L^{\infty }$ spaces, which we specify in the next section.

2.5. Exponential Convergence in Weighted $L^{\infty }$ Spaces

Another way of deriving exponential decay estimates for the semi-group ${(e^{t\mathcal{L}})}_{t\ge 0}$, which are sufficient to establish a central limit theorem for certain observables, is by means of well established Lyapunov techniques. These techniques have been formulated originally for discrete-time Markov processes/Markov chains [31,32,33] and have been subsequently extended to continuous time solutions processes of SDEs [22,23,34,35]. The function space on which decay estimates are shown in these references are Banach spaces of the form $(L_{\mathcal{K}}^{\infty },{\parallel ·\parallel }_{L_{\mathcal{K}}^{\infty }})$, with

$$L_{\mathcal{K}}^{\infty }:=\left\{\phi \mathrm{measurable}|\frac{\phi }{\mathcal{K}}\in L^{\infty }\left({\Omega }_{\mathit{x}}\right)\right\},$$

(20)

where $\mathcal{K}\in {\mathcal{C}}^2({\Omega }_x,[1,\infty ))$ is a positive function, and

$${\parallel \phi \parallel }_{L_{\mathcal{K}}^{\infty }}:={\parallel \frac{\phi }{\mathcal{K}}\parallel }_{L^{\infty }}.$$

(21)

Exponential convergence in the sense of

$${\parallel e^{t\mathcal{L}}\phi -\mu \left[\phi \right]\parallel }_{L_{\mathcal{K}}^{\infty }}\le C{e}^{-\kappa t}{\parallel \phi -\mu \left[\phi \right]\parallel }_{L_{\mathcal{K}}^{\infty }},$$

(22)

is typically shown provided that the following two assumptions are satisfied.

Assumption 1 (Infinitesimal Lyapunov condition).

There is a function $\mathcal{K}\in {\mathcal{C}}^2({\Omega }_x,[1,\infty ))$ with ${\lim }_{∥\mathit{x}∥\to \infty }\mathcal{K}\left(\mathit{x}\right)=\infty$, and real numbers $a\in (0,\infty ),b\in (0,\infty )$ such that,

$$\mathcal{L}\mathcal{K}\le -a\mathcal{K}+b.$$

(23)

Assumption 2.

For some $t>0$ there exists a constant $\eta \in (0,1)$ and a probability measure ν such that

$$\underset{\mathit{x}\in \mathcal{C}}{\inf }{e}^{t{\mathcal{L}}^{†}}{\delta }_{\mathit{x}}\left(\mathrm{d}\mathit{y}\right)\ge \eta \nu \left(\mathrm{d}\mathit{y}\right)$$

where $\mathcal{C}=\{\mathit{x}\in {\Omega }_x:\mathcal{K}\left(\mathit{x}\right)\le {\mathcal{K}}_{\max }\}$ for some ${\mathcal{K}}_{\max }>1+2b/a,$ where $a,b$ are the same constants as in (23).

As outlined in [23,35], Assumption 2 follows if $\mathcal{L}$ satisfies the parabolic Hörmander condition and the SDE (5) is controllable in the sense that there is a $t>0$ such that for any pair ${\mathit{x}}_{-},{\mathit{x}}_{+}\in \mathcal{C}$, there exists a continuous control $u\in L^1([0,t],{\Omega }_x)$, such that the solution $\tilde{\mathit{x}}$ of the differential equation

$$\dot{\tilde{\mathit{x}}}=\mathit{a}\left(\tilde{\mathit{x}}\right)+\mathit{B}\left(\tilde{\mathit{x}}\right)u,$$

satisfies $\tilde{\mathit{x}}\left(0\right)={\mathit{x}}_{-}$ and $\tilde{\mathit{x}}\left(t\right)={\mathit{x}}_{+}$. The following theorem is derived in [22] using results from [33]. Similar results can be found e.g., in [23,36].

Theorem 2.

Suppose that Assumptions 1 and 2 hold. Then, the solution of the SDE (5) admits a unique invariant probability measure μ, such that

$${\int }_{{\Omega }_x}\mathcal{K}\mathrm{d}\mu <\infty .$$

(24)

Moreover, there exist $C>0$ and $\kappa >0$ such that (22) holds for all $\phi \in L_{\mathcal{K}}^{\infty }\left({\Omega }_x\right)$ and any $t\in [0,\infty )$.

Finally, let us demonstrate how decay estimates in spaces $L_{\mathcal{K}}^{\infty }$ can be used to derive a central limit theorem for certain functions. Let $\mathcal{V}$ be a Lyapunov function such that the conditions for Theorem 2 are satisfied for $\mathcal{K}=\mathcal{V}$. Note that if the conditions of Theorem 2 are also satisfied for ${\mathcal{V}}^2$, then this implies that a central limit theorem holds for all observables $\phi \in L_{\mathcal{V}}^{\infty }$, since (24) being valid for $\mathcal{K}={\mathcal{V}}^2$ implies

$$L_{\mathcal{V}}^{\infty }\subset L^2\left(\mu \right).$$

Thus, the inequality (19) for $\mathcal{S}=L_{{\mathcal{V}}^2}^{\infty }$ again implies that the solution $\Phi$ of (16) is contained in $L^2\left(\mu \right)$ for $\phi \in L_{\mathcal{V}}^{\infty }$, so that by [28] indeed a central limit theorem of the form (15) holds for $\phi \in L_{\mathcal{V}}^{\infty }$. This motivates to show (22) for a wide class of Lyapunov functions. In the next section, we consider the case of the Langevin dynamics with position dependent coefficients (1) and (2) and show (22) for Lyapunov functions being of the form of polynomials of even but arbitrarily high degree.

3. Langevin Dynamics with Configuration-Dependent Diffusion

We now return to the examples from the introduction, specifically Langevin dynamics with space-dependent friction (1) and (2), and show that in case the variable friction tensor $\mathbf{\Gamma }$ is positive definite and the diffusion tensor $\mathbf{\Sigma }$ has full rank, at all points of the phase space, then the system satisfies the conditions described in the previous section for ergodicity and exponential decay estimates. As illustrations, we consider three applications: (i) a particle model with a temperature gradient; (ii) a simple 2-dimensional Langevin diffusion with a non-conservative force term; (iii) an illustrative Langevin dynamics model for of flocking/swarming, as commonly arises in studies of active particle systems.

3.1. Geometric Ergodicity of Langevin Dynamics With Space-Dependent Coefficients

The main theorem of this section generalizes well known results regarding the ergodicity properties and exponential decay estimates of the associated evolution operator for the underdamped Langevin with constant scalar friction and diffusion coefficient, such as presented in [23,35]. Specifically we extend existing results by

• allowing the systematic force $\mathit{F}$ to be non-conservative,
• explicitly considering the case of $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$ being matrix-valued functions of $\mathit{q}$.

In the current case, no assumption is made regarding a fluctuation-dissipation relation or the form of the invariant density. We will see that under certain conditions on the coefficients $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$ and the non-conservative force $\mathit{F}$, which are detailed below, the proof of the extended ergodicity criterion for Langevin dynamics follows in a straightforward way from the corresponding proof of the constant-coefficient result [35]. These generalization are nonetheless of high practical relevance and allow us to conclude ergodicity for a wide range of relevant modelling applications.

Let, for a square matrix $\mathbf{A}\in {\mathbb{R}}^{n\times n}$,

$$\sigma \left(\mathbf{A}\right):=\{\lambda \in \mathbb{R}:\exists \mathit{v}\in {\mathbb{R}}^n,\mathit{v}\ne 0,\lambda \mathit{v}=\mathbf{A}\mathit{v}\},$$

denote its spectrum. The following assumption on the spectrum of $\mathbf{\Gamma }$ and $\mathbf{\Sigma }{\mathbf{\Sigma }}^T$ ensures the existence of a suitable Lyapunov function in the case these matrices are not constant in $\mathit{q}$.

Assumption 3.

(i) 

The spectrum $\sigma (\mathbf{\Gamma }(\mathit{q}\left)\right)$ is uniformly bounded in $\mathit{q}\in {\Omega }_q$ from above and away from 0, i.e., there are positive constants ${\lambda }_{\min },{\lambda }_{\max }$, such that

$${\lambda }_{\min }\le \underset{\mathit{q}\in {\Omega }_q}{\inf }\min \sigma (\mathbf{\Gamma }\left(\mathit{q}\right)),$$

and

$$\underset{\mathit{q}\in {\Omega }_q}{\sup }\max \sigma (\mathbf{\Gamma }\left(\mathit{q}\right))\le {\lambda }_{\max }.$$

(ii) 

The diffusion matrix has full rank, i.e.,

$$\mathrm{rank}(\mathbf{\Sigma }(\mathit{q}\left)\right)=n,$$

for all $\mathit{q}\in {\Omega }_q$, and the spectrum of $\mathbf{\Sigma }{\mathbf{\Sigma }}^T$ is bounded in the sense that

$$\exists \overline{\sigma }>0:\underset{\mathit{q}\in {\Omega }_q}{\sup }\max \sigma (\mathbf{\Sigma }\left(\mathit{q}\right){\mathbf{\Sigma }}^T\left(\mathit{q}\right))\le \overline{\sigma }.$$

Obviously, Assumption 3 is automatically satisfied for ${\Omega }_q=L{\mathbb{T}}^n$ as long as the coefficients $\mathbf{\Gamma },\mathbf{\Sigma }$ are smooth and $\mathbf{\Gamma }\left(\mathit{q}\right)$ is positive definite and $\mathbf{\Sigma }\left(\mathit{q}\right)$ has full rank at every point $\mathit{q}\in {\Omega }_q$. The next assumption ensures the existence of a suitable Lyapunov function in the case of a non-conservative force. Again, it is trivially satisfied for ${\Omega }_q=L{\mathbb{T}}^n$ and $\mathit{F}\in {\mathcal{C}}^{\infty }({\mathbb{T}}^n,\mathbb{R})$.

Assumption 4.

There exists a potential function $U\in {\mathcal{C}}^2({\Omega }_q,\mathbb{R})$ with the following properties

(i) 

there exists $G\in \mathbb{R}$ such that

$$\forall \mathit{q}\in {\Omega }_q,\langle \mathit{q},\mathit{F}\left(\mathit{q}\right)\rangle \le -\langle \mathit{q},{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)\rangle +G.$$

for all $\mathit{q}\in {\Omega }_q$.

(ii) 

the potential function is bounded from below, i.e., there exists $u_{\min }>-\infty$ such that

$$\forall \mathit{q}\in {\Omega }_q,U\left(\mathit{q}\right)\ge u_{\min }.$$

(iii) 

there exist constants $D,E>0$ and $F\in \mathbb{R}$ such that

$$\forall \mathit{q}\in {\Omega }_q,\langle \mathit{q},{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)\rangle \ge DU\left(\mathit{q}\right)+E{\parallel \mathit{q}\parallel }_2^2+F.$$

(25)

We point out that if $\mathit{F}$ is a conservative force, then Assumption 4 reduces to the same asymptotic growth criteria commonly assumed in the derivation of geometric ergodicity of Langevin dynamics with constant coefficients $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$ on an unbounded configurational domain, ${\Omega }_q={\mathbb{R}}^n$ (See again e.g., [23,35]).

Theorem 3.

In (1) and (2), let the force, the friction and diffusion tensors be smooth functions, i.e., $\mathit{F}\in {\mathcal{C}}^{\infty }({\Omega }_q,{\mathbb{R}}^n)$, $\mathbf{\Gamma }\in {\mathcal{C}}^{\infty }({\Omega }_q,{\mathbb{R}}^{n\times n})$, $\mathbf{\Sigma }\in {\mathcal{C}}^{\infty }({\Omega }_q,{\mathbb{R}}^{n\times m})$ such that Assumptions 3 and 4 hold. There is a unique invariant probability measure

$$\tilde{\mu }(\mathrm{d}\mathit{q},\mathrm{d}\mathit{p}):=\tilde{\rho }(\mathit{q},\mathit{p})\mathrm{d}\mathit{q}\mathrm{d}\mathit{p},$$

with a $(\mathit{F},\mathbf{\Gamma },\mathbf{\Sigma })$-dependent density $\tilde{\rho }\in {\mathcal{C}}^{\infty }({\Omega }_q\times {\mathbb{R}}^n,{\mathbb{R}}_{+})$, and for all $l\in \mathbb{N}$ there are constants ${\kappa }_l,C_l>0$ so that

$${\parallel e^{t\mathcal{L}}\phi -\tilde{\mu }\left[\phi \right]\parallel }_{L_{{\mathcal{K}}_l}^{\infty }}\le C_l{e}^{-t{\kappa }_l}{\parallel \phi -\tilde{\mu }\left[\phi \right]\parallel }_{L_{{\mathcal{K}}_l}^{\infty }},$$

(26)

for all $\phi \in L_{{\mathcal{K}}_l}^{\infty }$ and $t\ge 0$, where

$${\mathcal{K}}_l(\mathit{q},\mathit{p})={(\langle \mathit{p},\mathit{p}\rangle +1)}^l,$$

(27)

in the case of ${\Omega }_q=L{\mathbb{T}}^n$ and

$${\mathcal{K}}_l(\mathit{q},\mathit{p})={\left(cU\left(\mathit{q}\right)+a\langle \mathit{q},\mathit{q}\rangle +2b\langle \mathit{q},\mathit{p}\rangle +\frac{c}{2}\langle \mathit{p},\mathit{p}\rangle \right)}^l,$$

(28)

with suitably chosen positive constants $a,b,c>0$ in the case of ${\Omega }_q={\mathbb{R}}^n$. Furthermore, the probability measure $\tilde{\mu }$ is such, that

$${\int }_{{\Omega }_q\times {\mathbb{R}}^n}{\mathcal{K}}_l\mathrm{d}\tilde{\mu }<\infty .$$

Proof. 

Lyapunov condition

We first show that Assumption 1 holds for ${\mathcal{K}}_l$ as defined in (27) and (28) in the respective setups. Note that

$$\mathcal{L}={\mathcal{L}}_H+{\mathcal{L}}_O,$$

with

$$\begin{array}{cc}\hfill {\mathcal{L}}_H& =\mathit{p}·{\nabla }_{\mathit{q}}-{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)·{\nabla }_{\mathit{p}},\hfill \\ \hfill {\mathcal{L}}_O& =-{\mathit{p}}^T\mathbf{\Gamma }\left(\mathit{q}\right)·{\nabla }_{\mathit{p}}+\frac{1}{2}{\mathbf{\Sigma }}^T\left(\mathit{q}\right)\mathbf{\Sigma }\left(\mathit{q}\right):{\nabla }_{\mathit{q}}^2.\hfill \end{array}$$

We first show the existence of constants $a_l,b_l$ such that Assumption 1 is satisfied for $a=a_l,b=b_l$ in the case $l=1$. For $l\ge 2$, the existence of suitable constants $a_l,b_l$ follows inductively.

By Assumption 4(i), it follows that

$$\begin{array}{cc}\hfill {\mathcal{L}}_H{\mathcal{K}}_1(\mathit{q},\mathit{p})& =2b\langle \mathit{q},\mathit{F}\left(\mathit{q}\right)\rangle +2a\langle \mathit{p},\mathit{q}\rangle +2b\langle \mathit{p},\mathit{q}\rangle \hfill \\ & \le -2b\langle \mathit{q},{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)\rangle +G+2a\langle \mathit{p},\mathit{q}\rangle +2b\langle \mathit{p},\mathit{q}\rangle .\hfill \end{array}$$

Similarly, by Assumption 3(ii), we find

$$\begin{array}{cc}\hfill {\mathcal{L}}_O{\mathcal{K}}_1(\mathit{q},\mathit{p})& =-2b\langle \mathit{q},\mathbf{\Gamma }\left(\mathit{q}\right)\mathit{p}\rangle -c\langle \mathit{p},\mathbf{\Gamma }\left(\mathit{q}\right)\mathit{p}\rangle +c\mathrm{tr}\left({\mathbf{\Sigma }}^T\left(\mathit{q}\right)\mathbf{\Sigma }\left(\mathit{q}\right)\right)\hfill \\ & \le -2b\langle \mathit{q},\mathbf{\Gamma }\left(\mathit{q}\right)\mathit{p}\langle -c\langle \mathit{p},\mathbf{\Gamma }\left(\mathit{q}\right)\mathit{p}\rangle +cn{\lambda }_{\max }.\hfill \end{array}$$

Putting both inequalities together thus yields

$$\begin{array}{cc}\hfill \mathcal{L}{\mathcal{K}}_1(\mathit{q},\mathit{p})& \le -2b\left(DU\left(\mathit{q}\right)+{E\parallel \mathit{q}\parallel }_2^2+F\right)+2(a+b)\langle \mathit{p},\mathit{q}\rangle -2b\langle \mathit{q},\mathbf{\Gamma }\left(\mathit{q}\right)\mathit{p}\rangle -c\langle \mathit{p},\mathbf{\Gamma }\left(\mathit{q}\right)\mathit{p}\rangle +cn{\lambda }_{\max }\hfill \\ & =-2bDU\left(\mathit{q}\right)-{\left(\begin{array}{c}\mathit{q}\\ \mathit{p}\end{array}\right)}^T\mathcal{M}\left(\begin{array}{c}\mathit{q}\\ \mathit{p}\end{array}\right)-2bF+cn{\lambda }_{\max },\hfill \end{array}$$

with

$$\mathcal{M}=\left(\begin{array}{cc}Eb{\mathit{I}}_n& b\mathbf{\Gamma }\left(\mathit{q}\right)-a{\mathit{I}}_n\\ b\mathbf{\Gamma }\left(\mathit{q}\right)-a{\mathit{I}}_n& -2b{\mathit{I}}_n+c\mathbf{\Gamma }\left(\mathit{q}\right)\end{array}\right),$$

where we used Assumption 4(iii) in the inequality. The existence of suitable constants a and b such that Assumption 1 is satisfied for the case $\mathcal{K}={\mathcal{K}}_1$, follows directly if the matrix $\mathcal{M}$ is positive definite. The positive definiteness of the block matrix $\mathcal{M}$ is implied if (See e.g., [37]), the matrices $-2b{\mathit{I}}_n+c\mathbf{\Gamma }\left(\mathit{q}\right)$ and the Schur complement

$$Eb{\mathit{I}}_n-{(b\mathbf{\Gamma }\left(\mathit{q}\right)-a{\mathit{I}}_n)}^T{\left(-2b{\mathit{I}}_n+c\mathbf{\Gamma }\left(\mathit{q}\right)\right)}^{-1}(b\mathbf{\Gamma }\left(\mathit{q}\right)-a{\mathit{I}}_n)$$

are both positive definite. Indeed, since the spectrum of $\mathbf{\Gamma }$ is uniformly bounded on ${\Omega }_q={\mathbb{R}}^n$ according to Assumption 3(i) the positive definiteness of both these matrices can be ensured by choosing a, b sufficiently small and c sufficiently large.

Minorization condition

A simple calculation shows that $\mathbf{\Sigma }\left(\mathit{q}\right)$ having rank n for all $\mathit{q}\in {\Omega }_q$ immediately implies that the SDE (1) and (2) satisfies the parabolic Hörmander condition. It therefore only remains to show that for arbitrary $T>0$ and $({\mathit{q}}^{-},{\mathit{p}}^{-}),({\mathit{q}}^{+},{\mathit{p}}^{+})\in {\mathbb{R}}^{2n}$, the control problem

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\mathit{q}}{\mathrm{d}t}& =\mathit{p},\hfill \\ \hfill \frac{\mathrm{d}\mathit{p}}{\mathrm{d}t}& =\mathit{F}\left(\mathit{q}\right)-\mathbf{\Gamma }\left(\mathit{q}\right)\mathit{p}+\mathbf{\Sigma }\left(\mathit{q}\right)\mathbf{u},\hfill \\ \hfill \mathrm{subject}\mathrm{to}\left(\mathit{q}\right(0),\mathit{p}(0\left)\right)& =({\mathit{q}}^{-},{\mathit{p}}^{-}),(\mathit{q}\left(T\right),\mathit{p}\left(T\right))=({\mathit{q}}^{+},{\mathit{p}}^{+}),\hfill \end{array}$$

(29)

has a continuous solution $u\in L^1([0,T],{\mathbb{R}}^m)$. It is easy to verify that there exists a smooth path $\tilde{\mathit{q}}\in {\mathcal{C}}^2([0,T],{\mathbb{R}}^n)$ such that

$$(\tilde{\mathit{q}}\left(0\right),\dot{\tilde{\mathit{q}}}\left(0\right))=({\mathit{q}}^{-},{\mathit{p}}^{-}),(\tilde{\mathit{q}}\left(T\right),\dot{\tilde{\mathit{q}}}\left(T\right))=({\mathit{q}}^{+},{\mathit{p}}^{+}).$$

Rewrite (29) as a second order differential equation:

$$\ddot{\mathit{q}}=-{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)-\mathbf{\Gamma }\left(\mathit{q}\right)\dot{\mathit{q}}+\mathbf{\Sigma }\left(\mathit{q}\right)\mathbf{u}.$$

Since $\mathrm{rank}(\mathbf{\Sigma }(\mathit{q}\left)\right)=n$ for all $\mathit{q}\in {\Omega }_q$, there exists ${\mathbf{\Sigma }}^g\in {\mathcal{C}}^{\infty }({\Omega }_q,{\mathbb{R}}^{m\times n})$ such that

$$\mathbf{\Sigma }\left(q\right){\mathbf{\Sigma }}^g\left(\mathit{q}\right)\mathit{p}=\mathit{p},$$

for all $\mathit{q},\mathit{p}\in {\mathbb{R}}^n$, thus,

$$\mathit{u}\left(t\right)={\mathbf{\Sigma }}^g\left(\tilde{\mathit{q}}\left(t\right)\right)\left[\ddot{\tilde{\mathit{q}}}+{\nabla }_{\mathit{q}}U\left(\tilde{\mathit{q}}\left(t\right)\right)+\mathbf{\Gamma }\left(\tilde{\mathit{q}}\left(t\right)\right)\dot{\tilde{\mathit{q}}}\left(t\right)\right]$$

(30)

is a solution of (29). ☐

Note that the case where $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$ are constant was already studied in Mattingly et al. [35] and the proof of the above theorem resembles the structure of the proof therein.

In what follows we provide three examples of models whose ergodic properties can be characterised by the above Theorem 3.

3.2. Single-Particle System with Non-Conservative Force

As an example of a system with a non-conservative force term, which satisfies the condition of the above theorem we consider a Langevin equation of the form (1) and (2) where the force term $\mathit{F}$ is of the form

$$\mathit{F}\left(\mathit{q}\right)=-{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)+\alpha \mathit{J}\mathit{q},$$

(31)

where $\mathit{J}\ne \mathbf{0}$ is a skew-symmetric matrix and $\alpha \in \mathbb{R}$. For $\alpha \ne 0$ the additional term $\alpha \mathit{J}\mathit{q}$ obviously does not correspond to the negative gradient of a smooth potential energy function, thus in this case the force (31) is indeed non-conservative. It is easy to verify that U satisfying the growth conditions (ii) and (iii) in Assumption 4 and

$$U\left(\mathit{q}\right)\in O\left({\parallel \mathit{q}\parallel }^{2+ϵ}\right),$$

with $ϵ>0$ as $\parallel \mathit{q}\parallel \to \infty$ implies (i) in Assumption 4. Therefore, as long as the remaining conditions on the coefficients $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$ in Theorem 3 are satisfied, it follows from the same theorem that the respective (non-equilibrium) dynamics possesses an invariant measure, ${\mu }_{\alpha }(\mathrm{d}\mathit{q},\mathrm{d}\mathit{p})={\rho }_{\alpha }(\mathit{q},\mathit{p})\mathrm{d}\mathit{q}\mathrm{d}\mathit{p}$, to which it converges exponentially fast in $L_{\mathcal{K}}^{\infty }$ as specified above.

3.3. Multi-Particle Systems

In the remainder of this section we consider the application of Theorem 3 to two different types of particle systems, which can be seen as instances of the underdamped Langevin equation with non-constant coefficients. With some abuse of notation whenever a particle model is considered, let ${\mathit{q}}_i\in L{\mathbb{T}}^d,i=1,\dots ,N$ where $L>0$ and ${\mathit{p}}_i\in {\mathbb{R}}^d,i=1,\dots ,N$ denote the position and momentum vectors of the particles, respectively. If, on the other hand, we want to refer to the i-th entry of the vector $\mathit{v}$ where $\mathit{v}=\mathit{q}$ or $\mathit{v}=\mathit{p}$, we write ${\Pi }_i\mathit{v}$, where ${\Pi }_i:{\mathbb{R}}^n\to \mathbb{R}$ with $i\in \{1,2,\dots ,n\},n=Nd$, denotes the operator which selects the i-th cartesian coordinate, i.e., ${\Pi }_i\mathit{x}={\mathbf{e}}_i·\mathit{x}$. Furthermore, if not stated otherwise, we will assume, that the force is conservative and corresponds to the gradient of a potential function $U\left(\mathit{q}\right)$, which is composed as the sum of smooth pair potentials ${\tilde{U}}_{i,j}\in {\mathcal{C}}^{\infty }(\mathbb{R},\mathbb{R}),1\le j<i\le N$., i.e.,

$$U\left(\mathit{q}\right)=\sum _{i=1}^N\sum _{j=1}^{i-1}{\tilde{U}}_{i,j}\left(\parallel {\mathit{q}}_i-{\mathit{q}}_j\parallel \right).$$

(32)

3.4. Particle System with Temperature Gradient

As a first application we consider an N-particle system with periodic boundary condition where the particles are coupled to a heat bath whose temperature varies depending on the position of the particle within the periodic simulation box. The system we consider in the following is of the form

$$\begin{array}{cc}\hfill {\dot{\mathit{q}}}_i& =-{\mathit{p}}_i,\hfill \\ \hfill {\dot{\mathit{p}}}_i& =-{\nabla }_{{\mathit{q}}_i}U\left(\mathit{q}\right)-\gamma {\mathit{p}}_i+\sqrt{2\gamma {\left(\beta \left({\mathit{q}}_i\right)\right)}^{-1}}{\dot{W}}_i,\hfill \end{array}$$

(33)

for $i=1,\dots ,N$, where $W_i,i=1,\dots ,N$, are independent white noise processes in ${\mathbb{R}}^d$. We further assume that the friction coefficient $\gamma >0$ is a positive constant and $\beta$ a smooth positive function on $L{\mathbb{T}}^d$, i.e., $\beta \in {\mathcal{C}}^{\infty }(L{\mathbb{T}}^d,{\mathbb{R}}_{+})$. In light of (29) this corresponds to a constant diagonal friction matrix

$$\begin{array}{cc}\hfill \mathbf{\Gamma }\left(\mathit{q}\right)& =\gamma {\mathit{I}}_n,\hfill \end{array}$$

and diffusion tensor of the form

$$\mathbf{\Sigma }\left(\mathit{q}\right)=\sqrt{2\gamma }{\mathit{\beta }}^{-\frac{1}{2}}\left(\mathit{q}\right)\otimes {\mathit{I}}_d,$$

where

$${\mathit{\beta }}^{-\frac{1}{2}}\left(\mathit{q}\right):=\mathrm{diag}{\left(\sqrt{\beta \left({\mathit{q}}_1\right)},\dots ,\sqrt{\beta \left({\mathit{q}}_N\right)}\right)}^{-1}.$$

The matrices $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$ are clearly positive definite and invertible, hence by Theorem 3, there exists a unique invariant measure ${\mu }_{\gamma ,\beta }$ to which the solution of (33) converges (exponentially fast) in law. Note that, even for this relatively simple generalization of the standard underdamped Langevin equation, the form of the invariant measure ${\mu }_{\gamma ,\beta }$ depends nontrivially on $\gamma$ and $\beta$, and one can in general not find an analytic solution of the corresponding stationary Fokker-Planck equation. As long as ${\mathit{\beta }}^{-1}$ is bounded from above and away from 0, and the potential U is modified so that Assumption 4 is satisfied (e.g., by adding a confining potential to U), it is easy to see that Theorem 3 applies also for the case ${\Omega }_q={\mathbb{R}}^n$.

3.5. Stochastic Cucker-Smale Model

As a second application of Theorem 3, we consider a stochastic model for flocking which is a variant of the model presented in [4]. In fact, the primary model considered in that paper replaced the deterministic Cucker-Smale system [38], in which a collection of active particles interact with a configuration-dependent friction, by one in which the particles were additionally perturbed by bounded noise. An SDE model was presented there without analysis. A subsequent paper [5] provided numerical evidence for flocking states. In these papers, the SDE approach consists of a Langevin-type model with coordinate-dependent friction and additive (typically uniform) white noise, that is they take the form

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}\mathit{q}}{\mathrm{d}t}& =& {\mathit{M}}^{-1}\mathit{p},\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}\mathit{p}}{\mathrm{d}t}& =& -{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)-\mathbf{\Gamma }\left(\mathit{q}\right)\mathit{p}+\sigma {\mathit{I}}_{Nd}\dot{\mathit{W}}\left(t\right).\hfill \end{array}$$

(35)

For a system of particles moving in one dimension, one has a friction matrix $\mathbf{\Gamma }\left(\mathit{q}\right)={\left[{\gamma }_{ij}\left(\mathit{q}\right)\right]}_{1\le i,j\le N}\otimes {\mathbf{I}}_d$, where

$${\gamma }_{ij}\left(\mathit{q}\right)=-\psi (\parallel {\mathit{q}}_j-{\mathit{q}}_i\parallel ),i\ne j$$

(36)

and

$${\gamma }_{ii}=-\sum _{j\ne i}{\gamma }_{ij},$$

(37)

for some given scalar kernel function $\psi \left(r\right)\ge 0$. That is, the matrix $\mathbf{\Gamma }\left(\mathit{q}\right)$ is a (weighted) graph Laplacian which reflects the interaction structure of the problem. If $\psi \left(r\right)>0$, this interaction graph is complete.

The inclusion of conservative forces derived from a potential energy function U is an addition to the models mentioned above and allows to incorporate direct attraction and repulsion effects. Note that since ${\left[{\gamma }_{ij}\right]}_{1\le i,j\le N}$ is of the form of a graph Laplacian, it is positive semi-definite with at least one eigenvalue being singular. If we assume that $\psi \left(r\right)>0$, then it follows that ${\left[{\gamma }_{ij}\right]}_{1\le i,j\le N}$ possesses exactly one singular eigenvector, $1_N\in {\mathbb{R}}^N$, i.e., the vector whose entries all are equal to 1. Consequently, $\mathbf{\Gamma }\left(\mathit{q}\right)$, has exactly d singular eigenvalues, each of them being of the form ${\mathbf{u}}_l=1_N\otimes {\mathit{e}}_l$, where ${\mathit{e}}_l$ denotes the l-th canonical vector in ${\mathbb{R}}^d$. These singular eigenvectors can be associated with the mean momentum

$$\overline{\mathit{p}}=N^{-1}\sum _{i=1}^N{\mathit{p}}_i,$$

of the collection of particles via the relation

$${\overline{\mathit{p}}}_l=N^{-1}{\mathbf{u}}_l^T\mathit{p},$$

for $l\in \{1,\dots ,d\}$. There are several issues regarding the model (34) and (35). Most importantly, while the diffusion matrix in (35) has full rank, the friction tensor $\mathbf{\Gamma }\left(\mathit{q}\right)$ is only of rank $(N-1)d$, which means that in directions of the singular vectors ${\mathbf{u}}_l,1\le l\le d$ there is no dissipation and therefore the kinetic energy of the system will be unbounded as $t\to \infty$, and hence one cannot expect (34) and (35) to possess an invariant measure. More specifically, if U is composed of pair potentials one can show that the mean momentum $\overline{\mathit{p}}$ is a Brownian motion,

$$\dot{\overline{\mathit{p}}}={\left[{\mathbf{u}}_l\right]}_{1\le l\le d}^T\dot{\mathit{W}}.$$

3.5.1. Regularized Stochastic Cucker-Smale Dynamics

A simple fix to the model (34) and (35), which ensures ergodicity of the dynamics, is to add additional dissipation, which is uniform in each component of $\mathit{p}$. This can be achieved by replacing the friction tensor in (35) by ${\mathbf{\Gamma }}_{ϵ}\left(\mathit{q}\right)$, which for $ϵ>0$ is defined such that

$${\mathbf{\Gamma }}_{ϵ}\left(\mathit{q}\right)=\mathbf{\Gamma }\left(\mathit{q}\right)+ϵ{\mathbf{I}}_{Nd}.$$

It follows directly from the Gershgorin circle theorem that ${\mathbf{\Gamma }}_{ϵ}\left(\mathit{q}\right)$ is positive definite for all $\mathit{q}\in {\Omega }_q$ and any choice of $ϵ>0$.

3.5.2. Modified Stochastic Cucker-Smale Dynamics

While the above described regularized stochastic Cucker-Smale dynamics is a valid extension of the original model which ensures that the conditions of Theorem 3 are satisfied and hence geometric ergodicity for (34) and (35) holds, the form of the corresponding invariant measure does depend in a non-trivial way on $\mathbf{\Gamma },\sigma$ and $ϵ$ and unless $\mathbf{\Gamma }$ is constant in $\mathit{q}$, one cannot easily find a closed form of the invariant measure. We therefore propose another novel modification of (34) and (35), which is geometrically ergodic with an invariant measure of closed form. We construct this model as a superposition, i.e., $\mathit{x}={\mathit{x}}_{\perp }+{\mathit{x}}_{\left|\right|}$, of two independent stochastic processes, ${\mathit{x}}_{\perp }=({\mathit{q}}_{\perp },{\mathit{p}}_{\perp })$ and ${\mathit{x}}_{\left|\right|}=({\mathit{q}}_{\left|\right|},{\mathit{p}}_{\left|\right|})$ taking both values in ${\Omega }_x\times {\mathbb{R}}^n$, respectively. We construct these processes such that the parametrization of the former process determines the statistical properties of the stochastic inter-particle interactions and the parametrization of the latter process the collective motion of the flock, i.e., the diffusive properties of the centre of mass. More specifically, we construct ${\mathit{x}}_{\perp }$ as the solution of the SDE

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathit{q}}_{\perp }}{\mathrm{d}t}& ={\mathit{p}}_{\perp },\hfill \\ \hfill \frac{\mathrm{d}{\mathit{p}}_{\perp }}{\mathrm{d}t}& =-{\nabla }_{\mathit{q}}U\left({\mathit{q}}_{\perp }\right)-{\gamma }_{\perp }{\mathbf{\Gamma }}_{\perp }\left({\mathit{q}}_{\perp }\right){\mathit{p}}_{\perp }+\sqrt{2{\gamma }_{\perp }{T}_{\perp }}{\mathbf{\Sigma }}_{\perp }\left({\mathit{q}}_{\perp }\right)\dot{\mathit{W}}\left(t\right).\hfill \end{array}$$

(38)

and ${\mathit{x}}_{\left|\right|}$ as the solution of

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathit{q}}_{\left|\right|}}{\mathrm{d}t}& ={\mathit{p}}_{\left|\right|},\hfill \\ \hfill \frac{\mathrm{d}{\mathit{p}}_{\left|\right|}}{\mathrm{d}t}& =-{\gamma }_{\left|\right|}{\mathbf{\Gamma }}_{\left|\right|}{\mathit{p}}_{\left|\right|}+\sqrt{2{\gamma }_{\left|\right|}{T}_{\left|\right|}}{\mathbf{\Sigma }}_{\left|\right|}\dot{\mathit{W}}\left(t\right).\hfill \end{array}$$

(39)

where ${\gamma }_{\perp },T_{\perp },{\gamma }_{\left|\right|},T_{\left|\right|}\ge 0$ and $\mathit{W}$ as specified in (3) and (4). We refer to (38) as the equation of the peculiar dynamics and to (39) as the equation of the consensus dynamics. This naming is motivated by the following choices of the respective friction tensor and diffusion matrix.

• Peculiar dynamics:
  As in (34) and (35), we choose

  $${\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)={\left[{\gamma }_{ij}\left(\mathit{q}\right)\right]}_{1\le i,j\le N}\otimes {\mathbf{I}}_d,$$

  (40)

  where ${\gamma }_{ij}$ are defined as above in (36) and (37). By choosing the diffusion matrix ${\mathbf{\Sigma }}_{\perp }$ such that

  $${\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)={\mathbf{\Sigma }}_{\perp }\left(\mathit{q}\right){\mathbf{\Sigma }}_{\perp }^T\left(\mathit{q}\right),$$

  we ensure that the total momentum in each dimension of the physical domain remains constant, i.e.,

  $$\frac{\mathrm{d}}{\mathrm{d}t}\sum _{j=1}{{\mathit{p}}_{\perp }}_j=\mathbf{0}.$$

  (41)
  This follows directly from the fact, that

  $${\mathbf{u}}_l\in \ker \left({\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)\right)=\ker \left({\mathbf{\Sigma }}_{\perp }\left(\mathit{q}\right)\right),l\in \{1,\dots ,d\}.$$

  as well as

  $${\mathbf{u}}_l·{\nabla }_{\mathit{q}}U\left({\mathit{q}}_{\perp }\right)=0,$$

  (42)

  (since U is composed of pair potentials) and therefore by (38),

  $${\mathbf{u}}_l·\frac{\mathrm{d}}{\mathrm{d}t}{\mathit{p}}_{\perp }=0,l\in \{0,\dots ,d\},$$

  hence

  $$\frac{\mathrm{d}}{\mathrm{d}t}\sum _{j=1}^N{{\mathit{p}}_{\perp }}_j=\sum _{l=1}^d{\mathit{e}}_l{\mathbf{u}}_l·\frac{\mathrm{d}}{\mathrm{d}t}{\mathit{p}}_{\perp }=\mathbf{0}.$$

  (43)
• Consensus dynamics:
  We construct the matrix ${\mathbf{\Gamma }}_{\left|\right|}$ such that the difference in the momenta of all particle pairs remains constant under the dynamics (39), i.e.,

  $$\frac{\mathrm{d}}{\mathrm{d}t}\left({{\mathit{p}}_{\left|\right|}}_j-{{\mathit{p}}_{\left|\right|}}_i\right)=\mathbf{0},$$

  for all $1\le i,j\le N$. We achieve that by choosing ${\mathbf{\Gamma }}_{\left|\right|}$ as

  $${\mathbf{\Gamma }}_{\left|\right|}={\mathbf{I}}_N\otimes {\mathbf{I}}_d={\left[(i+j)\mathrm{mod}d\right]}_{1\le i,j\le Nd},$$

  (44)

  and ${\mathbf{\Sigma }}_{\left|\right|}$ such that

  $${\mathbf{\Gamma }}_{\left|\right|}={\mathbf{\Sigma }}_{\left|\right|}{\mathbf{\Sigma }}_{\left|\right|}^T,$$

  i.e.,

  $${\mathbf{\Sigma }}_{\left|\right|}=N^{-1/2}{\mathbf{I}}_N\otimes {\mathbf{I}}_d.$$
• Combined dynamics:
  We first observe that although the processes ${\mathit{x}}_{\perp }$ and ${\mathit{x}}_{\left|\right|}$ are driven by the same Wiener process $\mathit{W}$ they are indeed independent. This follows since the column vectors of ${\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)$ are orthogonal to the column vectors of ${\mathbf{\Gamma }}_{\left|\right|}$ in the sense that

  $$\ker \left({\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)\right)=\mathrm{span}\left({\mathbf{\Gamma }}_{\left|\right|}\right),$$

  i.e.,

  $${\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right){\mathbf{\Gamma }}_{\left|\right|}^T={\mathbf{\Gamma }}_{\left|\right|}{\mathbf{\Gamma }}_{\perp }^T\left(\mathit{q}\right)=\mathbf{0},$$

  (45)

  for all $\mathit{q}\in {\Omega }_q$, and hence also

  $${\mathbf{\Sigma }}_{\perp }\left(\mathit{q}\right){\mathbf{\Sigma }}_{\left|\right|}^T={\mathbf{\Sigma }}_{\left|\right|}{\mathbf{\Sigma }}_{\perp }^T\left(\mathit{q}\right)=\mathbf{0},$$

  (46)

  so that ${\mathbf{\Sigma }}_{\perp }\left(\mathit{q}\right)\mathit{W}$ and ${\mathbf{\Sigma }}_{\left|\right|}\mathit{W}$ are independent processes, which implies that also the solution processes of the respective SDEs are independent. Moreover, since U is composed of pair-potentials, we have

  $$U({\mathit{q}}_{\perp }+{\mathit{q}}_{\left|\right|})=U\left({\mathit{q}}_{\perp }\right).$$

  (47)
  Using (45)–(47) it directly follows that $\mathit{x}={\mathit{x}}_{\perp }+{\mathit{x}}_{\left|\right|}$ can be identified with the solution of (1) and (2) with

  $$\mathbf{\Gamma }\left(\mathit{q}\right)={\gamma }_{\perp }{\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)+{\gamma }_{\left|\right|}{\mathbf{\Gamma }}_{\left|\right|},{\gamma }_{\perp }\ge 0,{\gamma }_{\left|\right|}\ge 0,$$

  (48)

  with

  $${\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)={\left[{\gamma }_{ij}\left(\mathit{q}\right)\right]}_{1\le i,j\le N}\otimes {\mathbf{I}}_d,$$

  (49)

  and

  $$\mathbf{\Sigma }\left(\mathit{q}\right)=\sqrt{2{\gamma }_{\perp }{T}_{\perp }}{\mathbf{\Sigma }}_{\perp }\left(\mathit{q}\right)+\sqrt{2{\gamma }_{\left|\right|}{T}_{\left|\right|}}{\mathbf{\Sigma }}_{\left|\right|},T_{\perp }>0,T_{\left|\right|}>0.$$

  (50)

  or more explicitly with the solution of

  $$\begin{array}{cc}\hfill \frac{\mathrm{d}\mathit{q}}{\mathrm{d}t}& =\mathit{p},\hfill \\ \hfill \frac{\mathrm{d}\mathit{p}}{\mathrm{d}t}& =-{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)-{\gamma }_{\perp }{\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)\mathit{p}-{\gamma }_{\left|\right|}{\mathbf{\Gamma }}_{\left|\right|}\mathit{p}+\left[\sqrt{2{\gamma }_{\perp }{T}_{\perp }}{\mathbf{\Sigma }}_{\perp }\left(\mathit{q}\right)+\sqrt{2{\gamma }_{\left|\right|}{T}_{\left|\right|}}{\mathbf{\Sigma }}_{\left|\right|}\right]\dot{\mathit{W}}\left(t\right).\hfill \end{array}$$

  (51)

Remark 1.

We note that the above choice of ${\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)$ and ${\mathbf{\Sigma }}_{\perp }\left(\mathit{q}\right)$ is very similar to the friction tensor and diffusion tensor in dissipative particle dynamics. In fact (38) would exactly correspond to a DPD system, if instead of (40) one constructs the friction tensor such that dissipation is aligned with the relative orientation of particle pairs

$${\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)={\left[{\gamma }_{ij}\left(\mathit{q}\right){\widehat{\mathit{q}}}_{i,j}\otimes {\widehat{\mathit{q}}}_{i,j}\right]}_{1\le i,j\le N}.$$

with

$${\widehat{\mathit{q}}}_{i,j}=({\mathit{q}}_j-{\mathit{q}}_i)/\parallel {\mathit{q}}_j-{\mathit{q}}_i\parallel .$$

Proposition 1.

Let $\mathbf{\Gamma }\left(\mathit{q}\right)$ and $\mathbf{\Sigma }\left(\mathit{q}\right)$ defined as above in (49) and (50), with $T_{\perp }>0,T_{\left|\right|}>0$. Let further ${\Omega }_q={\mathbb{T}}^{N\times d}$ and U be of the form (32), then the SDE (51) possesses an invariant measure ${\mu }_{T_{\perp },T_{\left|\right|}}\left(\mathrm{d}\mathit{q}\mathrm{d}\mathit{p}\right)={\rho }_{T_{\perp },T_{\left|\right|}}(\mathit{q},\mathit{p})\mathrm{d}\mathit{q}\mathrm{d}\mathit{p}$ of the form

$${\rho }_{T_{\perp },T_{\left|\right|}}(\mathit{q},\mathit{p})=\frac{1}{Z}{e}^{-\left[T_{\perp }^{-1}U\left(\mathit{q}\right)+\frac{1}{2}{\mathit{p}}^T{\mathbf{C}}^{-1}\mathit{p}\right]}$$

(52)

with covariance matrix

$$\mathbf{C}=\frac{1}{N}\left[T_{\left|\right|}{\mathbf{1}}_N+T_{\perp }\mathit{L}\right]\otimes {\mathbf{I}}_d$$

where ${\mathbf{1}}_N\in {\mathbb{R}}^{N\times N}$ is the all-ones matrix, i.e., every entry of ${\mathbf{1}}_N$ is equal to 1. The matrix $\mathit{L}\in {\mathbb{R}}^{N\times N}$ with

$${\mathit{L}}_{i,j}=\left\{\begin{array}{cc}N-1,\hfill & i=j,\hfill \\ -1,\hfill & i\ne j\hfill \end{array}\right.$$

denotes the graph Laplacian of a fully connected graph. In particular,

$${\mathbb{E}}_{\mu }\left[({\mathit{p}}_j-{\mathit{p}}_i)\otimes ({\mathit{p}}_j-{\mathit{p}}_i)\right]=2T_{\perp }{\mathbf{I}}_d$$

(53)

where $1\le i<j\le N$, and for $\overline{\mathit{p}}:=N^{-1}{\sum }_{i=1}^N{\mathit{p}}_i\in {\mathbb{R}}^d$,

$${\mathbb{E}}_{\mu }\left[\overline{\mathit{p}}\otimes \overline{\mathit{p}}\right]=T_{\left|\right|}{\mathbf{I}}_d$$

(54)

where $1\le j\le d$.

Proof. 

We show that (52) is a stationary solution of the Fokker-Planck equation associated with the SDE (51), i.e.,

$${\mathcal{L}}^{†}{\rho }_{T_{\perp },T_{\left|\right|}}=0,$$

(55)

with

$${\mathcal{L}}^{†}={\mathcal{L}}_H^{†}+{\gamma }_{\perp }{\mathcal{L}}_{O\perp }^{†}+{\gamma }_{\left|\right|}{\mathcal{L}}_{O\left|\right|},$$

where the action of each of these operators applied to $\rho \in {\mathcal{C}}^2({\Omega }_q,\mathbb{R})$ is given as

$$\begin{array}{cc}\hfill {\mathcal{L}}_H^{†}\rho & =\left[-\mathit{p}·{\nabla }_{\mathit{q}}+{\nabla }_{\mathit{p}}·{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)\right]\rho ,\hfill \\ \hfill {\mathcal{L}}_{O\perp }^{†}\rho & ={\nabla }_{\mathit{p}}·\left({\mathbf{\Gamma }}_{\perp }\mathit{p}\rho \right)+{\nabla }_{\mathit{p}}^2:T_{\perp }{\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)\rho ,\hfill \\ \hfill {\mathcal{L}}_{O\left|\right|}^{†}\rho & ={\nabla }_{\mathit{p}}·\left({\mathbf{\Gamma }}_{\left|\right|}\mathit{p}\rho \right)+{\nabla }_{\mathit{p}}^2:T_{\left|\right|}{\mathbf{\Gamma }}_{\left|\right|}\rho .\hfill \end{array}$$

Before we show (55), we first note that since

$$\ker \left({\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)\right)=\mathrm{span}\left({\mathbf{\Gamma }}_{\left|\right|}\right),$$

for all $\mathit{q}\in {\Omega }_q$, we have

$${\mathbf{\Sigma }}_{\perp }\left(\mathit{q}\right){\mathbf{\Sigma }}_{\left|\right|}^T={\mathbf{\Sigma }}_{\left|\right|}{\mathbf{\Sigma }}_{\perp }^T\left(\mathit{q}\right)=\mathbf{0},$$

hence

$$\begin{array}{cc}\hfill \mathbf{\Sigma }\left(\mathit{q}\right){\mathbf{\Sigma }}^T\left(\mathit{q}\right)& ={\gamma }_{\perp }{T}_{\perp }{\mathbf{\Sigma }}_{\perp }\left(\mathit{q}\right){\mathbf{\Sigma }}_{\perp }^T\left(\mathit{q}\right)+{\gamma }_{\left|\right|}{T}_{\left|\right|}{\mathbf{\Sigma }}_{\left|\right|}{\mathbf{\Sigma }}_{\left|\right|}^T\hfill \\ & ={\gamma }_{\perp }{T}_{\perp }{\mathbf{\Gamma }}_{\perp }\left(\mathit{q}\right)+{\gamma }_{\left|\right|}{T}_{\left|\right|}{\mathbf{\Gamma }}_{\left|\right|}.\hfill \end{array}$$

(56)

Furthermore,

$${\mathbf{C}}^{-1}=N^{-1}\left[T_{\left|\right|}^{-1}{\mathbf{1}}_N+T_{\perp }^{-1}\mathit{L}\right]\otimes {\mathbf{I}}_d,$$

which implies

(i)

${\mathbf{C}}^{-1}\mathit{p}=T_{\left|\right|}^{-1}\mathit{p},$ for $\mathit{p}\in \mathrm{span}\left({\mathbf{\Gamma }}_{\left|\right|}\right)$,

(ii)

${\mathbf{C}}^{-1}\mathit{p}=T_{\perp }^{-1}\mathit{p},$ for $\mathit{p}\in \mathrm{span}\left({\mathbf{\Gamma }}_{\perp }\right)$.

Finally, since $U\left(\mathit{q}\right)$ is composed of pair potential functions, we have

$$U\left(\mathit{q}\right)=U([c{\mathbf{1}}_N+N^{-1}\mathit{L}]\otimes {\mathbf{I}}_d\mathit{q}),$$

for all $c\in \mathbb{R}$, hence in particular

$$\begin{array}{cc}\hfill T_{\perp }^{-1}{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)& =T_{\perp }^{-1}{\nabla }_{\mathit{q}}U\left(\left[\frac{T_{\perp }}{T_{\left|\right|}}{\mathbf{1}}_N+N^{-1}\mathit{L}\right]\otimes {\mathbf{I}}_d\mathit{q}\right)\hfill \end{array}$$

(57)

Given the above identities, we conclude

$$\begin{array}{ccc}\hfill {\mathcal{L}}_H{\rho }_{T_{\perp },T_{\left|\right|}}& =\left[T_{\perp }^{-1}\mathit{p}·{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)-{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)·{\mathbf{C}}^{-1}\mathit{p}\right]{\rho }_{T_{\perp },T_{\left|\right|}}\hfill & \hfill =0\end{array}$$

due to (57), and

$$\begin{array}{cc}\hfill {\mathcal{L}}_{O\perp }^{†}{\rho }_{T_{\perp },T_{\left|\right|}}& =0,\hfill \\ \hfill {\mathcal{L}}_{O\left|\right|}^{†}{\rho }_{T_{\perp },T_{\left|\right|}}& =0,\hfill \end{array}$$

due to (ii) and (i), respectively. ☐

Corollary 1.

Let $\psi \left(r\right)>0$ in the definition (36) of ${\gamma }_{ij}$, furthermore ${\gamma }_{\perp }>0$ and ${\gamma }_{\left|\right|}>0$. The invariant measure ${\mu }_{T_{\perp },T_{\left|\right|}}$ specified in Proposition 1 is unique and the law of (51) converges exponentially fast towards ${\mu }_{T_{\perp },T_{\left|\right|}}$ in the sense of (26) with ${\mathcal{K}}_l,l\in \mathbb{N}$ as constructed in the proof of Theorem 3.

4. Numerical Discretization

We next describe the construction of numerical methods for the Langevin system. The discretization schemes described here are based on a general splitting framework as explained in the introduction, adapted for the variable coefficient structure.

Following [15], we break the Langevin system (1) and (2) into three parts: A and B corresponding to the deterministic flow:

$${\mathcal{L}}_{\mathrm{A}}f=\mathit{p}·{\mathit{M}}^{-1}{\nabla }_{\mathit{q}}f$$

and

$${\mathcal{L}}_{\mathrm{B}}f=\mathit{F}·{\nabla }_{\mathit{p}}f,$$

and O, which is associated to the isolated momenta diffusion process defined by an Ornstein-Uhlenbeck type equation, i.e., the stochastic system

$$\dot{\mathit{q}}=0,\dot{\mathit{p}}=-\mathbf{\Gamma }\left(\mathit{q}\right)\mathit{p}+\mathbf{\Sigma }\left(\mathit{q}\right)\dot{\mathit{W}}\left(t\right).$$

(58)

A stochastic splitting method is then obtained by concatenating the (stochastic) flow maps corresponding to the B, A and O part. Although other decompositions can be used for splitting, experience in the constant coefficient case has shown that it is beneficial to maintain the same balance of noise and dissipation of the original equations by keeping these terms together. Among stochastic splitting schemes based on this decomposition, the symmetric integration sequence BAOAB was observed to yield a lower discretization bias for ergodic averages in comparison to other integration sequences requiring only one evaluation of the gradient ${\nabla }_{\mathit{q}}U$ [16]; moreover, in the case of constant coefficients, this integration sequence has been shown to yield a superconvergent (4th order) error in configurational ($\mathit{q}$-dependent) observables [15].

In the case of non-constant coefficients $\mathbf{\Gamma }\left(\mathit{q}\right),\mathbf{\Sigma }\left(\mathit{q}\right)$, exact solution of the O-step can be computationally costly. More specifically, bearing in mind that $\mathit{q}\left(t\right)\equiv \mathit{q}$ is fixed during the isolated Ornstein Uhlenbeck process (58), a time-step $\Delta t>0$ can be written:

$$\mathit{p}(t+\Delta t)={\mathit{G}}_t\mathit{p}\left(t\right)+{\mathit{S}}_t\mathit{R},$$

(59)

where $\mathit{R}\sim \mathcal{N}(\mathbf{0},{\mathbf{I}}_{Nd})$, and

$${\mathit{G}}_t=e^{-\Delta t\mathbf{\Gamma }\left(\mathit{q}\right)}.$$

(60)

The matrix ${\mathit{S}}_t\in {\mathbb{R}}^{Nd\times Nd}$ is related to ${\mathit{G}}_t$ as

$${\mathit{S}}_t{\mathit{S}}_t^T={\mathbf{C}}_t-{\mathit{G}}_t{\mathbf{C}}_t{\mathit{G}}_t^T,$$

(61)

where ${\mathbf{C}}_t$ is the solution of the Lyapunov equation

$$\mathbf{\Gamma }\left(\mathit{q}\right){\mathbf{C}}_t+{\mathbf{C}}_t{\mathbf{\Gamma }}^T\left(\mathit{q}\right)=\mathbf{\Sigma }\left(\mathit{q}\right){\mathbf{\Sigma }}^T\left(\mathit{q}\right).$$

(62)

This means that in order to compute an exact solution of the O-part as an (59), first the Lyapunov Equation (62) has to solved, followed by the computation of the Cholesky decomposition (61) and matrix exponential (60). Each of these operations is without any additional assumptions on the structure of the matrices $\mathbf{\Gamma }\left(\mathit{q}\right)$ and $\mathbf{\Sigma }\left(\mathit{q}\right)$ of computational complexity $\mathcal{O}\left(N^3{d}^3\right)$. We circumvent these computations by instead integrating the O-part using a numerical method. We construct a symmetric splitting method based on the decomposition

$$\dot{\mathit{p}}=\underset{=:\mathrm{D}}{\underbrace{-\mathbf{\Gamma }\left(\mathit{q}\right)\mathit{p}}}+\underset{=:\mathrm{F}}{\underbrace{\mathbf{\Sigma }\left(\mathit{q}\right)\dot{\mathit{W}}}},$$

and use the integration sequence DFD, which results in an update of the form

$$\mathit{p}\leftarrow {\mathit{G}}^2\mathit{p}+\sqrt{\Delta t}\mathit{G}\mathbf{\Sigma }\left(\mathit{q}\right)\mathit{R},\mathit{R}\sim \mathcal{N}(\mathbf{0},{\mathbf{I}}_n),$$

(63)

with

$$\mathit{G}=\exp \left(-\frac{\Delta t}{2}\mathbf{\Gamma }\left(\mathit{q}\right)\right).$$

In order to avoid large numerical errors induced by this approximation we employ a multiple timestepping approach [17], meaning that instead of executing a single integration step of (63) for a full time step, we repeat for $K\in \mathbb{N}$, the integrations step (63) K-times using a step size $\Delta t/K$.

We refer this variant of the BAOAB integration scheme using a multiple time stepping solution of the O-part as “multiple time stepping BAOAB” (m-BAOAB). We provide the implementation of this method in Algorithm 1. The performance evaluation of schemes such as this can be demanding as there is a trade off between computational costs and accuracy of the approximation of the solution of the Ornstein-Uhlenbeck process, which in turn affects the accuracy of the invariant measure of the numerical method. The parameter K in the m-BAOAB algorithm allows control of the accuracy with which the OU process is resolved. The method was tested in the example of Section 5.3, below, but the detailed numerical analysis of this method will be taken up in a forthcoming article.

Algorithm 1 m-BAOAB

1: INPUT $(\mathit{q},\mathit{p}),\mathbf{\Gamma },\mathbf{\Sigma },\Delta t,K$ 2: $\mathit{q}\leftarrow \mathit{q}+(\Delta t/2)\mathit{p}$ 3: $\mathit{p}\leftarrow \mathit{p}+(\Delta t/2)\mathit{F}\left(\mathit{q}\right)$ 4: $\mathit{G}\leftarrow \exp (-\frac{\Delta t}{2K}\mathbf{\Gamma }\left(\mathit{q}\right))$ 5: for $i=1$ to K do 6:  $\mathit{p}\leftarrow {\mathit{G}}^2\mathit{p}+\mathit{G}\sqrt{\Delta t/K}\mathbf{\Sigma }\left(\mathit{q}\right){\mathit{R}}_i,{\mathit{R}}_i\sim \mathcal{N}(\mathbf{0},{\mathbf{I}}_n)$ 7: end for 8: $\mathit{p}\leftarrow \mathit{p}+(\Delta t/2)\mathit{F}\left(\mathit{q}\right)$ 9: $\mathit{q}\leftarrow \mathit{q}+(\Delta t/2)\mathit{p}$ 10: return $(\mathit{q},\mathit{p})$

5. Numerical Experiments

In this section we consider three systems that fit the form of (1) and (2). In the first, we consider a system of particles diffusing in a temperature gradient. The second model is a two-dimensional example incorporating a nonconservative stirring force. The third example consists of a particle flocking model of the stochastic Cucker-Smale type, as discussed in Section 3. In all of these examples, the model is sufficiently complicated that we do not possess an analytical solution for the nonequilibrium steady state. For this reason, we are only able to argue for the correctness of the numerical results based on qualitative features (or from our understanding of an unperturbed equilibrium model, in the second example). We are also able to explore the robustness of the numerical solution obtained with respect to variation of the initial conditions, which illustrates the ergodic property. Ultimately comparisons need to be performed with respect to laboratory experiments which can assess the efficacy of the modelling approach in particular situations.

The simulation of the stochastic Cucker-Smale model required the above introduced m-BAOAB scheme for efficient computation. The second model relies on constant (scalar) friction and diffusion coefficient. Due to the diagonal structure of $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$ in the first model, the additional computational costs incurred for the computation of the coefficients in (59) are minimal.

5.1. Particle System with Temperature Gradient

We first study a simple particle system with a position dependent temperature parameter as described in Section 3.4, comprising $N=64$ particles on a two dimension torus, i.e., ${\mathit{q}}_i\in L{\mathbb{T}}^2$ with $L=5$, which are heated by a source at the center of the simulation box. The heat source is modeled by choosing the position-dependent heat bath temperature as

$${\left(\beta \left({\mathit{q}}_i\right)\right)}^{-1}=\Psi \left(|{\mathit{q}}_i-c_L|\right),$$

where $c_L=\frac{1}{2}\left(L,L\right)$ is the center point of the simulation box and $\Psi$ is a smooth bump function of the form

$$\Psi \left(r\right)=\left\{\begin{array}{cc}{T}_{\min }+(T_{\max }-T_{\min })\exp (-\frac{1}{1-{(r/r_{\max })}^2}),\hfill & r\le r_{\max }\hfill \\ T_{\min },\hfill & r>r_{\max }\hfill \end{array}\right.$$

(64)

The constant $T_{\max }>0$ corresponds to the maximum heat bath temperature at the center point $c_L$ of the simulation box and $T_{\min }>0$ describes the heat bath temperature outside of the disk

$$B_{r_{\max }}\left(c_L\right):=\{\mathit{q}\in L{\mathbb{T}}^2:|\mathit{q}-c_L|\le r_{\max }\}.$$

The potential function U is modeled as the sum of pairwise potentials, i.e.,

$$U\left(\mathit{q}\right)=\sum _{i=1}^{N-1}\sum _{j=1}^{i-1}\omega \left(r_{i,j}\right)$$

with $r_{i,j}:=|{\mathit{q}}_j-{\mathit{q}}_i|$ and $\omega$ is a simple repulsive soft potential of the form

$$\omega \left(r_{i,j}\right)=\left\{\begin{array}{cc}\frac{k}{2}{(r_{i,j}-c_r)}^2,\hfill & r_{i,j}<c_r\hfill \\ 0,\hfill & r_{i,j}\ge c_r.\hfill \end{array}\right.$$

where k and $c_r$ are positive constants. The pair interaction described by $\omega$ corresponds to a harmonic spring of stiffness k and rest length $c_r$. Particle systems involving this type of pair potential are commonly used as benchmark systems in the context of dissipative particle dynamics [39,40]. Due to the isolated jump discontinuity in the second derivative of the potential Theorem 3 does not strictly speaking, apply in this case. Although it would be easy to modify the potential to have any desired level of smoothness, we expect, based on our experience with molecular dynamics problems where similar such issues arise, that the results will be very similar to those obtained with the potential given here. The simulation results reported in the remainder of this section are obtained for a parameterization of the model where $k=25,c_r=1$ and $r_{\max }=1$ with $T_{\min }=1/10$ and $T_{\max }=4$. The particle positions and momenta were initialized on an equidistant grid such that ${\mathit{q}}_i=\frac{L}{N}(\lfloor \frac{i}{N}\rfloor ,i\mathrm{mod}N)$ and ${\mathit{p}}_i=(0,0)$, respectively. $N_t={10}^5$ timesteps of stepsize $h=2\times {10}^{-2}$ were simulated with varying values $\gamma \in \{{10}^i,i=-1,0,1,2\}$ of the friction coefficient. Define by

$${\beta }_{\gamma }^{-1}\left(r\right):=\mathbb{E}\left({\mathit{p}}_i^2\right||{\mathit{q}}_i-c_L|=r)$$

the effective temperature of the system at distance $r\in [0,L/2)$ from the center $c_L$ of the simulation box. Figure 1B, shows estimates of this function for different values of $\gamma$ calculated at points $r_i={\Delta }_r/2+i{\Delta }_r,i=0,\dots ,\lfloor \frac{L}{2{\Delta }_r}\rfloor$, as

$${\widehat{\beta }}_{\gamma }^{-1}\left(r_i\right)=\frac{c_i}{N_i}\sum _{j=1}^N\sum _{k=1}^{N_t}{𝟙}_{[r_i-{\Delta }_r/2,r_i+{\Delta }_r/2)}(|{\mathit{q}}_j^{\left(k\right)}-c_L|)\frac{1}{2}\left({({\mathit{p}}_{j,1}^{\left(k\right)})}^2+{({\mathit{p}}_{j,2}^{\left(k\right)})}^2\right),$$

(65)

with

$$c_i:=\frac{1}{[{(r_i+{\Delta }_r/2)}^2-{(r_i-{\Delta }_r/2)}^2]}{\left(\sum _{i=1}^{\lfloor \frac{L}{2{\Delta }_r}\rfloor }\frac{1}{{(r_i+{\Delta }_r/2)}^2-{(r_i-{\Delta }_r/2)}^2}\right)}^{-1},$$

and

$$N_i=\sum _{j=1}^N\sum _{k=1}^{N_t}{𝟙}_{[r_i-{\Delta }_r/2,r_i+{\Delta }_r/2)}(|{\mathit{q}}_j^{\left(k\right)}-c_L|).$$

A larger value of $\gamma$ leads to a tighter local coupling of the particle system with the heat bath and hence for large friction $\gamma ={10}^2$, the estimates of the effective temperature coincide very well with the heat bath temperature given by $\Psi$ at the respective distances from the center of the simulation box. For smaller friction values the position dependence of the effective temperature is less apparent and for very small friction values, e.g., $\gamma ={10}^{-2}$, it is nearly lost entirely. At the same time the decay of the temperature gradient for smaller friction values leads to an increase in the temperature outside the heat source region. Indication of this effect can be also found in the estimated radial density function (see Figure 1A). For small values of the friction coefficient the increased temperature leads to a stronger fluctuation of inter-particle distances outside the heat source area, and hence minima and maxima in the density of radial distribution are less apparent in comparison to the radial distributions for systems with a high friction value.

Due to the relatively small size of the heat source area the number of particles in this area is small in comparison to the number of particles outside of the heat source area and hence the shape of the radial density is mainly determined by the statistical properties of the particles outside the heat source area. The empirical particle densities shown in Figure 2 provide further insight into how the value of the friction coefficient affects the properties of ${\mu }_{\gamma ,\beta }$. For a small friction coefficient, i.e., $\gamma ={10}^{-2}$ the particle density is very close to uniform on $L{\mathbb{T}}^2$. If the heat bath temperature was constant in $L{\mathbb{T}}^2$ one would expect a uniform density as the potential function U is solely composed of pair potentials and thus translation-invariant in the coordinates of the particle density. More precisely, let ${\mathit{q}}_{i,1}$ and ${\mathit{q}}_{i,2}$ denote the first and second coordinate value of the i-th particle, then

$$U\left(\mathit{q}\right)=U\left({({\mathit{q}}_{i,1}+a_1,{\mathit{q}}_{i,2}+a_2)}_{1\le i\le N}\right),$$

for all $a_1,a_2\in \mathbb{R}$. Therefore, the observed uniform particle density for $\gamma ={10}^{-2}$ is consistent with the observation that the position dependency of the effective temperature vanishes for this friction value. While for the considered trajectory length/sample size the plot of the invariant for a friction coefficient $\gamma ={10}^2$ is very noisy, it still strongly suggests that the translation invariance of the particle density, i.e., the uniformity of the particle density, is broken in this regime. Within the interior of the heat source region particles are distributed approximately uniformly whereas close to the boundary of the heat source region the particle density is increased. Outside the heat source region the particle density is concentrated around positions close to energy minima of the potential energy function U.

5.2. System with Non-Conservative Force

We next consider an instance of a non-equilibrium system with a non-conservative force of a form as outlined in Section 3.2. More specifically, we let ${\Omega }_q={\mathbb{R}}^n$ with $n=2$ and let the force $\mathit{F}$ be of the form (31), i.e.,

$$\mathit{F}\left(\mathit{q}\right)=-{\nabla }_{\mathit{q}}U\left(\mathit{q}\right)+\alpha \mathit{J}\mathit{q},$$

with

$$U\left(\mathit{q}\right)=\sum _{i=1}^2{({\mathit{q}}_i^2-1)}^2,$$

and

$$\mathit{J}:=\frac{3}{2}\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right).$$

We furthermore assume that the system is driven by a standard Langevin diffusion, i.e.,

$$\mathbf{\Gamma }=\gamma {\mathit{I}}_n,$$

and

$$\mathbf{\Sigma }=\sqrt{2\gamma {\beta }^{-1}}{\mathit{I}}_n,$$

where we choose $\gamma =1$, and $\beta =1$. This means, that without the non-conservative force part, i.e., in the case $\alpha =0$, the system considered resembles to a particle moving in a 4-well potential driven by a standard Langevin equation at equilibrium at unit temperature. The non-conservative force part $\alpha \mathit{J}\mathit{q}$, corresponds to a stirring force, which pushes the system radially and in clockwise direction around the origin. The effect of the stirring force can be seen in Figure 3. In the absence of the stirring force (see Figure 3A) the invariant distribution is exactly the canonical distribution, i.e.,

$${\rho }_0(\mathit{q},\mathit{p})\propto \exp \left(-U\left(\mathit{q}\right)-\frac{1}{2}\mathit{p}·\mathit{p}\right).$$

The modes of the marginal density in $\mathit{q}$ are exactly positioned at the energy minima of U at $(\pm 1,\pm 1)$. Estimates of the mean momentum vectorfield

$$\overline{\mathit{p}}\left(\mathit{q}\right):={\int }_{{\Omega }_q}\mathit{p}{\rho }_{\alpha }(\mathit{q},\mathit{p})\mathrm{d}\mathit{p},$$

(66)

vanish at each point in the configurational phase space. In the presence of the stirring force (see Figure 3B), the invariant density is rotated slightly and smeared out over the energy barriers. The mean momentum resembles a vector field spiralling clockwise around the origin.

5.3. Stochastic Cucker-Smale Model

In this section we present simulation results for a modified stochastic Cucker-Smale model as described in Section 3.5.2. In particular, we demonstrate that the peculiar and consensus temperature as well as the consensus diffusion rate can be controlled independently and match in simulation with the analytically derived values, if the friction tensor $\mathbf{\Gamma }$ and the diffusion tensor $\mathbf{\Sigma }$ are suitably parametrized. Apart from these macroscopic quantities, we also demonstrate the effect of different parameter values of ${\gamma }_{\perp },{\gamma }_{\left|\right|},T_{\perp },T_{\left|\right|}$ on other macroscopic quantities which in general do not possess a closed form, such as the radial density of particles and peculiar velocity autocorrelation function. Several animations of the particle motion, which effect of parameter choices, may be viewed at the following web location: https://github.com/MatthiasSachs/StochasticCuckerSmale.git.

5.3.1. Model Parametrization

If not stated otherwise the results presented in this section were obtained from a system of 64 particles in a periodic simulation box of edge length 5 and dimension 2, i.e., ${\Omega }_q=L{\mathbb{T}}^{Nd}$ with $d=2,L=5$ and $N=64$. We assume the particle pair-potentials ${\tilde{U}}_{i,j}$ in the definition of U to be identical Morse potentials, i.e., for $i\ne j$

$${\tilde{U}}_{i,j}\left(r\right)=D{(1-e^{-a(r-c_r)})}^2,r>0,$$

with $D=1,a=1,c_r=1/2$. Following [4] we choose the functions $\psi$ in the definition of ${\gamma }_{ij}$ in (36) to be of the form

$$\psi \left(r\right)=\frac{K}{1+{\parallel r\parallel }^{\alpha }},$$

with $K=1/10$, and $\alpha =6$.

5.3.2. Independent Control of Peculiar and Consensus Temperature

We first demonstrate that as stated in Proposition 1, we can indeed control the peculiar and the consensus temperature independently, such that these quantities coincide with the values of the model parameters $T_{\left|\right|}$ and $T_{\perp }$, respectively, i.e., we show that the identities (53) and (54) are reproduced in simulations. For this purpose we compute an estimate of the peculiar temperature at time step $n\in \mathbb{N}$ as

$${\phi }_{T_{\perp }}({\mathit{p}}^{\left(k\right)})=\frac{1}{N(N-1)d}\sum _{i=1}^{N-1}\sum _{j=i+1}^N\sum _{l=1}^d{\left({\mathit{p}}_{j,l}^{\left(k\right)}-{\mathit{p}}_{i,l}^{\left(k\right)}\right)}^2$$

(67)

and an estimate of the consensus temperature at time step $k\in \mathbb{N}$ as

$${\phi }_{T_{\left|\right|}}({\mathit{p}}^{\left(k\right)})=\frac{1}{d}\sum _{l=1}^d{\left({\overline{\mathit{p}}}_l^{\left(k\right)}\right)}^2.$$

(68)

In Figure 4 we show the results for a system with $T_{\perp }=1/2$ and $T_{\left|\right|}=5$. We see that the cumulative average of the respective estimates converges, after a short equilibration period, to the target values.

5.3.3. Properties of the Flock

We first demonstrate the effect of the parameterization of the peculiar dynamics (the values of the peculiar friction ${\gamma }_{\perp }$ and the peculiar temperature $T_{\perp }$) on the flock size, flock formation, and inter-flock diffusion. While we vary the values of the peculiar temperature parameter and peculiar friction parameter in order to demonstrate the effect of these parameters on the above quantities, we leave the parameter values of the consensus dynamics unchanged in all simulations as ${\gamma }_{\left|\right|}=1$ and $T_{\left|\right|}=1$. In the first series of simulations we vary the value of $T_{\perp }$, while fixing the value of the peculiar friction parameter to ${\gamma }_{\perp }=1$. Figure 5 shows the radial density for different values of $T_{\perp }$. As one would expect, the probability mass is concentrated around the mean rest-length of the Morse-potential for low values of $T_{\perp }$ and the distribution spreads out for higher values of $T_{\perp }$. If we measure flock size in terms of the mean distance between particles in the flock, this implies that the flock size grows for increased temperature.

We next explore the effects of the peculiar friction parameter ${\gamma }_{\perp }$. We study the effect of the value of ${\gamma }_{\perp }$ on the flock formation using the mean distance between particles,

$${\phi }_{\mathrm{md}}\left(\mathit{q}\right)=\frac{2}{N(N-1)}\sum _{i=1}^{N-1}\sum _{j=i+1}^N\parallel {\mathit{q}}_j-{\mathit{q}}_j\parallel ,$$

as a measure of the progression of the flock formation. We initialize the particle position out of equilibrium on a equidistant square grid covering a square with side length $5/2$. Figure 6A shows the time evolution of observed mean distance ${\phi }_{\mathrm{md}}\left(\mathit{q}\right)$ in simulations using different values of ${\gamma }_{\perp }$. We find that for small value of ${\gamma }_{\perp }$ flock formation comes with strong oscillation in the flock size (measured in terms of the mean particle distance), which only slowly decays. With increased values of ${\gamma }_{\perp }$ these oscillations are more strongly damped so that the flock size quickly approaches its equilibrium value.

We next explore the effect of the value of the model parameter ${\gamma }_{\perp }$ on the mobility of particles within the flock. In order to measure the strength of diffusion within the flock we consider the mean distance of the particle i to the center of mass, i.e.,

$${\phi }_{\mathrm{mc}}^i\left(\mathit{q}\right)=\parallel {\mathit{q}}_i-\overline{\mathit{q}}\parallel .$$

The autocorrelation of this observable can be used as a measure of mobility within the flock. Figure 6B, shows the autocorrelation function for ${\phi }_{\mathrm{mc}}^i$ calculated from a trajectory in equilibrium, i.e., after the initial transitional flock-formation phase described in the previous paragraph. We can see that the decay of the autocorrelation function becomes slower for increasing values of ${\gamma }_{\perp }$, which indicates that for an increased value of ${\gamma }_{\perp }$ the mobility of particles within the flock is reduced.

5.3.4. Collective Motion

We next explore the effects of the values of ${\gamma }_{\left|\right|}$ and $T_{\left|\right|}$ on the collective motion of the flock, i.e., the diffusive behaviour of the center of mass. Since the motion of the center of mass is described by the consensus dynamics, which is driven by an Ornstein-Uhlenbeck process, we find the diffusion constant D of the center of mass (viewed as a free particle in space) to be

$$D=\frac{T_{\left|\right|}}{N{\gamma }_{\left|\right|}}.$$

Table 1. Estimates of $ND$ for various values of ${\gamma }_{\left|\right|}$ and $T_{\left|\right|}$. $N=16$ denotes the number of particles and D the diffusion coefficient of the diffusive motion of the center of mass. Values in parentheses show the exact values.

${\mathit{T}}_{\left|\right|}$	1	10	100
${\mathit{\gamma }}_{\left|\right|}^{-1}$
1	0.987, (1)	9.775, (10)	101.8, (100)
3	3.059, (3)	29.445, (30)	308.1, (300)
9	8.476, (9)	91.994, (90)	861.4, (900)

shows the estimated diffusion coefficients for various values of ${\gamma }_{\left|\right|}$ and $T_{\left|\right|}$. We find a good match between theoretically predicted values and observed values.

6. Conclusions

In this article we have provided a general treatment of the convergence of Langevin dynamics to a stationary state, including for systems with configuration-dependent friction and noise, as well as nonconservative forces. We have demonstrated the concepts in applications to systems with temperature gradients and stochastic models of active particle systems. Our approach does not assume the usual fluctuation-dissipation rule, so it can be applied to a wide range of nonequilibrium molecular and particle systems where the form of the stationary distribution is a priori unknown. Future work might look at the use of configuration-dependent memory kernels (within a generalized Langevin equation setting) and proofs of ergodicity for degenerate thermostats, for example the pairwise adaptive thermostats [40], which provide an alternative method for controlling observables.