Kullback–Leibler Divergence of a Freely Cooling Granular Gas

Finding the proper entropy-like Lyapunov functional associated with the inelastic Boltzmann equation for an isolated freely cooling granular gas is a still unsolved challenge. The original H-theorem hypotheses do not fit here and the H-functional presents some additional measure problems that are solved by the Kullback–Leibler divergence (KLD) of a reference velocity distribution function from the actual distribution. The right choice of the reference distribution in the KLD is crucial for the latter to qualify or not as a Lyapunov functional, the asymptotic “homogeneous cooling state” (HCS) distribution being a potential candidate. Due to the lack of a formal proof far from the quasielastic limit, the aim of this work is to support this conjecture aided by molecular dynamics simulations of inelastic hard disks and spheres in a wide range of values for the coefficient of restitution (α) and for different initial conditions. Our results reject the Maxwellian distribution as a possible reference, whereas they reinforce the HCS one. Moreover, the KLD is used to measure the amount of information lost on using the former rather than the latter, revealing a non-monotonic dependence with α.


Introduction
Thermodynamics and information theory are clearly connected via the entropy concept. This idea allows physicists to understand plenty of details and consequences in the evolution and intrinsic behavior of physical systems. However, finding the entropy-like Lyapunov functional for a given problem is not an easy task. Thankfully, information theory provides tools that one can use in physics problems, usually proving a rewarding feedback. Along this paper, and as usually done in the context of information theory [1,2] and nonequilibrium statistical mechanics [3,4], we borrow from equilibrium statistical mechanics and thermodynamics the use of the term "entropy" in a broader sense. The same applies to the term "temperature," introduced in Equation (3) below.
In this work, we address the quest of finding the Lyapunov functional of an isolated freely cooling monodisperse granular gas, modeled by identical inelastic and smooth hard disks (d = 2) or hard spheres (d = 3) with constant coefficient of restitution (α). The interest of this study does not only reside in the mathematical challenge, but also in the physical consequences for granular matter. Typically, for a classical gas, Boltzmann's H-theorem provides the desired entropy-like Lyapunov functional [5,6]. Nevertheless, inelasticity plays a fundamental role in the dynamics, and the hypotheses of the latter theorem are not applicable. Previous works have proposed the Kullback-Leibler divergence (KLD) [7,8] as the proper alternative to the H-functional [9][10][11][12]. One of the aims of this paper is to explore with molecular dynamics (MD) simulations [13] the validity of the KLD as a Lyapunov functional in the whole range of definition of α and for both disks and spheres.
The freely cooling one-particle velocity distribution function (VDF) of our granular-gas model is expected to asymptotically reach a scaled form, the so-called "homogenous cooling state" (HCS), f HCS . Although its explicit form is unknown, there is a vast amount of literature about it [14][15][16][17][18][19][20][21][22] and recent experiments have demonstrated some of their properties [23]. While computational and experimental evidence supporting the HCS are overwhelming, a rigorous mathematical proof on its existence and long-time approach has only been achieved for inelastic Maxwell models described by the Boltzmann equation [24][25][26][27][28].
The HCS VDF f HCS is usually expressed as an infinite expansion around the Maxwellian VDF in terms of Sonine polynomials [5,14,15], even though the expansion may break down for large inelasticities [29,30]. Here, in order to provide a detailed description of the problem for both the stationary and transient regimes, we revisit some well-known results and also provide new simulation data and theoretical expressions obtained from a truncation in the Sonine expansion up to the sixth cumulant. In particular, our MD simulation results for the HCS fourth and sixth cumulants are compared with previous "direct simulation Monte Carlo" (DSMC) results [18,19] and a good agreement is found.
The paper is structured as follows. In Section 2, the Sonine expansion formalism is presented and simulation and theoretical results for the fourth and sixth cumulants are provided. The measure problem introduced by the original H-functional is established in Section 3 and the KLD for two different reference VDFs is studied and compared with MD simulation outcomes. Finally, in Section 4, some concluding remarks of this work are presented and discussed.

Boltzmann Equation and HCS
Consider a model of a monodisperse granular gas consisting of an isolated collection of inelastic hard d-spheres of mass m, diameter σ, and a constant coefficient of normal restitution α < 1. Under the molecular chaos ansatz (Stosszahlansatz), the free cooling of a homogeneous and isotropic gas can be described by the Boltzmann equation [14] ∂ t f (v 1 (1) where n is the number density, v 12 = v 1 − v 2 is the relative velocity of the two colliding particles, σ is a unit vector along the line of centers from particle 1 to particle 2, the subscript + in the integral over σ means the constraint v 12 · σ > 0, and are precollisional velocities. Note that we have defined the VDF with the normalization condition dv f (v; t) = 1. An important quantity is the granular temperature defined as Taking moments in Equation (1), one finds the cooling equation where the cooling rate is given by Let us introduce the thermal velocity v th (t) ≡ 2T(t)/m, which allows us to define the rescaled VDF φ(c; s) as where the variable s in φ(c; s) is a scaled time defined by Here, ν is the (nominal) collision frequency, so that s(t) represents the (nominal) accumulated average number of collisions per particle up to time t. In terms of these dimensionless quantities, the Boltzmann Equation (1) can be rewritten as where we have taken into account that Note that F 0 = F 2 = 0, since µ 0 = 0 and c 2 = d 2 . In the long-time limit, the free cooling is expected to reach an asymptotic regime (the HCS) in which the scaled VCF is stationary , i.e., φ(c; Henceforth, a subscript or superscript H on a quantity means that the quantity is evaluated in the HCS. Within that regime, Equation (5a) shows that ζ H (t)/ T H (t) = const, so that the solution to Equation (4) gives rise to the well-known cooling Haff's law [14,15,31] T H (t) = T H (t 0 ) where t 0 being an arbitrary time belonging to the HCS regime. Also in the HCS regime, µ 2 (s) → µ H 2 = const and thus Equation (4) becomes ∂ s T H (s) = −(4/κd)µ H 2 T H (s), whose solution is Therefore, in the HCS, the temperature decays exponentially with the average number of collisions per particle.

Sonine Expansion Formalism
The Maxwell-Boltzmann VDF φ M (c) = π −d/2 e −c 2 is not a solution of the HCS Boltzmann Equation (10). While its analytic form has not been found, the HCS solution is known to be rather close to φ M in the domain of thermal velocities (c∼1) [20]. Thus, it is convenient to represent the time-dependent VDF in terms of a Sonine polynomial expansion, where where are Sonine (or generalized Laguerre) polynomials, which satisfy the orthogonalization condition In Equation (13), the Sonine coefficient a k (s) is the 2k-th cumulant of the VDF at time s. According to Equation (15), In particular, a 0 (s) = 1, a 1 (s) = 0, and

Truncated Sonine Approximation
Thus far, all the results presented in Sections 2.1 and 2.2 are formally exact within the framework of the homogeneous Boltzmann Equation (1). However, in order to obtain explicit results, we need to resort to approximations.
As usual [15][16][17][18][19]29,32], we will start by neglecting the coefficients a k with k ≥ 4 in Equation (13), as well as the nonlinear terms a 2 2 , a 2 a 3 , and a 2 3 in the bilinear collision operator I[c|φ, φ]. Given a functional X[φ] of the scaled VDF φ(c), we will use the notation L 3 {X} to denote the result of that truncation and linearization procedure. Furthermore, if a 3 is also neglected, the corresponding approximation will be denoted by L 2 {X}. In particular, in the case of the collisional moments µ 2 , µ 4 , and µ 6 , one has (18) where the expressions for the coefficients A i , B i , and C i as functions of α and d can be found in Ref. [29] and in Appendix A of Ref. [19]. Obviously, L 2 {µ 2 }, L 2 {µ 4 }, and L 2 {µ 6 } are obtained by formally setting A 3 → 0, B 3 → 0, and C 3 → 0, respectively.
Let us first use the simple approximation L 2 to estimate a H 2 . From Equation (9), we have that F H 4 = 0. Thus, the obvious approximation [17] consists of where, in the last steps, use has been made of the explicit expressions of A 0 , A 2 , B 0 , and B 2 . However, this is not by any means the only possibility of estimating a H 2 [18,19,33]. In particular, one can start from the logarithmic time derivative of the fourth moment and then take Note that a H,a Both approximations (a H,a 2 and a H,b 2 ) are practically indistinguishable in the region 0.6 α < 1, but a H,b 2 is much more accurate than a H,a 2 for higher inelasticity [18,19]. Next, to estimate a H 3 , we start from the exact condition F H 6 = 0 and carry out either the linearization (22) or, alternatively, In Equations (22) and (23), a H 3 is expressed in terms of a H 2 . Using Equations (19) and (20), four possibilities in principle arise, namely have been obtained, we turn our attention to the evolution equations of a 2 (s) and a 3 (s). Approximating Equation (9) Its solution is Analogously, if Equation (9) with k = 6 is approximated as κ 2 ∂ s c 6 = L 3 {F 6 (s)}, the resulting evolution equation for a 3 is where Taking into account Equation (26), the solution to Equation (27) is where 2 F 1 (a, b; c; z) is the hypergeometric function [34]. As far as we know, Equations (26) and (29) had not been obtained before.

Comparison with MD Simulations
The approximate theoretical predictions for a H 2 and a H 3 were tested against results obtained from the DSMC simulation method in, for instance, refs. [18][19][20]. However, since the DSMC method is a stochastic scheme to numerically solve the Boltzmann equation [35], it does not prejudice by construction the hypotheses upon which the Boltzmann equation is derived, in particular the molecular chaos ansatz. Therefore, it seems important to validate the Sonine approximations for a H 2 and a H 3 by event-driven MD simulations as well. In addition, the theory allows us to solve the initial-value problem and predict the evolution of the fourth and sixth cumulants, as shown by Equations (26) and (29), and an assessment of those solutions is in order. In our MD simulations, we studied systems with densities nσ d ≈ 5 × 10 −4 and 2 × 10 −4 for disks and spheres, respectively. It is known that the HCS exhibits a shearing/clustering instability for sufficiently large systems [14,36]. To prevent this, the side length of the simulation box was chosen as L/σ ≈ 5 × 10 3 for disks and L/σ ≈ 4 × 10 2 for spheres (see Appendix A for technical details). These values are about 2 and 30 times smaller, respectively, than the critical values beyond which the HCS becomes unstable in the less favorable case considered (α = 0.1). Moreover, we have expressly verified that the systems remain stably homogeneous even for long times. Figures 1a,b show the α-dependence of a H 2 and a H 3 , respectively, for both hard disks (d = 2) and spheres (d = 3). An excellent agreement between the MD and DSMC simulation results for the whole range of α is observed. This means that the molecular chaos ansatz does not limit the applicability of the Boltzmann description, even for large inelasticities [14], at least for dilute granular gases. As for the approximate theoretical predictions, it is quite apparent that a H,b 2 (see Equation (20)) performs very well, even if the fourth cumulant is not small (e.g., a H 2 ∼0.2 at α = 0.1). The approximate sixth cumulant a H,ab 3 (see Equations (22) and (24)) is less accurate at a quantitative level, especially in the case of disks, but captures quite well the general influence of inelasticity. While a H 2 changes from negative to positive values at α 1/ √ 2 0.71, a H 3 is always negative. Note that, for large inelasticity, the cumulants a H 2 and a H 3 are comparable in magnitude. Given that the Sonine expansion (13) is only asymptotic [15,30], it is remarkable that a theoretical approach based on the assumptions |a H 3 | |a H 2 | 1 does such a good job for high inelasticity as observed in Figure 1.
Next, we study the evolution from a non-HCS state, as monitored by a 2 (s) and a 3 (s). We have chosen an initial state very far from the HCS: the particles are arranged in an ordered crystalized configuration and all have a common speed √ d/2v th (0) along uniformly randomized directions. Therefore, at s = 0, c k = (d/2) k/2 , so that a 2 (0) = − 2 d+2 and a 3 (0) = −  (26) and (29), respectively. Four representative values of the coefficient of restitution have been considered, namely α = 0.1 (very high inelasticity), 0.4 (high inelasticity), 0.87 (moderately small inelasticity), and 1 (elastic collisions); α = 0.87 has been included because it is practically at this value where a H 2 presents a local minimum, both for disks and spheres [see Figure 1a]. Note that, in the case of simulations, the quantity s represents the actual average number of collisions per particle and, consequently, is not strictly defined by Equation (7), in contrast to the case of theory. From Figure 2 we observe that, despite the large magnitude of the initial fourth cumulant (a 2 (0) = − 1 2 and − 2 5 for d = 2 and 3, respectively), the simple relaxation law (26) describes very well the full evolution of the cumulant. Discrepancies with the simulation results are visible only in the region (2 s 4) where the curves turn to their stationary values, especially in the case of disks. In what concerns the sixth cumulant, which also has a large initial magnitude (a 3 (0) = − 2 3 and − 16 35 for d = 2 and 3, respectively), the theoretical expression (29) is able to capture, at least, the main qualitative features, including the change from a non-monotonic (α = 0.1 and 0.4) to a monotonic (α = 0.87 and 1) evolution. Again, the agreement is better for spheres than for disks. Note also that the evolution curves for α = 0.87 and 1 are hardly distinguishable from each other. (see Equation (20)) and a H,ab 3 (see Equations (22) and (24)). The insets magnify the region 0.6 ≤ α ≤ 1. The error bars in the simulation data are smaller than the size of the symbols. In Figures 2 and 3, the initial values a 2 (0) and a 3 (0) are common to all the coefficients of restitution considered. In order to have a more complete picture, let us now fix the most inelastic systems (α = 0.1) and take five different initial conditions. The HCS values of the fourth and sixth cumulants at .150, −0.077} for d = 2 and d = 3, respectively. Thus, we have chosen the same initial distribution (hereafter labeled as δ) as in Figures 2 and 3 as a representative example of a 2 (0) < 0, the Maxwellian distribution (labeled as M) with a 2 (0) = 0, another one (labeled as I) with 0 < a 2 (0) < a H 2 , and two more (labeled as Γ and S) with a 2 (0) > a H 2 . The details of those five distributions can be found in Appendix B and the corresponding values of a 2 (0) and a 3 (0) are shown in Table A1. In the case of a 2 (s), Figure 4 shows again an excellent agreement between theory and simulation, except for the initial condition Γ and near the turning point already observed in Figure 2 for the initial condition δ. In what concerns a 3 (s), one can observe from Figure 5 that the performance of the approximation (29) is generally fair, especially for the initial conditions M and I. The limitations of Equation (26) for the initial condition Γ and of Equation (29) for the initial conditions Γ, S, and δ are due to the role played by higher-order cumulants in those cases.

