Relative entropy of freely cooling granular gases. A molecular dynamics study

Whereas the original Boltzmann's $H$-theorem applies to elastic collisions, its rigorous generalization to the inelastic case is still lacking. Nonetheless, it has been conjectured in the literature that the relative entropy of the velocity distribution function with respect to the homogeneous cooling state (HCS) represents an adequate nonequilibrium entropy-like functional for an isolated freely cooling granular gas. In this work, we present molecular dynamics results reinforcing this conjecture and rejecting the choice of the Maxwellian over the HCS as a reference distribution. These results are qualitatively predicted by a simplified theoretical toy model. Additionally, a Maxwell-demon-like velocity-inversion simulation experiment highlights the microscopic irreversibility of the granular gas dynamics, monitored by the relative entropy, where a short ``anti-kinetic'' transient regime appears for nearly elastic collisions only.


Introduction
Granular gases are modeled in their simplest form as inelastic hard spheres with a constant coefficient of restitution, α [1][2][3][4]. It is well known that granular gases are intrinsically out of equilibrium and that a description by means of kinetic theory is meaningful. In a kinetictheoretical description of a granular gas, one defines the granular temperature as the mean kinetic energy per particle, as an analogue to its definition for molecular gases. Even though this temperature is not a thermodynamic temperature, one can look for the nonequilibrium entropy-like functional of this system, i.e., a Lyapunov functional, in analogy with Boltzmann's H-functional and the celebrated H-theorem for elastic collisions The problem introduced and solved by Boltzmann in 1872 [5] is not easy to extend in the context of granular gases and the associated inelastic form of the renowned Boltzmann equation. The original H-functional is precisely the Shannon measure [6] of the one-particle velocity distribution function [7,8]. However, it is known that, in its continuous description, Shannon's entropy presents the so-called measure problem [9], i.e., it does not weigh properly the phase space. In the elastic case, this problem is easily solved by considering the relative entropy (or Kullback-Leibler divergence [10]) of the one-particle velocity distribution function with respect to the Maxwellian distribution, which becomes the original H-functional up to a constant in that case. Moreover, some relevant properties of the elastic-particle system, like collisional symmetry and reversibility, do not hold anymore in the inelas- * e-mail: albertom@unex.es * * e-mail: andres@unex.es tic scheme. Then, the proper entropy-like functional must solve these issues.
The quest of such a quantity in the homogeneous case has been addressed mathematically in Refs. [11][12][13] in the context of the inelastic Boltzmann equation, and in Ref. [14] from a stochastic point of view. Both approaches converge into a single functional, which is proved in the quasielastic limit, i.e., 1 − α ≪ 1, to be the entropy-like functional associated with this system. In the case of free cooling, the conjectured quantity is the relative entropy of the reduced velocity distribution function, φ, with respect to the homogeneous cooling state (HCS), φ H , chosen as the proper reference distribution. This conjecture was recently reinforced with computer simulations in the whole range of inelasticity [15].
In this work, we complement the study carried out in Ref. [15] with new simulations. First, we study the problem by means of a simplified toy model [15] and investigate how it highlights the possible Lyapunov character of the proposed functional for two different reference distributions, namely the Maxwellian and the HCS distributions. Next, molecular dynamics (MD) results are presented and compared with the predicted theoretical behavior, including three systems not considered in Ref. [15]. Finally, as a fully original contribution of this work, we report MD results for a sort of Maxwell-demon experiment where the irreversibility of the collisional process and the possibility of an anti-kinetic stage are discussed.

