Next Article in Journal
Generating Artificial Reverberation via Genetic Algorithms for Real-Time Applications
Next Article in Special Issue
Interpreting Social Accounting Matrix (SAM) as an Information Channel
Previous Article in Journal
Estimation of the Reliability of a Stress–Strength System from Poisson Half Logistic Distribution
Previous Article in Special Issue
Complexity as Causal Information Integration
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Kullback–Leibler Divergence of a Freely Cooling Granular Gas

1
Departamento de Física, Universidad de Extremadura, E-06006 Badajoz, Spain
2
Departamento de Física and Instituto de Computación Científica Avanzada (ICCAEx), Universidad de Extremadura, E-06006 Badajoz, Spain
*
Author to whom correspondence should be addressed.
Submission received: 7 October 2020 / Revised: 6 November 2020 / Accepted: 13 November 2020 / Published: 17 November 2020
(This article belongs to the Special Issue Entropy: The Scientific Tool of the 21st Century)

Abstract

:
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 α .

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), 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.

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]
t f ( v 1 ; t ) = n σ d 1 I [ v 1 | f , f ] n σ d 1 d v 2 + d σ ^ ( v 12 · σ ^ ) α 2 f ( v 1 ; t ) f ( v 2 ; t ) f ( v 1 ; t ) f ( v 2 ; t ) ,
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
v 1 = v 1 1 + α 2 α ( v 12 · σ ^ ) σ ^ , v 2 = v 2 + 1 + α 2 α ( v 12 · σ ^ ) σ ^
are precollisional velocities. Note that we have defined the VDF with the normalization condition d v f ( v ; t ) = 1 .
An important quantity is the granular temperature defined as
T ( t ) = m d v 2 , X ( v ) d v X ( v ) f ( v ; t ) .
Taking moments in Equation (1), one finds the cooling equation
t T ( t ) = ζ ( t ) T ( t ) ,
where the cooling rate is given by
ζ ( t ) = m n σ d 1 T ( t ) d d v v 2 I [ v | f , f ] = ( 1 α 2 ) m n σ d 1 T ( t ) π ( d 1 ) / 2 4 d Γ ( d + 3 2 ) v 12 3 ,
X ( v 1 , v 2 ) d v 1 d v 2 X ( v 1 , v 2 ) f ( v 1 ; t ) f ( v 2 ; t ) .
Let us introduce the thermal velocity v th ( t ) 2 T ( t ) / m , which allows us to define the rescaled VDF ϕ ( c ; s ) as
f ( v ; t ) = v th d ( t ) ϕ ( c ; s ) , c v v th ( t ) ,
where the variable s in ϕ ( c ; s ) is a scaled time defined by
s ( t ) = 1 2 0 t d t ν ( t ) , ν ( t ) κ n σ d 1 v th ( t ) , κ 2 π ( d 1 ) / 2 Γ ( d 2 ) .
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
κ 2 s ϕ ( c ; s ) + μ 2 ( s ) d c · c ϕ ( c ; s ) = I [ c | ϕ , ϕ ] , μ k ( s ) d c c k I [ c | ϕ , ϕ ] ,
where we have taken into account that ζ ( t ) / n σ d 1 v th ( t ) = 2 μ 2 ( s ) / d . The associated hierarchy of moment equations is
κ 2 s c k = F k ( s ) k μ 2 ( s ) d c k μ k ( s ) .
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 ; s ) ϕ H ( c ) , where ϕ H ( c ) satisfies the integrodifferential equation
μ 2 H d 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 ) / 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 ) 1 + 1 2 ζ H ( t 0 ) ( t t 0 ) 2 ,
where t 0 being an arbitrary time belonging to the HCS regime. Also in the HCS regime, μ 2 ( s ) μ 2 H = const and thus Equation (4) becomes s T H ( s ) = ( 4 / κ d ) μ 2 H T H ( s ) , whose solution is
T H ( s ) = T H ( s 0 ) e 4 μ 2 H ( s s 0 ) / κ 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 / 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,
ϕ ( c ; s ) = ϕ M ( c ) 1 + k = 2 a k ( s ) S k ( c 2 ) ,
where
S k ( x ) = L k ( d 2 1 ) ( x ) = j = 0 k ( 1 ) j Γ d 2 + k Γ d 2 + j ( k j ) ! j ! x j
where are Sonine (or generalized Laguerre) polynomials, which satisfy the orthogonalization condition
S k | S k d c ϕ M ( c ) S k ( c 2 ) S k ( c 2 ) = N k δ k , k , N k Γ d 2 + k Γ d 2 k ! .
In Equation (13), the Sonine coefficient a k ( s ) is the 2 k -th cumulant of the VDF at time s. According to Equation (15),
a k ( s ) = S k ( c 2 ) N k .
In particular, a 0 ( s ) = 1 , a 1 ( s ) = 0 , and
a 2 ( s ) = 4 d ( d + 2 ) c 4 1 , a 3 ( s ) = 1 + 3 a 2 8 d ( d + 2 ) ( d + 4 ) c 6 .

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 a k with k 4 in Equation (13), as well as the nonlinear terms a 2 2 , a 2 a 3 , and a 3 2 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
L 3 μ 2 = A 0 + A 2 a 2 + A 3 a 3 , L 3 μ 4 = B 0 + B 2 a 2 + B 3 a 3 , L 3 μ 6 = C 0 + C 2 a 2 + C 3 a 3 ,
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 2 H . From Equation (9), we have that F 4 H = 0 . Thus, the obvious approximation [17] consists of
L 2 F 4 H = 0 a 2 H , a = ( d + 2 ) A 0 B 0 B 2 ( d + 2 ) ( A 2 + A 0 ) = 16 ( 1 α ) ( 1 2 α 2 ) 9 + 24 d ( 41 8 d ) α + 30 ( 1 α ) α 2 ,
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 2 H [18,19,33]. In particular, one can start from the logarithmic time derivative of the fourth moment and then take
L 2 F 4 H c 4 H = 0 a 2 H , b = ( d + 2 ) A 0 B 0 B 2 B 0 ( d + 2 ) A 2 = 16 ( 1 α ) ( 1 2 α 2 ) 25 + 24 d ( 57 8 d ) α 2 ( 1 α ) α 2 .
Note that
a 2 H , a a 2 H , b = 1 + a 2 H , a = 1 1 a 2 H , b .
Both approximations ( a 2 H , a and a 2 H , b ) are practically indistinguishable in the region 0.6 α < 1 , but a 2 H , b is much more accurate than a 2 H , a for higher inelasticity [18,19].
Next, to estimate a 3 H , we start from the exact condition F 6 H = 0 and carry out either the linearization
L 3 F 6 H = 0 a 3 H , a = G a ( a 2 H ) C 0 3 4 ( d + 2 ) ( d + 4 ) A 0 + C 2 3 4 ( d + 2 ) ( d + 4 ) ( 3 A 0 + A 2 ) a 2 H 3 4 ( d + 2 ) ( d + 4 ) ( A 3 A 0 ) C 3
or, alternatively,
L 3 F 6 H c 6 H = 0 a 3 H , b = G b ( a 2 H ) C 0 3 4 ( d + 2 ) ( d + 4 ) A 0 + C 2 3 C 0 3 4 ( d + 2 ) ( d + 4 ) A 2 a 2 H 3 4 ( d + 2 ) ( d + 4 ) A 3 C 3 C 0 .
In Equations (22) and (23), a 3 H is expressed in terms of a 2 H . Using Equations (19) and (20), four possibilities in principle arise, namely
a 3 H , a a = G a ( a 2 H , a ) , a 3 H , a b = G a ( a 2 H , b ) , a 3 H , b a = G b ( a 2 H , a ) , a 3 H , b b = G b ( a 2 H , b ) .
Comparison with DSMC results shows that the best general estimates are provided by a 3 H , a a and a 3 H , a b . In what follows, we choose a 2 H , b for the fourth cumulant and, for the sake of consistency with that choice, we adopt a 3 H , a b for the sixth cumulant. To simplify the notation, we make a 2 H , b a 2 H and a 3 H , a b a 3 H .
Once the (approximate) HCS values a 2 H and a 3 H have been obtained, we turn our attention to the evolution equations of a 2 ( s ) and a 3 ( s ) . Approximating Equation (9) with k = 4 as κ 2 s ln c 4 = L 2 F 4 ( s ) / c 4 , one obtains
s a 2 ( s ) = K 2 1 + a 2 ( s ) a 2 ( s ) a 2 H , K 2 8 d ( d + 2 ) κ B 2 B 0 ( d + 2 ) A 2 .
Its solution is
a 2 ( s ) = a 2 H + 1 + a 2 H X 0 e γ s 1 , X 0 1 + a 2 ( 0 ) a 2 ( 0 ) a 2 H , γ 1 + a 2 H K 2 .
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
s a 3 ( s ) = 3 s a 2 ( s ) K 2 a 2 ( s ) a 2 H K 3 a 3 ( s ) a 3 H ,
where
K 2 16 d ( d + 2 ) ( d + 4 ) κ 3 4 ( d + 2 ) ( d + 4 ) ( A 2 + 3 A 0 ) C 2 ,
K 3 16 d ( d + 2 ) ( d + 4 ) κ 3 4 ( d + 2 ) ( d + 4 ) ( A 3 A 0 ) C 3 .
Taking into account Equation (26), the solution to Equation (27) is
a 3 ( s ) = a 3 H + Y 0 e K 3 s + 1 + a 2 H 3 X 0 e γ s 1 + K 2 K 3 + 3 2 F 1 1 , K 3 γ ; K 3 γ + 1 ; X 0 e γ s ,
Y 0 a 3 ( 0 ) a 3 H 1 + a 2 H 3 X 0 1 + K 2 K 3 + 3 2 F 1 1 , K 3 γ ; K 3 γ + 1 ; X 0 ,
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.