KLD as a Lyapunov Functional
In this section, we restrict ourselves to spatially homogeneous states.

Boltzmann's H-Functional
The introduction of the H-theorem by Ludwig Boltzmann [37] was a revolution in physics and became an inspiration for new mathematical and physical concepts. This theorem is a direct consequence of the Boltzmann kinetic equation for classical rarefied gases, derived under its molecular chaos assumption [5,6]. Beneath this hypothesis for a classical gas which evolves via elastic collisions, the H-functional defined as is proved to be a non-increasing quantity; in other words, S = −H, up to a constant, is a non-decreasing entropy-like Lyapunov functional for the assumed gaseous system. After almost a century, once Information Theory was developed, Boltzmann's H-functional was interpreted as Shannon's measure [1] for the one-particle VDF of a rarefied gas. Nonetheless, the model considered in this paper for a rarefied monocomponent granular gas (inelastic and smooth hard d-spheres with a constant coefficient of restitution) violates Boltzmann's hypothesis of elastic collisions. In fact, a key role in the demonstration of the H-theorem for elastic collisions is played by the condition of collisional symmetry [37]. Consider two colliding particles with precollision velocities {v 1 , v 2 } and a relative orientation characterized by the unit vector − σ (with v 12 · σ < 0). After collision, the velocities are, in agreement with Equation (2), given by Next, suppose 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, Comparison with Equation (2) shows that Thus, C σ C − σ {v 1 , v 2 } = {v 1 , v 2 } unless α = 1 and, therefore, the H-functional, as defined by Equation (30), is not ensured to be non-increasing anymore if α < 1. Furthermore, Boltzmann's H-functional for the model of inelastic particles presents the so-called measure problem [38]. Shannon's measure is invariant under unitary transformations, but not for rescaling. In fact, under the transformation (6), From Haff's law, Equation (12), it turns out that (in the HCS) H * H is stationary but H H (s) grows linearly with the average number of collisions s. Then, one could naively think that a possible candidate to the Lyapunov functional would be H * (s), but the latter is still non-invariant under a change of variables where J ≡ |∂ c/∂c| is the Jacobian of the invertible transformation c = w(c). As will be seen below, whereas Shannon's measure presents a problematic weighting of the phase space, the KLD solves this non-invariance issue.

