The random gas of hard spheres

The inconsistency between the time-reversible Liouville equation and time-irreversible Boltzmann equation has been pointed out long ago by Loschmidt. To avoid Loschmidt's objection, here we propose a new dynamical system to model the motion of atoms of gas, with their interactions triggered by a random point process. Despite being random, this model can approximate the collision dynamics of rigid spheres via adjustable parameters. We compute the exact statistical steady state of the system, and determine the form of its marginal distributions for a large number of spheres. We find that the Kullback-Leibler entropy (a generalization of the conventional Boltzmann entropy) of the full system of random gas spheres is a nonincreasing function of time. Unlike the conventional hard sphere model, the proposed random gas model results in a variant of the Enskog equation, which is known to be a more accurate model of dense gas than the Boltzmann equation. We examine the hydrodynamic limit of the derived Enskog equation for spheres of constant mass density, and find that the corresponding Enskog-Euler and Enskog-Navier-Stokes equations acquire additional effects in both the advective and viscous terms. In the dilute gas approximation, the Enskog equation simplifies to the Boltzmann equation, while the Enskog-Euler and Enskog-Navier-Stokes equations become the conventional Euler and Navier-Stokes equations.


Introduction
It is known that the atoms in an electrostatically neutral monatomic gas interact via the Lennard-Jones potential [32]. At short range, the Lennard-Jones potential is inversely proportional to the 12th power of the distance between the centers of two interacting atoms. Due to this high order singularity, the mean field approximation [44], which is often used to close the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [9,11,28] for the long-range (that is, electrostatic or gravitational) potentials, cannot be applied to the molecular collisions, since the relevant spatial integrals diverge as the distance between the atoms approaches zero. To work around this issue, in molecular kinetic theory the Lennard-Jones potential interaction is replaced with the hard sphere collision model [12-15, 21, 26]. According to this model, the gas molecules are presumed to be rigid impenetrable spheres, which interact instantaneously and elastically according to the mechanics of the collision of two hypothetical billiard balls, which have nonzero mass, but do not possess the moment of inertia.
In the conventional setting of the hard sphere gas dynamics, one starts with the Liouville equation [12][13][14]21] for the probability density function of the complete system of all participating spheres. The Liouville equation is a homogeneous transport equation, whose characteristics are the straight lines of free flight of the spheres. In order to describe the collisions, the Liouville equation is endowed with special deflection conditions at the collision surface (the set of points in the coordinate space where a pair of spheres is separated by their diameter). These deflection conditions restrict the solutions of the Liouville equation to those for which the probability of a pair of spheres entering the collision surface equals the probability of a pair of spheres simultaneously exiting the collision surface along the corresponding direction of deflection. Assuming that the Liouville equation is already solved with the solution being known, one then computes the corresponding BBGKY hierarchy [9,11,28]. The lowest-order identity in the BBGKY hierarchy is then converted to the Boltzmann equation [10, 12-15, 21, 24, 26] via an approximation known as the Boltzmann hierarchy, and the subsequent factorization of the joint probability density of two spheres into the product of single-sphere densities.
Upon the examination of the standard closure of hard sphere dynamics [13,14,21], we observe that the derivation of the Boltzmann equation from the leading order BBGKY identity violates the assumptions under which the BBGKY identity itself was derived. Also, the observed contradictions are similar to those pointed out by Loschmidt [34]. Yet, the Boltzmann equation is known in practice to yield an accurate approximation to the observable gas dynamics. Thus, we propose that the Boltzmann equation instead originates, in a consistent manner, from a different model of hard sphere dynamics, which, on one hand, does not violate Loschmidt's observation, but, on the other hand, is on par with the hard sphere model at approximating real molecular interactions.
To avoid Loschmidt's objection [34] and derive the Boltzmann equation in a consistent fashion, we propose a new random process to model the underlying hard sphere dynamics, where the changes in the velocities of the spheres still obey the mechanics of the rigid sphere collision, but the "collisions" themselves are triggered by a point process. This random dynamical system possesses the infinitesimal generator, so that the corresponding forward Kolmogorov equation for its probability density is readily available via the integration by parts. The proposed random hard sphere process can approximate the deterministic collision process by increasing the intensity of the triggering point process via a parameter. We compute some of the steady states of the proposed process, and find that, while the conventional Boltzmann entropy can both increase and decrease, depending on a solution, the Kullback-Leibler entropy [29] between a solution and the steady state is a nonincreasing function of time. We also examine the structure of marginal distributions of the steady state in the limit of infinitely many spheres.
Then, we compute the forward equation for a single-sphere marginal distribution, under the assumption that the corresponding multisphere probability density is invariant under an arbitrary reordering of spheres. The closure that we use is, however, different from conventional -rather than using a direct factorization of the probability state (used, for example, in the conventional transition from the Boltzmann hierarchy to the Boltzmann equation), we take into account the structure of the previously computed steady state of the full multisphere system. We then find that, in the limit as the intensity of the point process increases to infinity, a variant of the Enskog equation [8,19,22,31,39,43] emerges, which, of course, simplifies to the Boltzmann equation in the dilute gas approximation. We then examine the hydrodynamic limit of the Enskog equation for spheres of constant mass density, as it appears to be physically plausible for atoms of the noble gases [16]. In this limit, we find that the resulting Enskog-Euler and Enskog-Navier-Stokes equations acquire additional nonvanishing terms, which are not present in the conventional gas dynamics equations originating from the Boltzmann equation [12-14, 24, 26]. These additional effects disappear in the dilute gas approximation, which yields the usual Boltzmann, Euler, and Navier-Stokes equations, respectively.