2.4. Comparison with MD Simulations

The approximate theoretical predictions for a 2 H and a 3 H 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 2 H and a 3 H 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.
Figure 1a,b show the α -dependence of a 2 H and a 3 H , 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 2 H , b (see Equation (20)) performs very well, even if the fourth cumulant is not small (e.g., a 2 H ∼0.2 at α = 0.1 ). The approximate sixth cumulant a 3 H , a b (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 2 H changes from negative to positive values at α 1 / 2 0.71 , a 3 H is always negative. Note that, for large inelasticity, the cumulants a 2 H and a 3 H 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 3 H | | a 2 H | 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 / 2 v 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 ) = 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 a 2 H 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.
In Figure 2 and Figure 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 α = 0.1 are { a 2 H , a 3 H } = { 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 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 2 H , and two more (labeled as Γ and S) with a 2 ( 0 ) > a 2 H . 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.

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
H ( t ) = d v f ( v ; t ) ln f ( 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 { 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
C σ ^ { v 1 , v 2 } = { v 1 , v 2 } , v 1 , 2 = v 1 , 2 1 + α 2 ( v 12 · σ ^ ) σ ^ .
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,
C σ ^ { v 1 , v 2 } = C σ ^ C σ ^ { v 1 , v 2 } = { v 1 , v 2 } , v 1 , 2 = v 1 , 2 1 + α 2 ( v 12 · σ ^ ) σ ^ .
Comparison with Equation (2) shows that
v 1 , 2 = v 1 , 2 ± 1 α 2 2 α ( v 12 · σ ^ ) σ ^ , v 12 · σ ^ = α 2 v 12 · σ ^ .
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),
H ( s ) = d v f ( v , t ) ln f ( v , t ) = H * ( s ) d 2 ln 2 T ( s ) m , H * ( s ) d c ϕ ( c , s ) ln ϕ ( c , s ) .
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 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
D KL ( f g ) = X d x f ( x ) ln f ( x ) g ( x ) ,
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
D KL ( f f ref ) = D KL ( ϕ ϕ ref ) = d c ϕ ( 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 ) = 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., D KL ( ϕ ϕ ref ) = D KL ( ϕ ˜ ϕ ˜ ref ) for any invertible transformation c c ˜ = w ( c ) .
If D KL ( ϕ ϕ 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 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 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
D KL ( ϕ ϕ M ) = H * ( s ) + d 2 1 + ln π ,
where H * ( s ) is defined in Equation (34). Thus, D KL ( ϕ ϕ M ) differs from H * ( s ) by a constant, so that s D KL ( ϕ ϕ M ) = s H * ( s ) .
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 lim s D KL ( ϕ ϕ M ) = D KL ( ϕ H ϕ M ) > 0 , it is impossible that s D KL ( ϕ ϕ M ) 0 for all s. Nevertheless, in principle, it might happen that s D KL ( ϕ ϕ M ) 0 for the class of initial conditions such that D KL ( ϕ ϕ M ) s = 0 D KL ( ϕ H ϕ M ) , while s D KL ( ϕ ϕ M ) 0 for the complementary class of initial conditions such that D KL ( ϕ ϕ M ) s = 0 D KL ( ϕ H ϕ M ) . 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
D KL ( ϕ ϕ M ) = d c ϕ M ( c ) 1 + k = 2 a k ( s ) S k ( c 2 ) ln 1 + k = 2 a k ( s ) S k ( c 2 ) .
Now, in the spirit of the truncation approximation of Section 2.3, we can write the approximate expression
D KL ( ϕ ϕ M ) d c ϕ M ( c ) 1 + a 2 ( s ) S 2 ( c 2 ) + a 3 ( s ) S 3 ( c 2 ) ln 1 + a 2 ( s ) S 2 ( c 2 ) + a 3 ( s ) S 3 ( c 2 ) ,
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 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 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 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 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.

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 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 Figure 2, Figure 3, and Figure 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.,
D KL ( ϕ ϕ H ) d c ϕ M ( c ) 1 + a 2 ( s ) S 2 ( c 2 ) + a 3 ( s ) S 3 ( c 2 ) ln 1 + a 2 ( s ) S 2 ( c 2 ) + a 3 ( s ) S 3 ( c 2 ) 1 + a 2 H S 2 ( c 2 ) + a 3 H S 3 ( c 2 ) ,
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 Figure 4, Figure 5, and Figure 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 ϵ ,
ϕ ( c ; s ) ϕ M ( c ) ln ϕ ( c ; s ) ϕ H ( c ) = ϵ k = 2 a k ( s ) a k H S k ( c 2 ) + ϵ 2 2 k , k = 2 a k ( s ) a k H a k ( s ) a k H S k ( c 2 ) S k ( c 2 ) + O ( ϵ 3 ) .
Next, taking into account the orthogonality condition (15), we get
D KL ( ϕ ϕ H ) = ϵ 2 2 k N k a k ( s ) a k H 2 + O ( ϵ 3 ) ,
s D KL ( ϕ ϕ H ) = ϵ 2 k N k a k ( s ) a k H s a k ( s ) + O ( ϵ 3 ) .
Interestingly, this approximation preserves the positive-definiteness of the KLD. Note also that, to order ϵ 2 , D KL ( ϕ ϕ H ) is symmetric under the exchange ϕ ϕ H , i.e., D KL ( ϕ ϕ H ) D KL ( ϕ H ϕ ) = O ( ϵ 3 ) . Finally, consistent with the derivation of Equations (20) and (25), we neglect the cumulants a k with k 3 and apply Equation (25) to obtain
D KL ( ϕ ϕ H ) d ( d + 2 ) 16 a 2 ( s ) a 2 H 2 ,
s D KL ( ϕ ϕ H ) d ( d + 2 ) 8 K 2 1 + a 2 ( s ) a 2 ( s ) a 2 H 2 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 s D KL ( ϕ ϕ 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 a k H 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 2 H and a 3 ( s ) a 3 H . 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 D KL ( ϕ H ϕ M ) d = 2 < D KL ( ϕ H ϕ M ) d = 3 for very large inelasticity to D KL ( ϕ H ϕ M ) d = 2 > D KL ( ϕ 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 D KL ( ϕ H ϕ M ) as a function of α agrees with a toy model analogous to that of Equation (43a), namely D KL ( ϕ H ϕ M ) d ( d + 2 ) 16 a 2 H 2 .

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., 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 2 H ) and sixth ( a 3 H ) 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 2 H 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 2 H and a 3 H 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 2 H ] 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.

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:
DSMCDirect simulation Monte Carlo
HCSHomogenous cooling state
KLDKullback–Leibler divergence
MDMolecular dynamics
VDFVelocity 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 / σ ) = ( 10 4 , 4 641.58 ) and ( 1.35 × 10 4 , 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 L c / σ 8.3 × 10 3 (disks) and 1.3 × 10 4 (spheres) for the most demanding case ( α = 0.1 ). This represents ratios L c / L 1.8 (disks) L c / L 32.8 (spheres). Those ratios generally increase with decreasing inelasticity. For instance, at α = 0.87 , one finds L c / 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,
Φ ( c ; s ) = c d 1 d c ^ ϕ ( c ; s ) = Ω d c d 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
D KL ( ϕ ϕ ref ) = 0 d c Φ ( 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.,
Φ ( c i ; s ) N i ( s ) N Δ c , Φ H ( c i ) N i H N Δ c , c i = i 1 2 Δ c , i = 1 , 2 , , M .
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 i H 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
Φ M ( c i ) π d / 2 Ω d Δ c c i Δ c / 2 c i + Δ c / 2 d c c d 1 e c 2 = e c i Δ c 2 2 e c i + Δ c 2 2 Δ c , ( d = 2 ) , erf c i + Δ c 2 erf c i Δ c 2 Δ c + 2 π c i Δ c 2 e c i Δ c 2 2 c i + Δ c 2 e c i + Δ c 2 2 Δ c , ( d = 3 ) ,
where erf ( x ) = 2 π 0 x d t e t 2 is the error function.
Next, the KLD (A2) with ϕ ref ( c ) = ϕ M ( c ) and with ϕ ref ( c ) = ϕ H ( c ) are approximated in the simulations by
D KL ( ϕ ϕ M ) i = 1 M N i ( s ) N ln N i ( s ) / N Δ c Φ M ( c i ) , D KL ( ϕ ϕ H ) i = 1 M N i ( s ) N ln N i ( s ) N i H ,
where Φ M ( c i ) is given by Equation (A4). Analogously,
D KL ( ϕ H ϕ M ) i = 1 M N i H N ln N i H / N Δ c Φ M ( c i ) .
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., D KL ( ϕ ϕ ref ) D KL ( ϕ ϕ ref ) D KL ( ϕ H ϕ M ) α = 1 .
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 Figure 2, Figure 3, Figure 6, and Figure 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
Φ δ ( c ) = δ c d / 2 ,
which will be labeled with the Greek letter δ . The second initial distribution is just a Maxwellian (label M), i.e.,
Φ M ( c ) = 2 Γ ( d 2 ) c d 1 e c 2 .
Next, we choose the gamma distribution (label Γ ) normalized to c 2 = d 2 , namely
Φ Γ ( c ) = 2 θ d 2 θ Γ ( d 2 θ ) c d / θ 1 e c 2 / θ ,
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 θ ) 8 , so that a 2 = 2 ( θ 1 ) d + 2 and a 3 = 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 a 2 and a 3 corresponding to, in order of increasing a 2 , the initial states δ , M, I, Γ , and S.
Table A1. Values of the fourth and sixth cumulants for the initial distributions δ , M, I, Γ , and S (see text).
Table A1. Values of the fourth and sixth cumulants for the initial distributions δ , M, I, Γ , and S (see text).
δ MI Γ S
a 2 ( 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 )
a 3 ( 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 s D KL ( ϕ ϕ ref )

The aim of this appendix is to derive a formal expression for s D KL ( ϕ ϕ 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]
J [ ψ ] d c ψ ( c ) I [ c 1 | ϕ , ϕ ] = 1 2 d c 1 d c 2 + d σ ^ ( c 12 · σ ^ ) ϕ ( c 1 ) ϕ ( c 2 ) ψ ( c 1 ) + ψ ( c 2 ) ψ ( c 1 ) ψ ( c 2 ) .
Next, we perform the change of variables { c 1 , c 2 , σ ^ } { c 1 , c 2 , σ ^ } and take into account that d c 1 d c 2 = α d c 1 d c 2 and c 12 · σ ^ = α c 12 · σ ^ to obtain
J [ ψ ] = α 2 2 d c 1 d c 2 + d σ ^ ( c 12 · σ ^ ) ϕ ( c 1 ) ϕ ( c 2 ) ψ ( c 1 ) + ψ ( c 2 ) ψ ( c 1 ) ψ ( c 2 ) = α 2 2 d c 1 d c 2 + d σ ^ ( c 12 · σ ^ ) ϕ ( c 1 ) ϕ ( c 2 ) ψ ( c 1 ) + ψ ( c 2 ) ψ ( c 1 ) ψ ( c 2 ) ,
where in the second equality we have just renamed { c 1 , c 2 , c 1 , c 2 } { c 1 , c 2 , c 1 , c 2 } . Taking the average between Equations (A10) and (A11), we arrive at
J [ ψ ] = 1 4 d c 1 d c 2 + d σ ^ ( c 12 · σ ^ ) { ϕ ( c 1 ) ϕ ( c 2 ) ψ ( c 1 ) + ψ ( c 2 ) ψ ( c 1 ) ψ ( c 2 ) ϕ ( c 1 ) ϕ ( c 2 ) α 2 ψ ( c 1 ) + ψ ( c 2 ) ψ ( c 1 ) ψ ( c 2 ) } .
Now, we start from the KLD defined by Equation (36) and use the Boltzmann equation (8) to get
κ 2 s D KL ( ϕ ϕ ref ) = J ln ϕ ϕ ref μ 2 d d c ln ϕ ( c ) ϕ ref ( c ) c · c ϕ ( c ) .
where we have taken into account that ϕ ref ( c ) and d c ϕ ( c ) = 1 are independent of time. Integration by parts of the second term on the right-hand side of Equation (A13) yields
κ 2 s D KL ( ϕ ϕ ref ) = J ln ϕ ϕ ref μ 2 1 + 1 d d c ϕ ( c ) c · c ln ϕ ref ( c ) .
Finally, making use of Equation (A12) with ψ ( c ) = ln [ ϕ ( c ) / ϕ ref ( c ) ] , we obtain
κ 2 s D KL ( ϕ ϕ ref ) = 1 4 d c 1 d c 2 + d σ ^ ( c 12 · σ ^ ) [ ϕ ( c 1 ) ϕ ( c 2 ) ln ϕ ( c 1 ) ϕ ( c 2 ) ϕ ref ( c 1 ) ϕ ref ( c 2 ) ϕ ( c 1 ) ϕ ( c 2 ) ϕ ref ( c 1 ) ϕ ref ( c 2 ) ϕ ( c 1 ) ϕ ( c 2 ) α 2 ln ϕ ( c 1 ) ϕ ( c 2 ) ϕ ref ( c 1 ) ϕ ref ( c 2 ) ϕ ( c 1 ) ϕ ( c 2 ) ϕ ref ( c 1 ) ϕ ref ( c 2 ) ] μ 2 d d c ϕ ( c ) c · c ln ϕ ref ( c ) ϕ M ( c ) ,
where we have taken into account that d c ϕ ( c ) c · c ln ϕ M ( c ) = 2 d c c 2 ϕ ( 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 ) = c 2 + const and, therefore, J [ ln ϕ M ] = μ 2 [see Equation (8)]. As a consequence,
κ 2 s D KL ( ϕ ϕ M ) = 1 4 d c 1 d c 2 + d σ ^ ( c 12 · σ ^ ) [ ϕ ( c 1 ) ϕ ( c 2 ) ln ϕ ( c 1 ) ϕ ( c 2 ) ϕ ( c 1 ) ϕ ( c 2 ) ϕ ( c 1 ) ϕ ( c 2 ) α 2 ln ϕ ( c 1 ) ϕ ( c 2 ) ϕ ( c 1 ) ϕ ( c 2 ) ] μ 2 .
In the special case of elastic collisions ( α = 1 ), one has μ 2 = 0 and c i = c i , so that the standard H-theorem is recovered, namely
κ 2 s D KL ( ϕ ϕ M ) α = 1 = 1 4 d c 1 d c 2 + d σ ^ ( c 12 · σ ^ ) ϕ ( c 1 ) ϕ ( c 2 ) ϕ ( c 1 ) ϕ ( c 2 ) ln ϕ ( c 1 ) ϕ ( c 2 ) ϕ ( c 1 ) ϕ ( c 2 ) 0 .

References

  1. Shannon, C.E. A mathematical theory of communication. Bell Syst. Tech. J. 1948, 27, 379–423. [Google Scholar] [CrossRef] [Green Version]
  2. Gray, R.M. Entropy and Information Theory, 2nd ed.; Springer: New York, NY, USA, 2011. [Google Scholar]
  3. Brey, J.J.; Santos, A. Nonequilibrium entropy of a gas. Phys. Rev. A 1992, 45, 8566–8572. [Google Scholar] [CrossRef] [PubMed]
  4. Kremer, G.M. Thermodynamics and kinetic theory of granular materials. In Perspectives and Challenges in Statistical Physics and Complex Systems for the Next Decade; Viswanathan, G.M., Raposo, E.P., da Luz, M.G.E., Eds.; World Scientific: Singapore, 2014; pp. 287–299. [Google Scholar] [CrossRef]
  5. Chapman, S.; Cowling, T.G. The Mathematical Theory of Non-Uniform Gases, 3rd ed.; Cambridge University Press: Cambridge, UK, 1970. [Google Scholar]
  6. Garzó, V.; Santos, A. Kinetic Theory of Gases in Shear Flows: Nonlinear Transport; Fundamental Theories of Physics; Springer: Dordrecht, The Netherlands, 2003. [Google Scholar]
  7. Kullback, S.; Leibler, R.A. On Information and Sufficiency. Ann. Math. Statist. 1951, 22, 79–86. [Google Scholar] [CrossRef]
  8. Kullback, S. Information Theory and Statistics; Dover: New York, NY, USA, 1978. [Google Scholar]
  9. Santos, A.; Kremer, G.M. Relative Entropy of a Freely Cooling Granular Gas. AIP Conf. Proc. 2012, 1501, 1044–1050. [Google Scholar] [CrossRef] [Green Version]
  10. Bettolo Marconi, U.M.; Puglisi, A.; Vulpiani, A. About an H-theorem for systems with non-conservative interactions. J. Stat. Mech. 2013, P08003. [Google Scholar] [CrossRef]
  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. [Google Scholar] [CrossRef] [Green Version]
  12. Plata, C.A.; Prados, A. Global stability and H theorem in lattice models with nonconservative interactions. Phys. Rev. E 2017, 95, 052121. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  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. [Google Scholar] [CrossRef] [Green Version]
  14. Garzó, V. Granular Gaseous Flows. A Kinetic Theory Approach to Granular Gaseous Flows; Springer Nature: Cham, Switzerland, 2019. [Google Scholar]
  15. Brilliantov, N.V.; Pöschel, T. Kinetic Theory of Granular Gases; Oxford University Press: Oxford, UK, 2004. [Google Scholar]
  16. Brilliantov, N.; Pöschel, T. Deviation from Maxwell distribution in granular gases with constant restitution coefficient. Phys. Rev. E 2000, 61, 2809–2812. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef]
  18. Montanero, J.M.; Santos, A. Computer simulation of uniformly heated granular fluids. Granul. Matter 2000, 2, 53–64. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef]
  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. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  21. Ahmad, S.R.; Puri, S. Velocity distributions in a freely evolving granular gas. Europhys. Lett. 2006, 75, 56–62. [Google Scholar] [CrossRef]
  22. Ahmad, S.R.; Puri, S. Velocity distributions and aging in a cooling granular gas. Phys. Rev. E 2007, 75, 031302. [Google Scholar] [CrossRef]
  23. Yu, P.; Schröter, M.; Sperl, M. Velocity Distribution of a Homogeneously Cooling Granular Gas. Phys. Rev. Lett. 2020, 124, 208007. [Google Scholar] [CrossRef]
  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. [Google Scholar] [CrossRef]
  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. [Google Scholar] [CrossRef]
  26. Bolley, F.; Carrillo, J.A. Tanaka Theorem for Inelastic Maxwell Models. Commun. Math. Phys. 2007, 276, 287–314. [Google Scholar] [CrossRef] [Green Version]
  27. Carrillo, J.A.; Toscani, G. Contractive probability metrics and asymptotic behavior of dissipative kinetic equations. Riv. Mat. Univ. Parma 2007, 6, 75–198. [Google Scholar]
  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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef] [Green Version]
  30. Noskowicz, S.H.; Bar-Lev, O.; Serero, D.; Goldhirsch, I. Computer-aided kinetic theory and granular gases. EPL 2007, 79, 60001. [Google Scholar] [CrossRef]
  31. Brito, R.; Ernst, M.H. Extension of Haff’s cooling law in granular flows. Europhys. Lett. 1998, 43, 497–502. [Google Scholar] [CrossRef] [Green Version]
  32. Goldshtein, A.; Shapiro, M. Mechanics of collisional motion of granular materials. Part 1. General hydrodynamic equations. J. Fluid Mech. 1995, 282, 75–114. [Google Scholar] [CrossRef]
  33. Coppex, F.; Droz, M.; Piasecki, J.; Trizac, E. On the first Sonine correction for granular gases. Physica A 2003, 329, 114–126. [Google Scholar] [CrossRef] [Green Version]
  34. Abramowitz, M.; Stegun, I.A. (Eds.) Handbook of Mathematical Functions; Dover: New York, NY, USA, 1972. [Google Scholar]
  35. Bird, G.A. Molecular Gas Dynamics and the Direct Simulation of Gas Flows; Clarendon: Oxford, UK, 1994. [Google Scholar]
  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. [Google Scholar] [CrossRef] [Green Version]
  37. Boltzmann, L. Lectures on Gas Theory; Dover: New York, NY, USA, 1995. [Google Scholar]
  38. Maynar, P.; Trizac, E. Entropy of Continuous Mixtures and the Measure Problem. Phys. Rev. Lett. 2011, 106, 160603. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef]
  42. Pettersson, R. On Solutions to the Linear Boltzmann Equation for Granular Gases. Transp. Theory Stat. Phys. 2004, 33, 527–543. [Google Scholar] [CrossRef]
  43. Esipov, S.E.; Pöschel, T. The granular phase diagram. J. Stat. Phys. 1997, 86, 1385–1395. [Google Scholar] [CrossRef] [Green Version]
  44. Kuninaka, H.; Hayakawa, H. Anomalous Behavior of the Coefficient of Normal Restitution in Oblique Impact. Phys. Rev. Lett. 2004, 93, 154301. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. Plot of (a) the HCS fourth cumulant a 2 H and (b) the HCS sixth cumulant a 3 H 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 a 2 H , b (see Equation (20)) and a 3 H , a b (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 1. Plot of (a) the HCS fourth cumulant a 2 H and (b) the HCS sixth cumulant a 3 H 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 a 2 H , b (see Equation (20)) and a 3 H , a b (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.
Entropy 22 01308 g001
Figure 2. Evolution of the fourth cumulant a 2 ( 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 2. Evolution of the fourth cumulant a 2 ( 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.
Entropy 22 01308 g002
Figure 3. Evolution of the sixth cumulant a 3 ( 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 3. Evolution of the sixth cumulant a 3 ( 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 .
Entropy 22 01308 g003
Figure 4. Evolution of the fourth cumulant a 2 ( 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 4. Evolution of the fourth cumulant a 2 ( 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.
Entropy 22 01308 g004
Figure 5. Evolution of the sixth cumulant a 3 ( 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 5. Evolution of the sixth cumulant a 3 ( 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.
Entropy 22 01308 g005
Figure 6. Evolution of D KL ( ϕ ϕ 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 D KL ( ϕ ϕ M ) 10 4 for α = 1 .
Figure 6. Evolution of D KL ( ϕ ϕ 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 D KL ( ϕ ϕ M ) 10 4 for α = 1 .
Entropy 22 01308 g006
Figure 7. Evolution of D KL ( ϕ ϕ 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 D KL ( ϕ ϕ M ) 10 4 for the initial condition M.
Figure 7. Evolution of D KL ( ϕ ϕ 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 D KL ( ϕ ϕ M ) 10 4 for the initial condition M.
Entropy 22 01308 g007
Figure 8. Evolution of D KL ( ϕ ϕ 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 D KL ( ϕ ϕ M ) 10 4 .
Figure 8. Evolution of D KL ( ϕ ϕ 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 D KL ( ϕ ϕ M ) 10 4 .
Entropy 22 01308 g008
Figure 9. Evolution of D KL ( ϕ ϕ 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 D KL ( ϕ ϕ M ) 10 4 .
Figure 9. Evolution of D KL ( ϕ ϕ 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 D KL ( ϕ ϕ M ) 10 4 .
Entropy 22 01308 g009
Figure 10. Plot of D KL ( ϕ 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 a 2 ( s ) a 2 H and a 3 ( s ) a 3 H . The inset magnifies the region 0.6 α 1 . The error bars in the simulation data are smaller than the size of the symbols.
Figure 10. Plot of D KL ( ϕ 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 a 2 ( s ) a 2 H and a 3 ( s ) a 3 H . The inset magnifies the region 0.6 α 1 . The error bars in the simulation data are smaller than the size of the symbols.
Entropy 22 01308 g010
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Megías, A.; Santos, A. Kullback–Leibler Divergence of a Freely Cooling Granular Gas. Entropy 2020, 22, 1308. https://0-doi-org.brum.beds.ac.uk/10.3390/e22111308

AMA Style

Megías A, Santos A. Kullback–Leibler Divergence of a Freely Cooling Granular Gas. Entropy. 2020; 22(11):1308. https://0-doi-org.brum.beds.ac.uk/10.3390/e22111308

Chicago/Turabian Style

Megías, Alberto, and Andrés Santos. 2020. "Kullback–Leibler Divergence of a Freely Cooling Granular Gas" Entropy 22, no. 11: 1308. https://0-doi-org.brum.beds.ac.uk/10.3390/e22111308

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop