# Collisional statistics of the hard-sphere gas

Paolo Visco, Frédéric van Wijland, Emmanuel Trizac

► To cite this version:

Paolo Visco, Frédéric van Wijland, Emmanuel Trizac. Collisional statistics of the hard-sphere gas. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics [2001-2015], 2008, 77 (4), pp.041117. 10.1103/PhysRevE.77.041117. hal-00303734

HAL Id: hal-00303734

https://hal.science/hal-00303734v1

Submitted on 16 May 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 CC BY 4.0 - Attribution - International License
Collisional statistics of the hard-sphere gas

Paolo Visco,¹,² Frédéric van Wijland,³ and Emmanuel Trizac¹

¹Université Paris-Sud, LPTMS, UMR 8626, Orsay Cedex, F-91405 and CNRS, Orsay, F-91405

²SUPA, School of Physics, University of Edinburgh, Mayfield Road, Edinburgh, EH9 3JZ, UK

³Laboratoire Matière et Systèmes Complexes (CNRS UMR 7057), Université Denis Diderot (Paris VII),
10 rue Alice Domon et Léonie Duquet, 75205 Paris cedex 13, France

(Dated: October 22, 2018)

We investigate the probability distribution function of the free flight time and of the number of collisions in a hard sphere gas at equilibrium. At variance with naive expectation, the latter quantity does not follow Poissonian statistics, even in the dilute limit which is the focus of the present analysis. The corresponding deviations are addressed both numerically and analytically. In writing an equation for the generating function of the cumulants of the number of collisions, we came across a perfect mapping between our problem and a previously introduced model: the probabilistic ballistic annihilation process [Coppex et al., Phys. Rev. E 69 11303 (2004)]. We exploit this analogy to construct a Monte-Carlo like algorithm able to investigate the asymptotically large time behavior of the collisional statistics within a reasonable computational time. In addition, our predictions are confronted against the results of Molecular Dynamics simulations and Direct Simulation Monte Carlo technique. An excellent agreement is reported.

PACS numbers: 05.40.-a, 51.10.+y, 02.50.-r

# I. INTRODUCTION

Despite the great interest that the hard sphere gas has triggered since the early days of statistical physics, there are still, to this day, simple questions which have not been deeply investigated yet. In this paper we address one such question: what is the probability that a tagged particle suffers a given number of collisions in a time $t$? Apart from the simplicity of the question, relevant in its own right, the collisional statistics of the hard sphere gas turns out to be useful in several purposes, such as the estimation of transport coefficients or the characterization of the dissipation in the closely related granular gases.

In a low density hard sphere gas at equilibrium, the velocities of two colliding particles are uncorrelated just before the collision: this is the molecular chaos (strosszahlansatz) statement. This remark could naively lead to the conclusion that collisions are uncorrelated random events, implying therefore that the number of collisions is simply a Poisson random variable. However, even if molecular chaos is exactly verified, the collision number is *not a Poisson random variable*. The non-Poissonian nature of collisional statistics has already been noticed in the literature (see e.g. [1]) but to our knowledge, has resisted analytical investigations. It is our purpose here to fill this gap.

The reason for the non-Poissonian behavior alluded to above is that the probability for a collision to take place depends on the scattering cross section of the colliding pair (which depends on the relative velocity $g$ of the pair). Namely, for the hard sphere gas, the probability of a collision behaves roughly as $g \cdot \hat{\sigma}\Theta(g \cdot \hat{\sigma})$, where $\hat{\sigma}$ is a unitary vector joining the centers of mass of the two particles at contact (directed along the apse line), and $\Theta$ is the Heaviside step function. For a gas made up of particles interacting through a potential other than the hard-sphere potential, the probability of having a “collision” is of course different, which affects the distribution of the number of collisions. Of particular inter-

est is the gas of particle interacting through a pair potential $V(r) \sim 1/r^{2(d-1)}$, where $d$ is the space dimension. Such particles are known as Maxwell molecules [2] and lead to a velocity independent probability of having a collision. Within this model, it then appears that the collisions are truly uncorrelated random events, and hence that the number of collision is a Poisson random variable.

In the next section, the free flight time distribution of a hard sphere gas is addressed, together with the distribution of free path lengths. These results pertain to $N = 0$ collision properties. In section III, the analysis is extended to consider the full probability $P(N, t)$ encoding the number of collisions statistics. The large time behavior and a complementary perturbative treatment are worked out analytically. Explicit expressions are obtained for the cumulants of the number of collisions $N$. In order to put the theoretical predictions to the test, three types of numerical simulations are performed. The first two, Molecular Dynamics and Direct Simulation Monte Carlo (DSMC) are routinely employed in the field and beyond. The third type, of the Monte Carlo family, is discussed in IV, and specifically devised to solve the particular problem at hand. It relies on a reinterpretation of the eigenvalue equation for the generating function for the cumulants of $N$, in terms of a population dynamics with cloning and annihilation. A Markov chain with the desired properties in then constructed, which allows for a direct measure of several key quantities involved in the analytical treatment. In doing so, we uncover a fruitful mapping with the probabilistic ballistic annihilation model proposed in Ref [3]. The three numerical methods provide results that are in excellent agreement with the analytical predictions. Preliminary accounts of part of this work has appeared in [4], where the numerical aspect of the work was limited to its Molecular Dynamics content, and where the theory is restricted to the zeroth order of the treatment put forward here.
## II. FREE FLIGHTS TIME DISTRIBUTION

### A. General results

We consider a hard sphere gas in $d$ dimensions composed of $N$ particles of equal mass $m$ and equal diameter $\sigma$. The gas is thermalized at some temperature $T_0$ in a homogeneous state of (constant) density $\rho$, so that the one point distribution function of the system is a Gaussian:

$$ \phi(\mathbf{r}, \mathbf{v}) \equiv \phi(\mathbf{v}) = \rho(2\pi T_0)^{-d/2} \exp\left(-\frac{v^2}{2T_0}\right) . \quad (1) $$

If one decides to follow the evolution of one particle among the $N$, then the velocity probability distribution function (pdf) of that particle will evolve according to the homogeneous linear Boltzmann equation:

$$ \frac{\partial}{\partial t} f(\mathbf{v}_1, t) = \frac{1}{\ell} \int d\mathbf{v}_2 \int d\hat{\sigma} \, \mathbf{v}_{12} \cdot \hat{\sigma} \, [f_1^{**} \phi_2^{**} - f_1 \phi_2] , \quad (2) $$

where $\ell = (\sigma^{d-1}\rho\chi)^{-1}$ is a length proportional to the mean free path denoted $l$ below, and $\mathbf{v}_{12} = \mathbf{v}_1 - \mathbf{v}_2$. The factor $\chi$ is the pair correlation function at contact and is the so called Enskog correction factor [5]. In the dilute limit where the molecular chaos assumption is justified, it tends to unity. Moreover, the primed integral means that the integration has to be performed in the domain for which $\mathbf{v}_{12} \cdot \hat{\sigma} > 0$. Finally, we are employing the short-hand notation $f_i^{**} \equiv f(v_i^{**}, t)$, where the two star superscript refers to the precollisional velocity:

$$ \mathbf{v}_1^{**} = \mathbf{v}_1 - (\mathbf{v}_{12} \cdot \hat{\sigma})\hat{\sigma}, \quad \mathbf{v}_2^{**} = \mathbf{v}_2 + (\mathbf{v}_{12} \cdot \hat{\sigma})\hat{\sigma} . \quad (3) $$

In this framework the evolution of our tagged particle is exactly a Markov process. The probability of hopping from a state with velocity $\mathbf{v}$ to another state of velocity in a narrow interval $d\mathbf{v}'$ around $\mathbf{v}'$ and in a time interval $dt$ is given by $W(\mathbf{v}'|\mathbf{v})d\mathbf{v}'dt$, where

$$ W(\mathbf{v}'|\mathbf{v}) = \frac{1}{\ell} |\mathbf{v}' - \mathbf{v}|^{2-d} (2\pi T_0)^{-1/2} \exp \left\{ -\frac{1}{2T_0} \times \left[ v'^2 - v^2 + \mathbf{v} \cdot \left( \frac{\mathbf{v}' - \mathbf{v}}{|\mathbf{v}' - \mathbf{v}|} \right) \right] \right\} \quad (4) $$

is the transition rate density per unit time (see e.g. [6]). For any general Markov process, the probability of leaving a given configuration $\Gamma$ in a time interval $dt$ is $r(\Gamma)dt$, where

$$ r(\Gamma) = \int d\Gamma' W(\Gamma'|\Gamma) . \quad (5) $$

This is simply the loss term of the Linear Boltzmann equation (2), and it reads for the hard sphere gas:

$$ r(v) = \frac{\omega}{\sqrt{2}} \left( \frac{v^2}{dT_0} {}_1F_1 \left( \frac{1}{2}, 1 + \frac{d}{2}, -\frac{v^2}{2T_0} \right) + e^{-\frac{v^2}{2T_0}} {}_1F_1 \left( \frac{d-1}{2}, \frac{d}{2}, \frac{v^2}{2T_0} \right) \right) , \quad (6) $$

where ${}_1F_1$ is the confluent hyper-geometric function of the first kind [7], and $\omega$ is the collision frequency of the gas:

$$ \omega = \int d\mathbf{v} \phi(\mathbf{v}) r(v) = \frac{2\pi^{\frac{d-1}{2}} \sqrt{T_0}}{\Gamma(d/2)} . \quad (7) $$

From the Markovian property, it follows that the probability of having a given velocity $\mathbf{v}$ for a given time $t$ and making a collision exactly at time $t$ is exponential:

$$ P_{FFT}(t|\mathbf{v}) = r(v) \exp(-r(v)t) . \quad (8) $$

The probability density of the free flight time follows:

$$ P_{FFT}(t) = \int d\mathbf{v} \phi_{coll}(\mathbf{v}) P_{FFT}(t|\mathbf{v}) . \quad (9) $$

Here $\phi_{coll}(\mathbf{v})$ is the velocity pdf of the colliding particles (i.e. obtained sampling the velocity of the particle only on collision):

$$ \phi_{coll}(\mathbf{v}) = \frac{r(v)}{\omega} \phi(\mathbf{v}) , \quad (10) $$

where the normalization factor $\omega$ is exactly the collision frequency of the gas. The free flight time distribution for hard spheres was already investigated in [8], but the authors used the equilibrium velocity distribution $\phi$ instead of the on collision velocity distribution $\phi_{coll}$ as a weight in (9), and hence obtained an incorrect result (this issue is further commented in Appendix A). Two measurements of the free flight time distribution in an event driven molecular dynamic simulation of hard disks (see section IV for more details) and in a Direct Simulation Monte Carlo (DSMC) [9], are shown in Fig. 1, together with the expression of Eq. (9). The agreement is excellent. The molecular dynamics simulations are performed at a low density ($\rho\sigma^2 = 10^{-2}$) to ensure that molecular chaos holds, whereas the DSMC technique relies by construction on the molecular chaos approximation, and can therefore be considered as providing a $\rho \to 0$ benchmark. Figure 1 also provides a comparison with the exponential law that would hold for a constant $r(v)$, as is the case for Maxwell molecules. The latter expectation [$\omega \exp(-\omega t)$] is seen to hold at early times, but significantly fails at large $t$, a time regime that we investigate in more details in the next section.

### B. Large time and high dimension analysis

Apart from simplified models (see Appendix B), it does not seem possible to obtain a simple exact analytical expression for the free flight time distribution, mainly because of the particular behavior of the hyper-geometric function involved in the expression of the collision rate $r(v)$. In order to get a simpler approximated result, one has to look at limiting cases. For large times, the integral involved in Eq. (9) can be estimated through the saddle-point approximation (a.k.a. method of steepest descent, see e.g. [10]), yielding:

$$ P_{FFT}(t) \sim \exp\left(-\frac{\omega t}{\sqrt{2}}\right) \frac{\omega}{2} \left(1 - \frac{2}{d} + \frac{\omega t}{\sqrt{2}d}\right)^{-d/2} . \quad (11) $$
<figure><img src="image_2.png" /><figcaption>FIG. 1: (Color online) Free flights time distribution of a hard disc gas. The circles and squares correspond to the results of Molecular Dynamics (MD) and Direct Simulation Monte-Carlo (DSMC) methods. MD simulations in all this work are shown at density $\rho = 0.01\sigma^{-2}$. The full line is the numerical integration of Eq. (9). The dot-dashed line shows the asymptotic large time behavior of Eq. (11), and the dashed line is the exponential distribution, associated to the Poisson distribution. Time is measured in units of the mean free time ($\omega = 1$).</figcaption></figure>

The above function is plotted in Fig. 1, and seen to be in very good agreement with our numerical data provided $t$ is large enough. Note that expression (11) does not provide a normalized probability (a feature visible on Fig. 1), but only an asymptotic expansion.

Besides, a way to get an explicit expression for the collision rate $r(v)$ is to investigate the infinite dimension limit. In this case, the saddle-point approximation yields:

$$r(v) \approx \frac{\omega}{\sqrt{2}} \sqrt{1 + \frac{v^2}{2T_0d}}. \qquad (12)$$

A detailed derivation of this result is given in Appendix C.

### C. Free path length distribution

From the knowledge of the free flight time distribution, it is straightforward to compute the distribution of the free path length. In fact, for a particle with a velocity $v$, the length covered in a time $t$ is $x = |v|t = vt$. Therefore, the conditional distribution of free path length is obtained as:

$$P_{FPL}(x|\mathbf{v}) = \frac{1}{v} P_{FFT} \left( \frac{x}{v} \middle| \mathbf{v} \right) \\ = \frac{r(v)}{v} \exp \left( -\frac{r(v)}{v} x \right). \quad (13)$$

The free path length distribution can be obtained averaging the above expression with the weight $\phi_{coll}$:

$$P_{FPL}(x) = \int d\mathbf{v} \phi_{coll}(\mathbf{v}) P_{FPL}(x|\mathbf{v}). \quad (14)$$

<figure><img src="image_3.png" /><figcaption>FIG. 2: (Color online) Plot of the free path length distribution from MD simulations (at density $\rho = 0.01\sigma^{-2}$) and from a numerical integration of Eq. (14). In the inset, the dashed line is the dominant exponential term, Eq. (16).</figcaption></figure>

Although the conditional probability of free path length is very similar to the conditional probability of free flight time, the step of averaging over the collisional velocity distribution may drastically change the shape of the distribution. A case worthy of attention is discussed in Appendix D, for Maxwell molecules: the free flight time distribution is a pure exponential, while the free path length distribution has a stretched exponential decay. In the case of hard particles, the dominant large path behavior remains exponential, as for the free flight time. Nevertheless, the leading exponential prefactor changes. The free flight time distribution is dominated by $e^{-r(0)t}$, where $r(0)$ is the minimum of the function $r(v)$. In the case of free path length, the dominant exponential behavior is $e^{Ax}$, where $A$ is the minimum of $r(v)/v$, which is obtained in the limit $v \to \infty$:

$$A = \lim_{v \to \infty} \frac{r(v)}{v} = \frac{1}{l\sqrt{2}} = \frac{\pi^{\frac{d-1}{2}}}{\ell \Gamma(\frac{d+1}{2})}, \quad (15)$$

where $l = \omega/\langle v \rangle$ is the mean free path. This yields:

$$P_{FPL}(x) \sim e^{\frac{x}{\sqrt{2}l}}, \quad (16)$$

(here the $\sim$ should be understood as an equivalence between the logarithms in the large $x$ limit). Unfortunately, this expression only gives the behavior of the dominating exponential term, which becomes visible at a length scales significantly larger than the mean free path. Hence, if one analyzes the result of MD simulations as in Fig. 2, one sees that there are sub-leading (algebraic) terms in the asymptotic behavior, which still play a role and that are responsible for the mismatch between Eq. (16) and MD data in Fig. 2. A similar feature holds for the free flight time distribution: upon neglecting the algebraic sub-dominant correction in Eq. (11), the agreement displayed at late times in the inset of Fig. 1 would be spoiled.
### III. THE NUMBER OF COLLISIONS

#### A. Uncorrelated approach

After having shown how the probability of the time between two subsequent collisions behaves, our interest goes to the probability of the sum of many such times, and in particular, the probability of having $N$ collisions in a time $t$. The time interval $t$ in which the $N$ collisions take place can be decomposed in a sum of $N+1$ (correlated) times:

$$t = t_{in} + t_{1} + \cdots + t_{N-1} + t_{f}. \qquad (17)$$

We are considering that a given particle begins a trajectory with a given initial velocity, and then waits a time $t_{in}$ before it makes a collision and gets a new velocity. During a time $t_1$, it then flies freely, before colliding again etc. The successive behavior of the particle's trajectory will be a sequence of collisions, spaced by other free flight times $t_2, t_3, \dots, t_{N-1}$, until the last $N$-th collision, which will be followed by a final time interval $t_f$ which is not ended with a collision. The free flight times $t_1$ to $t_{N-1}$ are of course distributed following the probability $P_{FFT}$ given by Eq. (9). We will denote by $P_{in}$ and $P_f$ the pdf of $t_{in}$ and $t_f$, although their precise statistics are irrelevant for our purposes, as becomes clear below.

If one would assume these $N+1$ times to be uncorrelated, then the probability $P(N, t)$ of having $N$ collisions in a time $t$ could be deduced from the knowledge of the free flight time probability as follows. We first introduce the Laplace transform:

$$\tilde{P}(z) = \int dt \, e^{-zt} P(t), \qquad (18)$$

where $P$ can be either $P_{FFT}$, $P_{in}$, or $P_f$. The transform of $P(N, t)$ can hence be expressed as:

$$\tilde{P}(N, z) = \tilde{P}_{in}(z) \tilde{P}_{FFT}(z)^{N-1} \tilde{P}_{f}(z). \qquad (19)$$

Of course, since $P_{FFT}$ cannot be computed analytically, a closed form expression for $\tilde{P}(N, z)$ is not available. One interesting limit to analyze is that of large times. To get an asymptotic form of the following integral

$$P(N, t) = \int dz \, e^{zt} \tilde{P}_{in}(z) \tilde{P}_{FFT}(z)^{N-1} \tilde{P}_{f}(z), \qquad (20)$$

we make the remark that the number of collisions increases with time linearly on average, and we define a time intensive counterpart to $N$, a fluctuating collision rate, which we denote by $n \equiv N/t$. The integral of Eq. (20) then reads:

$$P(N, t) = \int dz \exp \left( t(z + n \ln(\tilde{P}_{FFT}(z))) + \mathcal{O}(1) \right). \qquad (21)$$

Finally, when $t \to \infty$, the above probability behaves as:

$$P(N, t) \sim e^{\pi(n)t}, \qquad (22)$$

where $\pi(n)$ is a large deviation function, and is related to $\tilde{P}_{FFT}$ through the following relation:

$$\pi(n) = \min_{z} (z + n \ln(\tilde{P}_{FFT}(z))). \qquad (23)$$

#### B. Large time behavior

The result of the previous subsection leads to an approximate evaluation of the large deviation function of the number of collisions that neglects temporal correlations. We will now investigate this large deviation function keeping into account these correlations. Let us define the joint probability $f(v, N, t)$ of having a velocity $v$ and having suffered $N$ collisions up until time $t$. In a homogeneous state, the time evolution of the above defined probability is governed by a slightly different form of the linear Boltzmann equation:

$$\partial_t f(\mathbf{v}_1, \mathcal{N}, t) = \int d\mathbf{v}_2 \int' d\hat{\sigma}(\mathbf{v}_{12} \cdot \hat{\sigma}) \times [f(\mathbf{v}_1^{**}, \mathcal{N}-1, t)\phi(\mathbf{v}_2^{**}) - f(\mathbf{v}_1, \mathcal{N}, t)\phi(\mathbf{v}_2)] . \qquad (24)$$

The function $\phi(v)$ still represents the one point velocity pdf of the gas, which is in equilibrium. Moreover, since the particle whose collisions we are counting has the same velocity pdf, we enforce that:

$$\sum_{\mathcal{N}=0}^{\infty} f(\mathbf{v}, \mathcal{N}) = \phi(\mathbf{v}). \qquad (25)$$

It is useful to introduce the generating function of $f$ as:

$$\hat{f}(\mathbf{v}, \lambda, t) = \sum_{\mathcal{N}=0}^{\infty} e^{-\lambda \mathcal{N}} f(\mathbf{v}, \mathcal{N}, t). \qquad (26)$$

It can be seen that $\hat{f}$ evolves according to:

$$\partial_t \hat{f}(\mathbf{v}_1, \lambda, t) = \int d\mathbf{v}_2 \int d\hat{\sigma} \Theta(\mathbf{v}_{12} \cdot \hat{\sigma})(\mathbf{v}_{12} \cdot \hat{\sigma}) \times [e^{-\lambda} \hat{f}(\mathbf{v}_1^{**}, \lambda, t) \phi(\mathbf{v}_2^{**}) - \hat{f}(\mathbf{v}_1, \lambda, t) \phi(\mathbf{v}_2)] . \qquad (27)$$

The large time behavior of the solution of the above equation is dominated by the largest eigenvalue $\mu(\lambda)$ of the evolution operator:

$$\mu(\lambda) \tilde{f}(\mathbf{v}_1, \lambda) = L_\lambda \tilde{f}(\mathbf{v}_1, \lambda) = \int d\mathbf{v}_2 \int d\hat{\sigma} \Theta(\mathbf{v}_{12} \cdot \hat{\sigma})(\mathbf{v}_{12} \cdot \hat{\sigma}) \times [e^{-\lambda} \tilde{f}(\mathbf{v}_1^{**}, \lambda) \phi(\mathbf{v}_2^{**}) - \tilde{f}(\mathbf{v}_1, \lambda) \phi(\mathbf{v}_2)], \qquad (28)$$

where $\tilde{f}$ denotes the eigenfunction of $L_\lambda$ associated with $\mu$:

$$\hat{f}(\mathbf{v}, \lambda, t) \propto e^{\mu(\lambda)t} \tilde{f}(\mathbf{v}, \lambda). \qquad (29)$$

Moreover, since $\hat{P}(\lambda) \sim e^{\mu(\lambda)t}$, one sees that $\mu(\lambda)$ is proportional to the cumulant generating function:

$$\langle \mathcal{N}^p \rangle_c = t(-1)^p \frac{\partial^p \mu}{\partial \lambda^p} \Big|_{\lambda=0} . \qquad (30)$$

Furthermore, $P(N) \sim \int d\lambda e^{\lambda N} \hat{P}(\lambda)$, and hence, applying the saddle point method for $t \to \infty$, one has that

$$P(N, t) \sim \exp(t\pi(n)), \qquad (31)$$

where $n = N/t$ and the large deviation function $\pi$ is related to $\mu$ through a Legendre transform:

$$\pi(n) = \min_{\lambda} (\mu(\lambda) + \lambda n). \qquad (32)$$
1. Large $\lambda$ behavior

Before trying to solve Eq. (28) extending the methods of kinetic theory, we shall first provide exact results, which can be extracted from the analysis of asymptotically large values of $\lambda$. The evolution operator appearing in Eq. (28) can be cast in the form:

$$L_{\lambda} = L_{0} + e^{-\lambda} L_{1}, \quad (33)$$

where

$$L_0 \tilde{f}(\mathbf{v}_1, \lambda) = - \int d\mathbf{v}_2 \int' d\hat{\sigma}(\mathbf{v}_{12} \cdot \hat{\sigma}) \tilde{f}(\mathbf{v}_1, \lambda) \phi(\mathbf{v}_2) = -r(v_1) \tilde{f}(\mathbf{v}_1, \lambda), \quad (34a)$$

$$L_1 \tilde{f}(\mathbf{v}_1, \lambda) = \int d\mathbf{v}_2 \int' d\hat{\sigma}(\mathbf{v}_{12} \cdot \hat{\sigma}) \tilde{f}(\mathbf{v}_1^{**}, \lambda) \phi(\mathbf{v}_2^{**}). \quad (34b)$$

For large values of $\lambda$, the coefficient $e^{-\lambda}$ plays the role of a small parameter, and therefore the eigenvalue equation (28) can be solved in perturbation theory. We will therefore try to get an expression of $\mu(\lambda)$ for large $\lambda$ of the form:

$$\mu(\lambda) = \mu^{(0)} + e^{-\lambda}\mu^{(1)} + \mathcal{O}(e^{-2\lambda}). \quad (35)$$

To zero-th order, the largest eigenvalue of $L_\lambda$ is given by the maximum of the function $-r(v)$, whose expression is written in Eq. (6). The maximum of this function occurs at $v = 0$, and

$$\mu^{(0)} = -r(0) = -\frac{\omega}{\sqrt{2}}. \quad (36)$$

The eigenfunction associated with this eigenvalue is indeed a delta function, centered in $\mathbf{v} = \mathbf{0}$. In order to get the first order correction to the eigenvalue (36) one has to project the zero-th order eigenfunction on the operator proportional to the small parameter:

$$\mu^{(1)} = \int d\mathbf{v}_1 L_1 \delta(\mathbf{v}_1) = \frac{\omega}{\sqrt{2}}. \quad (37)$$

Finally one finds that for large $\lambda$, the largest eigenvalue $\mu$ behaves like:

$$\mu(\lambda) = \frac{\omega}{\sqrt{2}}(e^{-\lambda} - 1) + \mathcal{O}(e^{-2\lambda}). \quad (38)$$

Hence, the probability of having $N$ collisions behaves, for values of $N$ small with respect to its average $\langle N \rangle = \omega t$, as a Poisson distribution with a frequency equal to $\omega/\sqrt{2}$:

$$P(N) \sim \frac{e^{-\frac{\omega t}{\sqrt{2}}}}{N!} \left(\frac{\omega t}{\sqrt{2}}\right)^N, \quad \text{for } N \ll \omega t. \quad (39)$$

For such a distribution, the large deviation function $\pi$ easily follows:

$$\pi(n) = n - n \log(n\sqrt{2}/\omega) - \omega/\sqrt{2} \quad (40)$$

The result obtained in Eq. (39) is compatible with the large time behavior embodied in Eq. (11): the time derivative of the probability of having $N=0$ collisions is of course (minus) the free flight time distribution function. We emphasize that both results (39) and (11) hold for large times. In addition, we note that the only space dimension dependence involved is through the collision frequency $\omega$ [see Eq. (7)].

2. An approximate solution

In order to get an approximate expression for the largest eigenvalue $\mu(\lambda)$, we will suppose the associated eigenvector $\tilde{f}(\mathbf{v}, \lambda)$ to be a Gaussian with a given temperature $T(\lambda)$. We expect that this approximation will provide accurate results for small values of $\lambda$, given that for $\lambda = 0$, $\tilde{f}(\mathbf{v}, 0)$ is exactly a Gaussian with a temperature $T_0$, but it is not a priori a systematic approximation for larger values of $\lambda$. Projecting the Boltzmann-like equation (28) onto the first two velocity moments of the eigenfunction we want to compute, we are left with two closed equations for $\mu(\lambda)$ and for $T(\lambda)$:

$$\mu(\lambda) = \nu_0, \quad \mu(\lambda)dT(\lambda) = \nu_2, \quad (41)$$

where the $\nu_n$ are the collisional moments (the expression of the first ones is given in Appendix E):

$$\nu_n = \int d\mathbf{v}_1 \int d\mathbf{v}_2 \int' d\hat{\sigma} v_1^n (\mathbf{v}_{12} \cdot \hat{\sigma}) \times [\tilde{f}(\mathbf{v}_1^{**}, \lambda)\phi(\mathbf{v}_2^{**}) - \tilde{f}(\mathbf{v}_1, \lambda)\phi(\mathbf{v}_2)]. \quad (42)$$

Solving simultaneously Eqs. (41), one obtains:

$$\mu(\lambda) = -\frac{\omega}{\sqrt{2}}(1 - e^{-\lambda})\sqrt{1 + \frac{T(\lambda)}{T_0}}, \quad (43a)$$

$$T(\lambda) = \frac{\sqrt{2}T_0}{\sqrt{1 + e^{\lambda}}}. \quad (43b)$$

One may notice that this result satisfies some of the previous requirements obtained from the asymptotic large $\lambda$ analysis. In particular $\mu(\infty) = -\omega/\sqrt{2}$, as in the previous subsection. Moreover, the fictitious temperature $T(\lambda)$ vanishes for infinite $\lambda$, meaning that the eigenfunction associated with $\mu$ indeed tend towards a delta function. Nevertheless, while these features concerning the zero-th order perturbation results of the previous section are fulfilled, the behavior of Eq. (43) is different at the next order. In fact, one can see from Eq. (43) that, for $\lambda \to \infty$:

$$\mu(\lambda) \sim -\frac{\omega}{\sqrt{2}} - \omega \frac{\sqrt{e^{-\lambda}}}{2} + \mathcal{O}(e^{-\lambda}), \quad (44)$$

to be compared with Eq. (38). This is a deficiency of the Gaussian approximation.

In order to get an expression of $\pi(n)$ one should compute the Legendre Transform of $\mu(\lambda)$. This seems not feasible analytically, but can be achieved numerically. The function $\pi(n)$
<figure><img src="image_4.png" /><figcaption>FIG. 3: (Color online) Plot of $\pi(n)$ obtained in several fashions ($\omega = 1$). The solid line has been obtained from the numerical Legendre Transform of $\mu(\lambda)$ from Eq. (43). The dashed line corresponds to the uncorrelated result obtained in Eq. (23), while the dot-dashed line is the large deviation function of the Poisson distribution of average 1. For completeness, we also show the large deviation function of the Poisson distribution with average $1/\sqrt{2}$ (dotted line).</figcaption></figure>

is plotted in Fig. 3, together with the result of uncorrelated estimation of Eq. (23), as well as two Poisson distributions, of average $\omega$ and $\omega/\sqrt{2}$. One can see that the two estimations carried out are in very close agreement, apart from a slight difference in the behavior of the right tail. Moreover, the two new results presented here are clearly different from the Poisson distribution with mean $\omega$ [denoted Poisson(1) since $\omega = 1$ in Fig. (3)], which is narrower, and underestimates extreme events, characterized by either very few or many collisions with respect to typical realizations.

### 3. Sonine perturbation and cumulants

We now further exploit the property that for $\lambda = 0$ the solution of the Boltzmann equation (28) is exactly a Gaussian. This incites us, for $\lambda$ close to 0, to search for a solution $\tilde{f}(\mathbf{v}, \lambda)$ as a small perturbation of a Gaussian distribution. One of the most useful expansions in kinetic theory is the Sonine polynomials expansion $^1$. In practice this expansion consists in looking for solutions expressed as a Gaussian times a series of Sonine polynomials, denoted $S_n(x)$:

$$ \tilde{f}(\mathbf{v}, \lambda) = (2\pi)^{-d/2} \exp\left(-\frac{v^2}{2T}\right) \sum_{n=0}^{\infty} a_n S_n\left(-\frac{v^2}{2T}\right). \quad (45) $$

The first Sonine polynomials are:

$$ S_0(x) = 1, \quad (46a) $$

$$ S_1(x) = -x + \frac{d}{2}, \quad (46b) $$

$$ S_2(x) = \frac{1}{2}x^2 - \frac{d+2}{2}x + \frac{d(d+2)}{8}. \quad (46c) $$

These polynomials have the property of being orthogonal with respect to a Gaussian measure in dimension $d$, and are therefore related to Laguerre polynomials. From this feature, it follows that the coefficient of the series (45) $a_0 = 1$ and that $a_1 = 0$. Hence, the first nontrivial correction to the Gaussian approximation comes from the term proportional to $a_2$ in the expansion (45). The procedure to get an estimate of the coefficient $a_2$ consists of solving a closed system of equations obtained projecting the equation (28) onto the first velocity moments:

$$ \mu(\lambda)m_n = \nu_n, \quad (47) $$

where:

$$ m_n(\lambda) = \int d\mathbf{v} \tilde{f}(\mathbf{v}, \lambda), \quad (48) $$

and $\nu_n(\lambda)$ denotes the collisional moment of order $n$ defined in Eq. (42). Truncating the expansions (45) up to the second Sonine polynomial, one gets for the moments $m_n$:

$$ m_0 = 1, \quad m_2 = dT(\lambda), \quad (49a) $$

$$ m_4 = d(d+2)T^2(\lambda) (1 + a_2(\lambda)). \quad (49b) $$

The expression of the first collisional moments in the Sonine approximation is given in the Appendix E. Hence, taking the moment equation (48) for $n = 0, 2$ and $4$ gives a closed system of equations for $\mu(\lambda)$, $T(\lambda)$ and $a_2(\lambda)$. This system can be solved perturbatively expanding its solutions in power series around $\lambda = 0$. The cumulants of $\mathcal{N}$ obtained by means of the latter expansion turn out to be remarkably accurate. The expansion of $\mu(\lambda)$ up to the 3rd order in $\lambda$ gives access to the first three cumulants, which read:

$$ \frac{\langle \mathcal{N} \rangle_c}{\omega t} = 1, \quad (50a) $$

$$ \frac{\langle \mathcal{N}^2 \rangle_c}{\omega t} = \frac{9}{64} \left( 8 + \frac{1}{4d+3} \right), \quad (50b) $$

$$ \frac{\langle \mathcal{N}^3 \rangle_c}{\omega t} = \frac{28d(64d(320d+729) + 35775) + 257391}{8192(4d+3)^3}, \quad (50c) $$

As can be noted from these values, the variance of the number of collisions already deviates by more than 12.5% from its Poisson value (equal to unity).

## IV. NUMERICAL RESULTS

In this section, we compare the theoretical results obtained in the previous sections against numerical simulations. In

<footnote n="1">¹ For a general treatment of Sonine polynomials in kinetic theory, see e.g. [11] and references therein.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 4: (Color online) The circles are the numerical measurements of $\mu(\lambda)$ for the Maxwell model. The solid line is the generating function of the Poisson distribution. The inset shows the same data, shifted vertically by one (in such a way that $\mu(\lambda) + 1$ is always positive), in semi-logarithmic scale.</figcaption></figure>

addition to the two aforementioned techniques of Molecular Dynamics and Direct Simulations Monte Carlo (see Fig. 1), we have shaped a third numerical tool, constructing a Monte Carlo algorithm in order to directly solve the eigenvalue equation (28). We start by setting the stage for the latter method, before briefly commenting on the Molecular Dynamics simulations used to measure the statistics of the number of collisions suffered by tagged particles.

### A. Monte Carlo approach

In order to derive an algorithm for solving Eq. (28), it is useful to rewrite this equation in the form:

$$ \mu(\lambda)\tilde{f}(\mathbf{v}, \lambda) = e^{-\lambda} I[\tilde{f}|\phi] + (1 - e^{-\lambda}) A[\tilde{f}|\phi], \quad (51) $$

where

$$ I[f|g] = \int d\mathbf{v} \int' d\hat{\sigma}(\mathbf{v}_{12} \cdot \hat{\sigma}) [f_1^{**} g_2^{**} - f_1 g_2], \quad (52) $$

is the collision integral describing the elastic collision between two particles having respectively a velocity pdf $f$ and $g$. We are using the short-hand notation $f_i^{**} \equiv f(v_i^{**})$. The functional $A[f|g]$ is the loss term of the above collision integral. It actually describes the statistics of hard spheres which annihilate after each collision:

$$ A[f|g] = - \int d\mathbf{v} \int' d\hat{\sigma}(\mathbf{v}_{12} \cdot \hat{\sigma}) f_1 g_2 . \quad (53) $$

In the context of the nonlinear Boltzmann Equation, i.e. when $\phi = \tilde{f}$ in Eq. (51), and for positive values of $\lambda$, the rhs of Eq. (51) exactly describes the time evolution of the velocity pdf of the probabilistic ballistic annihilation process. This model, introduced in [3], consists in a system of $N$ particles

<figure><img src="image_6.png" /><figcaption>FIG. 5: (Color online) Numerical measurement of $\mu(\lambda)$ for the hard-sphere model. The circles show the results of Monte Carlo simulations. The solid line is the theoretical prediction within the Gaussian approximation. The dashed line is the theoretical prediction (obtained solving numerically (47-48) in the framework of the first Sonine correction). Finally the dotted line is the generating function of a Poisson distribution of parameter $\omega/\sqrt{2}$, which should asymptotically dominate for large $\lambda$ (see Fig. 6). The inset shows a zoom near $\lambda = 0$, where the difference between Gaussian and Sonine orders can only be appreciated for the larger values of $\lambda$ displayed.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6: (Color online) Same data as in Fig. 5, divided by the generating function of a Poisson distribution of parameter $\omega/\sqrt{2}$.</figcaption></figure>

which move ballistically, and interact when at contact. The interaction may be an elastic collision, with probability $e^{-\lambda}$, or an annihilation (with probability $1 - e^{-\lambda}$). In the context of the linear Boltzmann equation, this process can be extended in the following way. Consider a set of $N$ independent systems, one of each is just a single particle, characterized by its velocity $v_i(t)$ ($i = 1, \dots, N$), assuming spatial homogeneity. If each of these systems (particles) evolves in a hard sphere gas (in the thermodynamic limit) thermalized at temperature $T_0$, with a one point velocity pdf $\phi$, and can both collide elastically (with probability $e^{-\lambda}$), or annihilate (with probability $1 - e^{-\lambda}$), then, the reduced one point velocity pdf $g(\mathbf{v}, t) = \langle \sum_i^N \delta(\mathbf{v} - \mathbf{v}_i(t)) \rangle$ will verify the Boltzmann Equation.
<figure><img src="image_8.png" /><figcaption>FIG. 7: (Color online) Temperature (defined as the variance) of the eigenfunction $\tilde{f}(\mathbf{v}, \lambda)$ measured in Monte Carlo simulations. The solid and dashed lines are the predictions of the Gaussian and Sonine approximations respectively.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8: (Color online) Coefficient $a_2$ ($\lambda$-dependent velocity kurtosis), measured in a Monte Carlo simulation. The dashed line corresponds to the numerical solution of the system (49), which is correct only for small values of $a_2$.</figcaption></figure>

tion

$$ \frac{\partial g(\mathbf{v}, t)}{\partial t} = L_{\lambda} g(\mathbf{v}, t). \qquad (54) $$

Since $N(t) = \int d\mathbf{v}g(\mathbf{v}, t)$ is the number of particles at time $t$, and since we know that for long times $g(\mathbf{v}, t) \sim e^{\mu t}$, it is clear that for long times one has that $N(t) \sim N_0e^{\mu t}$. Note that the present interpretation implicitly assumes positive values of $\lambda$. Moreover, the particles can only annihilate or collide; then the total number of particles can only decrease, and hence $\mu$ must be negative. The above observations provide a numerical tool for measuring $\mu(\lambda)$ as the decay rate of the total number of particles. The main difficulty, which we have successfully addressed, with the above algorithm, is that the number of particles constantly decreases, and there is no steady state but the trivial state $N = 0$. To circumvent this practical difficulty (that would lead to somewhat noisy statistics in the simulation), we have introduced an external source of particles (systems), acting in such a way that the total number of particles is conserved. If every time that an annihilation takes place, a new particle is inserted as the clone of one of the $N-1$ remaining particles (chosen uniformly among this population), then the evolution equation of $g(\mathbf{v}, t)$ reads:

$$ \partial_t g(\mathbf{v}, t) = L_{\lambda} g(\mathbf{v}, t) + s g(\mathbf{v}, t), \qquad (55) $$

where $s$ is a constant rate. At late times, if the above equation has a steady state, the largest eigenvalue of the operator $L_{\lambda}+s$ vanishes, and hence one has that $\mu(\lambda)+s=0$. Finally one can measure $\mu(\lambda)$ simply as (minus) the steady state average of the number of particles injected by the external source.

So far, we showed how to construct a Markov chain in order to simulate the eigenvalue equation (28) for positive values of $\lambda$. For negative values of $\lambda$ the procedure is almost identical. Introducing a new time scale $\tilde{t} = (2e^{-\lambda} - 1)t$, Eq. (54) can be rewritten as:

$$ \partial_t g(\mathbf{v}, \tilde{t}) = \frac{e^{-\lambda}}{2e^{-\lambda} - 1} I[g|\phi] - \frac{e^{-\lambda} - 1}{2e^{-\lambda} - 1} A[g|\phi]. \qquad (56) $$

All the previous physical interpretations and remarks still hold, except that now instead of having annihilation (with probability $(e^{-\lambda} - 1)/(2e^{-\lambda} - 1)$), one has duplication (or cloning). Hence $\mu(\lambda)$ will be positive, and one can add an external source in order to remove particles when new particles are created. Summarizing, the algorithm proceeds as follows:

(o) The velocity of the $N$ particles are stored in a $N \times d$ matrix. A scalar $s$ is set to 0. A (small) time step $dt$ is chosen.

(i) A particle is chosen randomly in the population with uniform probability. Its velocity is denoted $\mathbf{v}_1$.

(ii) An interaction is accepted with a probability $\propto \mathbf{v}_{12} \cdot \hat{\sigma}\Theta(\mathbf{v}_{12} \cdot \hat{\sigma}) \times dt$, where $\hat{\sigma}$ is a random direction in $d$ dimensions, and $\mathbf{v}_2$ a $d$-dimensional zero-mean Gaussian random variable of variance $T_0$. In practice, $dt$ has to be chosen in such a way that $|\mathbf{v}_{12} \cdot \hat{\sigma}|dt$ is always smaller than (or equal to) one.

(iii) When the interaction is accepted, if $\lambda > 0$ (resp. $\lambda < 0$), the particle will have a post-collisional velocity $\mathbf{v}_1^*$ with probability $e^{-\lambda} \left( \text{resp. } \frac{e^{-\lambda}}{2e^{-\lambda} - 1} \right)$ or will be removed (resp. duplicated) otherwise.

(iv) If the particle has been removed (resp. duplicated) in step (iii), one of the $N-1$ remaining particles is chosen randomly and uniformly, and is duplicated (resp. removed). $s$ is increased by 1.

(v) Time is increased by an amount $dt$ (resp. $\tilde{t}$).

(vi) Back to (i).

This algorithm bears some similarities with the approach proposed in [12], also intended to directly measure large deviation functions. Nonetheless, the version proposed here is
more inspired by some variants of the DSMC algorithm for systems which do not conserve the total number of particles [13, 14, 15]. In order to check the reliability of the algorithm, we have first performed simulations in the case of Maxwell molecules, where the number of collisions is exactly distributed following the Poisson distribution (cf. the Introduction). The measurements of $\mu(\lambda)$ for this particular model are shown in Fig. 4, together with the generating function of the Poisson distribution. The agreement is excellent. In the simulations the temperature scale is set by the temperature of the heat bath $T_0$, which we set to unity. The time scale is set by the mean free time, which we also set to unity. In the simulation data, the mean collision frequency is therefore $\omega = 1$.

In the case of the hard-sphere model, a measurement of the largest eigenvalue $\mu(\lambda)$ is shown in Fig. 5. One can see that the numerical results are in very good agreement both with the Gaussian and the Sonine approximations. Figure 6 shows the large $\lambda$ behavior of $\mu(\lambda)$ divided by the prediction (38). In this case, a good agreement between the numerical data and the theoretical predictions is found. In the framework of the above described Monte Carlo algorithm, it is also possible to measure the stationary velocity pdf $g(v, t = \infty, \lambda)$ which is equal to the eigenfunction $\tilde{f}(v, \lambda)$ associated with $\mu(\lambda)$. Hence one can also compare the analytical predictions for the temperature $T(\lambda)$ and for the Sonine coefficient $a_2(\lambda)$ with the Monte Carlo results, see Figs. 7 and 8. In Fig. 7, one can see the temperature $T(\lambda)$. Here the Gaussian approximation is already able to capture the general behavior of this “effective” temperature, and for small values of $\lambda$, the Sonine corrections compare very well with the simulation data. In Fig. 8, one can see the coefficient $a_2$ as a function of $\lambda$. As expected, for small values of $\lambda$ the expansion carried out in the previous section correctly describes the result of the simulations.

## B. Molecular Dynamics

The most direct way to measure the large deviation function of the number of collisions $N$ is of course to count it in a (numerical) experiment, and then construct its probability distribution function. To this end, we have performed Molecular Dynamics (MD) simulations of a two dimensional hard-disk gas of $N = 10^3$ particles of diameter $\sigma = 1$ at density $\rho = N/V = 10^{-3}\sigma^{-2}$ and $10^{-2}\sigma^{-2}$. The particles evolve in a square box with periodic boundary conditions, and the time is measured in mean free time units, in such a way that $\langle N \rangle = t$. We measured the statistics of the number of collisions $N$ suffered by each particle. The value of the first cumulants is reported in Table I, together with the Poisson prediction, as well as the results from the Gaussian and Sonine approximations. It seems that, when time increases, the cumulants converge towards a value which is close to the Sonine predictions. However, it must be noted that even if increasing time makes finite time corrections smaller, the statistics become poorer and poorer. The agreement with the results from the Sonine approximation is very good, while the Gaussian order already provides a reliable estimation. This is further shown in Figure 9, where two measurements of the large deviation function of the number of collisions for two different times, $t = 10$ and $t = 50$ (with $\sigma = 1$ and $T_0 = 1$), are compared to the numerical inverse Legendre Transform of $\mu(\lambda)$ at Gaussian order (given by Eq. (43)). Here again, the Poisson prediction is distinctly off.

<figure><img src="image_10.png" /><figcaption>FIG. 9: (Color online) Large deviation function $\pi(n)$ of the number of collisions $N$ suffered by each single particle. Symbols show the results of Molecular Dynamics simulations. The solid line is the result in the framework of the Gaussian approximation, and the dashed line is the large deviation function of the Poisson distribution, of average $1: \pi(n) = n - n \log n - 1$.</figcaption></figure>

<table><caption>TABLE I: Cumulants for the number of collisions $N$ from MD simulations, and comparison with the Gaussian and Sonine approximations (time is measured in units of the mean-free time).</caption><thead><tr><th></th><th>⟨N⟩<sub>c</sub>/t</th><th>⟨N<sup>2</sup>⟩<sub>c</sub>/t</th><th>⟨N<sup>3</sup>⟩<sub>c</sub>/t</th></tr></thead><tbody><tr><td>t = 10</td><td>1.</td><td>1.1228</td><td>1.1282</td></tr><tr><td>t = 50</td><td>1.</td><td>1.1354</td><td>1.1045</td></tr><tr><td>Poisson</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Gaussian</td><td>1</td><td>1.125</td><td>1.1289</td></tr><tr><td>Sonine</td><td>1</td><td>1.1377</td><td>1.1073</td></tr></tbody></table>

## V. CONCLUSION

We have shown that the collisional statistics of the hard sphere gas exhibits clear deviations from the Poisson distribution. These deviations have been consistently quantified both analytically and numerically. In the analytical treatment, the cumulant generating function $\mu(\lambda)$ plays a pivotal role. The eigenvalue equation defining this object can be fruitfully interpreted in terms of population dynamics with annihilation and cloning events, which eventually leads to an efficient algorithm allowing to compute the various quantities involved in the theoretical analysis. The corresponding numerical method, of Monte Carlo type, should not be confused with the more conventional Direct Simulation Monte Carlo technique, which we also implemented, and that is intended to solve a different kinetic equation (the Boltzmann equation).
Finally, a third numerical method (Molecular Dynamics) was employed. These three routes provide complementary and valuable results, that strongly support our predictions.

Interestingly, the present formalism can be extended to the study of out of equilibrium systems, such as granular gases, where the question of the collisional statistics has been the focus of recent interest [17, 18, 19, 20, 21]. In particular, the strong effect of dissipation on the distribution of free flight times reported in [19, 20] calls for further investigations. In addition, if one considers a gas of inelastic smooth hard spheres, kept in a steady state by a velocity independent force (as e.g. a vibrating wall of the container, or a stochastic force acting independently on each particle), then the phase space volume of the system has been reduced, after a time $t$, by a factor $(1 - \alpha)^{N(t)}$, where $\alpha$ is the coefficient of normal restitution [16]. Hence one sees that in this non-equilibrium system, the number of collisions can be exactly identified, up to a constant prefactor, with the integrated phase space contraction rate. This quantity has already been the subject of many works in non-equilibrium statistical mechanics, and is often intimately related with the irreversible entropy production (see e.g. [22] and references therein).

## Acknowledgments

The authors acknowledge useful discussions with J. Pasecki, J. M. J. van Leeuwen, M. H. Ernst, D. Frenkel, H. van Beijeren and P. Krapivsky. This work was supported by the French Ministry of Education through a ANR-05-JCJC-44482 grant.

## APPENDIX A: COMMENT ON THE FREE FLIGHT TIME DISTRIBUTION

We shall give here some additional arguments on the in-correctness of Eq. (9), when used with a Gaussian weight $\phi$, instead of the weight $\phi_{coll}$ defined in Eq.(10). We shall refer to this (erroneous) distribution as $P_{FFT}^W$:

$$P_{FFT}^{W}(t) = \int d\mathbf{v} \phi(\mathbf{v}) P(t|\mathbf{v}) . \quad (A1)$$

A first argument bears on the inconsistency between the above relation and the definition of the collision frequency $\omega = \langle r(v) \rangle$, where the brackets denote an average over a Gaussian weight. Indeed, the average time $\tau$ between two subsequent collisions (mean free time) of a given particle is equal to the inverse of the collision frequency: $\tau = 1/\langle r(v) \rangle$. Besides, from Eq. (A1) the mean free time $\tau$ is obtained as:

$$\tau = \int_{0}^{\infty} dt P_{FFT}^{W}(t) = \left\langle \frac{1}{r(v)} \right\rangle \neq \frac{1}{\langle r(v) \rangle} . \quad (A2)$$

Hence we see that expression (A1) is in contradiction with the definition of the collision frequency. On the other hand, the counterpart of Eq. (A2) with the distribution provided by (9)

and (10), yields

$$\tau = \int_{0}^{\infty} dt \int d\mathbf{v} \frac{r^2(v)}{\omega} e^{-r(v)t} \phi(\mathbf{v}) = \frac{1}{\langle r \rangle}, \quad (A3)$$

which is the required relation.

Second, expression (A1) explicitly slightly differs from the correct free flight time distribution in two limiting cases.

* At small values of $t$, $P_{FFT}^W$ behaves as:

$$P_{FFT}^{W}(t) \sim \langle r(v) \rangle - \langle r(v)^2 \rangle t + \mathcal{O}(t^2), \quad (A4)$$

while the true distribution behaves as

$$P_{FFT}(t) \sim \frac{\langle r(v)^2 \rangle}{\omega} - \frac{\langle r(v)^3 \rangle}{\omega^2} t + \mathcal{O}(t^2), \quad (A5)$$

so that even the value in $t = 0$ is different, although this difference is numerically small (see also Fig. 10). For instance, in $d = 2$ one has:

$$\frac{\langle r(v)^2 \rangle}{\omega^2} \approx 1.06354 \quad (A6)$$

and

$$\frac{\langle r(v)^3 \rangle}{\omega \langle r(v)^2 \rangle} \approx 1.13962 . \quad (A7)$$

* The large time behavior of the probability also slightly differs from Eq. (11):

$$P_{FFT}^{W} \sim \exp \left( -\frac{\omega t}{\sqrt{2}} \right) 2^{\frac{d-1}{2}} \omega \left( 2 - \frac{2}{d} + \frac{\sqrt{2}\omega t}{d} \right)^{-d/2} . \quad (A8)$$

In the above expression, the leading exponential term remains the same as in Eq. (11), since it is only determined by the minimum of the function $r(v)$. Nonetheless, the subleading prefactor is slightly different, as it can be appreciated in Fig. 10.

## APPENDIX B: FREE FLIGHTS TIME DISTRIBUTION FOR VERY HARD PARTICLES

For the sake of completeness, we report in this appendix the result for the free flight time distribution arising in the very hard particle (VHP) model. This framework allows for an explicit analytic computation of $P_{FFT}(t)$. The model, introduced by Ernst and Hendriks [23], consists in a choice of the collision rate proportional to the total energy of the system. In the case of a single tagged particle, the corresponding linear Boltzmann equation reads:

$$\frac{\partial}{\partial t} f(\mathbf{v}_1, t) = \frac{1}{\ell} \int d\mathbf{v}_2 \int d\hat{\sigma} \frac{(\mathbf{v}_{12} \cdot \hat{\sigma})^2}{\sqrt{T_0}} [f_1^{**} \phi_2^{**} - f_1 \phi_2] . \quad (B1)$$

The velocity dependent collision rate is thus quadratic in $v$:

$$\frac{r_{VHP}(\mathbf{v})}{\omega_{VHP}} = \frac{1}{2} + \frac{v^2}{2dT_0}, \quad (B2)$$
<figure><img src="image_11.png" /><figcaption>FIG. 10: (Color online) Free flights time distribution of a hard disc gas. The circles corresponds to the results of Molecular Dynamics (MD) simulations at density $\rho = 0.01\sigma^{-2}$. The full line is the numerical integration of Eq. (9). The dashed line is the numerical integration of $P_{FFT}^W$, and the dotted line its large time behavior described by Eq. (A8). The inset shows the same data in semi-logarithmic scale, where it is possible to note that at large times the prediction of $P_{FFT}^W$ is always above the true distribution. Time is measured in units of the mean free time ($\omega = 1$).</figcaption></figure>

where

$$ \omega_{VHP} = \frac{2\pi^{d/2}\sqrt{T_0}}{\ell\Gamma(d/2)} \qquad (B3) $$

is the collision frequency for VHP. Due to the simpler form of the collision rate, the free flight time distribution can be expressed analytically:

$$ P_{FFT}^{VHP}(t) = \frac{1}{4} d^{d/2} e^{-t/2} \times (d+t)^{-\frac{d}{2}-2} (4d^2 + (4t+2)d + t^2), \quad (B4) $$

where the time is here in units of the mean collision rate $\omega_{VHP}$. Note that here the large time behavior is $\sim e^{-\omega t/2}$, which is even slower than the Hard-Sphere gas behaviour (which is $\sim e^{-\omega t/\sqrt{2}}$) with respect to the Poissonian case ($\sim e^{-\omega t}$). A comparison between these three distributions is shown in Fig. 11. The observed large time behaviors suggest that Maxwell and VHP models provide a upper and lower bounds for the hard-sphere model, a phenomenon reminiscent to that observed in the long time behaviour of dynamics that do not conserve the density [15].

## APPENDIX C: SADDLE-POINT APPROXIMATION FOR THE COLLISION RATE

In this appendix we show how to recover Eq. (12) with the saddle point method. We first note that $r(v)$ can be expressed as the difference of two different integrals:

$$ r(v) = \Omega_{d-1}(I_1(v) - I_2(v)), \qquad (C1) $$

<figure><img src="image_12.png" /><figcaption>FIG. 11: (Color online) Free flights time distributions for three different interaction kernels, in two dimensions. The (red) dashed line shows the result for VHP (Eq. B4), the (blue) solid line and the (green) dotted line show respectively the results for Hard Spheres and Maxwell Molecules. Time is measured in units of the mean free time ($\omega = 1$).</figcaption></figure>

where

$$ I_1(v) = \int_{-\infty}^{+\infty} dv_{2\sigma} \int_0^{\pi} d\theta \sin \theta^{d-2} (v \cos \theta) \times \Theta(v \cos \theta - v_{2\sigma}) \phi(v_{2\sigma}), \quad (C2) $$

$$ I_2(v) = \int_{-\infty}^{+\infty} dv_{2\sigma} \int_0^{\pi} d\theta \sin \theta^{d-2} v_{2\sigma} \times \Theta(v \cos \theta - v_{2\sigma}) \phi(v_{2\sigma}). \quad (C3) $$

For the first integral it is more convenient to perform first the integration over the angle $\theta$, which leads to:

$$ I_1(v) = \frac{v}{d-1} \int_{-\infty}^{+\infty} dv_{2\sigma} \left(1 - \left(\frac{v_{2\sigma}}{v}\right)^2\right)^{\frac{d-1}{2}} \phi(v_{2\sigma}). \qquad (C4) $$

Then, defining the rescaled variables $\tilde{v} = v/\sqrt{2T_0d}$ and $\tilde{v}_{2\sigma} = v_{2\sigma}/\sqrt{2T_0d}$ one finds, to leading order in $d$:

$$ I_1(\tilde{v}) \approx \tilde{v} \sqrt{\frac{T_0}{\pi}} \int_{-\infty}^{+\infty} d\tilde{v}_{2\sigma} \times \exp \left[ d \left( -\tilde{v}_{2\sigma}^2 + \frac{1}{2} \ln \left( 1 - \left( \frac{\tilde{v}_{2\sigma}}{\tilde{v}} \right)^2 \right) \right) \right]. \quad (C5) $$

When $d$ is very large, the above integral is dominated by the maximum of the function inside the exponential, which turns out to be located in $\tilde{v}_{2\sigma} = 0$. One can then perform a series expansion around $\tilde{v}_{2\sigma} = 0$ up to the second order. This results in a Gaussian integral which is easily integrated, and yields:

$$ I_1(\tilde{v}) \approx \sqrt{\frac{T_0}{d}} \frac{2\tilde{v}^2}{\sqrt{1+2\tilde{v}^2}}. \qquad (C6) $$
<figure><img src="image_13.png" /><figcaption>FIG. 12: (Color online) Free path length distribution for Maxwell molecules in two dimensions. The circles are the result of a DSMC simulation, while the solid line is the numerical integration of Eq. (D3). The dashed line is the large length approximation of Eq. (D4). All data are reported in units of mean free path.</figcaption></figure>

As for the second integral involved in the expression of $r(v)$ the simplest is first to perform the integration over $v_{2\sigma}$; hence one obtains, to leading order in $d$:

$$I_2(\tilde{v}) = -\sqrt{\frac{T_0}{2\pi}} \int_0^\pi \exp \left[ d \left( -\tilde{v}^2 \cos^2 \theta + \ln \sin \theta \right) \right] . \quad (C7)$$

The only maximum of the function in the exponential between 0 and $\pi$ is in $\theta = \pi/2$. Then, expanding as usual this function around the maximum up to second order, and extending the range of integration from $-\infty$ to $+\infty$, one finds:

$$I_2(\tilde{v}) \approx \sqrt{\frac{T_0}{d}} \frac{1}{\sqrt{1 + 2\tilde{v}^2}} . \quad (C8)$$

Finally, noting that $\Gamma(d/2)/\Gamma((d-1)/2) \stackrel{d \to \infty}{\sim} \sqrt{d/2}$, Eq. (12) follows.

## APPENDIX D: FREE PATH LENGTH DISTRIBUTION FOR MAXWELL MOLECULES

In this appendix, we investigate the free path length distribution for Maxwell molecules. As already mentioned, in this case the collision rate is a constant, $\omega$, independent of the velocity of the particle. It follows then that the free flight time probability is exponential:

$$P_{FT}^{M}(t|\mathbf{v}) = P_{FT}^{M}(t) = \omega e^{-\omega t} . \quad (D1)$$

Besides, the free path length distribution for a given velocity $v$ reads:

$$P_{FPL}^{M}(x|\mathbf{v}) = \frac{\omega}{v} \exp\left(-\frac{\omega}{v}x\right), \quad (D2)$$

and the free path length distribution is simply the average of the above probability over a Gaussian weight (for Maxwell molecules the on collision distribution $\phi_{coll}$ is still a Gaussian):

$$P_{FPL}^{M}(x) = \int d\mathbf{v} \phi(\mathbf{v}) P_{FPL}^{M}(x|\mathbf{v}) . \quad (D3)$$

This last expression can be expressed analytically in terms of Meijer G-functions [7], but here we shall focus only on the large length behavior, for which simpler expressions are available. In particular, the saddle point approximation gives:

$$P_{FPL}^{M}(x) \sim \frac{\sqrt{2\omega\Gamma\left(\frac{d-1}{2}\right)}}{\sqrt{T} 3^{(d-1)/2} \Gamma\left(\frac{d}{2}\right)} \exp\left(-3\left(\frac{\omega x}{\sqrt{T}}\right)^{2/3}\right) \\ \times {}_1F_1\left(\frac{d-1}{2}, \frac{1}{2}, \frac{3}{2}; \left(\frac{\omega x}{\sqrt{T}}\right)^{2/3}\right), \quad (D4)$$

which leads to a stretched exponential behavior at large $x$ ($\sim \exp(-x^{2/3})$ (cf Fig. 12).

## APPENDIX E: COLLISIONAL MOMENTS

Here we provide the expressions of the first collisional moments, defined by Eq. (42), for $\tilde{f}$ a Gaussian, and a Gaussian multiplied by a Sonine Polynomial. We denote by

$$\nu_n^{(0)} = \int d\mathbf{v}_1 v_1^n L_\lambda \frac{e^{-\frac{v_1^2}{2T}}}{(2\pi T)^{d/2}}, \quad (E1)$$

$$\nu_n^{(1)} = \int d\mathbf{v}_1 v_1^n L_\lambda \frac{e^{-\frac{v_1^2}{2T}}}{(2\pi T)^{d/2}} S_2\left(\frac{v_1^2}{2T}\right) . \quad (E2)$$

where $L_\lambda$ is the operator appearing in (28). Hence in the Gaussian approximation one has $\nu_n \equiv \nu_n^{(0)}$, while in the Sonine approximation $\nu_n \equiv \nu_n^{(0)} + a_2\nu_n^{(1)}$. The expressions of the first $\nu_n^{(i)}$ are:

$$\nu_0^{(0)} = - \left( \frac{(e^{-\lambda} - 1) \sqrt{T + T_0}}{\sqrt{2\pi}} \right) \quad (E3)$$

$$\nu_2^{(0)} = \frac{T^2 + 2 T T_0 + 2 T_0^2 - e^{\lambda} T (3 T + 2 T_0)}{e^{\lambda} \sqrt{2\pi} \sqrt{T + T_0}} \quad (E4)$$
$$ \nu_4^{(0)} = \left\{ \sqrt{\pi} \left( - (e^\lambda T^2 (8T^2 + 12TT_0 + 3T_0^2)) + T_0^2 (3T^2 + 12TT_0 + 8T_0^2) \right) - 4T(T+T_0) \left( \frac{3(-1+e^\lambda)\sqrt{\pi}T(T+T_0)}{4} + \frac{\sqrt{\pi}(e^\lambda T(2T+T_0) - T_0(T+2T_0))}{2} \right) \right\} / \left\{ \sqrt{2} e^\lambda \pi (T+T_0)^{\frac{3}{2}} \right\} \quad (E5) $$

$$ \nu_0^{(1)} = \frac{(-1 + e^\lambda) T^2}{8 e^\lambda \sqrt{2\pi} (T + T_0)^{\frac{3}{2}}} \quad (E6) $$

$$ \nu_2^{(1)} = \frac{T^2 (-3(1 + e^\lambda + 2(-1 + e^\lambda)) T^2 - (-6 + 22e^\lambda) T T_0 - 2(-3 + 8e^\lambda) T_0^2)}{8 e^\lambda \sqrt{2\pi} (T + T_0)^{\frac{5}{2}}} \quad (E7) $$

$$ \nu_4^{(1)} = \left\{ T^2 (-45(-1 + 5e^\lambda)T^4 + 6(28 - 132e^\lambda)T^3 T_0 - (-228 + 1000e^\lambda)T^2 T_0^2 + 4(36 - 128e^\lambda)T T_0^3 + \right. \\ \left. \qquad + 8(3 - 8e^\lambda)T_0^4) \right\} / \left\{ 8e^\lambda \sqrt{2\pi} (T + T_0)^{\frac{7}{2}} \right\} \quad (E8) $$

[1] L. Lue, J. Chem. Phys. **122**, 044513 (2005).

[2] M. H. Ernst, Phys. Reports **78**, 1 (1981).

[3] F. Coppex, M. Droz, and E. Trizac, Phys. Rev. E **69**, 011303 (2004).

[4] P. Visco, F. van Wijland and E. Trizac, to appear in Journal of Physical Chemistry B, arXiv:0801.2728.

[5] D. Enskog, Physik Z. **12**, 56 (1911).

[6] A. Puglisi, P. Visco, E. Trizac, and F. van Wijland, Phys. Rev. E **73**, 021301 (2006).

[7] M. Abramowitz and I. A. Stegun, Handbook of mathematical functions with formulas, graphs and mathematical tables (Dover, New York, 1972).

[8] F. W. Wiegel and J. P. J. Michels, Chem. Phys. Lett. **40** (1976).

[9] G. A. Bird, *Molecular Gas Dynamics and the Direct Simulation of Gas Flows* (Clarendon, Oxford, 1994).

[10] P. M. Morse and H. Feshbach, *Methods of Theoretical Physics* (McGraw-Hill, New York, 1953).

[11] S. Chapman and T. G. Cowling, *The mathematical theory of nonuniform gases* (Cambridge University Press, London, 1960).

[12] C. Giardinà, J. Kurchan, and L. Peliti, Phys. Rev. Lett. **96**, 120603 (2006).

[13] E. Trizac, Phys. Rev. Lett. **88**, 160601 (2002).

[14] J. Piasecki, E. Trizac and M. Droz, Phys. Rev. E **66**, 066111 (2002).

[15] E. Trizac and P. L. Krapivsky, Phys. Rev. Lett. **91**, 218302 (2003).

[16] T. Pöschel and S. Luding editors, *Granular Gases*, Berlin, 2001, Springer, Lecture Notes in Physics 564.

[17] S. Luding, in: T.A.S.K. Quarterly, Scientific Bulletin of Academic Computer Centre of the Technical University of Gdansk, **2** (3), 417 (1998).

[18] I. Pagonabarraga, E. Trizac, T.P.C. van Noije and M.H. Ernst, Phys. Rev. E **65**, 011303 (2001).

[19] D. L. Blair and A. Kudrolli, Phys. Rev. E **67**, 041301 (2003).

[20] D. Paolotti, C. Cattuto, U. Marini Bettolo Marconi and A. Puglisi, Granular Matter **5**, 75 (2003).

[21] E. Falcon, S. Aumaître, P. Évesque, F. Palencia, C. Lecoutre-Chabot, S. Fauve, D. Beysens and Y. Garrabos, Europhys. Lett. **74**, 830 (2006).

[22] P. Visco, A. Puglisi, A. Barrat, E. Trizac and F. van Wijland, Europhys. Lett. **72**, 55 (2005).

[23] M. H. Ernst and E. M. Hendriks, Phys. Lett. A **70**, 183 (1979).