A toy model
Let us consider a granular-gas model of inelastic and smooth hard spheres with collisional rules for the precollisional velocities, where v 12 is the relative velocity, and σ the unit intercenter vector. We will assume that the system satisfies the inelastic Boltzmann equation, which in reduced units reads Here, κ ≡ 2 √ 2π is a constant, c = v/v th is the reduced velocity, v th (t) = √ 2T (t)/m is the thermal velocity, m is the mass of a particle, T = m v 2 /3 is the granular temperature, which decreases monotonically following Haff's law [3,4,16,17] is the collisional operator in reduced units, and The relative entropy, or Kullback-Leibler divergence, of a velocity distribution, φ, with respect to a reference distribution, φ R , is defined as This functional is convex, non-negative, and identically zero if and only if φ = φ R [10]. We assume that both φ and φ R are isotropic and can be expanded around the Maxwellian φ M (c) = π −3/2 e −c 2 in terms of Sonine polynomials, where S k is the k-th Sonine polynomial and a k is the 2k-th cumulant of the distribution, defined as a k = S k /N k with N k = (2k + 1)!!/2 k k!. By definition, a 0 = 1 and a 1 = 0, so that the first nontrivial coefficient is the fourth cumulant a 2 = 4 15 c 2 − 1. Let us now construct a toy model of D KL (φ φ R ) [15] for an arbitrary reference distribution. Imagine a perturbative parameter ε in front of the Sonine summation in Eq. (4). Expanding in powers of ε and keeping terms up to second order, we have where ∆ R a k (s) ≡ a k (s) − a R k , a R k being the Sonine coefficients for the reference distribution function φ R (c). Inserting this expression into Eq. (3) and using the orthogonality condition of the Sonine polynomials, one obtains The r.h.s of Eq. (6a) is non-negative, and it is zero if and only if a k = a R k ∀k ≥ 2, that is, φ = φ R , in accordance with the properties of the relative entropy. Next, we neglect terms of O(ε 3 ), formally take ε = 1, and, as usually done in the literature [3,4,[18][19][20][21][22][23][24], discard terms with k ≥ 3. The result is Here, in consistency with the neglect of a k (s) for k ≥ 3, we have used the evolution equation ∂ s a 2 (s) = −K[1 + a 2 (s)]∆ H a 2 (s) for the fourth cumulant [15], where ∆ H a 2 (s) ≡ a 2 (s)−a H 2 and K is a positive constant. Whereas the r.h.s of Eq. (7a) is non-negative, the sign of the r.h.s of Eq. (7b) is determined by the relative signs of ∆ R a 2 (s) and ∆ H a 2 (s).
Let us consider two different reference distributions: the Maxwellian velocity distribution function, φ M , and the HCS velocity distribution function, φ H . In the first case (R = M), one has ∆ M a 2 (s) = a 2 (s), and, therefore, Thus, our toy model shows that a monotonic relaxation of D KL (φ φ M ) is not guaranteed. Let us assume, for instance, that the initial value a 2 (0) is negative and α 0.71, so that the steadyvalue a H 2 is positive [15,19,[21][22][23][24]; due to Bolzano's theorem, during its evolution a 2 (s) must cross the zero value,

Molecular dynamics simulations
In order to check the predictions of the toy model for the two considered reference distributions, we have performed MD simulations using the DynamO software [25] for this model of granular gases. It is well known that the free cooling of granular gases presents long-wavelength instabilities [4,26]. In order to avoid them, we have simulated systems formed by N = 1.35 × 10 4 particles in a cubic box of side length L/σ = 407.16, which is at least 30 times smaller than the critical length for the development of instabilities, which are not observed.
Initially, all particles are arranged in an ordered crystalized configuration from which the system melts. The initial velocities are oriented along randomized directions with either a common magnitude (initial distribution δ) or with a magnitude drawn from a Gamma (Γ) distribution. The respective initial values of the fourth cumulant are a 2 (0) = −0.4 (δ distribution) and a 2 (0) = 0.4 (Γ distribution). Thus, according to the toy model, a nonmonotonic relaxation of D KL (φ φ M ) is expected for the initial distribution δ if α 0.71 and for the initial distribution Γ if α 0.71.
In Fig. 1 one can observe that, as predicted by the toy model, a local minimum is actually observed during the evolution of D KL (φ φ M ) for α = 0.1 and 0.4 when starting from the initial condition δ, and for α = 0.87 when starting from the initial condition Γ. In the other three cases, however, the evolution of D KL (φ φ M ) is monotonic. In contrast, the relative entropy D KL (φ φ H ) decays monotonically for the six cases, in qualitative agreement with the toy model.

Velocity-inversion experiment
A discussion about entropy is not complete if the issue of irreversibility is not included. In the case of elastic hard disks, a simulated velocity-inversion experiment (produced by a sort of Maxwell's demon) was proposed more than forty years ago [27][28][29], where schemes with "anti-kinetic" parts in the evolution were tested [30] and Loschmidt's paradox was discussed. In Orban and Bellemans' pioneering works [27,28], during the evolution toward equilibrium the velocities of all elastic disks (simulated by MD) were inverted at a given waiting time t w and Boltzmann's H-functional was analyzed and seen to revert its decay by retracing its past values (anti-kinetic stage), in agreement with the underlying reversibility of the equations of motion. However, the H-functional resumed its decay after time t = 2t w and, moreover, due to unavoidable error propagation [31], the initial value of H was not exactly recovered if the velocity inversion took place after a sufficiently long waiting time. In a study involving irreversible particle dynamics, Aharony [29] observed that the anti-kinetic stage was not symmetric, the system rapidly forgetting the correlations it had at t w , and thereafter continuing to approach equilibrium.
In this section we revisit the velocity-inversion experiment in a freely cooling granular gas, modeled as inelastic hard spheres, where the collisional rules are given by (1). In this system, collisional symmetry is broken down by the inelasticity of collisions, closely related to a violation of microscopic reversibility. Consider two colliding particles with precollision velocities {v 1 , v 2 } and a relative orientation characterized by the unit vector σ (with v 12 · σ > 0). In that case, the postcollisional velocities are Now, we invert the velocities {v ′ 1 , v ′ 2 } and obtain the subsequent postcollision velocities, Therefore, where I is the inversion-velocity operator) unless α = 1. We studied this effect from MD simulations in a computer experiment similar to those of the works discussed above [27][28][29]. A waiting time s w = 0.5 was chosen, several values of α were considered, and the evolution was monitored by D KL (φ φ H ), which plays the role of H in the elastic case.  Figure 2 shows the time evolution of D KL (φ φ H ) when starting from the δ initial condition and then applying the velocity inversion. The coefficients of restitution considered are α = 0.1, 0.4, 1/ √ 2, 0.87, 0.99, and 1. In the elastic case (α = 1), one recovers the results of Ref. [27], i.e., the system almost reaches the original configuration at s = 1 but afterwards it evolves toward equilibrium again. Whereas one expects that inelastic collisions erase completely the possibility of a reversible period, in the quasielastic case α = 0.99, although it is short, an antikinetic transient stage exists after the velocity inversion; this effect is translated into a small growth of D KL (φ φ H ). Of course, the duration of the anti-kinetic regime becomes longer as α comes closer to 1. On the other hand, as inelasticity increases (α ≤ 0.87), the influence of the velocity inversion is noticeable by a change of curvature only, and this short effect shrinks with increasing inelasticity, as expected. The effect of inelasticity on the microscopic irreversiblity reflected by the behavior of D KL (φ φ H ) is analogous to that observed by Aharony [29] for the conventional H-functional in the evolution toward equilibrium.

Concluding remarks
In this paper we have provided further evidence from MD simulations on the conjecture that the Kullback-Leibler divergence D KL (φ φ H ) is a possible entropy-like functional for the case of isolated freely cooling granular gases [14,15], even for strongly inelastic systems. Furthermore, this conjecture is supported by a simple toy model, which, on the other hand, predicts a nonmonotonic behavior of D KL (φ φ M ) if a 2 (0) and a H 2 have opposite signs. This theoretical expectation has been nicely confirmed by our simulations.
Finally, the classical velocity-inversion experiment [27][28][29][30], originally devised for systems relaxing to equilibrium, has been applied on granular gases relaxing to the HCS and monitored via D KL (φ φ H ). While, as expected, the initial configuration is almost perfectly recovered if the collisions are elastic (α = 1), microscopic reversibility is frustrated by inelasticity, no matter how small. In fact, a (short) anti-kinetic stage, where ∂ s D KL (φ φ H ) > 0, is only possible in the quasielastic regime (e.g., α = 0.99) and disappears for sufficiently high inelasticity (α 0.9).