# Dynamical mechanisms leading to equilibration in two-component gases

Stephan de Bievre, Carlos Mejia-Monasterio, Paul E. Parris

► To cite this version:

Stephan de Bievre, Carlos Mejia-Monasterio, Paul E. Parris. Dynamical mechanisms leading to equilibration in two-component gases. Physical Review E, 2016, 10.1103/PhysRevE.93.050103. hal-01408060

HAL Id: hal-01408060

https://inria.hal.science/hal-01408060

Submitted on 22 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
This is the accepted manuscript made available via CHORUS. The article has been published as:

# Dynamical mechanisms leading to equilibration in two-component gases

Stephan De Bièvre, Carlos Mejía-Monasterio, and Paul E. Parris

Phys. Rev. E 93, 050103 — Published 23 May 2016

DOI: 10.1103/PhysRevE.93.050103
Dynamical mechanisms leading to equilibration in two-component gases

Stephan De Bièvre*

Laboratoire Paul Painlevé, CNRS, et UFR de Mathématiques,
Université Lille 1; Équipe-Projet Mephysto, INRIA Lille-Nord Europe, France

Carlos Mejía-Monasterio†

Laboratory of Physical Properties, Technical University of Madrid, Av. Complutense s/n 28040 Madrid, Spain

Paul E. Parris‡

Department of Physics, Missouri University of Science & Technology, Rolla, MO 65409, USA

Demonstrating how microscopic dynamics cause large systems to approach thermal equilibrium remains an elusive, longstanding, and actively-pursued goal of statistical mechanics. We identify here a dynamical mechanism for thermalization in a general class of two-component dynamical Lorentz gases, and prove that each component, even when maintained in a non-equilibrium state itself, can drive the other to a thermal state with a well-defined effective temperature.

PACS numbers: 05.70.-a, 02.50.Ey, 05.20.-y

That isolated systems with many degrees of freedom evolve asymptotically in time towards thermal equilibrium lies at the heart of classical thermodynamics. Statistical mechanics teaches that for systems described by a Hamiltonian $H$, the thermal states are those described by the canonical Boltzmann relation $\rho = Z^{-1} \exp(-\beta H)$. This follows from original arguments of Maxwell, marginal distributions that arise from microcanonical ensembles, and the properties of the maximum entropy states to which systems thermodynamically tend [1].

While such statistical arguments identify the thermal state, they provide no insight into the problem of how the microscopic dynamics of diverse large systems each lead towards equilibrium from an arbitrary initial state.

Considerable progress has been made in understanding *return to equilibrium* [2] in which a small system coupled to a large thermal reservoir thermalizes to the temperature of the latter. The more general problem of *approach to equilibrium*, in which many mutually interacting elements of an isolated system are initially out of equilibrium, is less understood [3]. In spite of recent progress, including identification of the thermalization mechanism in the Fermi-Pasta-Ulam problem [4], there is continuing discussion and even some controversy regarding the role played by various properties of the dynamics, such as chaos, mixing, and resonances [5, 6], as well as on the emergence of non-equilibrium effective temperatures [7].

As an important step in this longstanding problem, we identify here a dynamical mechanism for thermalization in a general class of two component systems [8–16], in which non-interacting point particles move freely through a spatially fixed array of isolated dynamical “scatterer”

particles, each of which has a few rotationally-invariant degrees of freedom with which the itinerant point particles interact and exchange energy when they are within a fixed range.

In particular, we prove here that, starting from an arbitrary initial state, an ensemble of probe subsystems (either the itinerant point particles or the “scatterer” particles) when subjected to repeated weak interactions with members of the complementary component gas, acting as a homogeneous and stationary energy reservoir, approaches a Boltzmann state at a well-defined temperature, provided merely that the reservoir is stationary–it need not be in thermal equilibrium.

While providing an understanding of the approach to equilibrium of the entire system as a whole, the result also provides insight into and allows a justification of the concept of non-equilibrium effective temperatures.

In systems of this type studied previously, the scatterers were taken to be disks [8–11, 17] or needles in 2D [12] that rotate about fixed centers, and harmonic oscillators [13–16]. In each of these systems, any particle or scatterer in the gas follows a non-interacting Newtonian evolution during time intervals separating the collisions it experiences. In each collision an energy conserving interaction occurs between the particle and scatterer involved. Such systems, in which the scatterers themselves possess an internal dynamics, have been referred to as dynamical Lorentz gases [13], since they provide a natural generalization of the two-component gases originally introduced by Lorentz [18], and subsequently studied by many others [19], that feature *inert* (i.e., non-dynamical) scatterers, and which for over a century have been an essential tool for understanding diffusion, related equilibrium and non-equilibrium statistical properties, and the so-called Boltzmann-Grad limit.