The hard sphere collision model and the Boltzmann equation
For the sake of clarity of the exposition, we start by citing the standard derivation of the Boltzmann equation from the collision mechanics of hard spheres [12][13][14]21]. For simplicity, we consider only two identical hard spheres, each of diameter σ. We denote by x and v the coordinate of the center and velocity of the first sphere, respectively, and by y and w the corresponding coordinate of the center and velocity of the second sphere. In the absence of contact, the spheres maintain their constant velocities v and w, while their respective coordinates x and y are given via Whenever the distance x − y between the centers of the two spheres equals their diameter σ, their velocities are changed, at that instance of time, to where v and w are the new values of velocities. Such a transformation preserves the momentum and kinetic energy of the system of the two spheres: Here and below, we assume that the total momentum of the system (the sum of velocities of the spheres) is zero without loss of generality, as otherwise the momentum can be set to zero via a suitable Galilean shift of the reference frame. The transformation above in (1.2) is fully symmetric; indeed, subtracting the first relation from the second, we obtain x − y x − y 2 . Scalar-multiplying by (x − y)/ x − y on both sides, we further obtain (1.5) ( Substituting the above expression into (1.2), we arrive at It is also easy to see that the Jacobian of the change of variables (v, w) → (v , w ) is unity; indeed, observe that, for (x − y) taken as a fixed parameter, which can be verified via the rank-one update lemma for determinants.
1.1. The Liouville problem for two spheres. Let F(t, x, y, v, w) denote the probability density function of the system of two spheres. Then, F satisfies the following conditions: • If the distance between the spheres is less than their diameter (that is, the spheres are overlapped), then the density is set to zero: This condition ascertains that the spheres are rigid and may not overlap. • If the distance between the spheres is greater than their diameter (that is, the spheres are separated), then F obeys the Liouville equation [12][13][14]: • The Liouville equation (1.9) is solved in the open set x − y > σ -therefore, a boundary condition is needed at the collision surface x − y = σ. This boundary condition is given via (1.10) F(x, y, v , w ) = F(x, y, v, w), where v and w are given in (1.2) as functions of x, y, v and w. • Observe that in the open set x − y < σ, F = 0 also solves the Liouville equation in (1.9) -albeit with zero initial and boundary conditions. This leads to a possible discontinuity of F at the collision surface x − y = σ, which will be taken into account below. • Assuming that there is a solution of (1.9) in both open regions x − y < σ and x − y > σ (with boundary conditions given via zero and (1.10), respectively), on the collision surface x − y = σ itself we assign F to be equal to the "outer" boundary condition (1.10). This definition of F at the collision surface means that F is continuous as x − y ↓ σ, and possibly discontinuous as In what follows, we assume for convenience that the spatial part of the domain has a finite volume, but no boundary effects other than those in (1.10). For example, one can assume that the space of coordinates is periodic (that is, if a sphere leaves the coordinate "box" through a wall, it immediately re-enters from the opposite wall).
According to (1.10), the probability that a pair of spheres enters the collision surface at the point (x, y) and time t with velocities (v, w) is the same as that of a pair of spheres exiting at the same point and time, but with the velocities (v , w ). The condition in (1.10) preserves the normalization of F separately in the region x − y > σ.
Note that the formulation of the problem in (1.8)-(1.10) is not fully equivalent to the actual dynamics in (1.1) and (1.2). For the formulation (1.8)-(1.10) to correspond precisely to (1.1) and (1.2), one needs the characteristic curves of the Liouville equation in (1.9) to be the trajectories of the system in (1.1) and (1.2). However, this is not the case in (1.8)-(1.10) -the characteristic curves are straight lines given via constant parameters v and w, which pierce the collision surface x − y = σ. Instead, the "collisions" in (1.8)-(1.10) are implemented via reassigning the values of F between different characteristics according to (1.10).
1.2. The BBGKY identity for two spheres. Here we follow [14] and derive the BBGKY hierarchy (which, for two spheres, consists of a single identity) and, subsequently, the Boltzmann equation. In what follows, we take advantage of the fact that F = 0 inside the overlapped region x − y < σ also satisfies the Liouville equation in (1.9) with zero initial and boundary conditions.
Let us assume that a solution F is computed for the Liouville problem (1.8)-(1.10) for both overlapped and non-overlapped regions, so that the relation in (1.9) is no longer an equation, but rather an identity. In what follows, we assume that F is symmetric under the permutations of two spheres -that is, F(t, x, y, v, w) = F(t, y, x, w, v). Our goal is to manipulate the Liouville identity in (1.9) so as to obtain an appropriate identity for the marginal distribution of a single sphere First, we integrate the identity in (1.9) in dydw over the whole space, which includes both overlapped (that is, x − y < σ) and non-overlapped (that is, x − y > σ) regions: Next, we exchange the order of differentiation and integration above. For the t-derivative term, the exchange is done in a direct manner, since the dydw-integration is unrelated to t: However, due to the discontinuity in F at the collision surface x − y = σ, one cannot directly swap the spatial derivatives and the spatial integration. For the term with yderivative, we use the Gauss theorem: Above, dS x (y) is the surface area element at the point y of the sphere of radius σ centered at x, the unit normal vector n x (y), which originates at x in the direction of y, is given via Also, F (1.8) denotes the boundary value for x − y < σ (which is zero), F (1.10) denotes the boundary value for x − y > σ, given via (1.10), and we replace F (1.10) with F in the last identity, as it was agreed above that F at the collision surface x − y = σ is assigned the boundary value given via (1.10). For the term with the x-derivative, we need to use the generalized Leibniz rule. First, we split Then, for the first term in the right-hand side we write For the second term, we repeat the calculations in a similar way: Adding the two terms, recalling that F (1.8) is zero, F (1.10) is F itself on x − y = σ, and using (1.11), we arrive at the identity Fv · n x (y)dS x (y).
Combining the above expression with (1.13) and (1.14), and substituting into (1.12), we arrive at the identity For further convenience, we can rename the dummy variables of integration so that the surface integration is computed over a unit sphere (that is, over dn). Namely, since y follows the surface of the sphere of radius σ centered at x, we denote In the above variables, the identity in (1.2) becomes where we note that there is no dependence on σ. The integral over the collision surface, written in the new variables, yields the following identity: In order to transform the surface integral in the right-hand side above into the Boltzmann collision integral, one further needs to rearrange the surface integration with the help of (1.10). Consider the distance y(t) − x(t) between the spheres as a function of time.
The sign of its time derivative indicates whether the spheres are approaching or escaping each other: Clearly, whenever the time derivative above is positive (that is, the spheres escape each other), the dot-product (w − v) · n is positive, and vice versa. Following [13,14,21], we use the condition in (1.10) for F and rearrange the collision integral in (1.23), with help of the Heaviside step-function Θ(x), as follows: Above, we split the collision integral into two hemispheres -one where (w − v) · n > 0 (that is, the spheres are escaping each other), and another one where the spheres are approaching each other. In the first hemisphere, we replace F(v, w) with F(v , w ), with v = v (v, w, n) and w = w (v, w, n) given via (1.22). This is a valid rearrangement since the condition (1.10) requires F(v, w) and F(v , w ) to be equal anywhere on the surface of integration. In the second hemisphere, we change the sign of n to the opposite. Then, we merge the integrals. The resulting BBGKY identity is given via Clearly, the BBGKY identity above is the same identity as in (1.23) [13,14,21], assuming that F is symmetric under an arbitrary permutation of the spheres. This results in the multiple BBGKY identities for marginal distributions of various orders, chain-linked to each other starting with the highest-order. The lowest-order identity in the BBGKY hierarchy for K spheres is given via where F (2) is the two-sphere marginal distribution of the full K-sphere density F: The Boltzmann equation [10, 12-14, 24, 26] is obtained from the BBGKY identity (1.27) in two following steps. First, one assumes that which is known as the Boltzmann-Grad limit [26]. As σ → 0, the following assumption is made [14,21]: This assumption transforms the BBGKY hierarchy into what is known as the Boltzmann hierarchy [14,21], whose lowest-order relation is given via Next, the joint two-sphere marginal density F (2) above is approximated as follows: Substituting the approximations in (1.32) into (1.31), one arrives at The relation above is known as the Boltzmann equation [10, 12-14, 24, 26]. Unlike the BBGKY identity in (1.27), built upon a (presumably known) solution of the Liouville problem for K spheres, (1.33) is treated as a closed, self-contained equation for the marginal distribution f . The factor (K − 1)σ 2 in front of the integral above can be changed to Kσ 2 , since the Boltzmann equation is the result of the Boltzmann-Grad limit in (1.29).

Inconsistencies in the derivation of the Boltzmann equation
Here we point out some contradictions in the conventional derivation of the Boltzmann equation (1.33) from the Liouville problem (1.8)-(1.10), presented in the previous section. We start with a brief summary of the derivation.
(1) The Liouville equation (1.9) by itself does not have any collision effects. The effect of collision is imposed separately for every pair of spheres with coordinates x and y, on the surface x − y = σ. (2) Due to the effect on the collision surface, the probability density of states F is discontinuous on this surface -it is zero for all x − y < σ according to (1.8) (as the spheres are impenetrable) and is generally nonzero otherwise.
(3) Due to the discontinuity of F on the collision surface, the collision surface integrals emerge according to the Gauss theorem and the Leibniz rule, when the BBGKY hierarchy is constructed. The resulting identity is given by (1.23) -observe that it is valid for any F, which satisfies (1.8), (1.9), and is discontinuous on the collision surface x − y = σ. (4) In addition to the discontinuity, the velocity deflection condition in (1.10) is imposed on F. This condition allows to rewrite the collision integral in (1.23) in the form (1.26). Note that it is the same exact integral as obtained originally in (1.23) via Gauss theorem and Leibniz rule, only written in a different manner thanks to (1.10). Similarly, one can obtain the K-sphere BBGKY identity in (1.27) [14,21].
After the BBGKY identity is obtained in (1.26) (or its K-sphere version (1.27)), the following additional steps are taken to obtain the Boltzmann equation: Remarkably, none of the steps in (a)-(b) above are consistent with the conditions under which steps in (1)-(4) were valid. Below we elaborate on the inconsistencies between (a)-(b) and the conditions under which (1)-(4) were obtained.

A contradiction between the Liouville problem and the Boltzmann hierarchy.
Observe that if F is a distribution of K rigid spheres, then x, x, v, w) = 0, for all x, v, w and for any σ.
The reason for this is the following. Observe that F (2) with the arguments given above is the two-sphere marginal distribution of the full probability density F, where the latter is computed at the point where the coordinates of the first and second spheres are identical.
In such a state, the first and second spheres are fully overlapped for any value of the diameter σ. Due to the impenetrability requirement in (1.8) (which can be extended onto K spheres in a direct manner), the value of F at this state is guaranteed to be zero irrespectively of the states of all other spheres: The identity above, in turn, means that the marginal integral (1.28) of such an F is also zero, which leads to (2.1). Now, observe that the marginal distribution F (2) , computed for the states of the form in (2.1), is used in the collision integral of the lowest-order equation of the Boltzmann hierarchy in (1.31), which was obtained from the BBGKY equation in (1.27) via (1.30). Comparing the right-hand side of (1.31) with (2.1), we find that the collision integral in (1.31) must always be zero, if the impenetrability condition in (1.8) is indeed satisfied. Therefore, if the approximation of the BBGKY equation (1.27) via the Boltzmann hierarchy equation (1.31) is formally valid, this could only mean that any physical effects from the collision integral of the BBGKY equation in (1.27) somehow vanish in the Boltzmann-Grad limit (1.29). However, even if that is indeed the case, it is not the dynamical regime which is physically interesting or relevant -in typical applications of gas dynamics, the molecular collisions and the associated effects (such as viscosity and heat conductivity) are usually important.