KLD
In general, given two distribution functions f (x) and g(x), one defines the KLD from g to f (or relative entropy of f with respect to g) as [7,8], as where x is a random vector variable defined on the set X. The quantity D KL ( f g) is convex and non-negative, being identically zero if and only if f = g. While it is not a distance or metric function (it does not obey either symmetry or triangle inequality properties), D KL ( f g) somehow measures how much a reference distribution g diverges from the actual distribution f or, equivalently, the amount of information lost when g is used to approximate f . Therefore, it seems convenient to define the KLD as the entropy-like Lyapunov functional for our problem, where the (stationary) reference function φ ref must be an attractor to ensure the Lyapunov-functional condition. Thus, if we choose φ ref (c) = lim s→∞ φ(c; s), assuming that this limit exists, it will minimize the KLD for asymptotically long times. In addition, the definition (36) solves the measure problem posed above, i.e., is indeed the Lyapunov functional of our problem, the natural conjecture is that φ ref (c) = φ H (c) [11]. As a consequence, the challenge is to prove that ∂ s D KL (φ φ H ) ≤ 0 (see Appendix C for a formal expression of ∂ s D KL (φ φ ref ) in the context of the inelastic Boltzmann equation). While in this paper we do not intend to address such a proof from a mathematical point of view, we will provide support by means of MD simulations (see Appendix A for technical details). Before doing that, and in order to put the problem in a proper context, we consider the alternative

is chosen in Equation (36), one simply has
where H * (s) is defined in Equation (34).
Note that ∂ s D KL (φ φ M ) cannot be semi-definite negative for arbitrary initial conditions. For instance, if the initial condition is a Maxwellian, i.e., φ(c; 0) = φ M (c), then it is obvious that D KL (φ φ M )| s=0 = 0 and, given that If that were the case, one could say that the quantity [D KL (φ φ M ) − D KL (φ H φ M )] 2 would always decrease for every initial condition, thus qualifying as a Lyapunov functional. As will be seen below, this expectation is frustrated by our simulation results.
From the formal Sonine expansion (13), one has Now, in the spirit of the truncation approximation of Section 2.3, we can write the approximate expression where a 2 (s) and a 3 (s) are given by Equations (26) and (29), respectively. Since the truncated Sonine approximation is not positive definite, we will take the real part of the right-hand side of Equation (39) for times such that 1 + a 2 (s)S 2 (c 2 ) + a 3 (s)S 3 (c 2 ) < 0 for a certain range of velocities. Figure 6 shows the evolution of D KL (φ φ M ) for the same initial conditions and the same values of α as in Figures 2 and 3, as obtained from our MD simulations (for details, see Appendix A) and from the crude approximation (39). For that initial condition, one clearly has D KL (φ φ M )| s=0 > D KL (φ H φ M ). A monotonic behavior ∂ s D KL (φ φ M ) ≤ 0 is observed only in the cases of small or vanishing inelasticity. For α = 0.1 and 0.4, however, D KL (φ φ M ) does not present a monotonic decay but tends to its asymptotic value D KL (φ H φ M ) from below, there existing a time (s∼2) at which D KL (φ φ M ) exhibits a local minimum. This non-monotonic behavior is certainly exaggerated by the truncated Sonine approximation (39), but it is clearly confirmed by our MD simulations, especially in the case of spheres. Therefore, it is quite obvious that, not unexpectedly, both D KL (φ φ M ) and [D KL (φ φ M ) − D KL (φ H φ M )] 2 must be discarded as a Lyapunov functional for the free cooling of granular gases.
In order to examine how generic the non-monotonic behavior of D KL (φ φ M ) is for high inelasticity, we have taken the case α = 0.1 and considered the same five different initial conditions as in Figures 4 and 5 (see Appendix B). The results are displayed in Figure 7, where we can observe that only the initial condition δ exhibits a non-monotonic behavior, whereas D KL (φ φ M ) decays (grows) monotonically in the cases of the initial conditions Γ and S (M and I). This shows that the nonmoniticity in the time evolution of D KL (φ φ M ) is a rather subtle effect requiring high inelasticity and special initial conditions.

