# Stochastic motion in phase space on a surface of constant energy

Tânia Tomé and Mário J. de Oliveira  Universidade de São Paulo, Instituto de
Física, Rua do Matão, 1371, 05508-090 São Paulo, SP, Brasil.

###### Abstract

We study closed systems of particles that are subject to stochastic forces in
addition to the conservative forces. The stochastic equations of motion are
set up in such a way that the energy is strictly conserved at all times. To
ensure this conservation law, the evolution equation for the probability
density is derived using an appropriate interpretation of the stochastic
equation of motion that is not the Itô nor the Stratonovic interpretation. The
trajectories in phase space are restricted to the surface of constant energy.
Despite this restriction, the entropy is shown to increase with time,
expressing irreversible behavior and relaxation to equilibrium. This main
result of the present approach contrasts with that given by the Liouville
equation, which also describes closed systems, but does not show
irreversibility.

##  I Introduction

One of the fundamental problems of non-equilibrium statistical mechanics is
the explanation of the irreversible decay to equilibrium of closed systems,
which is expressed by the increase of entropy. If we consider a closed system
of particles interacting through conservative forces, then we are naturally
led to a description of the system through the Liouville equation gibbs1902 ;
borel1925 ; tolman1938 ; khinchin1949 ; landau1958 ; prigogine1962 ;
hobson1971 ; zubarev1974 ; arnold1978 ; reichl1980 ; salinas2001 . However, a
direct result of this equation is the invariance of entropy with time and as a
consequence the equation is unable to describe the desired irreversible decay
to equilibrium.

The Liouville equation describes appropriately systems in equilibrium. Any
probability density $\rho$ which depends on the dynamic variables through the
Hamiltonian ${\cal H}$ is a stationary solution of this equation. This is the
case of the microcanonical Gibbs distribution

$\rho_{e}=\frac{1}{\Omega}\delta(E-{\cal H}),$ (1)

which describes a system in thermodynamic equilibrium with a fixed energy $E$
gibbs1902 . A relevant property of the Liouville equation is the conservation
of the energy. If the initial probability density is defined on a surface of a
constant given energy, then the density will remain forever defined on this
surface.

The Liouville equation predicts that the entropy remains invariant in time, a
result which is a direct consequence of the constance of the probability
density along a trajectory in phase space. If we start with a probability
density with an entropy distinct from that of the Gibbs equilibrium
distribution, we conclude that this distribution will never be reached and
that the Liouville equation does not account for the irreversibility. This
inconvenience may be circumvent by introducing random forces that changes the
Hamiltonian trajectory while keeping the conservation of energy and momentum.

The Boltzmann kinetc equation huang1963 ; rumer1980 ; cercignani1988 ;
kremer2010 ; oliveira2019  is in fact a realization of this proposal. It
incorporates random forces and as a consequence the entropy increases while
the energy remains strictly constant. Here we derive a stochastic equation
that is similar to Boltzmann equation in the sense that the entropy also
increases at strictly constant energy. To this end we set up a stochastic
equation of motion, or Langevin equations, that strictly conserves the energy,
that is the trajectories in phase space lie on the surface of constant energy.
From the stochastic equation of motion one derives the associated evolution
equation for the probability density, which turns up to be an equation of the
Fokker-Planck type kampen1981 ; risken1989 ; gardiner2009 ; tome2015 , and it
is not an integro-differential equation as is the case of the Boltzmann
equation.

The time evolution equation for the probability density, which incorporates
stochastic forces, is a type of equations that is used in approaches to
stochastic thermodynamics tome2010 ; vandebroeck2010 ; spinney2012 ;
seifert2012 ; santillan2013 ; luposchainsky2013 ; tome2015a ; oliveira2020 ;
oliveira2020a , which describe systems in contact with heat reservoirs,
including the exchange of heat and entropy as well as the production of
entropy. The main difference of these approaches from ours is the strictly
conservation of energy observed in the present approach. In this sense we may
say that the present approach provides a stochastic thermodynamics for closed
systems.

To derive the evolution equation, it is necessary to give an interpretation to
the stochastic equation of motion because the noise is multiplicative
kampen1981 . Usually, the interpretation is that given by Stratonovich and not
that given by Itô. However, we find that the proper interpretation necessary
to preserve energy at all times is none of them.

The origin of stochastic trajectories may be attributed to forces that are of
the hyperbolic type such as those occurring when two hard spheres collide with
each other. Given the initial positions and velocities of two hard spheres,
then the positions and velocities after the collision will be uniquely
determined by the Newton equations of motion. However, if we consider two very
similar initial conditions, that differ only slightly by the positions, the
velocities being the same, then the velocities after the collision will not
differ slightly but can be very different. Since the initial condition is
almost the same, this result looks as if more than one trajectory is possible
from a single initial condition.

We will not pursue in more detail the problem of how stochastic trajectories
emerge from pure mechanics, that is, from the Newton equations of motion
applied to a system acted by internal conservative forces. Here we adopt the
point of view that the stochastic trajectories or rather the stochastic forces
are an assumption of the present approach.

The stochastic equation of motion that we have introduced can be used as a
numerical method to simulate a system of interacting particles at constant
energy, as is the case of the method of molecular dynamics allen1987 ;
tuckerman2010 . In equilibrium, the static properties will be the same as
those obtained from the molecular dynamics because the stochastic equations of
motion lead to the Gibbs distribution. However, the two-time correlations will
be distinct due to presence of the stochastic forces.

##  II Impulsive stochastic force

###  II.1 Stochastic equation of motion

We consider two particles interacting through random forces that act on a
short period of time. We denote by ${\bf p}_{1}=m_{1}{\bf v}_{1}$ and ${\bf
p}_{2}=m_{2}{\bf v}_{2}$ the momenta of the particles, where $m_{1}$ and
$m_{2}$ are their masses and $\bf v_{1}$ and $\bf v_{2}$ their velocities. The
random force acting on particle 1 due to the particle 2 is denoted by $\xi$.
The reaction on particle 2 is $-\mbox{\boldmath$\xi$}$ and the equations of
motion are

$\frac{d{\bf p}_{1}}{dt}={\mbox{\boldmath$\xi$}},\qquad\frac{d{\bf
p}_{2}}{dt}=-{\mbox{\boldmath$\xi$}}.$ (2)

The total momentum ${\bf P}=\bf p_{1}+\bf p_{2}$ is conserved because from the
stochastic equation of motion it follows that $d{\bf P}/dt=0$.

