==== Front Entropy (Basel) Entropy (Basel) entropy Entropy 1099-4300 MDPI 33287073 10.3390/e22111308 entropy-22-01308 Article Kullback–Leibler Divergence of a Freely Cooling Granular Gas https://orcid.org/0000-0002-7889-1312Megías Alberto 1 https://orcid.org/0000-0002-9564-5180Santos Andrés 2* 1 Departamento de Física, Universidad de Extremadura, E-06006 Badajoz, Spain; albertom@unex.es 2 Departamento de Física and Instituto de Computación Científica Avanzada (ICCAEx), Universidad de Extremadura, E-06006 Badajoz, Spain * Correspondence: andres@unex.es; Tel.: +34-924-289-651 17 11 2020 11 2020 22 11 130807 10 2020 13 11 2020 © 2020 by the authors.2020Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).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 α. Kullback–Leibler divergencegranular gaseskinetic theorymolecular dynamics ==== Body 1. 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), fHCS. 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 fHCS 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. 2. Free Cooling Evolution of Velocity Cumulants 2.1. 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] (1) ∂tf(v1;t)=nσd−1I[v1|f,f]≡nσd−1∫dv2∫+dσ^(v12·σ^)α−2f(v1″;t)f(v2″;t)−f(v1;t)f(v2;t), where n is the number density, v12=v1−v2 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 v12·σ^>0, and (2) v1″=v1−1+α2α(v12·σ^)σ^,v2″=v2+1+α2α(v12·σ^)σ^ are precollisional velocities. Note that we have defined the VDF with the normalization condition ∫dvf(v;t)=1. An important quantity is the granular temperature defined as (3) T(t)=md〈v2〉,〈X(v)〉≡∫dvX(v)f(v;t). Taking moments in Equation (1), one finds the cooling equation (4) ∂tT(t)=−ζ(t)T(t), where the cooling rate is given by (5a) ζ(t)=−mnσd−1T(t)d∫dvv2I[v|f,f]=(1−α2)mnσd−1T(t)π(d−1)/24dΓ(d+32)〈〈v123〉〉, (5b) 〈〈X(v1,v2)〉〉≡∫dv1∫dv2X(v1,v2)f(v1;t)f(v2;t). Let us introduce the thermal velocity vth(t)≡2T(t)/m, which allows us to define the rescaled VDF ϕ(c;s) as (6) f(v;t)=vth−d(t)ϕ(c;s),c≡vvth(t), where the variable s in ϕ(c;s) is a scaled time defined by (7) s(t)=12∫0tdt′ν(t′),ν(t)≡κnσd−1vth(t),κ≡2π(d−1)/2Γ(d2). 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 (8) κ2∂sϕ(c;s)+μ2(s)d∂∂c·cϕ(c;s)=I[c|ϕ,ϕ],μk(s)≡−∫dcckI[c|ϕ,ϕ], where we have taken into account that ζ(t)/nσd−1vth(t)=2μ2(s)/d. The associated hierarchy of moment equations is (9) κ2∂s〈ck〉=Fk(s)≡kμ2(s)d〈ck〉−μk(s). Note that F0=F2=0, since μ0=0 and 〈c2〉=d2. 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;s)→ϕH(c), where ϕH(c) satisfies the integrodifferential equation (10) μ2Hd∂∂c·cϕH(c)=I[c|ϕH,ϕH]. 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)/TH(t)=const, so that the solution to Equation (4) gives rise to the well-known cooling Haff’s law [14,15,31] (11) TH(t)=TH(t0)1+12ζH(t0)(t−t0)2, where t0 being an arbitrary time belonging to the HCS regime. Also in the HCS regime, μ2(s)→μ2H=const and thus Equation (4) becomes ∂sTH(s)=−(4/κd)μ2HTH(s), whose solution is (12) TH(s)=TH(s0)e−4μ2H(s−s0)/κd. Therefore, in the HCS, the temperature decays exponentially with the average number of collisions per particle. 2.2. Sonine Expansion Formalism The Maxwell–Boltzmann VDF ϕM(c)=π−d/2e−c2 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, (13) ϕ(c;s)=ϕM(c)1+∑k=2∞ak(s)Sk(c2), where (14) Sk(x)=Lk(d2−1)(x)=∑j=0k(−1)jΓd2+kΓd2+j(k−j)!j!xj where are Sonine (or generalized Laguerre) polynomials, which satisfy the orthogonalization condition (15) 〈Sk|Sk′〉≡∫dcϕM(c)Sk(c2)Sk′(c2)=Nkδk,k′,Nk≡Γd2+kΓd2k!. In Equation (13), the Sonine coefficient ak(s) is the 2k-th cumulant of the VDF at time s. According to Equation (15), (16) ak(s)=〈Sk(c2)〉Nk. In particular, a0(s)=1, a1(s)=0, and (17) a2(s)=4d(d+2)〈c4〉−1,a3(s)=1+3a2−8d(d+2)(d+4)〈c6〉. 2.3. Truncated Sonine Approximation Thus far, all the results presented in Section 2.1 and Section 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 ak with k≥4 in Equation (13), as well as the nonlinear terms a22, a2a3, and a32 in the bilinear collision operator I[c|ϕ,ϕ]. Given a functional X[ϕ] of the scaled VDF ϕ(c), we will use the notation L3X to denote the result of that truncation and linearization procedure. Furthermore, if a3 is also neglected, the corresponding approximation will be denoted by L2X. In particular, in the case of the collisional moments μ2, μ4, and μ6, one has (18) L3μ2=A0+A2a2+A3a3,L3μ4=B0+B2a2+B3a3,L3μ6=C0+C2a2+C3a3, where the expressions for the coefficients Ai, Bi, and Ci as functions of α and d can be found in Ref. [29] and in Appendix A of Ref. [19]. Obviously, L2μ2, L2μ4, and L2μ6 are obtained by formally setting A3→0, B3→0, and C3→0, respectively. Let us first use the simple approximation L2 to estimate a2H. From Equation (9), we have that F4H=0. Thus, the obvious approximation [17] consists of (19) L2F4H=0⇒a2H,a=(d+2)A0−B0B2−(d+2)(A2+A0)=16(1−α)(1−2α2)9+24d−(41−8d)α+30(1−α)α2, where, in the last steps, use has been made of the explicit expressions of A0, A2, B0, and B2. However, this is not by any means the only possibility of estimating a2H [18,19,33]. In particular, one can start from the logarithmic time derivative of the fourth moment and then take (20) L2F4H〈c4〉H=0⇒a2H,b=(d+2)A0−B0B2−B0−(d+2)A2=16(1−α)(1−2α2)25+24d−(57−8d)α−2(1−α)α2. Note that (21) a2H,aa2H,b=1+a2H,a=11−a2H,b. Both approximations (a2H,a and a2H,b) are practically indistinguishable in the region 0.6≲α<1, but a2H,b is much more accurate than a2H,a for higher inelasticity [18,19]. Next, to estimate a3H, we start from the exact condition F6H=0 and carry out either the linearization (22) L3F6H=0⇒a3H,a=Ga(a2H)≡C0−34(d+2)(d+4)A0+C2−34(d+2)(d+4)(3A0+A2)a2H34(d+2)(d+4)(A3−A0)−C3 or, alternatively, (23) L3F6H〈c6〉H=0⇒a3H,b=Gb(a2H)≡C0−34(d+2)(d+4)A0+C2−3C0−34(d+2)(d+4)A2a2H34(d+2)(d+4)A3−C3−C0. In Equations (22) and (23), a3H is expressed in terms of a2H. Using Equations (19) and (20), four possibilities in principle arise, namely (24) a3H,aa=Ga(a2H,a),a3H,ab=Ga(a2H,b),a3H,ba=Gb(a2H,a),a3H,bb=Gb(a2H,b). Comparison with DSMC results shows that the best general estimates are provided by a3H,aa and a3H,ab. In what follows, we choose a2H,b for the fourth cumulant and, for the sake of consistency with that choice, we adopt a3H,ab for the sixth cumulant. To simplify the notation, we make a2H,b→a2H and a3H,ab→a3H. Once the (approximate) HCS values a2H and a3H have been obtained, we turn our attention to the evolution equations of a2(s) and a3(s). Approximating Equation (9) with k=4 as κ2∂sln〈c4〉=L2F4(s)/〈c4〉, one obtains (25) ∂sa2(s)=−K21+a2(s)a2(s)−a2H,K2≡8d(d+2)κB2−B0−(d+2)A2. Its solution is (26) a2(s)=a2H+1+a2HX0eγs−1,X0≡1+a2(0)a2(0)−a2H,γ≡1+a2HK2. Analogously, if Equation (9) with k=6 is approximated as κ2∂s〈c6〉=L3F6(s), the resulting evolution equation for a3 is (27) ∂sa3(s)=3∂sa2(s)−K2′a2(s)−a2H−K3a3(s)−a3H, where (28a) K2′≡16d(d+2)(d+4)κ34(d+2)(d+4)(A2+3A0)−C2, (28b) K3≡16d(d+2)(d+4)κ34(d+2)(d+4)(A3−A0)−C3. Taking into account Equation (26), the solution to Equation (27) is (29a) a3(s)=a3H+Y0e−K3s+1+a2H3X0eγs−1+K2′K3+32F11,K3γ;K3γ+1;X0eγs, (29b) Y0≡a3(0)−a3H−1+a2H3X0−1+K2′K3+32F11,K3γ;K3γ+1;X0, where 2F1(a,b;c;z) is the hypergeometric function [34]. As far as we know, Equations (26) and (29) had not been obtained before. 2.4. Comparison with MD Simulations The approximate theoretical predictions for a2H and a3H 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 a2H and a3H 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×103 for disks and L/σ≈4×102 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. Figure 1a,b show the α-dependence of a2H and a3H, 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 a2H,b (see Equation (20)) performs very well, even if the fourth cumulant is not small (e.g., a2H∼0.2 at α=0.1). The approximate sixth cumulant a3H,ab (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 a2H changes from negative to positive values at α≃1/2≃0.71, a3H is always negative. Note that, for large inelasticity, the cumulants a2H and a3H 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 |a3H|≪|a2H|≪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 a2(s) and a3(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/2vth(0) along uniformly randomized directions. Therefore, at s=0, 〈ck〉=(d/2)k/2, so that a2(0)=−2d+2 and a3(0)=−16(d+2)(d+4). Figure 2 and Figure 3 compare our MD results with the theoretical predictions (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 a2H 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 (a2(0)=−12 and −25 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 (a3(0)=−23 and −1635 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. In Figure 2 and Figure 3, the initial values a2(0) and a3(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 α=0.1 are {a2H,a3H}={0.206,−0.143} and {0.150,−0.077} for d=2 and d=3, respectively. Thus, we have chosen the same initial distribution (hereafter labeled as δ) as in Figure 2 and Figure 3 as a representative example of a2(0)<0, the Maxwellian distribution (labeled as M) with a2(0)=0, another one (labeled as I) with 0a2H. The details of those five distributions can be found in Appendix B and the corresponding values of a2(0) and a3(0) are shown in Table A1. In the case of a2(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 a3(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. 3. KLD as a Lyapunov Functional In this section, we restrict ourselves to spatially homogeneous states. 3.1. 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 (30) H(t)=∫dvf(v;t)lnf(v;t) 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 {v1″,v2″} and a relative orientation characterized by the unit vector −σ^ (with v12″·σ^<0). After collision, the velocities are, in agreement with Equation (2), given by (31) C−σ^{v1″,v2″}={v1,v2},v1,2=v1,2″∓1+α2(v12″·σ^)σ^. Next, suppose two colliding particles with precollision velocities {v1,v2} and a relative orientation characterized by the unit vector σ^ (with v12·σ^>0). In that case, (32) Cσ^{v1,v2}=Cσ^C−σ^{v1″,v2″}={v1′,v2′},v1,2′=v1,2∓1+α2(v12·σ^)σ^. Comparison with Equation (2) shows that (33) v1,2′=v1,2″±1−α22α(v12·σ^)σ^,v12′·σ^=α2v12″·σ^. Thus, Cσ^C−σ^{v1″,v2″}≠{v1″,v2″} 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), (34) H(s)=∫dvf(v,t)lnf(v,t)=H*(s)−d2ln2T(s)m,H*(s)≡∫dcϕ(c,s)lnϕ(c,s). From Haff’s law, Equation (12), it turns out that (in the HCS) HH* is stationary but HH(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 c→c˜=w(c), ϕ(c,s)→ϕ˜(c˜,s)=J−1ϕ(c,s), 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. 3.2. 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 (35) DKL(f∥g)=∫Xdxf(x)lnf(x)g(x), where x is a random vector variable defined on the set X. The quantity DKL(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), DKL(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 (36) DKL(f∥fref)=DKL(ϕ∥ϕref)=∫dcϕ(c;s)lnϕ(c;s)ϕref(c) 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)=lims→∞ϕ(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., DKL(ϕ∥ϕref)=DKL(ϕ˜∥ϕ˜ref) for any invertible transformation c→c˜=w(c). If DKL(ϕ∥ϕref) 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 ∂sDKL(ϕ∥ϕH)≤0 (see Appendix C for a formal expression of ∂sDKL(ϕ∥ϕ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 choice ϕref=ϕM. 3.3. MD Simulations 3.3.1. Maxwellian Distribution as a Reference (ϕref=ϕM) If ϕref=ϕM is chosen in Equation (36), one simply has (37) DKL(ϕ∥ϕM)=H*(s)+d21+lnπ, where H*(s) is defined in Equation (34). Thus, DKL(ϕ∥ϕM) differs from H*(s) by a constant, so that ∂sDKL(ϕ∥ϕM)=∂sH*(s). Note that ∂sDKL(ϕ∥ϕ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 DKL(ϕ∥ϕM)s=0=0 and, given that lims→∞DKL(ϕ∥ϕM)=DKL(ϕH∥ϕM)>0, it is impossible that ∂sDKL(ϕ∥ϕM)≤0 for all s. Nevertheless, in principle, it might happen that ∂sDKL(ϕ∥ϕM)≤0 for the class of initial conditions such that DKL(ϕ∥ϕM)s=0≥DKL(ϕH∥ϕM), while ∂sDKL(ϕ∥ϕM)≥0 for the complementary class of initial conditions such that DKL(ϕ∥ϕM)s=0≤DKL(ϕH∥ϕM). If that were the case, one could say that the quantity DKL(ϕ∥ϕM)−DKL(ϕ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 (38) DKL(ϕ∥ϕM)=∫dcϕM(c)1+∑k=2∞ak(s)Sk(c2)ln1+∑k=2∞ak(s)Sk(c2). Now, in the spirit of the truncation approximation of Section 2.3, we can write the approximate expression (39) DKL(ϕ∥ϕM)≈∫dcϕM(c)1+a2(s)S2(c2)+a3(s)S3(c2)ln1+a2(s)S2(c2)+a3(s)S3(c2), where a2(s) and a3(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+a2(s)S2(c2)+a3(s)S3(c2)<0 for a certain range of velocities. Figure 6 shows the evolution of DKL(ϕ∥ϕM) for the same initial conditions and the same values of α as in Figure 2 and Figure 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 DKL(ϕ∥ϕM)s=0>DKL(ϕH∥ϕM). A monotonic behavior ∂sDKL(ϕ∥ϕM)≤0 is observed only in the cases of small or vanishing inelasticity. For α=0.1 and 0.4, however, DKL(ϕ∥ϕM) does not present a monotonic decay but tends to its asymptotic value DKL(ϕH∥ϕM) from below, there existing a time (s∼2) at which DKL(ϕ∥ϕ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 DKL(ϕ∥ϕM) and DKL(ϕ∥ϕM)−DKL(ϕ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 DKL(ϕ∥ϕM) is for high inelasticity, we have taken the case α=0.1 and considered the same five different initial conditions as in Figure 4 and Figure 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 DKL(ϕ∥ϕ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 DKL(ϕ∥ϕM) is a rather subtle effect requiring high inelasticity and special initial conditions. 3.3.2. 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 ∂sDKL(ϕ∥ϕ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 ∂sDKL(ϕ∥ϕH)≤0. As an illustration, Figure 8 shows the evolution of DKL(ϕ∥ϕH) for α=0.1, 0.4, 0.87, and 1, starting from the same initial states as in Figure 2, Figure 3, and Figure 6. In the evaluation of DKL(ϕ∥ϕ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., (40) DKL(ϕ∥ϕH)≈∫dcϕM(c)1+a2(s)S2(c2)+a3(s)S3(c2)ln1+a2(s)S2(c2)+a3(s)S3(c2)1+a2HS2(c2)+a3HS3(c2), where again the real part of the right-hand side is taken if 1+a2(s)S2(c2)+a3(s)S3(c2)<0 for a certain range of velocities. The results (both from MD and from the approximate theory) displayed in Figure 8 show that DKL(ϕ∥ϕ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 DKL(ϕ∥ϕ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 Figure 4, Figure 5, and Figure 7 (see Appendix B). Figure 9 shows that the evolution of DKL(ϕ∥ϕ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 ∂sDKL(ϕ∥ϕ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 ϵ, (41) ϕ(c;s)ϕM(c)lnϕ(c;s)ϕH(c)=ϵ∑k=2∞ak(s)−akHSk(c2)+ϵ22∑k,k′=2∞ak(s)−akHak′(s)−ak′HSk(c2)Sk′(c2)+O(ϵ3). Next, taking into account the orthogonality condition (15), we get (42a) DKL(ϕ∥ϕH)=ϵ22∑k∞Nkak(s)−akH2+O(ϵ3), (42b) ∂sDKL(ϕ∥ϕH)=ϵ2∑k∞Nkak(s)−akH∂sak(s)+O(ϵ3). Interestingly, this approximation preserves the positive-definiteness of the KLD. Note also that, to order ϵ2, DKL(ϕ∥ϕH) is symmetric under the exchange ϕ↔ϕH, i.e., DKL(ϕ∥ϕH)−DKL(ϕH∥ϕ)=O(ϵ3). Finally, consistent with the derivation of Equations (20) and (25), we neglect the cumulants ak with k≥3 and apply Equation (25) to obtain (43a) DKL(ϕ∥ϕH)≈d(d+2)16a2(s)−a2H2, (43b) ∂sDKL(ϕ∥ϕH)≈−d(d+2)8K21+a2(s)a2(s)−a2H2≤0, 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 ∂sDKL(ϕ∥ϕH)≤0. 3.3.3. 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 akH 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)∼−c2 [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., DKL(ϕ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 DKL(ϕH∥ϕM) for both disks and spheres, as obtained from our MD simulations (see again Appendix A) and from the simple estimate (39) with a2(s)→a2H and a3(s)→a3H. 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 DKL(ϕH∥ϕM)d=2DKL(ϕH∥ϕM)d=3 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 DKL(ϕH∥ϕM) as a function of α agrees with a toy model analogous to that of Equation (43a), namely DKL(ϕH∥ϕM)≈d(d+2)16a2H2. 4. 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., DKL(ϕ∥ϕ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 (a2H) and sixth (a3H) cumulants, and have derived explicit time-dependent solutions, a2(s) and a3(s), for arbitrary (homogeneous) initial conditions. Comparison with our MD results shows an excellent general performance of a2H and a2(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 a2H and a3H 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 a2(s) and a3(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 DKL(ϕ∥ϕ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 DKL(ϕ∥ϕH) is monotonic for a wide spectrum of inelasticities and initial conditions. This is further supported by a simplified toy model, according to which ∂sDKL(ϕ∥ϕH)∼−[a2(s)−a2H]2≤0. While simulation results supporting the conjecture ∂sDKL(ϕ∥ϕ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 DKL(ϕ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 DKL(ϕ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 DKL(ϕ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 DKL(ϕ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. Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Author Contributions A.S. proposed the idea and A.M. carried out the simulations. Both authors participated in the analysis and discussion of the results and worked on the revision and writing of the final manuscript. All authors have read and agreed to the published version of the manuscript. Funding The authors acknowledge financial support from the Spanish Agencia Estatal de Investigación through Grant No. FIS2016-76359-P and the Junta de Extremadura (Spain) through Grant No. GR18079, both partially financed by Fondo Europeo de Desarrollo Regional funds. A.M. is grateful to the Spanish Ministerio de Ciencia, Innovación y Universidades for a predoctoral fellowship FPU2018-3503. Conflicts of Interest The authors declare no conflict of interest. Abbreviations The following abbreviations are used in this manuscript: DSMC Direct simulation Monte Carlo HCS Homogenous cooling state KLD Kullback–Leibler divergence MD Molecular dynamics VDF Velocity distribution function Appendix A. Simulation and Numerical Details Event-driven MD simulations were carried out using the DynamO software [13] on N particles in a d-dimensional cubic box of side L with periodic boundary conditions. We chose (N,L/σ)=(104,4641.58) and (1.35×104,407.16) for disks and spheres, respectively. Thus, the associated number densities were nσ2=4.64×10−4 (disks) and nσ3=2.00×10−4 (spheres). The critical lengths for the development of instabilities at those densities [14] are estimated to be Lc/σ≈8.3×103 (disks) and 1.3×104 (spheres) for the most demanding case (α=0.1). This represents ratios Lc/L≈1.8 (disks) Lc/L≈32.8 (spheres). Those ratios generally increase with decreasing inelasticity. For instance, at α=0.87, one finds Lc/L≈3.4 (disks) and 62.5 (spheres). Therefore, the simulations are performed in the region of parameters where the systems are stable. 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, (A1) Φ(c;s)=cd−1∫dc^ϕ(c;s)=Ωdcd−1ϕ(c;s),Ωd≡2πd/2Γ(d/2), 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 (A2) DKL(ϕ∥ϕref)=∫0∞dcΦ(c;s)lnΦ(c;s)Φref(c). The functions Φ(c;s) and ΦH(c) are numerically approximated by a discrete histogram, with a certain constant bin width Δc, i.e., (A3) Φ(ci;s)≈Ni(s)NΔc,ΦH(ci)≈NiHNΔc,ci=i−12Δc,i=1,2,…,M. Here, Ni(s) is the number of particles with a speed c inside the interval ci−Δc/2≤c0 can be freely chosen. The fourth- and sixth-order moments are 〈c4〉=d(d+2θ)4 and 〈c6〉=d(d+2θ)(d+4θ)8, so that a2=2(θ−1)d+2 and a3=−8(θ−1)(θ−2)(d+2)(d+4). Here, we have taken θ=2.16 and 2.45 for d=2 and 3, respectively. 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 a2 and a3 corresponding to, in order of increasing a2, the initial states δ, M, I, Γ, and S. entropy-22-01308-t0A1_Table A1Table A1 Values of the fourth and sixth cumulants for the initial distributions δ, M, I, Γ, and S (see text). δ M I Γ S a2(0) −0.500(d=2)−0.400(d=3) 0 0.151(d=2)0.111(d=3) 0.580(d=2)0.580(d=3) 0.885(d=2)0.792(d=3) a3(0) −0.667(d=2)−0.457(d=3) 0 −0.080(d=2)−0.046(d=3) −0.062(d=2)−0.149(d=3) −4.733(d=2)−2.219(d=3) Appendix C. Formal Expression for ∂sDKL(ϕ∥ϕref) The aim of this appendix is to derive a formal expression for ∂sDKL(ϕ∥ϕref) by following the same steps as in the proof of the conventional H-theorem [6]. Let us consider a generic test function ψ(c). By standard steps, one can easily obtain [14] (A10) J[ψ]≡∫dcψ(c)I[c1|ϕ,ϕ]=12∫dc1∫dc2∫+dσ^(c12·σ^)ϕ(c1)ϕ(c2)ψ(c1′)+ψ(c2′)−ψ(c1)−ψ(c2). Next, we perform the change of variables {c1,c2,σ^}→{c1′,c2′,−σ^} and take into account that dc1′dc2′=αdc1dc2 and c12′·σ^=−αc12·σ^ to obtain (A11) J[ψ]=α−22∫dc1′∫dc2′∫+dσ^(c12′·σ^)ϕ(c1)ϕ(c2)ψ(c1′)+ψ(c2′)−ψ(c1)−ψ(c2)=α−22∫dc1∫dc2∫+dσ^(c12·σ^)ϕ(c1′′)ϕ(c2′′)ψ(c1)+ψ(c2)−ψ(c1′′)−ψ(c2′′), where in the second equality we have just renamed {c1′,c2′,c1,c2}→{c1,c2,c1″,c2″}. Taking the average between Equations (A10) and (A11), we arrive at (A12) J[ψ]=14∫dc1∫dc2∫+dσ^(c12·σ^){ϕ(c1)ϕ(c2)ψ(c1′)+ψ(c2′)−ψ(c1)−ψ(c2)−ϕ(c1″)ϕ(c2″)α2ψ(c1′′)+ψ(c2′′)−ψ(c1)−ψ(c2)}. Now, we start from the KLD defined by Equation (36) and use the Boltzmann equation (8) to get (13) κ2∂sDKL(ϕ∥ϕref)=Jlnϕϕref−μ2d∫dclnϕ(c)ϕref(c)∂∂c·cϕ(c). where we have taken into account that ϕref(c) and ∫dcϕ(c)=1 are independent of time. Integration by parts of the second term on the right-hand side of Equation (A13) yields (A14) κ2∂sDKL(ϕ∥ϕref)=Jlnϕϕref−μ21+1d∫dcϕ(c)c·∂∂clnϕref(c). Finally, making use of Equation (A12) with ψ(c)=ln[ϕ(c)/ϕref(c)], we obtain (A15) κ2∂sDKL(ϕ∥ϕref)=14∫dc1∫dc2∫+dσ^(c12·σ^)[ϕ(c1)ϕ(c2)lnϕ(c1′)ϕ(c2′)ϕref(c1)ϕref(c2)ϕ(c1)ϕ(c2)ϕref(c1′)ϕref(c2′)−ϕ(c1″)ϕ(c2″)α2lnϕ(c1″)ϕ(c2″)ϕref(c1)ϕref(c2)ϕ(c1)ϕ(c2)ϕref(c1″)ϕref(c2″)]−μ2d∫dcϕ(c)c·∂∂clnϕref(c)ϕM(c), where we have taken into account that −∫dcϕ(c)c·∂∂clnϕM(c)=2∫dcc2ϕ(c)=d. Equation (A15) does not particularly simplify if ϕref=ϕH. However, in the case ϕref=ϕM, a somewhat simpler expression can be found. First, the last term on the right-hand side of Equation (A15) vanishes if ϕref=ϕM. Second, we can use the decomposition J[ln(ϕ/ϕM)]=J[lnϕ]−J[lnϕM] and take into account that lnϕM(c)=−c2+const and, therefore, J[lnϕM]=μ2 [see Equation (8)]. As a consequence, (A16) κ2∂sDKL(ϕ∥ϕM)=14∫dc1∫dc2∫+dσ^(c12·σ^)[ϕ(c1)ϕ(c2)lnϕ(c1′)ϕ(c2′)ϕ(c1)ϕ(c2)−ϕ(c1″)ϕ(c2″)α2lnϕ(c1″)ϕ(c2″)ϕ(c1)ϕ(c2)]−μ2. In the special case of elastic collisions (α=1), one has μ2=0 and ci″=ci′, so that the standard H-theorem is recovered, namely (A17) κ2∂sDKL(ϕ∥ϕM)α=1=−14∫dc1∫dc2∫+dσ^(c12·σ^)ϕ(c1′)ϕ(c2′)−ϕ(c1)ϕ(c2)lnϕ(c1′)ϕ(c2′)ϕ(c1)ϕ(c2)≤0. Figure 1 Plot of (a) the HCS fourth cumulant a2H and (b) the HCS sixth cumulant a3H versus the coefficient of restitution α. Symbols represent simulation results: MD (this work) for disks (∘) and spheres (∆), and DSMC [18,19,29] for disks (×) and spheres (□). The lines are the theoretical predictions a2H,b (see Equation (20)) and a3H,ab (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. Figure 2 Evolution of the fourth cumulant a2(s) as a function of the average number of collisions per particle for (a) disks and (b) spheres. Symbols represent MD simulation results, while the lines correspond to the theoretical prediction (26). The values of the coefficient of restitution are (from top to bottom) α=0.1 (□), 0.4 (×), 1 (∘), and 0.87 (∆). The error bars in the simulation data are smaller than the size of the symbols. Figure 3 Evolution of the sixth cumulant a3(s) as a function of the average number of collisions per particle for (a) disks and (b) spheres. Symbols represent MD simulation results, while the lines correspond to the theoretical prediction (29). The values of the coefficient of restitution are (from bottom to top on the right side) α=0.1 (□), 0.4 (×), 0.87 (∆), and 1 (∘). The error bars in the simulation data are smaller than the size of the symbols, except in the stationary regime for α=0.1. Figure 4 Evolution of the fourth cumulant a2(s) for a coefficient of restitution α=0.1 as a function of the average number of collisions per particle for (a) disks and (b) spheres. Symbols represent MD simulation results, while the lines correspond to the theoretical prediction (26). Five different initial conditions are considered (see Appendix B): δ (⋄), M (∘), Γ (×), I (□), and S (∆). The error bars in the simulation data are smaller than the size of the symbols, except in the early stage for the initial condition S. Figure 5 Evolution of the sixth cumulant a3(s) for a coefficient of restitution α=0.1 as a function of the average number of collisions per particle for (a) disks and (b) spheres. Symbols represent MD simulation results, while the lines correspond to the theoretical prediction (29). Five different initial conditions are considered (see Appendix B): δ (⋄), M (∘), Γ (×), I (□), and S (∆). The error bars in the simulation data are smaller than the size of the symbols, except in the early stage for the initial condition S. Figure 6 Evolution of DKL(ϕ∥ϕM) (in logarithmic scale) as a function of the average number of collisions per particle for (a) disks and (b) spheres. Symbols represent MD simulation results, while the lines correspond to the theoretical approximation (39) (the thin dashed lines for the first stage of the evolution mean that it was necessary to take the real part). The values of the coefficient of restitution are (from top to bottom on the right side) α=0.1 (□), 0.4 (×), 0.87 (∆), and 1 (∘). The error bars in the simulation data are smaller than the size of the symbols, except when DKL(ϕ∥ϕM)≲10−4 for α=1. Figure 7 Evolution of DKL(ϕ∥ϕM) (in logarithmic scale) for a coefficient of restitution α=0.1 as a function of the average number of collisions per particle for hard (a) disks and (b) spheres. Symbols represent MD simulation results. Five different initial conditions are considered (see Appendix B): δ (⋄), M (∘), Γ (×), I (□), and S (∆). The error bars are smaller than the size of the symbols, except when DKL(ϕ∥ϕM)≲10−4 for the initial condition M. Figure 8 Evolution of DKL(ϕ∥ϕH) (in logarithmic scale) as a function of the average number of collisions per particle for (a) disks and (b) spheres. Symbols represent MD simulation results, while the lines correspond to the theoretical prediction (40) (the thin dashed lines for the first stage of the evolution meaning that it was necessary to take the real part). The values of the coefficient of restitution are α=0.1 (□), 0.4 (×), 0.87 (∆), and 1 (∘) The error bars in the simulation data are smaller than the size of the symbols, except when DKL(ϕ∥ϕM)≲10−4. Figure 9 Evolution of DKL(ϕ∥ϕH) (in logarithmic scale) for a coefficient of restitution α=0.1 as a function of the average number of collisions per particle for hard (a) disks and (b) spheres. Symbols represent MD simulation results. Five different initial conditions are considered (see Appendix B): δ (⋄), M (∘), Γ (×), I (□), and S (∆). The error bars are smaller than the size of the symbols, except when DKL(ϕ∥ϕM)≲10−4. Figure 10 Plot of DKL(ϕH∥ϕM) as a function of the coefficient of restitution α for disks (– –, ∘) and spheres (—, ∆). Symbols represent MD simulation results, while the lines correspond to the theoretical prediction provided by Equation (39) with a2(s)→a2H and a3(s)→a3H. The inset magnifies the region 0.6≤α≤1. The error bars in the simulation data are smaller than the size of the symbols. ==== Refs References 1. Shannon C.E. A mathematical theory of communication Bell Syst. Tech. J. 1948 27 379 423 10.1002/j.1538-7305.1948.tb01338.x 2. Gray R.M. Entropy and Information Theory 2nd ed. Springer New York, NY, USA 2011 3. Brey J.J. Santos A. Nonequilibrium entropy of a gas Phys. Rev. A 1992 45 8566 8572 10.1103/PhysRevA.45.8566 9906956 4. Kremer G.M. Thermodynamics and kinetic theory of granular materials Perspectives and Challenges in Statistical Physics and Complex Systems for the Next Decade Viswanathan G.M. Raposo E.P. da Luz M.G.E. World Scientific Singapore 2014 287 299 10.1142/9789814590143_0016 5. Chapman S. Cowling T.G. The Mathematical Theory of Non-Uniform Gases 3rd ed. Cambridge University Press Cambridge, UK 1970 6. Garzó V. Santos A. Kinetic Theory of Gases in Shear Flows: Nonlinear Transport Fundamental Theories of Physics Springer Dordrecht, The Netherlands 2003 7. Kullback S. Leibler R.A. On Information and Sufficiency Ann. Math. Statist. 1951 22 79 86 10.1214/aoms/1177729694 8. Kullback S. Information Theory and Statistics Dover New York, NY, USA 1978 9. Santos A. Kremer G.M. Relative Entropy of a Freely Cooling Granular Gas AIP Conf. Proc. 2012 1501 1044 1050 10.1063/1.4769657 10. Bettolo Marconi U.M. Puglisi A. Vulpiani A. About an H-theorem for systems with non-conservative interactions J. Stat. Mech. 2013 P08003 10.1088/1742-5468/2013/08/P08003 11. de Soria M.I.G. Maynar P. Mischler S. Mouhot C. Rey T. Trizac E. Towards an H-theorem for granular gases J. Stat. Mech. 2015 P11009 10.1088/1742-5468/2015/11/P11009 12. Plata C.A. Prados A. Global stability and H theorem in lattice models with nonconservative interactions Phys. Rev. E 2017 95 052121 10.1103/PhysRevE.95.052121 28618629 13. Bannerman M.N. Sargant R. Lue L. DynamO: A Free O (N ) General Event-Driven Molecular Dynamics Simulator J. Comput. Chem. 2011 32 3329 3338 10.1002/jcc.21915 21953566 14. Garzó V. Granular Gaseous Flows. A Kinetic Theory Approach to Granular Gaseous Flows Springer Nature Cham, Switzerland 2019 15. Brilliantov N.V. Pöschel T. Kinetic Theory of Granular Gases Oxford University Press Oxford, UK 2004 16. Brilliantov N. Pöschel T. Deviation from Maxwell distribution in granular gases with constant restitution coefficient Phys. Rev. E 2000 61 2809 2812 10.1103/PhysRevE.61.2809 17. van Noije T.P.C. Ernst M.H. Velocity distributions in homogeneous granular fluids: The free and the heated case Granul. Matter 1998 1 57 64 10.1007/s100350050009 18. Montanero J.M. Santos A. Computer simulation of uniformly heated granular fluids Granul. Matter 2000 2 53 64 10.1007/s100350050035 19. Santos A. Montanero J.M. The second and third Sonine coefficients of a freely cooling granular gas revisited Granul. Matter 2009 11 157 168 10.1007/s10035-009-0132-8 20. Brey J.J. Ruiz-Montero M.J. Cubero D. Homogeneous cooling state of a low-density granular flow Phys. Rev. E 1996 54 3664 10.1103/PhysRevE.54.3664 9965515 21. Ahmad S.R. Puri S. Velocity distributions in a freely evolving granular gas Europhys. Lett. 2006 75 56 62 10.1209/epl/i2006-10071-3 22. Ahmad S.R. Puri S. Velocity distributions and aging in a cooling granular gas Phys. Rev. E 2007 75 031302 10.1103/PhysRevE.75.031302 23. Yu P. Schröter M. Sperl M. Velocity Distribution of a Homogeneously Cooling Granular Gas Phys. Rev. Lett. 2020 124 208007 10.1103/PhysRevLett.124.208007 32501095 24. Bobylev A.V. Cercignani C. Toscani G. Proof of an asymptotic property of self-similar solutions of the Boltzmann equation for granular materials J. Stat. Phys. 2003 111 403 417 10.1023/A:1022273528296 25. Bisi M. Carrillo J.A. Toscani G. Decay Rates in Probability Metrics Towards Homogeneous Cooling States for the Inelastic Maxwell Model J. Stat. Phys. 2006 124 625 653 10.1007/s10955-006-9035-9 26. Bolley F. Carrillo J.A. Tanaka Theorem for Inelastic Maxwell Models Commun. Math. Phys. 2007 276 287 314 10.1007/s00220-007-0336-x 27. Carrillo J.A. Toscani G. Contractive probability metrics and asymptotic behavior of dissipative kinetic equations Riv. Mat. Univ. Parma 2007 6 75 198 28. Carlen E.A. Carrillo J.A. Carvalho M.C. Strong Convergence towards homogeneous cooling states for dissipative Maxwell models Ann. I. H. Poincaré 2009 26 167 1700 10.1016/j.anihpc.2008.10.005 29. Brilliantov N. Pöschel T. Breakdown of the Sonine expansion for the velocity distribution of granular gases Europhys. Lett. 2006 74 424 430 Erratum: 2006 , 75 , 188, doi:10.1209/epl/i2006-10099-3 10.1209/epl/i2005-10555-6 30. Noskowicz S.H. Bar-Lev O. Serero D. Goldhirsch I. Computer-aided kinetic theory and granular gases EPL 2007 79 60001 10.1209/0295-5075/79/60001 31. Brito R. Ernst M.H. Extension of Haff’s cooling law in granular flows Europhys. Lett. 1998 43 497 502 10.1209/epl/i1998-00388-9 32. Goldshtein A. Shapiro M. Mechanics of collisional motion of granular materials. Part 1. General hydrodynamic equations J. Fluid Mech. 1995 282 75 114 10.1017/S0022112095000048 33. Coppex F. Droz M. Piasecki J. Trizac E. On the first Sonine correction for granular gases Physica A 2003 329 114 126 10.1016/S0378-4371(03)00593-4 34. Abramowitz M. Stegun I.A. Handbook of Mathematical Functions Dover New York, NY, USA 1972 35. Bird G.A. Molecular Gas Dynamics and the Direct Simulation of Gas Flows Clarendon Oxford, UK 1994 36. Brey J.J. Dufty J.W. Kim C.S. Santos A. Hydrodynamics for granular flow at low density Phys. Rev. E 1998 58 4638 4653 10.1103/PhysRevE.58.4638 37. Boltzmann L. Lectures on Gas Theory Dover New York, NY, USA 1995 38. Maynar P. Trizac E. Entropy of Continuous Mixtures and the Measure Problem Phys. Rev. Lett. 2011 106 160603 10.1103/PhysRevLett.106.160603 21599346 39. Mischler S. Mouhot C. Rodriguez Ricard M. Cooling Process for Inelastic Boltzmann Equations for Hard Spheres, Part I: The Cauchy Problem J. Stat. Phys. 2006 124 655 702 10.1007/s10955-006-9096-9 40. Mischler S. Mouhot C. Cooling Process for Inelastic Boltzmann Equations for Hard Spheres, Part II: Self-Similar Solutions and Tail Behavior J. Stat. Phys. 2006 124 703 746 10.1007/s10955-006-9097-8 41. Mischler S. Mouhot C. Stability, Convergence to Self-Similarity and Elastic Limit for the Boltzmann Equation for Inelastic Hard Spheres Commun. Math. Phys. 2009 288 431 502 10.1007/s00220-009-0773-9 42. Pettersson R. On Solutions to the Linear Boltzmann Equation for Granular Gases Transp. Theory Stat. Phys. 2004 33 527 543 10.1081/TT-200053937 43. Esipov S.E. Pöschel T. The granular phase diagram J. Stat. Phys. 1997 86 1385 1395 10.1007/BF02183630 44. Kuninaka H. Hayakawa H. Anomalous Behavior of the Coefficient of Normal Restitution in Oblique Impact Phys. Rev. Lett. 2004 93 154301 10.1103/PhysRevLett.93.154301 15524884