HCS Distribution as a Reference (φ ref = φ H )
By using formal arguments from refs. [39][40][41], García de Soria et al. [11] proved by means of a perturbation analysis around α = 1 that φ H is a unique local minimizer of the entropy production, implying that ∂ s D KL (φ φ H ) ≤ 0, in the quasielastic limit. Those authors also conjectured that this result keeps being valid in the whole inelasticity regime, this conjecture being supported by simulations for α ≥ 0.8 in the freely cooling case.
By performing MD simulations for a wide range of inelasticities (α = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1/ √ 2, 0.8, 0.87, 0.95, and 0.99), we have found further support for the inequality ∂ s D KL (φ φ H ) ≤ 0. As an illustration, Figure 8 shows the evolution of D KL (φ φ H ) for α = 0.1, 0.4, 0.87, and 1, starting from the same initial states as in Figures 2, 3, and 6. In the evaluation of D KL (φ φ H ), we have used the simulation results for both the transient distribution φ(c; s) and the asymptotic HCS distribution φ H (c) (see Appendix A). Our MD results are compared with a theoretical approximation similar to that of Equation (39), i.e., where again the real part of the right-hand side is taken if 1 + a 2 (s)S 2 (c 2 ) + a 3 (s)S 3 (c 2 ) < 0 for a certain range of velocities. The results (both from MD and from the approximate theory) displayed in Figure 8 show that D KL (φ φ H ) indeed decays monotonically to 0, even for very strong inelasticity, thus supporting its status as a very sound candidate of Lyapunov functional. It is also interesting to note that the characteristic relaxation time is generally shorter for disks than for spheres and tends to decrease with increasing inelasticity.  In order to reinforce the monotonic decay of D KL (φ φ H ) observed in Figure 8 for several representative values of the coefficient of restitution, let us now take the most demanding case (α = 0.1) and choose the five initial conditions already considered in Figures 4, 5, and 7 (see Appendix B). Figure 9 shows that the evolution of D KL (φ φ H ) keeps being monotonic for this wide spectrum of representative initial conditions, the relaxation to the HCS being again faster for disks than for spheres. It is also interesting to comment that, although the largest initial divergence corresponds to the initial distribution δ, this divergence decays more rapidly than the other four ones, and even seems to overtake the divergence associated with the initial condition Γ. While a rigorous mathematical proof of ∂ s D KL (φ φ H ) ≤ 0 is still lacking (see, however, Ref. [42] for the sketch of a proof in the context of the linear Boltzmann equation), we will now prove this inequality by using a simplified toy model. We start from the infinite series expansion (13) and imagine a formal bookkeeping parameter in front of the Sonine summation. Then, to the second order in , Next, taking into account the orthogonality condition (15), we get Interestingly, this approximation preserves the positive-definiteness of the KLD. Note also that, to order 2 , . Finally, consistent with the derivation of Equations (20) and (25), we neglect the cumulants a k with k ≥ 3 and apply Equation (25) to where we have formally set = 1. Although a certain number of approximations have been done to derive the toy model (43), it undoubtedly provides further support to the conjecture ∂ s D KL (φ φ H ) ≤ 0.