We wish to determine the properties of the stochastic vector variable $\xi$
that makes the energy of the two particles constant. As the forces are suppose
to be impulsive, which means that they act during a short period of time, the
conservation of energy means the conservation of the kinetic energy

${\cal K}=\frac{p_{1}^{2}}{2m_{1}}+\frac{p_{2}^{2}}{2m_{2}}.$ (3)

The derivative of the kinetic energy is

$\frac{d{\cal K}}{dt}=({\bf v}_{1}-{\bf v}_{2})\cdot{\mbox{\boldmath$\xi$}},$
(4)

Defining the relative velocity ${\bf v}={\bf v}_{1}-{\bf v}_{2}$, we see that
the condition

${\bf v}\cdot{\mbox{\boldmath$\xi$}}=0$ (5)

makes the derivative of the kinetic energy to vanish. This condition of
orthogonality implies that $\xi$ depend on the relative velocity $\bf v$.

To proceed further we use a more convenient procedure. As $\bf p_{1}$ and $\bf
p_{2}$ are connected by the conservation of the total momentum, if suffices to
use just one independent variable which we choose to be ${\bf p}=m{\bf v}$
where $m$ is the reduced mass. The two equations of motion are then reduced to
one,

$\frac{d{\bf p}}{dt}={\mbox{\boldmath$\xi$}}.\qquad$ (6)

In terms of the new variable ${\bf p}$, the kinetic energy, apart from a
constant, is given by

${\cal K}=\frac{p^{2}}{2m}.$ (7)

The stochastic variable $\xi$ is understood as follows. If we discretize the
time in intervals equal to $\tau$ the stochastic variable is proportional to
$1/\sqrt{\tau}$. Thus in a discretize formulation of the equation of motion
(6), we replace $\xi$ by ${\mbox{\boldmath$\xi$}}/\sqrt{\tau}$ with the
understanding that $\xi$ has a finite variance. The equation (6) becomes

$\Delta{\bf p}=\sqrt{\tau}\mbox{\boldmath$\xi$}.$ (8)

A finite increment of the kinetic energy is

$\Delta{\cal K}={\bf v}\cdot\Delta{\bf p}+\frac{1}{2m}(\Delta{\bf p})^{2}.$
(9)

Replacing (8) in this equation, we find

$\Delta{\cal K}=\sqrt{\tau}{\bf
v\cdot\mbox{\boldmath$\xi$}}+\frac{\tau}{2m}{\mbox{\boldmath$\xi$}}\cdot{\mbox{\boldmath$\xi$}},$
(10)

Using the condition (5) on $\xi$ established above, ${\bf
v}\cdot{\mbox{\boldmath$\xi$}}=0$, we see that $\Delta{\cal K}$ vanishes up to
order $\sqrt{\tau}$, but not up to order $\tau$. To overcome this
inconvenience, we observe that (8) is not the only possible interpretation of
the stochastic equation of motion (6). We may add any random variable
proportional to $1/\tau$ which means to interpret the stochastic equation of
motion as

$\Delta{\bf
p}=\sqrt{\tau}\mbox{\boldmath$\xi$}+\tau{\mbox{\boldmath$\omega$}}.$ (11)

where $\omega$ is a random variable with a finite variance. Replacing this
expression in (9), then up to terms of order $\tau$ we find

$\Delta{\cal K}=\sqrt{\tau}{\bf v}\cdot\mbox{\boldmath$\xi$}+\tau({\bf
v}\cdot{\mbox{\boldmath$\omega$}}+\frac{1}{2m}{\mbox{\boldmath$\xi$}}\cdot{\mbox{\boldmath$\xi$}}).$
(12)

The first term vanishes due to the condition (5). If we wish $\Delta{\cal K}$
to vanish to order $\tau$, the last term must vanish, which gives the
following relation

${\mbox{\boldmath$\omega$}}\cdot{\bf
p}+\frac{1}{2}{\mbox{\boldmath$\xi$}}\cdot{\mbox{\boldmath$\xi$}}=0,$ (13)

between $\omega$ and $\xi$. To find $\omega$ in terms of $\xi$, we
differentiate the condition (5) with respect to ${\bf p}$, and perform a dot
product with $\xi$ to obtain

${\mbox{\boldmath$\xi$}}\cdot(\frac{\partial}{\partial{\bf
p}}{\mbox{\boldmath$\xi$}})\cdot{\bf
p}+{\mbox{\boldmath$\xi$}}\cdot{\mbox{\boldmath$\xi$}}=0.$ (14)

Comparing the last two equation, we see that the expression

${\mbox{\boldmath$\omega$}}^{\prime}=\frac{1}{2}{\mbox{\boldmath$\xi$}}\cdot\frac{\partial}{\partial{\bf
p}}{\mbox{\boldmath$\xi$}},$ (15)

gives the desired solution. However, this is not the unique solution because
we may add to this expression any term orthogonal to ${\bf p}$ such as $\xi$
itself. At this point we postpone the problem of determining $\omega$ in terms
of $\xi$.

###  II.2 Evolution equation

To derive the Fokker-Planck equation associated to the Langevin equation (11),
which is the evolution equation for the probability density $\rho({\bf p})$,
we use a method that consists in determining the evolution of the
characteristic function $\Phi$, which is the Fourier transform of $\rho$,

$\Phi=\int e^{i\Omega}\rho\,d{\bf p},$ (16)

where $\Omega=i{\bf k}\cdot{\bf p}$ and ${\bf k}$ is a vector in the Fourier
space. The characteristic function is also the average

$\Phi=\langle e^{i\Omega}\rangle.$ (17)

Denoting by $\Delta\Phi$ the difference of the characteristic functions
between times $t+\tau$ and $t$, then

$\Delta\Phi=\langle e^{i\Omega}(e^{i{\bf k}\cdot\Delta{\bf p}}-1)\rangle,$
(18)

where the average is to be calculated using the probability of the variables
$\bf p$ and $\xi$ considered to be independent of each other. Next we replace
$\Delta{\bf p}$, given by (11) in this expression and expand $\Delta\Phi$ up
to terms of order $\tau$. Dividing the result by $\tau$, we find the time
derivative of the characteristic function in the form

$\frac{\partial\Phi}{\partial t}=\langle e^{i\Omega}i{\bf
k}\cdot{\mbox{\boldmath$\omega$}}\rangle+\frac{1}{2}\langle e^{i\Omega}i{\bf
k}\cdot{\mbox{\boldmath$\xi$}}{\mbox{\boldmath$\xi$}}\cdot i{\bf k}\rangle.$
(19)