A contradiction between the
, for all x, v, w, and all unit vectors n, where v = v (v, w, n) and w = w (v, w, n) are given via (1.22). However, solutions of the Boltzmann equation in (1.33) which are not simple traveling waves clearly may not have such a property, since (2.3) automatically turns the Boltzmann collision integral in the right-hand side of (1.33) into zero. Additionally, observe that, since (2.3) automatically sets the Boltzmann collision integral in (1.33) to zero, any solution, which is not a traveling wave, is required to violate (1.10). In other words, not only the violation of the deflection condition in (1.10) is "overlooked" in the Boltzmann equation, but it is, in fact, used as the quintessential device to create its nontrivial solutions. This happens despite the fact that (1.10) is the fundamental property of all solutions of the Liouville problem, whether stationary or not, which, together with (1.8), leads to the BBGKY identity in (1.26), from which the Boltzmann collision integral and the Boltzmann equation are subsequently derived.
Therefore, we conclude that what we observe in the conventional derivation of the Boltzmann equation is a logical fallacy -first, a chain of identities is established under the assumption that certain conditions hold; then, the resultant identity is replaced with an equation whose nontrivial solutions must violate requisite conditions under which the identity was derived to begin with.
2.3. Reversibility and Loschmidt's objection. One may try to "escape" the above reasoning by assuming that the transitions from BBGKY to the Boltzmann hierarchy in (1.30) and further to the factorization in (1.32) are only valid for incident velocities, but not recedent. More precisely, one may put forth an ansatz that, for recedent directions where v and w are functions of v, w and n, given in (1.22). In such a case, one can claim that the factorization in (1.32) "is not required" to satisfy (1.10), so that the relation (2.3) does not have to apply between the incident and recedent directions.
First, we observe that the ansatz in (2.4) is undoubtedly correct, and not just for the recedent, but also for incident directions, since F (2) (t, x, x, v, w) = 0 for all x, v, w and σ, as we already pointed out above in (2.1). However, even if one somehow "overlooks" this fact, still such a hypothetical dichotomy between the incident and recedent directions is not only unfounded, but is also clearly in contradiction with the original formulation of the Liouville problem (1.8)-(1.10). The reason is that in the Liouville problem all collisions are exactly reversible, and all states of the system reproduce themselves in reverse order if one uses the terminal coordinates and negatives of terminal velocities in place of an initial condition.
In the case of such a reversal, the incident density states become recedent and vice versa, however, the formulation of the Liouville problem in (1.8)-(1.10) remains invariant. Thus, the same exact reasoning as above involving the BBGKY hierarchy, the Boltzmann hierarchy in (1.30) and factorization (1.32) can subsequently be applied to the reversed states as well, in which case the formerly recedent states will be expressed via (1.30) and (1.32). Or, conversely, if the ansatz in (2.4) indeed holds for recedent states, then it also automatically holds for the incident states, in which case the transition from (1.27) to (1.33) cannot be carried out via (1.30) and (1.32) (which is in fact the case due to (2.1)).
A similar objection to the form of the Boltzmann collision integral has been posed by Loschmidt [34], who pointed out that, due to the entropy inequality (also known as the H-theorem [10]), the Boltzmann equation is time-irreversible, while the underlying dynamics of hard spheres, from which the collision integral is derived, are time-reversible. It is commonly known as Loschmidt's paradox, even though Loschmidt's observation is entirely valid and does not constitute a "paradox" by itself. Instead, the "paradox" appears due to the logical fallacy hidden between the formulation of the Liouville problem for hard spheres in (1.8)-(1.10) and the resultant Boltzmann collision integral in (1.33), as we exposed above.
2.4. Our proposal to amend the situation. At the same time, there is no doubt that the Boltzmann equation (1.33) is a very accurate model of a dilute gas, which is confirmed by numerous observations and experiments. This means that the collision integral in the right-hand side of (1.33) is a valid approximation of the statistical effects of molecular interactions in practical scenarios, even if there is a conceptual flaw in its conventional derivation [13,14,21]. In what follows, we propose a different model of hard sphere interactions in which such a problem does not manifest, and which also leads to the Boltzmann equation in the dilute gas approximation -albeit in a different way, via the Enskog equation [8,19,22,31,39,43].
From what is presented above, it is clear that the formulation of the dynamics of the spheres needs to be such that the collision integrals in the BBGKY equation (1.27) possess their form independently of what their integrand is. For that, we need to disentangle the instantaneous changes of sphere velocities from the properties of F in the Liouville equation, and arrange for them to be the property of the equation itself. This is possible to do if the underlying dynamical process possesses an infinitesimal generator. Generally, let z(t) be a Markov process in a d-dimensional Euclidean space. Assuming that z(t) = z is known at the time t, let us consider the conditional expectation E[ψ(z(t + ε))] of a suitable function ψ : R d → R, for some ε > 0. Assume that the following limit exists: where the (generally, integro-differential) linear operator L, which is independent of ψ, is called the infinitesimal generator. Clearly, L describes the underlying dynamics of the process -if z(t) is the phase state of the system of particles/spheres containing their coordinates and velocities, then L describes the interactions between the spheres. For example, for the point particles interacting via a long-range potential [44], L includes the gradient of the potential function of the force field. With L specified explicitly, it is easy to see that the forward equation for the probability density F(t, z) can be derived in a straightforward fashion [6,23] with the help of the adjoint operator L † , assuming that the latter can also be explicitly obtained. Integrating the above identity over R d against F(t, z)dz, we have First, it can be shown [23] that Next, with help of the adjoint operator L † , we write Combining the last two identities and stripping the ψ-integral, we arrive at the general form of the forward Kolmogorov equation [23,37] (also known as the Fokker-Planck equation [40]) for F alone: Since the infinitesimal generator L is independent of F, its adjoint L † is also independent of F. Thus, the form of the collision integral in the corresponding BBGKY hierarchy will be defined entirely by L † . This, in turn, will allow to replace F with a suitable closure without violating any prior assumptions on the integrand of the collision integral.
Below we introduce the infinitesimal generator into the hard sphere collision dynamics by randomizing the times at which the velocity jumps occur. More specifically, the velocity jumps will still occur according to (1.2), except that, rather then lying precisely at the collision surface, the coordinate (x, y) of the jump will be selected at random along the trajectory from within a small "tolerance" interval around the collision surface. In such a case, the conditional expectation E t [ψ] will be a differentiable function of t, even though the velocities of the spheres will still change instantaneously. Since a possibility arises that the spheres become overlapped, we will make provisions for their unimpeded subsequent separation. The aforementioned tolerance interval can be made as small as necessary, and the probability of the jump can be made as large as necessary, so as to mimic the deterministic collisions with prescribed accuracy. Also, Loschmidt's objection [34] will be avoided since the underlying dynamics will be inherently irreversible.
One can argue that the process we are to propose is a pure abstraction -clearly, the actual atoms in a gas are not known to interact randomly. However, note that the "hard sphere collision" is also an abstraction -in reality, no observable physical object can change its velocity instantaneously. In particular, the atoms in a real gas do not collide with each other instantaneously, but instead interact via the Lennard-Jones potential [32]. Therefore, neither the deterministic hard spheres, nor our random process are completely accurate models of molecular interaction on a microscopic level. What defines the quality of a model here is its ability to describe macroscopic effects in a gas while being mathematically and logically consistent with its own abstract foundation.