We focus here on thermalization in a specific two-dimensional system of this type, introduced in [8] as the rotating Lorentz gas model (RLG), in which the

<footnote>* stephan.de-bievre@univ-lille1.fr</footnote>

<footnote>† carlos.mejia@upm.es</footnote>

<footnote>‡ parris@mst.edu</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) A single rotating disk in a gas of point particles, and (b) a single point particle passing through a spatially fixed array of rotating disks.</figcaption></figure>

scatterers are rotating disks. The full many particle RLG exhibits realistic equilibrium and non-equilibrium behaviour [9], and has been used to study thermal rectification [10] and thermoelectricity [11]. Its time-reversible dynamics preserves phase space volume.

To understand equilibration of such RLGs as a whole, we investigate here the effect of repeated interactions on the phase space distribution of a single member of each of the two components making up the RLG. In the context of non-equilibrium temperatures, one might think of an ensemble of such single members as a “thermometer” locally probing the phase space distribution of the other component.

Thus, we first consider a single hollow disk that rotates freely about its center immersed in a gas of (and subjected to repeated impacts by) non-interacting point particles (Fig. 1a), with the particles of the gas reservoir being “probed” drawn independently from a stationary, but not necessarily thermal distribution.

This obviously ignores recollisions, which would modify the momentum distribution of gas particles in the full RLG as they repeatedly encounter different scatterers. In the second part of our analysis, therefore, we follow a single particle as it passes through (i.e., probes) a spatially-fixed gas of rotating disks (Fig. 1b), where now it is the disks of the scatterer reservoir that are drawn independently from an arbitrary stationary distribution.

As we prove, in either of these situations, the probe dynamics reduce to a Markov chain that leads in the limit of small average energy exchange to an approach of the corresponding probe energy distribution to a Boltzmann state with a well-defined effective temperature.

In either case, each step of the Markov chain involves a collision (Fig. 2) between a particle of mass $m = 1$ and initial momentum $p$ impinging with impact parameter $b$ upon a hollow disk of unit radius and mass $\mu = M/m$, rotating about its fixed center with initial angular velocity $\omega$. During the collision, particle and disk obtain new momentum and angular velocity [8]

$$p' = p - 2(p \cdot u)u - \frac{2\mu}{1 + \mu} (p_t - \omega) u_{\perp}, \quad (1)$$

$$\omega' = \omega + \frac{2}{1 + \mu} (p_t - \omega). \quad (2)$$

Here, $\mu$ represents both the mass of the disk and its mo-

<figure><img src="image_3.png" /><figcaption>FIG. 2. Collision geometry for a point particle of unit mass and initial momentum $p$ impinging on a freely rotating disk of mass $\mu$ and initial angular velocity $\omega$.</figcaption></figure>

ment of inertia, the unit vector $u$ links the disk center to the collision point, $u_{\perp}$ is obtained by rotating $u$ through $\pi/2$, and $p_t = p \cdot u_{\perp} = b\|p\|$. Under these rules, particle and disk exchange energy and angular momentum, conserving both.

**One rotating disk in a gas of particles** – Consider a single disk (the “thermometer”) with initial angular velocity $\omega_0$, subject to impacts by point particles (Fig. 1a) with i.i.d. momenta $p_n$ drawn from a rotationally invariant distribution $\rho_{\text{eff}}(|p|)$, with i.i.d. impact parameters $b_n$ uniform in $[-1, 1]$, at a Poissonian sequence $t_n$ of impact times. Note, the marginal distribution $\rho_{\text{eff}}(|p|) \propto \rho_0(|p|)|p|$ of particles striking the disk, related to the “effusive” flux $j(p) = \rho_0(|p|)p$ of particles on the disk, is different from the bulk gas distribution $\rho_0(|p|)$.

Under these conditions $\langle p_t \rangle = 0$. After $n$ collisions (2) gives

$$\omega_n = \gamma^n \omega_0 + \eta \sum_{k=0}^{n-1} \gamma^k p_{t,n-1-k}, \quad (3)$$

where $\eta = 2/(1+\mu)$ and $\gamma = 1-\eta \in (-1, 1)$. As $n \to \infty$, the first term in (3) vanishes, leaving a convergent sum

$$\omega = \lim_{n \to +\infty} \eta \sum_{k=0}^{n-1} \gamma^k p_{t,n-1-k} \quad (4)$$

of geometrically weighted i.i.d. random variables $p_t$, ensuring that the angular velocity distribution approaches a stationary limit $\rho(\omega)$.

When the particle gas is in thermal equilibrium, the $p_t$ in (4) are Gaussian. Hence, independent of the value $\mu$, the distribution $\rho(\omega)$ is thermal and has the same temperature as the bath (see below). Thus, in contact with thermal particles, a probe disk approaches a common thermal equilibrium with the particle bath.

Less obvious and more interesting is when the particle reservoir is not thermal [12]. Here, we characterize the limiting distribution by its low order moments. Assuming $\langle \omega_0 \rangle = 0$, then $\langle \omega_n \rangle = 0$, and

$$\langle \omega_n^2 \rangle = \gamma^{2n} \langle \omega_0^2 \rangle + \eta^2 \langle p_t^2 \rangle \left( \frac{1 - \gamma^{2n}}{1 - \gamma^2} \right). \quad (5)$$

As $n \to \infty$, the average disk energy then approaches

$$\varepsilon_d = \frac{\mu}{2} \langle \omega^2 \rangle = \frac{1}{2} \langle p_t^2 \rangle = \frac{1}{6} \langle p^2 \rangle_{\text{eff}}. \quad (6)$$
<figure><img src="image_4.png" /><figcaption>FIG. 3. Limiting value of $\kappa_\omega$ of a disk of mass $\mu$ subjected to repeated impacts by particles drawn from thermal (x), uniform (□), and microcanonical (+) distributions, and a superposition of two thermal distributions (○) with temperatures in ratio $T_1/T_2 = 5$. Curves are theoretical predictions, symbols are simulation results. Inset: Limiting value of kinetic energy ratio $\epsilon = \epsilon_d/\epsilon_p$ vs $\mu$. Theoretical predictions are $\epsilon = 1/2, \epsilon = 2/5$, and $\epsilon = 1/3$.</figcaption></figure>

On the right, the average $\langle b^2 \rangle = 1/3$ has been performed, and the remaining average is over $\rho_{\text{eff}}(|p|)$. To relate this to the particle energy $\epsilon_p = \langle p^2 \rangle_0/2$ [where $\langle \dots \rangle_0$ indicates averages over $\rho_0(|p|)$], we introduce an energy partitioning ratio $\epsilon = \epsilon_d/\epsilon_p$, which equals $1/2$ in equilibrium, satisfying the equipartition theorem. For a disk immersed in a non-equilibrium particle gas, however, $\epsilon = \langle p^2 \rangle_{\text{eff}}/\langle 3p^2 \rangle_0$ approaches a non-universal value that depends on the particle distribution, but is independent of the mass of the disk (see Fig. 3). For a thermal particle bath, $\epsilon = 1/2$ and (6) reduces to the equilibrium result $\mu\langle\omega^2\rangle/2 = k_B T/2$.

To study the shape of $\rho(\omega)$ we compute the asymptotic value of the excess disk kurtosis

$$ \kappa_{\omega} = \beta_2 - 3 = \eta \frac{1 + \gamma}{1 + \gamma^2} \left( \frac{\langle p_t^4 \rangle}{\langle p_t^2 \rangle^2} - 3 \right), \quad (7) $$

in which $\beta_2 = \langle \omega^4 \rangle / \langle \omega^2 \rangle^2$, and which vanishes when the limiting disk distribution is thermal. This is the case for a thermal particle bath, when the excess particle kurtosis $\kappa_p = (\langle p_t^4 \rangle / \langle p_t^2 \rangle^2 - 3)$ in (7) vanishes.

But $\rho(\omega)$ also reduces to a thermal distribution when the particle reservoir is not thermal, whenever the coupling is sufficiently weak. This occurs when the disk is very heavy ($\mu \gg 1, \eta \to 0, \gamma \to 1$), and very light ($\mu \ll 1, \eta \to 2, \gamma \to -1$). In these two limits (7) vanishes and $\rho(\omega)$ becomes Gaussian [20]. In weak coupling, this thermalization is universal, i.e., independent of the form of $\rho_{\text{eff}}(|p|)$.

Thus, when placed weakly in contact with a stationary, non-equilibrium particle gas, the probe disk “thermalizes to” (or “measures”) an apparent temperature such that

<figure><img src="image_5.png" /><figcaption>FIG. 4. Limiting values of the kinetic energy ratio $\epsilon = \epsilon_d/\epsilon_p$ and excess kurtosis $\kappa_p$ for a particle colliding with rotating disks of mass $\mu$ and angular velocities drawn from thermal (x), uniform (□), and microcanonical (+) distributions, and a superposition of two thermal distributions (○) with $T_1/T_2 = 5$.</figcaption></figure>

$k_B T/2 = \epsilon \langle p^2 \rangle_0 / 2$, where $\epsilon = \langle p^2 \rangle_{\text{eff}} / \langle 3p^2 \rangle_0$ depends on the particle distribution and is generally not equal to $1/2$. This thermalization for large and small $\mu$ is clearly seen in Fig. 3, which displays a numerical computation of the values of $\epsilon$ and $\kappa_\omega$, for $10^5$ disks, each subjected to $10^4$ repeated collisions for various particle distributions.

From Fig. 3 and Eq. (7) it is clear that for intermediate and strong coupling, with particle reservoirs for which $\kappa_p \neq 0$, the disk does not approach a thermal state. Instead, repeated strong interactions drive it to a non-thermal state with $\kappa_\omega \neq 0$.

**One particle in a gas of disks** – We now consider a single particle that collides with a sequence of rotating disks (Fig. 1b) whose angular velocities are i.i.d. variables drawn from a stationary distribution $\rho(|\omega|)$, again ignoring recollisions. Denote by $p_s$ the particle’s momentum before collision $s$, when it impinges with impact parameter $b_s$ on a disk with angular velocity $\omega_s$. According to (1), after this collision

$$ p_{s+1} = p_s - 2(p_s \cdot u_s)u_s - \frac{2\mu}{1+\mu} (p_{s,t} - \omega_s) u_{s,\perp}. \quad (8) $$

Collisions occur at times $t_{s+1} = t_s + \ell/p_s$, where $\ell$ is the collision mean free path. We assume $\langle p_0 \rangle = 0$.

The Markov chain (8) for $p$ is less tractable than for a single disk in a gas of particles. In weak coupling, however, analysis of the Markov chain to determine the limiting distribution is straightforward.

Figures 4 and 5 display simulation results in which $10^5$ particles each undergo a sequence of collisions with rotating disks initialized as described. In Fig. 4, values of $\epsilon$ and $\kappa_p$ for the limiting particle distribution, recorded
<figure><img src="image_6.png" /><figcaption>FIG. 5. Distribution function for velocity component $p_x$ for a particle undergoing repeated collisions with rotating disks of mass $\mu = 0.005$ in a microcanonical ensemble with $|\omega| = \sqrt{20}$, recorded (from broad to narrow) at times $t = 60, 300, 600, 2000$. Dashed curve corresponds to the limiting thermal distribution. Inset: convergence of the particle kurtosis $\beta_2^{(p)}(\sigma)$ to its limiting value 3 (dashed line), and kinetic energy ratio $\epsilon = \epsilon_d/\epsilon_p$ (□) to its limiting value 1/2 (solid line).</figcaption></figure>

at long fixed simulation time, are plotted for different disk distributions $\rho(|\omega|)$.

For thermal disks, the particle equilibrates to a thermal state with the same temperature. Moreover, when $\mu$ is very large or very small, independent of $\rho(\omega)$, the particle is also driven to a thermal state (see, e.g., Fig. 5), with a vanishing $\kappa_p$. Thus, we again observe “thermalization” of a small system in weak contact with a non-thermal reservoir. However, unlike what happens to a disk in a gas of particles, the limiting particle temperature always obeys equipartition, since $\epsilon \to 1/2$ in this regime. For intermediate $\mu$, provided $\kappa_\omega \neq 0$, the particle does not generally thermalize.

To analytically demonstrate the thermalization observed at large and small $\mu$, we express using (8) the change in particle energy

$$E_{s+1} - E_s = -\frac{4\mu}{(1+\mu)^2} E_s b_s^2 + \frac{2\mu}{(1+\mu)^2} \mu \omega_s^2 + 2\sqrt{\mu} \frac{1-\mu}{(1+\mu)^2} \sqrt{2E_s} \sqrt{\mu \omega_s} b_s \quad (9)$$

during collision $s$, in terms of its energy immediately before. The coefficients multiplying the dynamical variables in (9) set the scale for the energy change in any collision, and are small for $\mu \ll 1$ and $\mu \gg 1$. In the latter case the equation for large $\mu$ follows from that for small $\mu$ by replacing $\mu$ with $\mu^{-1}$.

It suffices to study this weak-coupling, small-step limit for $\mu \ll 1$. In this limit, (9) becomes

$$\Delta E_s = -4\mu E_s b_s^2 + 4\mu \epsilon_d k_s^2 + 2\sqrt{4\mu} \sqrt{\epsilon_d E_s} k_s b_s, \quad (10)$$

where

$$\epsilon_d = \frac{1}{2}\mu\langle\omega_s^2\rangle, \quad k_s = \sqrt{\frac{2\mu}{\epsilon_d}}\omega_s, \quad \langle k_s^2 \rangle = 1. \quad (11)$$

Averages here are over $\rho(|\omega|)$. The first two terms on the right hand side of (10) constitute a source of “dynamical friction” [21, 22] that counterbalances the stochastic acceleration caused by the fluctuating last term [13, 16, 21, 22]. This competition thus leads to a kind of fluctuation-dissipation like mechanism that naturally emerges from the deterministic dynamics.

Inspection of (10) suggests $\xi = 4\mu$ as a small parameter. Introducing (scaled) continuous collision number $\sigma = \xi s$, Eq. (10) can be described by the stochastic differential equation

$$dE(\sigma) = -\alpha(E)d\sigma + \lambda(E)dw(\sigma)$$

in which

$$\alpha(E) = \frac{1}{3}E(\sigma) - \epsilon_d, \quad \lambda(E) = \frac{2}{\sqrt{3}}\sqrt{\epsilon_d E(\sigma)}.$$

Here $dw$ is white noise, $\langle dw \rangle = 0$, and $\langle dw(\sigma)dw(\sigma') \rangle = \delta(\sigma - \sigma')$. Clearly, the results will depend on only the first two moments of $\rho(|\omega|)$.

From the corresponding Fokker-Planck equation, one finds the stationary distribution