Taking into account that $i{\bf k}\cdot e^{i\Omega}=\partial
e^{i\Omega}/\partial{\bf p}$, the first term on the right-hand side of (19) is
written as

$\int(\frac{\partial e^{i\Omega}}{\partial{\bf
p}})\cdot\langle{\mbox{\boldmath$\omega$}}\rangle\,\rho\,d{\bf p}=-\int
e^{i\Omega}\frac{\partial}{\partial{\bf
p}}\cdot\langle{\mbox{\boldmath$\omega$}}\rangle\,\rho\,d{\bf p},$ (20)

where we have performed an integration by parts, and now the average is on
$\xi$ only. This expression is the Fourier transform of

$-\frac{\partial}{\partial{\bf
p}}\cdot\langle{\mbox{\boldmath$\omega$}}\rangle\rho.$ (21)

In this and in other integration by parts, we are considering that the
integrated part disappears by the assumption of a rapid vanishing of $\rho$ at
the boundary of integration.

Using the same procedure, the second term on the right-hand side of (19) is
written as

$\frac{1}{2}\int e^{i\Omega}\frac{\partial}{\partial{\bf
p}}\cdot\frac{\partial}{\partial{\bf
p}}\cdot\langle{\mbox{\boldmath$\xi$}\mbox{\boldmath$\xi$}}\rangle\rho\,d{\bf
p},$ (22)

where we have performed two successive integrations by parts. This expression
is the Fourier transform of

$\frac{1}{2}\frac{\partial}{\partial{\bf p}}\cdot\frac{\partial}{\partial{\bf
p}}\cdot\langle{\mbox{\boldmath$\xi$}\mbox{\boldmath$\xi$}}\rangle\rho.$ (23)

To determine $\partial\rho/\partial t$, we bear in mind that
$\partial\Phi/\partial t$, given by (19), is its Fourier transform. Therefore,
the time derivative of $\rho$ is obtained by adding up the expressions (21)
and (23). The result is

$\frac{\partial\rho}{\partial t}=\frac{1}{2}\frac{\partial}{\partial{\bf
p}}\cdot\left(\frac{\partial}{\partial{\bf
p}}\cdot\langle{\mbox{\boldmath$\xi$}\mbox{\boldmath$\xi$}}\rangle\rho-2\langle{\mbox{\boldmath$\omega$}}\rangle\rho\right).$
(24)

It remains yet to determine $\omega$ in terms of $\xi$.

###  II.3 Lagrange multiplier

We have seen that the discretized equation of motion (11) preserves the
kinetic energy up to time $\tau$. It is not guaranteed that the kinetic energy
${\cal K}$ will be preserved forever. To ensure that this will happen for all
times, we impose the vanishing of the probability of a trajectory that do not
lie on a surface of constant kinetic energy. This amounts to say that if
initially $\rho$ is nonzero only on the surface of constant ${\cal K}$, then
this property should be preserved as $\rho$ evolves in time. Let us denote by
$\rho^{*}$ a probability density that is zero outside this surface. Then the
variation of $\rho^{*}$ to a point outside the surface vanishes, that is
$d\rho^{*}=0$. If $\rho$ is a generic solution of (24), then $\rho^{*}$ can be
found from $\rho$ by using the method of the Lagrange multipliers, which means

$d\rho^{*}=d\rho+\lambda d{\cal K},$ (25)

where $\lambda$ is a Lagrange multiplier. This equation is equivalent to

$\rho^{*}=\rho+\lambda{\cal K}.$ (26)

Replacing this expression in the evolution equation (24), written in the
simplified form

$\frac{\partial\rho}{\partial t}={\mathfrak{F}}\rho,$ (27)

where ${\mathfrak{F}}$ is a linear differential operator, we find that
$\rho^{*}$ is a solution of the evolution equation as long as

${\mathfrak{F}}{\cal K}=0.$ (28)

The condition (28) becomes fulfilled if the expression between parentheses in
(24) vanishes when $\rho$ is replaced by $K=p^{2}/2m$, that is,

$\frac{\partial}{\partial{\bf
p}}\cdot{\mbox{\boldmath$\xi$}\mbox{\boldmath$\xi$}}p^{2}-2{\mbox{\boldmath$\omega$}}p^{2}=0.$
(29)

Solving for $\omega$, we find

${\mbox{\boldmath$\omega$}}=\frac{1}{2p^{2}}\frac{\partial}{\partial{\bf
p}}\cdot{\mbox{\boldmath$\xi$}\mbox{\boldmath$\xi$}}p^{2},$ (30)

which, by the use of the relation ${\bf p}\cdot{\mbox{\boldmath$\xi$}}=0$, can
be written in the more simplified form

${\mbox{\boldmath$\omega$}}=\frac{1}{2}\frac{\partial}{\partial{\bf
p}}\cdot{\mbox{\boldmath$\xi$}\mbox{\boldmath$\xi$}}.$ (31)

which is the sought relation between $\omega$ and $\xi$.

Replacing (31) in (24) we obtain the evolution equation in the simple form

$\frac{\partial\rho}{\partial t}=\frac{1}{2}\frac{\partial}{\partial{\bf
p}}\cdot\langle{\mbox{\boldmath$\xi$}\mbox{\boldmath$\xi$}}\rangle\cdot\frac{\partial}{\partial{\bf
p}}\rho.$ (32)

Using (31), the discretized stochastic equation of motion (11) becomes

$\Delta{\bf
p}=\sqrt{\tau}\mbox{\boldmath$\xi$}+\frac{\tau}{2}\frac{\partial}{\partial{\bf
p}}\cdot{\mbox{\boldmath$\xi$}\mbox{\boldmath$\xi$}},$ (33)

which is the desired interpretation of the stochastic equation of motion
leading to the conservation of kinetic energy for all times. It is worth
mentioning that this equation does not correspond to the Stratonovich
interpretation of (6), which is

$\Delta{\bf
p}=\sqrt{\tau}\mbox{\boldmath$\xi$}+\frac{\tau}{2}{\mbox{\boldmath$\xi$}}\cdot\frac{\partial}{\partial{\bf
p}}{\mbox{\boldmath$\xi$}},$ (34)

and corresponds to use ${\mbox{\boldmath$\omega$}}^{\prime}$, given by (15).
Of course, it does not correspond either to the Itô interpretation which is
simply that given by (8).

###  II.4 Expression of $\xi$

It remains now to determine the explicit expression for the stochastic force
$\xi$. The only restriction that has to be fulfilled is the orthogonal
condition (5), that is,

${\bf v}\cdot{\mbox{\boldmath$\xi$}}=0.$ (35)

To meet this condition, we choose two vectors ${\bf a}$ and $\bf b$ that are
orthogonal to $\bf v$, and to each other,

${\bf v}\cdot{\bf a}=0,\qquad{\bf v}\cdot{\bf b}=0,\qquad{\bf a}\cdot{\bf
b}=0,$ (36)

and write

${\mbox{\boldmath$\xi$}}={\bf a}\sigma+{\bf b}\eta,$ (37)

where $\sigma$ and $\eta$ are independent stochastic variable with zero mean
and variance equal to $2\gamma$.

From these definitions we determine the average
$\langle{\mbox{\boldmath$\xi$}\mbox{\boldmath$\xi$}}\rangle$ appearing in the
evolution equation (32)

$\langle{\mbox{\boldmath$\xi$}\mbox{\boldmath$\xi$}}\rangle=2\gamma({\bf
a}{\bf a}+{\bf b}{\bf b}),$ (38)

which replaced in (24), gives

$\frac{\partial\rho}{\partial t}=\gamma\frac{\partial}{\partial{\bf
p}}\cdot({\bf a}{\bf a}+{\bf b}{\bf b})\cdot\frac{\partial}{\partial{\bf
p}}\rho.$ (39)

We write the Cartesian coordinates of the relative velocity $\bf v$ using
spherical coordinates as

${\bf v}=(v\sin\theta\cos\varphi,\,\,v\sin\theta\sin\varphi,\,\,v\cos\theta),$
(40)

where $v$ is the absolute value of the velocity, $\theta$ is the polar angle
and $\varphi$ is the azimuthal angle. The vectors $\bf a$ and $\bf b$ are
chosen to be unit vector with Cartesian coordinates given by

${\bf a}=(\cos\theta\cos\varphi,\,\,\cos\theta\sin\varphi,\,\,-\sin\theta),$
(41) ${\bf b}=(-\sin\theta\sin\varphi,\,\,\sin\theta\cos\varphi,\,\,0),$ (42)

It is easily seen that the conditions (36) are satisfied.

Up to now we have treated a system in three dimensions. To treat a two-
dimensional system, we have to consider that the vector quantities have two
Cartesian components. In this case, instead of two vectors orthogonal to ${\bf
v}$, just one is possible. Writing

${\bf v}=(v\cos\phi,\,\,v\sin\phi).$ (43)

this vector is

${\bf b}=(-\sin\phi,\,\,\cos\phi),$ (44)

In this case there is just one random variable $\eta$ and

${\mbox{\boldmath$\xi$}}={\bf b}\eta.$ (45)

If we wish to treat a system in one dimension, we see that this is
unattainable because it is not possible to meet the condition of orthogonality
(36). This reflects the following result concerning the collision of two
particles in one dimension. Given the velocities of the particles before the
collision, they are uniquely determined after the collision, if the energy and
momentum are conserved. Thus in one dimension there is no room for a
stochastic motion that conserves both the energy and momentum.

###  II.5 Entropy production

Although the energy is strictly conserved, thist is not the case of the
entropy. The entropy is defined by

$S=-k\int\rho\ln\rho\,d{\bf p},$ (46)

and its time variation can be obtained from the the evolution equation (39).
Deriving this expression with respect to time, we obtain

$\frac{dS}{dt}=-k\int\frac{\partial\rho}{\partial t}\ln\rho\,d{\bf p}.$ (47)

Replacing the derivative of $\rho$, given by the evolution equation (39), and
after an integration by parts, we reach the result

$\frac{dS}{dt}=k\gamma\int\frac{1}{\rho}\frac{\partial\rho}{\partial{\bf
p}}\cdot({\bf a}{\bf a}+{\bf b}{\bf b})\cdot\frac{\partial\rho}{\partial{\bf
p}}d{\bf p},$ (48)

which can be written in the form

$\frac{dS}{dt}=k\gamma\int\frac{1}{\rho}(A^{2}+B^{2})d{\bf p},$ (49)

where

$A={\bf a}\cdot\frac{\partial\rho}{\partial{\bf p}},\qquad B={\bf
b}\cdot\frac{\partial\rho}{\partial{\bf p}}.$ (50)

As the integral is positive definite, $dS\geq 0$, the entropy is a monotonic
increasing function of time. The right-hand side of (49) is understood as the
rate of entropy production. In the stationary state $\rho_{e}$ will be a
function of $K$ and $\partial\rho/\partial{\bf p}$ will be proportional to
${\bf p}$ which is orthogonal to $\bf a$ and $\bf b$. Thus $A$ and $B$ vanish,
the entropy production vanishes and the entropy reaches is maximum value.

##  III System of interacting particles

###  III.1 Stochastic equations of motion

Our attention is now directed toward a system of several interacting
particles. The position and the momentum of particle $i$ are denoted by ${\bf
x}_{i}$ and by ${\bf p}_{i}$, and the system is described by the Hamiltonian
function

${\cal H}=\sum_{i}\frac{p_{i}^{2}}{2m_{i}}+{\cal V},$ (51)

where $p_{1}^{2}=p_{xi}^{2}+p_{yi}^{2}+p_{zi}^{2}$, and $m_{i}$ is the mass of
the particle $i$, and ${\cal V}$ is a function of the coordinates and
represent the potential energy of the particles. The stochastic equations of
motion, or the Langevin equations, are

$\frac{d{\bf x}_{i}}{dt}={\bf v}_{i},\qquad\qquad\frac{d{\bf p}_{i}}{dt}={\bf
F}_{i}+{\mbox{\boldmath$\zeta$}}_{i},$ (52)

where ${\bf v}_{i}$ and ${\bf F}_{i}$ are the velocity and the conservative
force, respectively, associated to the particle $i$, and are given by

${\bf v}_{i}=\frac{\partial{\cal H}}{\partial{\bf p}_{i}},\qquad\qquad{\bf
F}_{i}=-\frac{\partial{\cal H}}{\partial{\bf x}_{i}}.$ (53)

We wish to describe a system by forces that strictly conserves the total
momentum of the collection of particles, given by

${\bf P}=\sum_{i}{\bf p}_{i},$ (54)

and the total energy, given by (51). The particles are subject only to
internal forces which means that for each force ${\bf F}_{i}$ there is a
reaction force with the opposite sign so that the sum of the conservative
forces vanish,

$\sum_{i}{\bf F}_{i}=0.$ (55)

Deriving ${\bf P}$ with respect to time and using (55), we get

$\frac{d{\bf P}}{dt}=\sum_{i}{\mbox{\boldmath$\zeta$}}_{i}.$ (56)

Therefore the conservation of momentum requires that the sum of the stochastic
forces vanishes.

Deriving the function ${\cal H}$ with respect to time, we find

$\frac{d{\cal H}}{dt}=\sum_{i}{\bf v}_{i}\cdot\mbox{\boldmath$\zeta$}_{i}.$
(57)

There is no contribution coming from the conservative forces. The right-hand
side of this equation is the total power of the stochastic force, which should
vanish.

The stochastic forces are chosen so that the right-hand sides of the equations
(56) and (57) vanish identically. To meet the first requirement, we choose
${\mbox{\boldmath$\zeta$}}_{i}$ as a sum of independent stochastic vector
variables ${\mbox{\boldmath$\xi$}}_{ij}$,

${\mbox{\boldmath$\zeta$}}_{i}=\sum_{j}{\mbox{\boldmath$\xi$}}_{ij}.$ (58)

with the properties

${\mbox{\boldmath$\xi$}}_{ji}=-{\mbox{\boldmath$\xi$}}_{ij},$ (59)

and ${\mbox{\boldmath$\xi$}}_{ii}=0$. We see that the right-hand side of (56)
vanishes identically, and the total momentum is strictly conserved. The vector
${\mbox{\boldmath$\xi$}}_{ij}$ is understood as the random force acting on
particle $i$ due to particle $j$ and ${\mbox{\boldmath$\xi$}}_{ji}$ as the
random force on particle $j$ due particle $i$. Since they were chosen to
differ only by their signs, they are interpreted as action and reaction,
leading to the conservation of the total momentum. We now replace (58) in the
right-hand side of equation (57) to find

$\frac{d{\cal H}}{dt}=\frac{1}{2}\sum_{ij}({\bf v}_{i}-{\bf
v}_{j})\cdot{\mbox{\boldmath$\xi$}}_{ij},$ (60)

where we have used the property (59). Requiring that

${\bf v}_{ij}\cdot{\mbox{\boldmath$\xi$}}_{ij}=0,$ (61)

for each pair $ij$, where ${\bf v}_{ij}={\bf v}_{i}-{\bf v}_{j}$ is the
relative velocity of particles $i$ and $j$, we see that the right-hand side of
(60) vanishes identically and $d{\cal H}/dt=0$. The random force $\xi$ij
acting on the particles $i$ and $j$ must be orthogonal to their relative
velocities ${\bf v}_{ij}$ and, therefore, performs no work.

To meet the condition (61), we choose two unit vectors ${\bf a}_{ij}$ and
${\bf b}_{ij}$, to be orthogonal to ${\bf v}_{ij}$, and to each other,

${\bf v}_{ij}\cdot{\bf a}_{ij}=0,\qquad{\bf v}_{ij}\cdot{\bf
b}_{ij}=0,\qquad{\bf v}_{ij}\cdot{\bf b}_{ij}=0$ (62)

We also introduce two new scalar random variables $\sigma_{ij}$ and
$\eta_{ij}$, and write

${\mbox{\boldmath$\xi$}}_{ij}={\bf a}_{ij}\sigma_{ij}+{\bf b}_{ij}\eta_{ij}.$
(63)

Due to the orthogonal property (62) of ${\bf a}_{ij}$ and ${\bf b}_{ij}$, we
see that the condition (61) is fulfilled.

The vectors ${\bf a}_{ij}$ and ${\bf b}_{ij}$ are chose to be given by the
expressions (41) and (42). We remark that when we interchange $i$ and $j$, the
vector ${\bf a}_{ij}$ preserves its sign and ${\bf b}_{ij}$ changes sign. As
${\mbox{\boldmath$\xi$}}_{ij}$ must change sign then the stochastic variable
$\sigma_{ij}$ should change its sign and $\eta_{ij}$ should preserve its sign.

The random variables $\sigma_{ij}$ and $\eta_{ij}$ are chosen to have zero
mean and the same variance $2\gamma_{ij}$, which represents the strength of
the random forces and might depend on the positions of the particles $i$ and
$j$. It is reasonable to assume that $\gamma_{ij}$ is nonzero only when the
particles $i$ and $j$ are close to each other.

The variance of ${\mbox{\boldmath$\xi$}}_{ij}$ is
$\langle{\mbox{\boldmath$\xi$}}_{ij}{\mbox{\boldmath$\xi$}}_{ij}\rangle=2{\mathbb{K}}_{ij}$,
where ${\mathbb{K}}_{ij}$ is a $3\times 3$ symmetric matrix given by

${\mathbb{K}}_{ij}=\gamma_{ij}({\bf a}_{ij}{\bf a}_{ij}+{\bf b}_{ij}{\bf
b}_{ij}),$ (64)

and ${\mathbb{K}}_{ii}=0$. The covariances of the random vector
$\mbox{\boldmath$\zeta$}_{i}$ are obtained from (58) and is related to this
matrix by

$\langle{\mbox{\boldmath$\zeta$}}_{i}{\mbox{\boldmath$\zeta$}}_{j}\rangle=-2{\mathbb{K}}_{ij},\qquad
i\neq j,$ (65)
$\langle{\mbox{\boldmath$\zeta$}}_{i}{\mbox{\boldmath$\zeta$}}_{i}\rangle=2\sum_{j}{\mathbb{K}}_{ij}.$
(66)

###  III.2 Evolution equation

To derive the equation that gives the time evolution of the probability
density function $\rho$ of the dynamic variables ${\bf x}_{i}$ and ${\bf
p}_{i}$, we assume the following discretized equations of motion

$\Delta{\bf x}_{i}=\tau{\bf v}_{i},$ (67) $\Delta{\bf p}_{i}=\tau{\bf
F}_{i}+\sqrt{\tau}{\mbox{\boldmath$\zeta$}}_{i}+\frac{\tau}{4}\frac{\partial}{\partial{\bf
p}_{i}}\cdot{\mbox{\boldmath$\zeta$}}_{i}{\mbox{\boldmath$\zeta$}}_{i},$ (68)

where ${\mbox{\boldmath$\zeta$}}_{i}$ is the sum of the independent stochastic
variables ${\mbox{\boldmath$\xi$}}_{ij}$, each one given by (63).

We derive the evolution equation by the method that we have used above. To
this end we define the characteristic function $\Phi$ by

$\Phi=\int e^{i\Omega}\rho\,dqdp,$ (69)

where the integration is over the phase space and

$\Omega=\sum_{j}({\bf k}_{j}\cdot{\bf x}_{j}+{\bf q}_{j}\cdot{\bf p}_{j}),$
(70)

where ${\bf k}_{i}$ and ${\bf q}_{i}$ are vectors in the Fourier space. The
characteristic function is also the average over $\rho$,

$\Phi=\langle e^{i\Omega}\rangle.$ (71)

The finite variation $\Phi$ is given by

$\Delta\Phi=\langle e^{i\Omega}(e^{i\Delta{\Omega}}-1)\rangle.$ (72)

Expanding the right-hand side of this equation up to terms of the order
$\tau$, we find the time derivative of $\Phi$ as

$\frac{\partial\Phi}{\partial t}=\sum_{j}\langle e^{i\Omega}(i{\bf
k}_{j}\cdot{\bf v}_{j}+i{\bf q}_{j}\cdot{\bf F}_{j})\rangle+$
$+\frac{1}{2}\sum_{j}\langle e^{i\Omega}(\frac{\partial}{\partial{\bf
p}_{j}}\cdot{\mbox{\boldmath$\zeta$}}_{j}{\mbox{\boldmath$\zeta$}}_{j}+\frac{1}{2}\sum_{\ell}i{\bf
q}_{\ell}\cdot{\mbox{\boldmath$\zeta$}}_{\ell}{\mbox{\boldmath$\zeta$}}_{j})\cdot
i{\bf q}_{j}\rangle.$ (73)

To obtain the evolution equation, it suffices to take the inverse Fourier
transform of the equation (73). To this end we observe that the first term of
(73) is the Fourier transform of the Poisson brackets

$\\{{\cal H},\rho\\}=\sum_{j}\left(\frac{\partial{\cal H}}{\partial{\bf
x}_{j}}\cdot\frac{\partial\rho}{\partial{\bf p}_{j}}-\frac{\partial{\cal
H}}{\partial{\bf p}_{j}}\frac{\partial\rho}{\partial{\bf x}_{j}}\right),$ (74)

a result which is reached by using the relations

$i{\bf k}_{j}e^{i\Omega}=\frac{\partial e^{i\Omega}}{\partial{\bf
x}_{j}},\qquad i{\bf q}_{j}e^{i\Omega}=\frac{\partial
e^{i\Omega}}{\partial{\bf p}_{j}}.$ (75)

Using (65) and (66), the second term of (73) can be written in the form

$\frac{1}{2}\sum_{j\neq\ell}\langle e^{i\Omega}(\frac{\partial}{\partial{\bf
p}_{j}}\cdot{\mathbb{K}}_{j\ell}+\frac{i}{2}({\bf q}_{j}-{\bf
q}_{\ell})\cdot{\mathbb{K}}_{j\ell})\cdot i{\bf q}_{j}\rangle.$ (76)

Taking into account that ${\mathbb{K}}_{j\ell}$ depends on ${\bf p}_{j}$
through the difference ${\bf p}_{j}-{\bf p}_{\ell}$, we may write

$\frac{1}{2}\sum_{j\neq\ell}\langle e^{i\Omega}({\bf
D}_{j\ell}\cdot{\mathbb{K}}_{j\ell}+i({\bf q}_{j}-{\bf
q}_{\ell})\cdot{\mathbb{K}}_{j\ell})\cdot i{\bf q}_{j}\rangle,$ (77)

where

${\bf D}_{j\ell}=\frac{\partial}{\partial{\bf
p}_{j}}-\frac{\partial}{\partial{\bf p}_{\ell}}.$ (78)

This expression is equal to the integral

$-\frac{1}{2}\sum_{j\neq\ell}\int e^{i\Omega}i{\bf
q}_{j}\cdot{\mathbb{K}}_{j\ell}\cdot{\bf D}_{j\ell}\rho\,dxdp,$ (79)

which was obtained by an integration by parts and using the result

${\bf D}_{j\ell}e^{i\Omega}=i({\bf q}_{j}-{\bf q}_{\ell})e^{i\Omega}.$ (80)

Employing again the second equality of (75) and performing another integration
by parts, we find a result which is the Fourier transform of

$\frac{1}{2}\sum_{j\neq\ell}\frac{\partial}{\partial{\bf
p}_{j}}\cdot{\mathbb{K}}_{j\ell}\cdot{\bf D}_{j\ell}\rho.$ (81)

The evolution equation is obtained by observing that $\partial\Phi/\partial t$
is the Fourier transform of $\partial\rho/\partial t$. Therefore, to reach the
equation it suffices to add up the results (74) and (81),

$\frac{\partial\rho}{\partial t}=\\{{\cal
H},\rho\\}+\frac{1}{2}\sum_{(ij)}{\bf D}_{ij}\cdot{\mathbb{K}}_{ij}\cdot{\bf
D}_{ij}\rho,$ (82)

Replacing ${\mathbb{K}}_{ij}$ in this expression, the evolution equation
acquires the form

$\frac{\partial\rho}{\partial t}=\\{{\cal
H},\rho\\}+\frac{1}{2}\sum_{(ij)}\gamma_{ij}{\bf D}_{ij}\cdot({\bf
a}_{ij}A_{ij}+{\bf b}_{ij}B_{ij}),$ (83)

where

$A_{ij}={\bf a}_{ij}\cdot{\bf D}_{ij}\rho,\qquad B_{ij}={\bf b}_{ij}\cdot{\bf
D}_{ij}\rho.$ (84)

We have argued above in the analysis of the stochastic motion of two particles
that the kinetic energy is a stationary solution of the evolution equation.
The reasoning based on the Lagrange multipliers can be generalized leading to
the rule that a conserved quantity is a stationary solution of the evolution
equation. This is indeed the case of ${\cal H}$ and ${\bf P}$ in relation to
the evolution equation (83), which can be verified by inspection. The Poisson
brackets vanish trivially when $\rho$ is replaced by ${\cal H}$. When it is
replaced by ${\bf P}$ it vanish as well if we use the condition (55). The
quantities $A_{ij}$ and $B_{ij}$ given by (84) vanish when $\rho$ is replaced
by ${\cal H}$ because ${\bf D}_{ij}{\cal H}$ equals ${\bf v}_{ij}$ which is
orthogonal to ${\bf a}_{ij}$ and ${\bf b}_{ij}$. The quantities $A_{ij}$ and
$B_{ij}$ also vanish when $\rho$ is replaced by ${\bf P}$ because ${\bf
D}_{ij}\cdot{\bf P}=0$.