Random dynamics of hard spheres
Here we present the details of the random dynamical system which models collisions of hard spheres, and at the same time possesses the infinitesimal generator. As before, first we consider the dynamics of two spheres with coordinates x and y, and velocities v and w, respectively. The dynamics of coordinates are, of course, given via (1.1). For the dynamics of velocities, we consider two separate configurations: (1) Collision configuration. The collision configuration is given by the following two criteria, which must hold concurrently: (a) The distance x − y between the centers of spheres satisfies where 0 < α 1 is a constant parameter. The condition above signifies that x − y /σ ≈ 1, within ±α-tolerance (that is, the spheres are separated approximately by their diameter). We will say that two spheres are in the "contact zone" whenever the condition in (3.1) holds for the coordinates of their centers. (b) The distance x − y between the centers of spheres also diminishes in time, that is, This condition signifies that the spheres are approaching each other. In the collision configuration, the velocities (v, w) may be randomly transformed according to (1.2), with a specified probability. For now, we write informally that in the collision configuration the velocities evolve according to (3.3) dv(t) = −dw(t) = random jump process, which changes the velocities instantaneously according to (1.2). The jump process must be random for the expectation of a jump to be a continuous function of t, despite the fact that the velocity jumps themselves remain instantaneous. The continuity of the expectation is the key property which allows the process to possess the infinitesimal generator. The exact representation of the requisite jump process will be provided below. (2) Free-flight configuration. If the two conditions above in (3.1) and (3.2) do not hold, then the spheres are in the free-flight configuration. In this case, the velocities are constant in time: At this point, we need to choose a appropriate type of the random jump process, which will be used as a "trigger" to change the velocities instantaneously. The simplest process which is suitable for the given task is the point process [18] -that is, a scalar, piecewiseconstant random process which starts at 0 and occasionally increments itself by 1, randomly and independently of previous increments. Before describing the random hard sphere process in detail, we formulate a general dynamical process driven by a point process, and compute the explicit form of its infinitesimal generator.

A dynamical system driven by an inhomogeneous point process.
Here we follow the theory in Section 7.4 of [18] (the original results are presented in [38]). Let (Ω, F , P) be the probability space, equipped with a filtration {F t }, t ∈ R ≥0 , such that for 0 ≤ t 1 ≤ t 2 < ∞, the corresponding sigma-algebras are nested as F t 1 ⊆ F t 2 ⊆ F . Let h : R ≥0 → R >0 be a bounded, strictly positive random variable adapted to {F t }. Let m : R ≥0 → Z ≥0 be a {F t }-adapted inhomogeneous point process with the conditional intensity h(t), and the compensator τ(t), given via By choosing the random variable h(t) appropriately, we can regulate the "temporal density" of jump points of m(t); indeed, depending on the magnitude of h(t), the jump points of m(t) can either arrive in a statistically rapid succession, or, to the contrary, disperse farther away from each other. For what is to follow, it is a necessary requirement that the intensity h(t) is a random variable -both m(t) and h(t) are, however, adapted to the same filtration {F t }, so that if the sequence of values of h up to t is given, so is the sequence of values m. We now denote the corresponding jump process of m(t) via ∆m(t): Above, the notation "t−" denotes the left-limit at t. Let M(t, ·) be the corresponding random measure of ∆m(t): As ∆m(t) only assumes the values of either 0 or 1, M(t, ·) is concentrated at 1.
Let us now define a stochastic process z(t) on a Euclidean space R d as follows: where f , g : R ≥0 → R d are suitable (for the purpose of stochastic integration above) random variables, adapted to {F t }. Clearly, z(t) is also {F t }-adapted by construction. Our task here is to compute the infinitesimal generator of z(t), that is, for a test function ψ(z), we would like to compute To compute the expectation above, we need to adapt the integral form of z in (3.8) to the Itô formula for Lévy-type stochastic integrals (see Chapter 4 of [6]). However, the problem is that the random measure in the right-hand side of (3.8) is not that of the standard Poisson point process with constant intensity, but that of the point process with random intensity h(t), defined above. Our first step here is, therefore, the transformation of the stochastic integral in the right-hand side of (3.8) to an integral against the random measure of a standard Poisson point process. According to Theorem 7.4.I of [18] (also see [38]), the random point process n : is the standard Poisson point process with intensity 1, where τ(t) is the random, albeit {F t }-adapted, compensator process of m(t), defined above in (3.5). Subsequently, we can write the stochastic integral in (3.8) from t to t + ε as with N(t, ·) being the random measure of the standard Poisson point process with intensity 1; in particular, its intensity measure EN [6] is given via Substituting the above integral into (3.8), we write Next, recalling the Itô formula for Lévy-type stochastic integrals in Chapter 4 of [6], we write Applying the conditional expectation on both sides, we obtain, for the left-hand side and the first term in the right-hand side, where in the second expression the conditional expectation disappears in the leading order term because the z(t) and f (t) are both {F t }-adapted. The expectation of the stochastic integral is somewhat more complicated. First, observe that , where in the right-hand side both leading order terms are {F t }-adapted. The integral then can be expressed as where in the second identity the dummy variable of integration s was replaced with its starting value τ(t).
At this point, observe that in the leading order integral above, the limits of integration, as well as the integrand, are {F t }-adapted. Thus, applying the conditional expectation to the leading order integral yields Assembling the terms together, we obtain the infinitesimal generator in the form where in the last identity we denote z = z(t), for brevity. It is worth noting that the form of the infinitesimal generator above extends naturally onto the free-flight configuration of the dynamics -it suffices to set the variable intensity h(t) = 0 above. In this case, z(t) in (3.8) will be driven solely by the integral over the vector field f (t) alone.

Random dynamics of two spheres.
To adapt the general stochastic process in (3.8) to the dynamics of spheres, we need to relate z(t), f (t), g(t) and h(t) to the variables of the dynamics. Obviously, z(t) is the state vector of the system, and thus it is going to incorporate the coordinates and velocities of both spheres: Subsequently, f (t) is related to the deterministic component of the dynamics, which is the evolution of the coordinates x(t) and y(t) for given velocities v and w according to (1.1): To specify g(t), we observe that the instantaneous change of velocities in (1.2) can be written, with help of the jump process ∆m(t), as where x = x(t), y = y(t). Therefore, we can define g(t) via Then we write the process in (1.1)+(3.22) as the following stochastic differential equation [6]: with z, f (z) and g(z) given via (3.20), (3.21) and (3.23), respectively. It remains to specify the variable intensity h(t), which should activate the point process m(t) when both (3.1) and (3.2) hold concurrently, and be zero in the free-flight configuration. Here, we define h(t) as Above, λ > 0, 0 < α 1 are constant parameters, and δ α (x) is the standard mollifier of the delta-function δ(x), given via where the constant parameter c ensures the proper normalization. For a function ψ(x), which is continuous at zero, we thus have that is, δ α can serve as the delta-function in the limit α → 0, while remaining smooth for finitely small α. We will denote the anti-derivative of δ α as Θ α : Clearly, as α → 0, Θ α (x) becomes the usual Heaviside step-function.
Observe that, in the collision configuration (3.1)-(3.2), the variable intensity of the point process in (3.25), as required in [18,38], is indeed F t -adapted, strictly positive and bounded, since, first, z(t) is F t -adapted by construction, and, second, whenever (3.1) and (3.2) hold, we have where E is the constant energy of the system of two spheres. The dynamics in (3.24)-(3.25) function as follows: • If x − y is away from σ, or if x − y is growing in time (that is, the spheres are escaping each other), then the triggering point process is not present, and the spheres are moving with constant velocities according to (3.40). Accordingly, the intensity h(z(t)) in (3.25) is zero in the infinitesimal generator of the process, and thus the generator consists of the free-flight term only. • Once x − y is close enough to σ and, at the same time x − y is decreasing in time, the spheres enter the contact zone (both conditions (3.1) and (3.2) are satisfied). In this case, the collision-triggering point process becomes present, with the intensity h(z(t)) being strictly greater than zero. Then, there are two possibilities: -A jump in the point process arrives so that the spheres "collide" according to (1.2). In this case, the spheres start escaping each other (so that (3.2) no longer holds) and the triggering point process is no longer present. In the infinitesimal generator, the Heaviside function becomes zero and so does the intensity h(z(t)). -A jump does not arrive, so that eventually δ ασ ( x − y − σ) decays back to zero together with the intensity h(z(t)) of the point process; in this case, the spheres pass through each other without interaction. In either scenario, the point process m(t) becomes dormant until the spheres approach each other again.
The existence of strong solutions to (3.24)- (3.25) in the collision configuration (3.1)-(3.2) is a subject that merits a separate discussion. For the purpose of this work, we will assume that bounded strong solutions are sufficiently generic for typical initial conditions, so that the corresponding statistical formulation (the forward Kolmogorov equation) of the dynamics is reliable enough for the description of large ensembles of solutions.
The definition of the jump intensity above in (3.25) also indicates that the probability that the jump in the point process does not arrive during the collision window is e −λ , regardless of the values of σ or α. Indeed, observe that one can write which means that if the contact zone is traversed completely (that is, the jump has not arrived), then the compensator τ(t) in (3.5) is incremented by λ. Clearly, to mimic the collisions of hard deterministic spheres in Section 1, one eventually needs to take α → 0 and λ → ∞, so that, first, the contact zone (3.1) becomes infinitely thin, and, second, the jump arrives with probability 1 whenever the spheres are in the collision configuration.
The corresponding infinitesimal generator of (3.24)-(3.25) is given via Changing back to the original variables x, y, v and w, we write the infinitesimal generator of (3.24) in the form where v and w are the functions of x, y, v and w given in (1.2). The jump portion of the generator above in (3.32) is not translationally invariant, and thus the process in (3.24) is not a Lévy process. However, it is a Lévy-type Feller process [6,20], whose infinitesimal generator can be reformulated in the Courrège form [17] via an appropriate change of variables. The next step is to obtain the corresponding forward Kolmogorov equation [23,37,40] for the probability density of states of the system, which is easily achieved via the integration by parts. Let F(t, x, y, v, w) be the corresponding probability distribution of the random process above. We can then integrate (3.32) against F and obtain, with help of (2.7), where dV 2 is the volume element of the coordinate space, and dS 2 is the area element of the sphere of zero momentum and constant energy (the subscript denotes the number of spheres in the system).
Above, the terms with spatial derivatives in x and y can be integrated by parts, with the condition that the boundary effects are not present. For the part with ψ(v , w ) we can write, for fixed x and y, (3.34) ψ where we used (1.5), (1.6) and (1.7) (note that v and w remain on the same zero momentum -constant energy sphere), and in the last identity renamed v → v, w → w and vice versa, since the integral occurs over the same velocity sphere. As a result, we can recombine the terms as Assuming that ψ is arbitrary, we can strip the integral over ψ and obtain the equation for F alone: Unlike the Liouville problem (1.8)-(1.10), here observe that the effect of collisions is present in the equation itself, and is not contingent upon additional properties imposed on F.