$$f_{\infty}(E) = \tilde{Z}^{-1} \lambda^{-2}(E) \exp\left[ - \int 2\alpha(E) \lambda^{-2}(E) dE \right] \\ = \frac{1}{Z'} E^{1/2} \exp\left(-\frac{E}{2\epsilon_d}\right)$$

for $E(\sigma)$ at large $\sigma$. Physically, $f(E, \sigma)dE$ is the fraction of particles after collision $\sigma$ with energy between $E$ and $E + dE$. Such particles stay in that state for time $\Delta t = l/\sqrt{2E}$. Thus, at long times, the distribution of particles with energy $E$ becomes thermal

$$f(E) = \frac{1}{Z} \exp\left(-\frac{E}{2\epsilon_d}\right) \quad (12)$$

as in a 2D non-interacting particle gas with apparent temperature given by $k_B T = 2\epsilon_d$. Thus, rotating disks not in equilibrium drive the particle to a thermal state compatible with equipartition.

**On the approach to equilibrium** – In weak coupling, as we have shown, both components of the full RLG are driven to a Gaussian thermal state. From (8), and (10), it is clear that equilibration is not a simple consequence of the central limit theorem. Indeed, the momentum or energy increments that particles experience are not i.i.d. random variables; they are steps in an associated Markov chain, with the size and variance of each being a function of the dynamical variables at each step. Thus, the limiting exponential distribution for the particle energy in this case can be understood as arising from a competition between impulsive fluctuating forces of zero mean, which tend to heat the particle, and a dynamical friction that
leads higher energy particles to lose energy to the disks. Thus, a fluctuation-dissipation mechanism emerges naturally from the deterministic dynamics associated with the interaction between the particles and the rotating disks.

The fact that in the present case both components of the full RLG equilibrate in the presence of a non-thermal bath allows us to identify the mechanism of approach to thermal equilibrium of the RLG as a whole.

At low particle density, e.g., each particle will scatter off many disks before any disk is likely to have interacted with more than a few particles. In weak coupling, therefore, the particle gas will equilibrate well before the gas of disks, to a thermal distribution in which equipartition of the particle and disk energies is obtained. Each disk will then be in contact with a thermal distribution of particles, and will only need therefore to undergo a *return* to the appropriate limiting thermal distribution. Beyond weak coupling, the situation becomes more complicated [12]. Moreover, extensive numerical simulations on large but finite systems with many particles and disks, to be reported elsewhere [12], confirm this basic picture of a

two-stage equilibration process of the whole system in the weak coupling regime.

In summary, we have identified a thermalization mechanism in a rotating disk, weak-coupling version of a two-component dynamical Lorentz gas. We expect the basic underlying dynamical friction mechanism to be effective more generally in systems in which the individual components undergo repeated scattering events (See, for example Ref. 13)

## ACKNOWLEDGMENTS

Part of this work was performed while C.M.-M. and P.E.P. visited the Université Lille 1 and the Labex CEMPI (ANR-11-LABX-0007-01). They thank those institutions for their hospitality. C.M.-M. acknowledges partial financial support from the Spanish MICINN grants MTM2012-39101-C02-01 and MTM2015-63914-P, and from ONRG Grant N62909-15-1-C076.

[1] A.I. Khinchin, *Mathematical Foundations of Statistical Mechanics*, (Dover, New York, 1949); F. Reif, *Fundamentals of Statistical and Thermal Physics*, (McGraw Hill, New York, 1965).

[2] V. Bach, J. Fröhlich and I. Sigal, J. Math. Phys. **41**, 3985-4060 (2000); V. Jakšić and C. Pillet, Acta Math. **181**, 245-282 (1998).

[3] E. Fermi, J. Pasta, and S. Ulam, *Studies of nonlinear problems I*, Los Alamos Scientific Laboratory Report, Los Alamos Scientific Laboratories, Los Alamos, New Mexico (1955); G.P. Berman and F.D. Izrailev, Chaos **15**, 015104 (2005).

[4] M. Onorato, L. Vozella, D. Proment, and Y. V. Lvov, Proc. Nat. Acad. Sci. U.S.A., **112**, 42084213 (2015).

[5] J. Casas-Vázquez and D. Jou, Rep. Prog. Phys., **66**, 1937-2023 (2003).

[6] A. Polkovnikov, K. Sengupta, A. Silva, M. Vengalattore Rev. Mod. Phys., **83**, 863 (2011).

[7] E. Dieterich et al., Nature Physics, **11**, 971-977 (2015).

[8] C. Mejía-Monasterio, H. Larralde, and F. Leyvraz, Phys. Rev. Lett. **86**, 5417 (2001).

[9] H. Larralde, F. Leyvraz, and C. Mejía-Monasterio, J. Stat. Phys. **113**, 197 (2003).

[10] J.P. Eckmann and C. Mejía-Monasterio, Phys. Rev. Lett. **97**, 094301 (2006).

[11] G. Casati, C. Mejía-Monasterio, and T. Prosen, Phys. Rev. Lett. **101**, 016601 (2008).

[12] S. De Bièvre, C. Mejía-Monasterio, and P.E. Parris, to be published.

[13] S. De Bièvre and P.E. Parris, J. Stat. Phys. **142**, 356385 (2011).

[14] A.A. Silvius, P.E. Parris, and S. De Bièvre, Phys. Rev. B **73**, 014304 (2006).

[15] S. De Bièvre, P.E. Parris, and P. Lafitte, J. Stat. Phys. **132**, 863-879 (2008).

[16] B. Aguer, S. De Bièvre, P. Lafitte, and P.E. Parris, J. Stat. Phys. **138**, 780814 (2010).

[17] R. Klages, K. Rateitschak, and G. Nicolis, J. Stat. Phys. **99**, 1339 (2000).

[18] A. Lorentz, in *KNAW, Proceedings*, **7**, 438 (Amsterdam, 1905).

[19] Carl P. Dettmann, Commun. Theor. Phys. **62**, 521-540 (2014).

[20] Li-Xin Zhang, Ann. Probab. **25**, 1621-1635 (1997).

[21] S. Chandrasekhar, Astrophys. J. **97**, 255262,(1943).

[22] S. Chandrasekhar, Astrophys. J. **97**, 263273, (1943).