From the results just obtained, it follows that any function of the conserved
quantities will be a stationary solution of the evolution equation. In
particular, the microcanonical Gibbs distribution (1) is a stationary solution
of the evolution equation, and represents the state of thermodynamic
equilibrium of the system. If we start with a probability distribution $\rho$
defined on a surface of constant energy, then it will remain on this surface.
If the energy is equal to $E$, the probability density will eventually reach
the Gibbs distribution (1). In other words the system will relax to
equilibrium. This result can be demonstrated by showing that the entropy is a
nondecreasing function of time and that its maximum value is the one
corresponding to the Gibbs distribution (1). This will shown in the following.

###  III.3 Entropy production

Let us determine the time evolution of the average $U$ of ${\cal H}$.
Multiplying equation (82) by ${\cal H}$ and integrating in the phase space, we
arrive at the equation

$\frac{dU}{dt}=-\frac{1}{2}\sum_{(ij)}\int{\bf
v}_{ij}\cdot{\mathbb{K}}_{ij}\cdot{\bf D}_{ij}\rho\,dxdp,$ (85)

where we have performed an integration by parts and used the result ${\bf
D}_{ij}{\cal H}={\bf v}_{ij}$. Taking into account the orthogonality (62), it
follows that the right-hand side of (85) vanishes and $dU/dt=0$. That is, $U$
remains constant in time, which is expected because ${\cal H}$ is conserved by
the evolution equation.

We determine now the time evolution of the entropy $S$, defined by

$S=-k\int\rho\ln\rho\,dxdp.$ (86)

Multiplying the evolution equation (82) by $-\ln\rho$ and integrating in the
phase space, we find

$\frac{dS}{dt}=\frac{k}{2}\sum_{(ij)}\int\frac{1}{\rho}({\bf
D}_{ij}\rho)\cdot{\mathbb{K}}_{ij}\cdot({\bf D}_{ij}\rho),$ (87)

where we have performed an integration by parts. Replacing ${\mathbb{K}}_{ij}$
given by (64), we reach the following expression

$\frac{dS}{dt}=\frac{k}{2}\sum_{(ij)}\int\frac{\gamma_{ij}}{\rho}(A_{ij}^{2}+B_{ij}^{2})dxdp.$
(88)

We see that the right-hand side of (88), which we denote by $\Pi$, is positive
definite and should be interpreted as the rate of entropy production since
there is no exchange of energy with the environment, in fact no exchange of
heat since there is no external work involved. We remark that $\Pi$ vanishes
in equilibrium because each one of the terms $A_{ij}$ and $B_{ij}$ vanish when
$\rho$ is replaced by the equilibrium distribution $\rho_{e}$, given by (1).
Since $S$ is monotonically increasing in time because $dS/dt=\Pi\geq 0$, it
reaches a maximum when $\Pi=0$, that is, when the distribution is $\rho_{e}$.

##  IV Liouville equation

If we set $\gamma_{ij}$ equal to zero in the Fokker-Planck equation (82),
which means that no stochastic forces are involved, it reduces to the form

$\frac{\partial\rho}{\partial t}+\\{\rho,{\cal H}\\}=0,$ (89)

which is the Liouville equation. Here we present the usual derivation of the
Liouville equation, which is carried out by the use of the Liouville theorem.
We also present some of its properties in order to compare with those of the
Fokker-Planck equation (82). We consider a system described by a Hamiltonian
${\cal H}$, and the equations of motion are

$\frac{dx_{i}}{dt}=\frac{\partial{\cal H}}{\partial
p_{i}},\qquad\frac{dp_{i}}{dt}=-\frac{\partial{\cal H}}{\partial x_{i}}.$ (90)

If at time $t$ we consider a region ${\cal R}$ in phase space which is
transformed, through the Hamiltonian motion, in a region ${\cal R}^{\prime}$
at time $t^{\prime}$, then the Liouville theorem states that the volume of
these regions are equal, that is,

$\int_{\cal R}dxdp=\int_{\cal R^{\prime}}dx^{\prime}dp^{\prime}.$ (91)

We may also write this result as $dxdp=dx^{\prime}dp^{\prime}$ which means
that the Jacobian of the transformation $(x,p)\to(x^{\prime},p^{\prime})$
equals the unity.

In the usual derivation, the Liouville equation is obtained by postulating
that the probability of the region ${\cal R}$ at time $t$ is equal to that of
the region ${\cal R^{\prime}}$ at time $t^{\prime}$, that is,

$\int_{\cal R}\rho(x,p,t)dxdp=\int_{\cal
R^{\prime}}\rho(x^{\prime},p^{\prime},t^{\prime})dx^{\prime}dp^{\prime},$ (92)

or
$\rho(x,p,t)dxdp=\rho(x^{\prime},p^{\prime},t^{\prime})dx^{\prime}dp^{\prime}$.
From this equation and the Liouville theorem, we find

$\rho(x,p,t)=\rho(x^{\prime},p^{\prime},t^{\prime}),$ (93)

and the probability density is invariant along a trajectory. This is
equivalent to say that the material time derivative of $\rho$ vanishes, that
is,

$\frac{\partial\rho}{\partial t}+\sum_{i}\left(\frac{\partial\rho}{\partial
x_{i}}\frac{\partial x_{i}}{dt}+\frac{\partial\rho}{\partial
p_{i}}\frac{\partial p_{i}}{dt}\right)=0.$ (94)

Taking into account the Hamilton equations of motion, the second term becomes
the Poisson brackets $\\{\rho,{\cal H}\\}$ and we arrive at the Liouville
equation (89).

The most relevant property of the Liouville equation is that the
microcanonical Gibbs distribution (1) is a stationary distribution of the
equation. However, not all initial distribution will relax to this
distribution. To discuss this point we use the property coming from the
Liouville equation that the entropy is invariant in time.

The entropy $S(t)$ at time $t$ is defined by (86). Using the invariance of the
probability density along a trajectory, given by equation (93), we reach the
result $S(t)=S(t^{\prime})$, at distinct instants of time. Alternatively, we
may employ the procedure used above in the stochastic approach. We multiply
the Liouville equation (89) by $k\ln\rho$ and integrate in the phase space.
After an integration by parts we find that $dS/dt=0$, and the entropy is
invariant.

To proceed in our analysis it is convenient to define the measure of a surface
A of constant energy $E$ in phase space. We consider another surface A′ of
constant energy $E^{\prime}$, with $E^{\prime}$ differing little from $E$. The
measure of A, which we denote by $\Omega_{\rm A}$, is the volume between the
two surfaces A′ and A divided by $E^{\prime}-E$. Thus $\Omega$ appearing in
equation (1) is the measure of the whole surface of constant energy $E$.