Extension to many spheres.
Here we extend the previously formulated dynamics onto K spheres, with the corresponding coordinates x i (t) and velocities v i (t), 1 ≤ i ≤ K. Observe that we have K(K − 1)/2 possible pairs of spheres. In order to define their random interactions, we introduce K(K − 1)/2 independent instances of the point process, each assigned to the pair of i-th and j-th spheres.
where the two nonzero entries in the G ij -vector above are in the (K + i)-th and (K + j)-th slots. Let m ij (t) be the set of K(K − 1)/2 independent inhomogeneous point processes with conditional intensities h ij (t), given via where h(z) is defined in (3.25). Let M ij (t, ·) be the set of corresponding random measures for m ij (t). Then, the K-sphere dynamics is defined via the following system of stochastic differential equations: The process above in (3.40) is also a Lévy-type Feller process [6,20], which lives on the sphere of zero momentum and constant energy As above for two spheres, here we assume that the total momentum of the system is zero without loss of generality. The infinitesimal generator of such process is, apparently, given via In the x i and v i variables, this translates into where the notation ψ(v i , v j ) specifies that all velocity arguments in ψ are set to the corresponding velocities v k , except for i-th and j-th, which are set to v i and v j via (1.2).
Observe that the dynamics in (3.40) is a direct extension of the dynamics of two spheres in (3.24) onto multiple spheres -the evolution of the coordinates is governed by the same equations, and the velocities of each sphere are coupled to all other spheres via the independent point processes m ij (t). In particular, there is no provision for a collision of more than two spheres at once (which is often discussed in the literature [12,14]); however, given the fact that the collisions in the hard sphere model are instantaneous, we will assume that the event of a three-sphere collision is improbable for a "generic" initial condition.
It is interesting that the properties of an infinitesimal generator similar to (3.43) were studied in [39] in the same context (that is, a system of K particles interacting according to (1.2)). However, the collision part of the infinitesimal generator in [39] was scaled differently, as if the intensity parameter λ in (3.43) was set to (ασ) 3 . Thus, as α → 0, the particles described by the generator in [39] ceased colliding upon contact. It is, however, unclear where the infinitesimal generator in [39] comes from; the generator in [39] appears to be postulated, rather than derived from an underlying SDE.
To obtain the corresponding forward equation, we follow the same principle as for the two spheres above. First, we integrate (3.43) against the probability density F and obtain where dV K is the volume element of the coordinate space of the K spheres, dS K is the area element of the corresponding velocity sphere of zero momentum and constant energy, and the term with the spatial derivatives is integrated by parts assuming that the boundary terms vanish. Now, for all terms with v i and v j we have, for fixed coordinates, where we used (1.5), (1.6) and (1.7), and observed that, for fixed coordinates, the variables v i and v j sample the same zero momentum and constant energy sphere as do v i and v j . Finally, stripping the integral over ψdV K dS K , we arrive at the forward Kolmogorov equation in the form (3.46) ∂F ∂t Note that the equation above in (3.46) admits solutions which are symmetric under the reordering of the spheres. These solutions are "physical", that is, they correspond to real-world scenarios where it is impossible to statistically tell the spheres apart. We solve (3.36) using the method of characteristics; namely, we treat the advection term of (3.36) as the ordinary, scalar spatial derivative in the direction of (1, v, w), with the directional parameter denoted as s. With this, t, x and y are given via the straight line (4.1) (t(s), x(s), y(s)) = (0, x 0 , y 0 ) + s (1, v, w).
On this straight line, (3.36) becomes Here, the two terms in the right-hand side are never nonzero simultaneously, due to the Heaviside step-functions. We assume that the initial condition (x 0 , y 0 ) satisfies x 0 − y 0 > σ(1 + α), and (x 0 − y 0 ) · (v − w) < 0, so that the spheres are not overlapping and on the "collision course". Then, as the characteristic traverses the contact zone for the first time, we have which results in the solution In the original variables, the solution translates into Sending α → 0 has no effect other than making the contact zone infinitely thin, and, subsequently, the transition between F 0 and e −λ F 0 instantaneous. As we continue along the characteristic, it eventually traverses the contact zone again, and exits the overlapped state. In this configuration, we have Here, the right-hand side contains the values of F(v , w ), which are carried along the characteristic curves with directions given via v (s) = v (x(s), y(s), v, w), w (s) = w (x(s), y(s), v, w), and should be treated as an external forcing. However, we already know from what we found above that we can express, for some t * , Without loss of generality, we can take t * large enough so that the second term under the exponent is λ (that is, the point back in time t * is sufficiently far away from the collision zone), which gives The equation for F thus becomes Assuming that the initial condition F 0 is continuous, and sending α → 0, we find, upon traversing the collision zone outward, where F 0 denotes the initial value of the solution that is carried from the outside to the point of the collision zone where our characteristic exits the overlapped region. Assembling the pieces together, we finally arrive at the following result in the limit as α → 0: after the characteristic traverses the overlapped region completely, the solution is given via Subsequently, sending λ → ∞ leads to F 0 being completely replaced by F 0 once the characteristic fully traverses the overlapped region x − y < σ. This is, obviously, the same solution as the one of the Liouville problem for hard spheres (1.8)-(1.10). Thus, as α → 0 and λ → ∞, the solution of (3.24)-(3.25) approximates the solution of (1.1)-(1.2) in the sense that, for the same initial condition, the corresponding probability densities of (1.8)-(1.10) and (3.36) converge to each other on the same characteristic curve.
The corresponding differential equation along the straight line is subsequently given via (4.14) dF ds + To compute steady solutions, it suffices to look for solutions with the property F(v i , v j ) = F, for all index pairs (i, j). In such a situation, the Heaviside step-functions in (3.46) are multiplied by identical F's, and thus coalesce into 1. The equation for F on the characteristic (4.12) thus becomes The equation above in (4.15) can obviously be integrated via separation of variables, but before we proceed with that, let us recall (3.30) and observe that the coefficient in front of F in (4.15) is by itself the time derivative of λΘ ασ , for each index pair (i, j). Then, the separation of variables yields the equation with the following solution on the straight line (4.12): Above, G is an arbitrary function of the starting point (x 10 , . . . , x K0 ) and the fixed velocity vector (v 1 , . . . , v K ).
To substitute the obtained F into (3.46), we change the notations back to the original variables: where the subscript refers to the number of spheres in the system, and we denote Above, S K is the area of the sphere of zero momentum and constant energy for K spheres. Observe thatF K is normalized to 1 and is a probability density by itself.