Relative Entropy of φ H with Respect to φ M
It is well known that, in a freely cooling granular gas, the HCS VDF is generally close (at least within the range of thermal velocities) to a Maxwellian. In particular, the cumulants a H k are rather small in magnitude, except at large inelasticity (see Figure 1). On the other hand, the HCS VDF exhibits an exponential high-velocity tail, ln φ H (c)∼−c, with respect to the Maxwellian behavior, ln φ M (c)∼−c 2 [17,23,43].
Here, we have one more tool to measure how far φ M (c) is from φ H (c), namely the KLD from φ M to φ H (or relative entropy of φ H with respect to φ M ), i.e., D KL (φ H φ M ). Note, however, that, as said at the beginning of this section, the KLD is not a real metric since it does not fulfill either symmetry or triangle inequality properties of a distance. Figure 10 displays the α-dependence of D KL (φ H φ M ) for both disks and spheres, as obtained from our MD simulations (see again Appendix A) and from the simple estimate (39) with a 2 (s) → a H 2 and a 3 (s) → a H 3 . We can observe that the theoretical truncated approach successfully captures (i) a weak influence of dimensionality (in contrast to the fourth and sixth cumulants plotted in Figure 1), (ii) a crossover from for smaller inelasticity, and (ii) a non-monotonic dependence on α, with a (small but nonzero) local minimum at about α = 1/ √ 2 0.71 and a local maximum at about α = 0.87. The latter property implies that, in the region 0.6 α < 1, three systems differing in the value of α may share the same divergence of φ M from φ H . The qualitative shape of D KL (φ H φ M ) as a function of α agrees with a toy model analogous to that of Equation (43a), namely

Summary and Conclusions
In this work, we have mainly focused on the role as a potential entropy-like Lyapunov functional played by the KLD of a reference VDF (φ ref ) with respect to the spatially homogeneous time-dependent VDF (φ), i.e., D KL (φ φ ref ), as supported by MD simulations in a freely cooling granular-gas model.
First, we have revisited the problem of obtaining, by kinetic theory methods, simple approximations for the HCS fourth (a H 2 ) and sixth (a H 3 ) cumulants, and have derived explicit time-dependent solutions, a 2 (s) and a 3 (s), for arbitrary (homogeneous) initial conditions. Comparison with our MD results shows an excellent general performance of a H 2 and a 2 (s) for values of the coefficient of restitution as low as α = 0.1 and for a variety of initial conditions. In the case of the sixth cumulant, however, the agreement is mainly semi-quantitative. In any case, our MD data for a H 2 and a H 3 agree very well with previous simulations of the inelastic Boltzmann equation [18][19][20]29], thus validating the applicability of kinetic theory (including the Stosszahlansatz) even for high inelasticity. We emphasize that, to the best of our knowledge, such a comprehensive MD analysis of the fourth and sixth cumulants had not been carried out before. We are not aware either of a previous (approximate) theoretical derivation of the time-dependent quantities a 2 (s) and a 3 (s).
As a first candidate to a Lyapunov functional, we have considered the KLD with a Maxwellian reference VDF (φ ref = φ M ). However, this possibility is clearly discarded as both simulation and a simple theoretical approach show that D KL (φ φ M ) does not relax monotonically for highly inelastic systems and certain initial conditions. On the other hand, when the asymptotic HCS VDF is chosen as a reference (φ ref = φ H ), the results show that the relaxation of D KL (φ φ H ) is monotonic for a wide spectrum of inelasticities and initial conditions. This is further supported by a simplified toy model, according to which ∂ s D KL (φ φ H )∼−[a 2 (s) − a H 2 ] 2 ≤ 0. While simulation results supporting the conjecture ∂ s D KL (φ φ H ) ≤ 0 had been presented before [11], it is subjected here to more stringent tests by considering highly dissipative collisions (α = 0.1 and 0.4) and a repertoire of different initial conditions. In fact, it is only under those more extreme conditions when one can reject the Maxwellian as a proper candidate for the reference VDF.
We have also used D KL (φ H φ M ) to characterize the departure of the Maxwellian distribution as an approximation to the actual HCS distribution. Interestingly, we found a non-monotonic influence of the coefficient of restitution on D KL (φ H φ M ), with a (nonzero) local minimum at α 1/ √ 2 0.71 and a (small) local maximum at α 0.87. This non-monotonicity implies a degeneracy of D KL (φ H φ M ) in the sense that three different coefficients of restitution (within the region 0.6 α < 1) may share a common value of the KLD from φ M to φ H . The analysis of D KL (φ H φ M ) is an additional asset of our work.
We expect that the results presented in this paper may stimulate further studies on the quest of proving (or disproving, if a counterexample is found) the extension of Boltzmann's celebrated H-theorem to the realm of dissipative inelastic collisions in homogeneous states. In this respect, it must be remarked that, since the simulation results we have presented are obtained from the MD technique (which numerically solves Newton's equations of motion) and not from the DSMC method (which numerically solves the Boltzmann equation), it is not obvious from a strict mathematical point of view that the obtained results imply the decay of the KLD in the context of the Boltzmann equation. On the other hand, on physical grounds, it is expected that such an implication holds.
As a final remark, it is worth emphasizing that, even if some kind of generalized H-theorem could be proved for homogeneous states, its extension to inhomogeneous situations would be far from trivial since the HCS is unstable under long-wavelength perturbations.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:
Since the DynamO code is designed for three-dimensional setups, we used it for the two-dimensional case by imposing a coordinate z = 0 to every particle and carefully avoiding any overlap in the initial ordered arrangement. The system melted very quickly and no inhomogeneities were observed thereafter. A velocity rescaling was done periodically in order to avoid numerical errors due to the cooling process and extremely small numbers.
To represent the VDF and the KLD in simulations, let us first introduce the probability distribution function of the velocity modulus, where in the second step we have assumed that the VDF φ(c; s) is isotropic and Ω d is the d-dimensional solid angle. Thus, Equation (36) can be rewritten as The functions Φ(c; s) and Φ H (c) are numerically approximated by a discrete histogram, with a certain constant bin width ∆c, i.e., Here, N i (s) is the number of particles with a speed c inside the interval c i − ∆c/2 ≤ c < c i + ∆c/2, N H i is evaluated by averaging N i (s) between s = 10 to s = 40 with a timestep δs = 0.2, and M is the total number of bins considered. In consistency with Equation (A3), the Maxwellian VDF is also discretized as where erf(x) = 2 √ π x 0 dt e −t 2 is the error function.
where Φ M (c i ) is given by Equation (A4). Analogously, A comment is now in order. In the case of elastic collisions (α = 1), one obviously should have Φ H (c i ) = Φ M (c i ) and hence D KL (φ H φ M )| α=1 = 0. However, since Φ H (c i ) is evaluated in simulations by Equation (A3) for any α, the equality Φ H (c i ) = Φ M (c i ) for α = 1 is not identically verified bin to bin due to fluctuations. As a consequence, in the simulations, D KL (φ H φ M )| α=1 ∼ 10 −5 = 0. This is an unavoidable background noise that was subtracted from the KLD obtained by simulations, i.e., We have chosen the values ∆c = 0.03 and M = 200. The results presented in the main text for any given quantity are obtained by averaging over 50 independent realizations.