Let us consider a probability distribution, denoted by $\rho_{\rm A}$, which
is nonzero only in a surface $A$ which is a subset of the whole surface of
constant energy $E$. Within $A$ it is given by

$\rho_{\rm A}=\frac{1}{\Omega}_{\rm A}\delta(E-{\cal H}).$ (95)

where $\Omega_{\rm A}$ is the measure of A. The corresponding entropy is
$S_{\rm A}=k\ln\Omega_{\rm A}$.

Suppose that the initial condition to the Liouville equation is $\rho_{\rm
A}$. As time evolves, the surface $A$ changes to another surface A′ and become
intertwined with the complementary surface to the whole surface. However, by
the Liouville theorem, the measure of A′ is the same as that of A, and the
entropy remains constant and equal to $S_{\rm A}$. Therefore the Gibbs
distribution will never be reached because the entropy of this distribution is
$S_{e}=k\ln\Omega$ which is larger than $S_{\rm A}$.

We present now an example where the initial distribution is nonzero and given
by

$\rho=a\rho_{\rm A}+b\rho_{\rm B},$ (96)

where $a+b=1$, and the surfaces A and B are complementary surfaces and make up
the whole surface of constant energy $E$. The corresponding entropy is

$S=ka\ln\frac{\Omega_{\rm A}}{a}+kb\ln\frac{\Omega_{\rm B}}{b},$ (97)

and it differs from $S_{e}=k\ln\Omega$, where $\Omega=\Omega_{\rm
A}+\Omega_{\rm B}$. In fact we can show that $S_{e}\geq S$, the equality
occurring when $\rho$ is the Gibbs distribution in which case $a=\Omega_{\rm
A}/\Omega$ and $b=\Omega_{\rm B}/\Omega$. Therefore, if the initial
distribution is of the type (96) but is not the Gibbs distribution, this Gibbs
distribution will never be reached.

##  V Conclusion

We have proposed a stochastic equation of motion that emerges as a consequence
of stochastic forces acting on each pair of particles. The random forces
acting on a pair of particles are action and reaction leading to the
conservation of momentum. These forces are perpendicular to the relative
velocity between the particles, performing therefore zero work on the two
particles, and thus preserving their kinetic energy.

These stochastic forces and the conservative forces leads to the conservation
of energy along a trajectory in phase space. The trajectories in phase space
are restricted to the surface of constant energy and is thus similar in this
aspect to the Hamiltonian flow. From the stochastic motion, we have derived
the evolution equation for the density which turns out to be a Fokker-Planck
equation, but distinct from the usual Fokker-Planck equation that describes
the contact of system with a heat reservoir.

To derive the evolution equation we introduced an appropriate interpretation
of the stochastic equation of motion that does not correspond to the
interpretation proposed by Itô nor that proposed by Stratonovich. The
interpretation proposed here makes the energy to be conserved at all times, a
result that was demonstrated by showing that the Hamiltonian function is a
stationary solution of the evolution equation.

In contrast to the Liouville equation, which also describes a closed system
and conserves the energy, the evolution equation that we have set up predicts
the increase of entropy and the relaxation to equilibrium. In other terms, the
equation describes an irreversible decay to equilibrium but we cannot acertain
whether the present approach could describe correctly the actual decay
observed in real closed systems.

## References

  * (1) J. W. Gibbs, Elementary Principles in Statistical Mechanics (Scribner, New York, 1902).
  * (2) E. Borel, Mécanique Statistique Classique (Gauthier-Villars, Paris, 1925).
  * (3) R. C. Tolman, The Principles of Statistical Mechanics (Oxford University Press, Oxford, 1938).
  * (4) A. I. Khinchin, Mathematical Foundations of Statistical Mechanics (Dover, New York, 1949).
  * (5) L.D. Landau and E.M. Lifshitz, Statistical Physics (Pergamon Press, London, 1958).
  * (6) I. Prigogine, Non-Equilibrium statistical (Interscience, New York, 1962).
  * (7) A. Hobson, Concepts in Statistical Mechanics (Gordon and Breach, New York, 1971).
  * (8) D. N. Zubarev, Nonequilibrium Statistical Thermodynamics (Consultants Bureau, New York, 1974).
  * (9) V. I. Arnold, Mathematical Methods of Classical Mechanics (Springer, 1978).
  * (10) L. E. Reichl, A Modern Course in Statistical Physics (University of Texas Press, Austin, 1980).
  * (11) S. R. Salinas, Introduction to Statistical Physics (Springer, New York, 2001).
  * (12) K. Huang, Statistical Mechanics (Wiley, New York, 1963).
  * (13) Y. B. Rumer and M. S. Ryvkin, Thermodynamics, Statistical Physics and Kinetics (Mir, Moscow, 1980).
  * (14) C. Cercignani, The Boltzmann Equation and its Applications (Springer, New York, 1988).
  * (15) G. M. Kremer, An Introduction to the Boltzmann Equation and Transport Processes in Gases (Springer, Berlin, 2010).
  * (16) M. J. de Oliveira, Phys. Rev. E 99, 052138 (2019).
  * (17) N. G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam, 1981).
  * (18) H. Risken, The Fokker-Planck Equation (Springer, Berlin, 1989).
  * (19) C. Gardiner, Stochastic Methods (Springer, Berlin, 2009).
  * (20) T. Tomé and M. J. de Oliveira, Stochastic Dynamics and Irreversibility (Springer, Cham, 2015).
  * (21) T. Tomé and M. J. de Oliveira, Phys. Rev. E 82, 021120 (2010).
  * (22) C. Van de Broeck and M. Esposito, Phys. Rev. E 82, 011144 (2010).
  * (23) R. E. Spinney and I. J. Ford, Phys. Rev. E 85, 051113 (2012)
  * (24) U. Seifert, Rep. Prog. Phys. 75, 126001 (2012).
  * (25) M. Santillan and H. Qian, Physica A 392, 123 (2013).
  * (26) D. Luposchainsky and H. Hinrichsen, J. Stat. Phys. 153, 828 (2013).
  * (27) T. Tomé, M. J. de Oliveira, Phys. Rev. E 91, 042140 (2015).
  * (28) M. J. de Oliveira, Phys. Rev. E 102, 032114 (2020).
  * (29) M. J. de Oliveira, Revista Brasileira de Ensino de Física 42, e20200210 (2020).
  * (30) P. M. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1987).
  * (31) M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation (Oxford University Press, Oxford, 2010).