Now, in order to satisfy the condition
preserves the momentum and energy in (3.41), we have to express G as a function of these two quantities, at the same time preserving the form in (4.20). Apparently, the only form of G which is consistent with both conditions is given via where the division by K is included for convenience. It is not difficult to see that At this point we recall that v 1 , . . . , v K belong to the set of zero momentum in (3.41). This forces G to lose the dependence on t and on the second argument, while the third argument becomes the constant energy E: Thus, we found a family of solutions F with F(v i , v j ) = F(v i , v j ) for all pairs (i, j), which are the steady states of (3.46), uniform on the velocity sphere of zero momentum and constant energy. Apparently, there are infinitely many such steady states for a given system of hard spheres. However, note that G above in (4.22) is a function of the center of mass of the system of spheres, which is invariant under the dynamics of (3.40). As it will become important later, observe that there is a family of initial conditions of (3.46) which has a uniform distribution of the center of mass of the spheres. For such a family of initial conditions, G = 1, and we arrive at the steady solution which consists purely ofF K in (4.15). Qualitatively,F K (x 1 , . . . , x K ) behaves as follows: •F K (x 1 , . . . , x K ) is constant outside contact zones (that is, the regions in the coordinate space where the mollifier δ ασ ( x i − x j − σ) > 0), and its transitions within contact zones are given via the exponents of anti-derivatives of the mollifier in (3.26). • Outside contact zones, the following condition holds: (4.23)F overlap = e −λ(n−1)F non-overlap , where n is the number of the simultaneously overlapping spheres at a given point (x 1 , . . . , x K ). Observe that the preceding discussion implicitly assumes that the steady state of the form in (4.19) and (4.22) is attracting for a given initial condition. This is certainly not the case for all initial conditions -there exist other steady states for which F(v i , v j ) = F(v i , v j ); for example, those which are not supported in the contact zones, such that the spheres do not interact at all. However, below we will argue that the stateF K in (4.19) is the most likely candidate for a "physical" (that is, statistically most common) steady state for (3.46).

Entropy inequality.
For the K-sphere dynamics in (3.46), the conventional Boltzmann (also Shannon [41]) entropy is given via Here we show that, while the Boltzmann entropy E can both increase and decrease in time, its appropriate modification, known as the Kullback-Leibler entropy [4,5,27,29,35], is a nonincreasing function of time. Also, unlike the Boltzmann entropy, the Kullback-Leibler entropy P is invariant under arbitrary changes of variables, and is also nonnegative, that is, for two probability densities F 1 and F 2 , where the equality is achieved only when F 1 = F 2 (or the difference between the two has zero volume measure on V K × S K ). Let ψ : R → R be a differentiable function. Then, multiplying both sides of (3.46) by its derivative ψ with F used for the argument, we obtain (4.27) ∂ ∂t Integrating over dV K dS K , we obtain where we denote F = F(v , w ). Via (1.5), (1.6) and (1.7), we write, for any i and j, the term with F(v i , v j ) in the right-hand side as which leads to For the entropy in (4.24), we substitute , and arrive at the following equation for E : Next, recalling the inequality x ln x ≥ x − 1, observe that, for two probability densities F 1 and F 2 , we have (4.33) which, upon substitution into (4.32), yields the following inequality: For the part with F(v i , v j ), we observe that, for any i and j, we have which allows to coalesce the two Heaviside functions into 1 and results in The right-hand side of (4.36) above can, in general, be negative, and, therefore, it is possible for the Boltzmann entropy of the system to decrease in general. As an example, consider the uniform initial condition in both coordinates and velocities, which, under the normalization constraint, maximizes the entropy over all possible states. However, this state is not a steady state of (3.46); indeed, the states corresponding to the overlapped spheres will be "drained", the solution will become non-uniform in coordinates, and the entropy will decrease as a result.
To examine the Kullback-Leibler entropy in (4.25), let us look at the expression where GF K is the steady state from (4.19) and (4.22). We rearrange the terms above as Using (3.46) and the fact that GF K is a steady state, we write (4.39) and subsequently obtain Upon the integration over dV K dS K , the terms with ln(GF K ) in the right-hand side above disappear, and we arrive at (4.41) Adding (4.41) to (4.36) and changing the sign on both sides, we arrive at

The steady solution for a system with independently distributed initial states.
Recall that the factor G in (4.22) is determined by the distribution of the center of mass of the spheres, which, in turn, is defined entirely by the corresponding initial condition for the forward equation (3.46). So far, we did not elaborate on the choice of the initial conditions, and thus G was presumed to be largely arbitrary. However, from the perspective of physics, some of the states are more realistic, while others are notably less so. As an example of a physically unrealistic state, one can arrange the spheres initially so they fly in straight lines without interacting. However, in real-world systems, such as gases and liquids, these states are not encountered in practice, and instead a strongly chaotic motion of frequently colliding molecules is observed.
Here we propose the "most common" steady state for a system with large number of spheres, based on a "physicist's reasoning". Namely, in what follows, we assume that the spheres are initially distributed independently of each other. We compute the distribution of their center of mass, and, observing that it is invariant under the dynamics, conclude that the relevant steady state must have the same distribution of the center of mass. Note that the independence of distributions for each sphere suggests that the initial measures of overlapped states may be nonzero, however, in our random gas model, it is permissible for the spheres to overlap.
Let the distribution F 0 of the K spheres in the system at the initial moment of time be given via the product of independent identical distributions f 0 for each sphere: Let the position of the center of mass of the system be given via y: Let us express the first coordinate, x 1 , through the position of the center of mass y, and the rest of the coordinates: (4.45) Then, the distribution of y is given via where f x 0 denotes the x-marginal of f 0 : In what follows, we will assume that f x 0 is nondegenerate, that is, its support is of the same dimension as that of the whole x-space. Letf denote the characteristic function of f x 0 , Clearly,f (0) = 1, since f x 0 is a probability density. Additionally, due to the fact that f x 0 is nondegenerate, |f (k)| < 1 for k = 0 [42]. Next, let us express g(y) via the characteristic function of f x 0 from (4.48): Recalling that, in the generalized sense, we obtain Observing that |e −ik·yf (k)| < 1 for k = 0, we find that, in the limit K → ∞, the integrand approaches zero for all k = 0. Thus, as K → ∞, g(y) loses its dependence on y and becomes the uniform distribution. We thus conclude that the uniform distribution of the center of mass of the system of spheres is the most physical one. The corresponding most physical steady state is thus given via G = 1 in (4.22), as, underF K , any given state of the spheres is exactly as likely as its arbitrary parallel coordinate translation, and thus the center of mass ofF K alone is distributed uniformly.
Observe that the preceding derivation can be modified so that the spheres are distributed independently, but not identically (that is, any i-th sphere could be distributed independently, but with its own distribution f 0i ). If so, for the same reasoning to hold, we need all corresponding characteristic functionsf i of each x-marginal to be bounded above by the same h(k), that is, the spheres should initially be distributed independently and "similarly", which is a reasonable requirement from the physical standpoint.
4.5. The structure of marginal distributions of the physical steady state. As it becomes important below, we need to discuss the structure of the marginal distributions of the physical steady stateF K (x 1 , . . . , x K ) in (4.19), that is, the functions of the form where S n+1...K (v 1 , . . . , v n ) is the surface area of the subset of the constant energy sphere which corresponds to fixed velocities v 1 , . . . , v n , and dV n+1...K is the volume element of the subset of the volume which corresponds to fixed coordinates x 1 , . . . , x n . The integral over dV n+1...K in the right-hand side of (4.53) can be expressed in the form Above, we introduce the spatial correlation function R (n) K involves integration over all factors e −λΘ ασ (σ− x i −x j ) which are part of F K , but notF n . Observe that the following normalization conditions hold concurrently: (4.56a) S n F n (x 1 , . . . , x n )dV n = 1, K does not vary "too much" throughout the domain, its magnitude should generally be of order 1. In situations where the overlapped states take a very small part of the total volume (which implies that the volume of the domain greatly exceeds the total volume of the spheres in it, the so-called "dilute gas"), the integrand of (4.55) equals 1 almost everywhere, and, therefore, R (n) For the reasons which will become clear below, of particular importance are the special cases with n = 1 and n = 2. For n = 1, we claim (4.57)F (1) where V is the volume of the coordinate space of a single sphere, and the independence on x follows from the fact thatF K depends only on the distances between the spheres, and not on their absolute positions. For the special case n = 2, we claim (4.58) R K (x, y) = R K ( x − y ), which follows from the symmetry and isotropy considerations. Indeed, observe that, sinceF K depends only on distances between different spheres, R K can depend neither on the absolute positions of its two spheres, nor on the direction of their mutual orientation, which leaves only the distance as a possible dependence. Taking into account (4.19), for the two-sphere marginalF (2) K we, therefore, arrive at (4.59)F Combining (4.57) and (4.59), we can express (4.60)