Appendix B. Initial Conditions
For the analysis of the evolution of a 2 (s), a 3 (s), D KL (φ φ M ), and D KL (φ φ H ) with α = 0.1, we have chosen five different initial conditions. The first one is the same as considered in Figures 2, 3, 6, and 8, i.e., an ordered crystalized configuration with isotropic velocities of a common magnitude. In terms of the distribution defined by Equation (A1), this initial condition reads which will be labeled with the Greek letter δ. The second initial distribution is just a Maxwellian (label M), i.e., Next, we choose the gamma distribution (label Γ) normalized to c 2 = d 2 , namely where θ > 0 can be freely chosen. The fourth-and sixth-order moments are c 4 = d(d+2θ) 4 and c 6 = d(d+2θ)(d+4θ)
The remaining two initial conditions are prepared by applying a coefficient of normal restitution α 0 and allowing the system to reach the corresponding steady state (in the scaled quantities). Then, at s = 0, the coefficient of restitution is abruptly changed to α = 0.1 and the evolution toward the corresponding HCS is monitored. We have taken two classes of values of α 0 : (a) α 0 < 1, corresponding to dissipative inelastic collisions (label I), and (b) α 0 > 1 [44], corresponding to "super-elastic" collisions (label S). More specifically, for the preparation of the initial state I, we have chosen α 0 = 0.29 and 0.27 for d = 2 and 3, respectively; the state S has been prepared with α 0 = 1.29 and 1.47 for d = 2 and 3, respectively. Table A1 displays the values of a 2 and a 3 corresponding to, in order of increasing a 2 , the initial states δ, M, I, Γ, and S.