The marginal distributions in the limit of infinitely many spheres.
In what follows, we are interested in situations where the number of spheres is large, and thus we need to to examine the limit of (4.60) as K → ∞. In such a limit, first we observe that This happens due to the fact that the diameter σ must decay as K −1/3 , in order for the total volume of the spheres to be bounded (to fit into the domain without overlapping). Therefore, the volume of the set of points where the integrand of Z 2 in (4.19) equals e −λ is proportional to K −1 . This means that the ratio above behaves as V 2 /Z 2 = 1 + O(K −1 ), and approaches 1 as K → ∞. Second, it can be shown via geometric arguments (see [1,2,35] and references therein) that, as K → ∞, where θ is the temperature of the system of the spheres, given via The relations in (4.62) lead to Combining (4.60), (4.61) and (4.64), we arrive, as K → ∞, to where we dropped the subscript K fromF (1) ,F (2) and R (2) to signify that these quantities are related to infinitely many spheres.

The forward equation for the marginal distribution of a single sphere
Let us consider an unsteady solution F of the forward equation in (3.46) which is invariant under the permutations of the spheres, and let us denote the two-sphere and one-sphere marginal distributions of F as where the integral in w occurs over the remaining area of S K for a fixed v. Then, integrating (3.46) over all coordinate-velocity pairs (x i , v i ) but one, we arrive at To arrive at the above result in the collision term, observe that, for fixed x and y, we, with the help of (1.5) and (1.7) have, for i, j > 1, i = j: which cancels out all terms in the collision which do not involve that particular sphere over which the integration does not occur. Next, we switch to the spherical coordinate system and replace where r is the nondimensional distance, n is the unit vector, dn is the area element of the unit sphere, "+" is used in the term with F (2) (v , w ), while "−" is used in the term with F (2) (v, w). In the new variables, the identity in (1.2) becomes (1.22), and (5.2) becomes the leading order identity of the corresponding BBGKY hierarchy [9,11,28] for (3.46): At this point, it might be tempting to assume that α is small enough so that one can neglect the variations in F (2) within the contact zone and integrate in dr, however, it is clearly not the case due to e.g. (4.65). Thus, in order to proceed further, we first need to find an appropriate closure to F (2) above in (5.5) in terms of the one-sphere marginal f .

5.1.
Approximating the two-sphere marginal via one-sphere marginals. Above in (5.5) we arrived at the ubiquitous closure problem of molecular dynamics -the two-sphere marginal F (2) must be approximated via the one-sphere marginal f in order to be able to transform the leading order BBGKY identity in (5.5) into a standalone closed equation for f . In the mean field approximation [44] for long-range potentials, the joint density of two particles is approximated via the product of single-particle distributions. However, in our situation, the relation in (4.65) for the steady state marginals does not permit such a direct factorization. Instead, here we choose to relate the unsteady marginals f and F (2) in the same way their steady counterparts are related in (4.65) in the limit of infinitely many spheres: where R (2) depends parametrically on α and λ: R (2) = R (2) λ,α . Observe that the closure in (5.6) becomes the exact solution if F (2) is the marginal distribution of the steady state as K → ∞. If F (2) is not a marginal distribution of the steady state, but sufficiently close to it, then we can estimate the improvement of (5.6) over (4.65) via the following appropriately generalized version of the marginal formula for the Kullback-Leibler divergence [3,36].
Let F(z 1 , . . . , z K ) be a probability density, ψ(z 1 , . . . , z K ) > 0 be an integrable function, and p 1 (z 1 ), . . . , p K (z K ) be a set of probability densities for z i ∈ R d , d > 0. Then, the following identity holds for the Kullback-Leibler divergence P: whereupon (5.7) follows from the first and second pair of terms in the integrand. The second term in the right-hand side of (5.7) is always nonnegative, which bounds the left-hand side below as: We can apply the formula above to our set-up by setting z 1 = (x, v), z 2 = (y, w), F = F (2) , ψ = e −λΘ ασ (σ− x−y ) R (2) ( x − y ), and p 1 = p 2 =F (1) from (4.65). One should note, however, that the right-hand side of (5.9) is not necessarily the Kullback-Leibler divergence due to the fact that (5.6) (which is the second argument of P) is not necessarily normalized to 1 for large σ. However, for a dilute gas (that is, when σ is much smaller than the average distance between the spheres), the second argument of P becomes normalized to 1, so that the right-hand side of the Kullback-Leibler bound in (5.9) indeed involves the probability densities in such a case. Additionally, the closure in (5.6) preserves the proportionality relations between the velocity moments in v or w variables separately. Now, taking into account that x − y = σr above in (5.5), we have and, with the approximation in (5.6), the forward equation for the marginals in (5.5) becomes

Thin contact zone and impenetrable spheres.
Above in (5.11), we can formally assume that the contact zone is "thin", that is, α → 0 so that, for the values of r for λ,0 (σ). In such a case, the integral over dr involves only the mollifier δ α (r − 1) and its antiderivative, and thus can be integrated across the thin contact zone exactly: This simplifies the forward equation in (5.11) to Now, observe that the factor (1 − e −λ ) in front of the equation above is the probability that at least one jump arrives in the point process during the collision. Sending the intensity of the point process λ → ∞ sets this probability to 1, which means that the spheres can be considered impenetrable. The resulting forward equation for a single impenetrable sphere is given via This is a variant of the Enskog equation for hard spheres [8,19,22,31,43], which we henceforth refer to as "the Enskog equation", even though it is not identical to what was originally proposed by Enskog. Apparently, the collision integral closure above in the Enskog equation (5.14) becomes exact if f is the single-sphere marginal distribution of the steady stateF K in (4.19), in the limit as the number of spheres K → ∞, the width of the mollifier α → 0, and λ → ∞.
Observe that if the gas is dilute (that is, the diameter σ of a sphere is much smaller than the average distance between the spheres), then we can neglect the higher order terms in σ in comparison to the leading order. This sets R , v), and the Enskog equation in (5.14) becomes the Boltzmann equation in (1.33).
6. The fluid dynamics of the Enskog equation in a physical hydrodynamic limit Above, we introduced a new, random model for hard sphere collision which avoids Loschmidt's objection [34], and where the form of the collision integral is a property of the forward Kolmogorov equation of the full multisphere dynamics. We have shown that, under suitable assumptions, this new model leads to the Boltzmann equation in (1.33). However, one of the key differences between our derivation of the Boltzmann equation and its conventional derivation [14,21] is that, as a precursor to the Boltzmann equation, we also obtain the Enskog equation (5.14) in a systematic fashion, by reducing and simplifying the forward Kolmogorov equation of the multisphere dynamics under suitable assumptions.
The hydrodynamic limit of the Boltzmann equation is a well-researched area (see, for example, [24] and references therein), and it is known that the equations for the velocity moments of the Boltzmann equation lead to the Euler [7] and Navier-Stokes [15,26] equations. However, the hydrodynamic limit of the Enskog equation is a less researched topic, although there are works on this subject as well (see [31] and references therein).
Below, we derive the fluid dynamics equations of the Enskog equation in (5.14) in the hydrodynamic limit (that is, as the diameter σ → 0) which, on one hand, is physically plausible, and, on the other hand, produces additional nonvanishing terms in the resulting fluid dynamics equations.
6.1. The mass-weighted equation and the hydrodynamic limit. Here we endow each sphere with a mass density ρ sp , which is a constant parameter. Above, f is the distribution density of a single sphere in the K-sphere system, and thus is normalized to 1. However, for the subsequent derivation of the fluid dynamics equations, it is more convenient to normalize the density by the total mass of the system. With this, we define the mass density g via where the factor in front of f above is the total mass of the system of K spheres, each of diameter σ and density ρ sp . The corresponding mass-weighted form of the Enskog equation for K spheres in (5.14) is, therefore, Here, observe that, as K → ∞, the ratio (K − 1)/K → 1, and we need to decide how the remaining coefficient in front of the integral behaves as σ → 0. Here, we argue that the density ρ sp of the spheres should be kept constant; in Figure 1, we plot the cubic root of the atomic mass versus the atomic diameter for noble gases, so that the straight line of the least squares fit approximates the hypothetical situation with constant density of the spheres ρ sp . Observe that, with the exception of neon, all atoms align on the least squares fit line, which indicates that the density of atoms of the actual noble gases tends to be constant. The resulting Enskog equation is given via where the factor in front of the collision integral above behaves as ∼ σ −1 in the constantdensity hydrodynamic limit K → ∞, σ → 0, ρ sp = const.
6.2. The Euler equations. With σ being treated as a small parameter, a standard way to look for a solution of (6.3) is to expand g in powers of σ [31]: and then recover the expansion terms sequentially starting from the leading order. We also note that R ∞,0 (0) = 1, and thus As a result, in the leading order of σ, we find which is the usual Boltzmann collision integral [14,15,24,26]. This means that the leading order expansion term g 0 is the standard Maxwell-Boltzmann thermodynamic equilibrium state, where the mass density ρ(t, x), average velocity u(t, x) and kinetic temperature θ(t, x) are given by the following mass-weighted moments of molecular velocity: Above, for further convenience, we additionally denoted the average kinetic energy via (6.9) For the next order in σ, we have Observe that the form of the equation above generally does not guarantee that g 0 retains its Gaussian form in (6.7), due to the additive collision integral in the right-hand side of (6.10). Thus, we will have to resort to a weaker formulation of (6.10) in the form of velocity moments of 1, v, and v 2 , in the same manner as in [26,31], which allows to exclude the terms with g 1 completely and close the equation with respect to ρ, u, and θ. First, for two arbitrary functions ψ(v) and h(v, w), the latter with the symmetry property h(v, w) = h(w, v), observe the following identities (assuming that the integrals are bounded): Above, with help of the symmetry of h(v, w), we renamed v as w and vice versa, and also changed the sign of n.
Second, observe the following chain of identities: Above, in the first identity we used (1.5); in the second identity we used (1.7) and also replaced n with −n, using the fact that v and w in (1.22) are invariant with respect to the sign of n; in the last identity we renamed v as v, w as w, and vice versa. Third, with help of (6.11) and (6.12), we can write For a constant ψ, the integral above is, obviously, zero regardless of the choice of h(v, w). Additionally, recall that v and w in (1.22) are chosen so that the momentum and energy of any two spheres are preserved during the collision. Thus, if ψ(v) is set to v or v 2 (or any linear combination thereof, including an additive constant), then the integral above in (6.13) is also zero irrespective of the choice of h(v, w).
Computing the velocity moments in (6.8) of both sides of the equation in (6.10), and setting h(v, w) = g 0 (v)g 1 (w) + g 1 (v)g 0 (w), we arrive at (6.14a) ∂ρ ∂t where by C we denote the corresponding collision integrals Next, recalling from (1.22) that we further obtain With g 0 provided explicitly in (6.7), the expressions above are integrated exactly into Substituting (6.18) into (6.14), we arrive at (6.19a) ∂ρ ∂t The equations above were also obtained in [31], where they were referred to as the "Enskog-Euler" equations. With ρ/ρ sp 1 (dilute gas), the equations in (6.19) become the conventional Euler equations [7,24]. However, observe that in the physical constantdensity limit, considered here, the additional terms are formally nonvanishing despite the fact that the diameter of the sphere σ → 0.

The Newton and Fourier laws.
Here we examine the next-order term g 1 from (6.4). With g 0 already computed above in (6.7), we can use (6.10) to determine the properties of g 1 . However, observe that, while the form of g 0 is given uniquely by the leading order identity (6.6), such a strict constraint on g 1 is not present -indeed, the next-order equation in (6.10) does not call for any particular restriction on g 1 . Therefore, the form of g 1 in the expansion (6.4) can be chosen as necessary.
Below, we are going to follow [26] and assume that the full density g can be expressed in the form of the Hermite polynomial expansion around the Gaussian density g 0 . This form is chosen solely due to convenience -as shown below, all subsequent moment integrals that arise are expressed in terms of elementary functions. However, this is not a unique choice for g; for an alternative, see, for example, [33], where the form of density was chosen so as to maximize the Boltzmann entropy under given moment constraints. The downside of the latter approach is that the moment integrals of the solutions of constrained maximum entropy optimization problems cannot, in general, be expressed in terms of elementary functions.
Substituting the identities in (6.28) into (6.27), we express the stress S and heat flux q in terms of ρ, u and θ: where µ is the usual viscosity for the hard sphere gas [26], and [ρS] B , [ρq] B denote the conventional Newton and Fourier law expressions for a dilute gas (ρ/ρ sp 1), respectively, derived purely from the Boltzmann equation [15,26]. Observe that the viscosity µ is O(σ). Therefore, in order to include viscous effects into the fluid dynamics equations, we will have to retain the O(σ)-term in the expansion (6.4).
6.4. The Navier-Stokes equations. Above, we computed the expansion terms g 0 and g 1 of (6.4). The g 0 is given via the Maxwell-Boltzmann state (6.7), while the first-order correction g 1 is given via the Hermite polynomial (6.22), with the parameters S and q provided via the Newton and Fourier laws in (6.29). Now, instead of proceeding with the formal hydrodynamic limit as the sphere diameter σ → 0, we will assume that σ is a finitely small constant. In such a case, we can truncate (6.4) to the first two leading order terms (assuming that all O(σ 2 ) terms are small enough so that they can be neglected). The result is known as the Grad state [26]: (6.30) g Grad = g 0 + σg 1 .
Then, we follow the same procedure for g Grad as above for g 0 ; namely, we substitute g Grad into (6.10), and compute the transport equations for the mass, momentum and energy moments. Since g Grad contains O(σ) terms, the terms of the same order must be retained in the collision term of (6.10). The resulting transport equation for the density is unchanged, however, the equations for the velocity and energy are now given via (6.31a) ∂(ρu) ∂t πρ sp (n · (w + v))(n · (w − v)) 2 Θ n · (w − v) −g 1 (v)n · ∂g 0 (w) ∂x − g 0 (v)n · ∂g 1 (w) ∂x + 1 2 g 0 (v)n T ∂ 2 g 0 (w) ∂x 2 n dndwdv.
Above, for convenience, the constant R * denotes the derivative of R (2) ∞,0 at zero: Due to the fact that both g 0 and g 1 are given explicitly via (6.7) and (6.22), respectively, all integrals in the right-hand side of (6.31) are computable in terms of elementary functions.
Above in (6.37), observe that the presence of viscous terms is entirely due to the fact that we treat the sphere diameter σ as being finitely small; indeed, if we formally send σ → 0 as in the constant-density hydrodynamic limit above, the viscous terms in (6.37) will disappear, and the Enskog-Euler equations (6.19) will emerge instead.

Summary
In the current work, we start by examining the conventional derivation of the Boltzmann equation [10,24,26] from the Liouville equation for the hard sphere dynamics [13,14,21]. The conventional derivation consists, roughly, of two parts -the first part is the derivation of the BBGKY hierarchy [9,11,28] from the Liouville equation under the assumption that the solution of the Liouville equation is known, while the second part is the derivation of the Boltzmann equation from the leading order BBGKY identity by, first, altering the integrand under the collision integral (Boltzmann hierarchy), and second, factorizing the joint probability distribution of two spheres into the product of the single-sphere distributions. We observe that the second part of the derivation contradicts the assumptions under which the first part of the derivation was carried out. We also note that some of the observed contradictions are similar to the objection pointed out by Loschmidt [34].
At the same time, it is known from observations and experiments that the Boltzmann equation is de facto an excellent model for a dilute gas in practical scenarios. We, therefore, propose that there should exist a different underlying dynamical model of the hard sphere gas, which does not violate Loschmidt's objection and, at the same time, leads to the Boltzmann equation in the dilute gas approximation. We subsequently formulate a random process with an infinitesimal generator, which triggers the necessary velocity jumps via a suitable point process [18,38] whenever the collision condition is detected. This process is a Lévy-type Feller process [6,20], and the general form of its characteristic function is given by Courrège's theorem [17]. We subsequently formulate the forward Kolmogorov equation [23] for the probability density of the random dynamics, compute those of its steady states which are uniform in the velocity variables, and show that the corresponding Kullback-Leibler entropy [29] of the complete system of spheres is a nonincreasing function of time.
We find that, in the case of many spheres, which are distributed independently and identically at the initial time, the corresponding steady state is uniform not only in velocities, but also in the coordinates of the spheres, except for the "contact zones" (that is, the sets of coordinates which satisfy the collision condition). With help of the computed marginal distributions of these steady states, we derive the forward equation for the dynamics of a single sphere under the assumption that the distribution of the full system is invariant under the reordering of the spheres.
As the total number of spheres in the system becomes large, we find that, for the limiting dynamics of thin contact zones and impenetrable spheres, the forward equation becomes a variant of the Enskog equation [8,19,22,31,43]. Further, assuming that the gas is dilute, we arrive at the Boltzmann equation [10,[12][13][14]24,26]. Finally, we find that, in the hydrodynamic limit of constant-density spheres, the resulting Enskog-Euler and Enskog-Navier-Stokes equations acquire additional effects in both the advective and viscous terms, which are absent in the conventional Euler and Navier-Stokes equations derived from the Boltzmann equation.