Temperature distribution in driven granular mixtures does not depend on mechanism of energy dissipation

We study analytically and numerically the distribution of granular temperatures in granular mixtures for different dissipation mechanisms of inelastic inter-particle collisions. Both driven and force-free systems are analyzed. We demonstrate that the simplified model of a constant restitution coefficient fails to predict even qualitatively a granular temperature distribution in a homogeneous cooling state. At the same time we reveal for driven systems a stunning result – the distribution of temperatures in granular mixtures is universal. That is, it does not depend on a particular dissipation mechanism of inter-particles collisions, provided the size distributions of particles is steep enough. The results of the analytic theory are compared with simulation results obtained by the direct simulation Monte Carlo (DSMC). The agreement between the theory and simulations is perfect. The reported results may have important consequences for fundamental science as well as for numerous application, e.g. for the experimental modelling in a lab of natural processes.


INTRODUCTION
Mixtures of granular particles of different size are ubiquitous in nature and technology [1][2][3][4].The examples in nature range from pebbles and sands to dust on the Earth.Many extraterrestrial objects are comprised of granular mixtures: One can mention interstellar dust clouds, protoplanetary discs [5] and planetary rings.Dense Saturn's rings demonstrate a tremendous size polydispersity, with the particles size ranging from centimeters up to a few meters [9,10].Granular dust also covers the surface of the Moon [6], Mars [7] and possibly other planets and satellites.Industrial granular materials, besides of pebbles and sands in building industry, are represented by powders in chemical and cosmetic production, as well as salt, sugar and cereals in food industry.
Granular materials demonstrate very rich behavior -depending on the applied load they can be in solid liquid or gaseous phase [2].If the applied load is small, granular materials resist the external force and keep their shape and volume as solids.With increasing external load they start to flow like fluids.Such unusual properties of granular material stem from the dissipative nature of the inter-particle collisions, which are quantified by the so-called restitution coefficient ε, see e.g.[2,14]: Here v ′ ki = v ′ k − v ′ i and v ki = v k − v i are the relative velocities of particles of masses m k and m i after and before a collision, and e is a unit vector directed along the inter-center vector at the collision instant.The post-collision velocities v ′ k and v ′ i are related to the pre-collision velocities v k and v i as follows, e.g.[14]: Here m eff = m 1 m 2 / (m 1 + m 2 ) is the effective mass of colliding particles.The restitution coefficient 0 ≤ ε < 1 shows that the after-collisional relative velocity is smaller than the pre-collisional one, since the mechanical energy is transformed into the internal degrees of freedom of the particles.Due to a permanent loss of the kinetic energy of particles in the collisions, a steady supply of energy is required to keep the system in liquid or gaseous phase (unless the system is in a force-free state, where it undergoes a homogeneous cooling).The nature of the driving forces that fluidize granular matter may be very different.These may be gravitational forces, as in the case of astrophysical objects (dust clouds, protoplanetary disks and planetary rings) or avalanches in mountains [11].It may be wind of atmospheric gases, initiating the motion of sand grains, which results in dune formation on the Earth [8] or trigger dust storms on Mars [7].The fluidization of granular materials in industry may be caused by the vibration of a container, or by moving parts of a system, like e.g.blades or a piston.The transport properties of granular fluids crucially depend on the mean kinetic energy of the grains, which is also termed as "granular temperature".Due to dissipative nature of the inter-particles collisions, the energy equipartition, valid for equilibrium molecular systems, does not hold for granular mixtures, where each species has its own temperature.For a mixture of i = 1, 2, . . .N species the granular temperature of k-th species is defined as follows 3 2 Here m k is the mass of the according granular species, v k -its velocity, f (v k , t) -the velocity distribution function, which quantifies the number of particles in the system of the kind k with the velocity v k at time t and n k is the number density of the k-th species of the granular fluid, In what follows we consider granular fluids with a low density, which are termed as "granular gases".It is expected that the mixture behaves as a gas, when the total packing fraction of all components does not exceed about 20%.The violation of the energy equipartition in granular mixtures has been recognized almost two decades ago.It was predicted theoretically, confirmed in computer simulations [14][15][16][17][18] and observed experimentally [19,20].An impressive natural example of a granular mixture with the broken energy equipartition is Saturn rings, where all granular species demonstrate different temperatures [21][22][23][24].The polydispersity in the rings arises due to coagulation and fragmentation of granular particles [10,25].Although the effect of broken equipartition is known for a long time, still the physical laws that determine the distribution of granular temperatures in granular mixtures are not known.Such laws should predict the granular temperature for each species as a function of (i) size and mass distribution of granular particles, (ii) of the dissipative mechanism of the particles collisions and (iii) of the driving mechanism, applied to the system to keep it fluidized.Force-free granular mixtures can exist in a gaseous state in the regime of homogeneous cooling; here the temperature distribution should be determined by the items (i) and (ii) above.The temperature distribution in granular mixtures with the simplified model of a constant restitution coefficient, ε = const, which is equal for all inter-particle collisions, has been addressed in our previous study [13]; some universality of the granular temperature distribution was reported [13].Although the assumption of a constant restitution coefficient drastically simplifies the analysis and is widely used, see e.g.[26][27][28][29][30][31][32][33][34], it contradicts, the experimental observations [35][36][37], as well as basic mechanical laws [12,38].The latter indicate that ε does depend on the impact velocity [12,37,[39][40][41].This dependence may be obtained by solving the equations of motion for colliding particles with the explicit account for the dissipative forces acting between the grains.The simplest, but still rigorous, first-principle model of inelastic collisions takes into account the viscoelastic properties of particles' material.This results in the corresponding inter-particle forces [39,46] and eventually, in the restitution coefficient for viscoelastic particles [12,41,42]: where h k are numerical coefficients, κ and A characterize respectively the elastic and dissipative properties of the particles material (see Methods).Viscoelastic model agrees well with the experimental data when the impact velocity is not very large [37,[43][44][45].If the dissipative mechanism is caused by plastic deformation of particles, one obtains the following expression for the restitution coefficient [47]: where V yield is the yield velocity.Dissipative mechanism associated with plastic deformation corresponds to rather high impact velocities [44,45].A phenomenological exponential model for the velocity-dependent restitution coefficient has been also employed for the description of experimental data in Ref. [48], where δ is a dimensionless parameter, ρ is the density of the particle material and Y is the Young modulus.
It is well known that the velocity dependence of the restitution coefficient may drastically change the qualitative behavior of granular systems.For instance, it changes the cooling law in a homogeneous cooling state [12,52], the velocity distribution function [52,54], the diffusion of granular particles [49] and even pattern formation, which becomes a transient process for the case of velocity-dependent ε [50].Therefore, to formulate the laws for the temperature distribution in a granular mixture one needs to consider in detail the dissipation mechanism of particles collision.This is done in the present study.We analyze the distribution of granular temperatures in mixtures of granular gases for different dissipative mechanisms and driving models.We show that the simplified model of a constant restitution coefficient fails to predict even qualitatively the granular temperature distribution in a homogeneous cooling state.At the same time for driven granular systems we arrive at an astonishing result -the distribution of temperatures in granular mixtures is universal, that is, it does not depend on a particular dissipation mechanism of particles collisions.This conclusion holds true for steep distributions of particles size.The results of the analytic theory of the present study are compared with simulation results obtained by the direct simulation Monte Carlo (DSMC).The agreement between the theory and simulations is perfect.

A. Model
We consider a polydisperse granular system with discrete distribution of masses of particles.Let the smallest particle mass be m 1 = 0.01 and masses of other particles read, m k = km 1 , where k = 1, 2, ...N are integers and N is the total number of different species in the system.The system is spatially uniform and dilute enough so that only pairwise collisions take place in the system and multiple collisions of the particles may be safely neglected.The mass-velocity distribution function f k (v k , t) gives concentrations of particles of mass m k with the velocity v k at time t.Since the deviation of the velocity distribution function from the Maxwellian distribution is relatively small [14], we assume for simplicity that f k (v k , t) is Maxwellian.This function evolves according to the Boltzmann equation, applicable for granular gases, where correlations of velocities of colliding particles may be neglected [14], In Eq. ( 7) I coll k is the Boltzmann collision integral [14]: where σ ki = (σ k + σ i ) /2, with σ k = (6m k /(πρ)) 1/3 being the diameter of particles of mass m k , ρ is the mass density of the particle material.The summation is performed over all species in the system.v ′′ k and v ′′ i are pre-collision velocities in the so-called inverse collision, resulting in the post-collision velocities v k and v i .The Heaviside function Θ(−v ki • e) selects the approaching particles and the factor χ equals the product of the Jacobian of the transformation ) and the ratio of the lengths of the collision cylinders of the inverse and the direct collisions [14]: In the case of a constant restitution coefficient χ = 1/ε 2 for viscoelastic particles it has a more complicated form [14].
The second term I heat k describes the driving of the system.It quantifies the energy injection into a granular gas to compensate its losses in dissipative collisions; it is zero for a gas in a homogeneous cooling state (HCS).Here we consider a uniform heating -the case when the grains suffer small random uncorrelated kicks throughout the volume [55,56].To mimic the external driving forces a few types of thermostat have been proposed [56].For a thermostat with a Gaussian white noise, the heating term has the form [56,57]: Here the constant Γ k characterizes the strength of the driving force.It may vary for different species, depending on the type of driving [58].When all species are supplied with the same energy, we have a driving equipartition, Γ k = Γ 1 = const.In the case of the force controlled driving Γ k ∝ 1/m k , while in the case of the velocity controlled driving Γ k ∝ m k [58].In our study we analyze a more general case of a power-law dependence of Γ k on a particle mass, namely, Γ k = Γ 1 k γ .The driving may also depend on the local velocity of granular particles [59], however we neglect this effect in the present study.
Multiplying the Boltzmann equation ( 7) by m k v 2 k /2 for k = 1...N and performing the integration over v k , we get the following system of equations for evolution of the granular temperatures of species of different masses: In a homogeneous cooling state Γ k = 0 and the granular system permanently cools down.Driven granular systems, that is, systems with a thermostat rapidly settle into a non-equilibrium steady state and all granular temperatures attain after a short time, some constant values, so that dT k /dt = 0 and the above system (11) turns into a set of algebraic equations, For the constant restitution coefficient the cooling rate ξ ki , quantifying the decrease of granular temperature of species of mass m k due to collisions with species of mass m i is given by the following expression [13]: In Ref. [13] it was assumed that ε ki = ε.In the case of viscoelastic particles the cooling rates may be also computed and the result reads (see Methods for detail): with Γ (x) being the Gamma-function.
While it is possible to obtain the explicit expressions for the cooling coefficients ξ ki for the viscoelastic dissipative model, it is not the case for the elasto-plastic model.Therefore it is practical to exploit a notion of a "quasi-constant" restitution coefficient, which corresponds to the effective restitution coefficient averaged over all collisions; it depends on the current temperature of granular fluid [54], but not on the impact velocity.

B. Effective restitution coefficient of colliding particles in a granular mixture
The equations ( 42), ( 5) and ( 6) for the restitution coefficient ε give this quantity for a collision with a particular impact velocity.Since a wide range of the impact velocities is observed, one deals with a wide range of restitution coefficients.Naturally, this is much more complicated than to deal with a single number, of a simplified model of a constant restitution coefficient.Therefore for practical reasons it is worth to define a restitution coefficient ε ki , averaged over all possible collisions.Such effective "quasi-constant" restitution coefficient may be defined as follows [54], Here v n = (v ki • e) is the normal component of the relative velocity of particles and |v n | is its collisional average: The Heaviside function Θ (−v n ) selects the approaching particles.This quantity has been introduced and tested in [54] for a uniform one-component granular gas and demonstrated its adequacy.Here we generalize this concept for a mixture of granular particles of different masses m i and granular temperatures T i .Using again the Maxwellian approximation for the velocity distribution function, the collisional average yields for the effective coefficient of the viscoelastic particles (see Methods for detail): In contrast to the case of one-component granular gas [54], the effective restitution coefficient depends on the granular temperatures of both components T k and T i , but not on the impact velocity, since it characterizes the collisions in average.The "quasi-constant" restitution coefficient (18) may be now used in Eq. ( 13) for the cooling rates in place of the constant restitution coefficients; that is, the following substitute may be used ε ki → ε ki .Although Eq. ( 18) refers to viscoelastic particles, one apply the collisional averaging (16) for other dissipation models.
C. Homogeneous cooling state: The failure of simplified model of constant ε In order to calculate the granular temperatures of particles of different masses in the homogeneous cooling state we have solved numerically the system of differential equations (11) with zero heating rate, Γ k = 0 and cooling coefficients ξ ki (Eqs.14), corresponding to viscoelastic particles.We exploit the size distributions of particles, which is steep enough, n k ≃ k −θ with θ = 3.The evolution of granular temperatures T k (t) is shown in Fig. 1a.The behavior of FIG. 1. a) Evolution of granular temperatures T1 and T200 in the granular gas in a homogeneous cooling statefor the velocitydependent restitution coefficient, Eq. ( 42).The dotted line shows the asymptotics ∼ t −5/3 .Symbols show the results of the DSMC simulations.b) Dependence of the granular temperatures T k in the HCS on the reduced mass of the particle k = m k /m1 at different times.Solid lines correspond to viscoelastic particles (solution of system of equations (11) with ξ ik in the form Eq. ( 14) with Aκ 2/5 = 0.441.Symbols show the results of the DSMC simulations.The dotted line shows the slope of the temperature distribution for the case of a constant ε. a granular gas of viscoelastic particles is drastically different as compared to the granular gas of particles colliding with a constant restitution coefficient.While in the case of constant restitution coefficient the granular temperature decrease with the same rate, corresponding to the Haff's law [14], the evolution of granular temperatures of viscoelastic particles is rather complicated.The temperature of monomers T 1 of mass m 1 cools down according to t −5/3 (the generalized Haff's law [14]), while the temperature of more massive particles decrease slower at the initial state of cooling (retarded cooling) and faster at the later state of cooling (accelerated cooling) (Fig. 1a).The difference becomes more pronounced with increasing mass, as it is shown in Fig. 2a, where the evolution of the ratios of granular temperatures is shown.In the case of a constant restitution coefficient the ratios of all granular temperatures tend to the steady state [54], while for viscoelastic particles the ratio of granular temperatures, T k /T 1 , first grows, then reaches its maximal value and then decreases with time, tending to unity, that is, tending to the equipartition (Fig. 2a).The larger the mass m i of the particle, the larger the granular temperature and at the later time the maximum of T k /T 1 is achieved.The temperature distribution T k /T 1 evolves with time and changes its form, see Fig. 1b.It does not correspond to the distribution T k ∼ k 1.85 , observed for the system with a constant restitution coefficient [13].Hence the behavior of granular temperatures with a realistic restitution coefficient qualitatively differs from that predicted for a model with a constant ε.In other words the simplified model of a constant restitution coefficient, widely used in the scientific literature, fails to describe a complicated behavior of a granular mixture.At the same time, as it is follows from Fig. 2a, the application of the corresponding quasi-constant, temperature-dependent restitution coefficient (18) allows to model a granular mixture with an acceptable accuracy.Figs. 1 and 2 also demonstrate that the theoretical results obtained by the solution of the rate equations (11) are in a perfect agreement with the numerical simulations by the DSMC (see Methods for more detail).

D. Driven granular mixtures: Universality of temperature distribution for all dissipative mechanisms
To study the evolution of a driven granular mixture we solve the system of equations (11) with the mass-dependent heating rates Γ k = Γ 1 k γ .We start from the case of viscoelastic particles and use the cooling coefficients (14) and then the cooling coefficients (13) for the impact-velocity independent restitution coefficients ε ki , with the use of the effective quasi-constant coefficients, ε ki from Eq. ( 18) in the place of ε ki .The evolution of granular temperatures in a heated gas is shown at Fig. 2b.As one can see from the figure, all temperatures relax to the steady-state values.The steady state temperatures T k form a stationary distribution, that behaves for large k as a power law, see Fig. 3.  11) with ξ ki in the form of (14).Dashed lines show the evolution of T k /T1 of particles, colliding with an effective restitution coefficient (18) with Aκ 2/5 = 0.063.b) In a uniformly heated granular gas with Γ k = Γ1 = 1 (all species are supplied with the same energy) for Aκ 2/5 = 0.063.The notations are the same as in the panel (a).
To understand this behavior theoretically, we assume the power-law distribution for the steady state temperatures, provided k are not small.Approximating for N ≫ 1 the summation by integration in Eq. ( 12) with the cooling rates, Eqs. ( 14), we get: If the particle size-distribution n i = n i (i) is steep enough the main contribution to the integral in Eq. ( 20) comes from i ≪ k.Expanding the integrand in Eq. ( 20) with respect to (i/k) ≪ 1 and keeping only the leading terms in the expansion we arrive at (see Methods for more detail): Here we exclude α < 0, since it may yield for i ≪ k a negative sign for the factor in the square brackets of the integrand in Eq. (20).The result of the integration, however, should be positive, as it gives the cooling rate.For steep distributions n i one can approximate N in the upper limits of the integrals in Eq. ( 22) by the infinity,  13) with the restitution coefficient ε ki , Eq. ( 18)), while the full solution of the system is given by the thick lines (cooling rates ξ ki , Eq. ( 14)).Two solutions are practically indistinguishable.Symbols show the results of the DSMC simulations.The viscoelastic parameter is Aκ 2/5 = 0.063.b) For γ = 0 (equal distribution for the energy input for all species) and negative values of the exponent γ: γ = −0.33,−0.47, −1.
The notations are the same as for the panel (a).The viscoelastic parameter is Aκ 2/5 = 0.441.
where p = 1 for α > 1 and p = (α+ 1)/2 for 0 < α < 1 (see Eq. ( 22)), so that the sum in (12) does not (asymptotically, for N ≫ 1) depend on N .Taking into account that the exponents of k in the both sides of Eq. ( 12) must be equal, we finally arrive at: Surprisingly, this result exactly coincides with the one for the velocity-independent restitution coefficient ε ik = ε = const.of Ref. [13].
In Fig. 3 the theoretical predictions for the temperature distribution are compared with the results of the numerical solution of equations (11) with cooling rates of viscoelastic particles, with the cooling rates for the effective quasiconstant restitution coefficient, as well as with the DSMC results (see Methods for the application detail of the DSMC).The results of all three approaches perfectly agree with each other as well as with the theoretical result, Eq. (24).For instance T k ∼ k 11/9 for γ = 1 and T k ∼ k 17/9 for γ = 2, which corresponds to the case of γ ≥ 2/3 [see Eq. ( 24)].Similarly, T k ∼ k 5/6 for γ = 0.5, corresponding to −1/3 ≤ γ < 2/3, see Fig. 3a.All these temperature distributions are the same as for the case of a velocity-independent restitution coefficient, with the same driving coefficient γ.
Interestingly, for the equal distribution of the external energy supply for all species (γ = 0) the energy equipartition does not hold for particles of different sizes: T k ∼ k 1/3 which is confirmed both by the scaling prediction and Monte Carlo simulations (Fig. 3b).The reason is that the losses of kinetic energy in collisions of smaller particles is larger than of bigger ones.In order to compensate this effect the input of the external energy should be larger for smaller particles, which can be observed for negative values of γ.For the system of viscoelastic particles the equipartition of energy takes place for γ ≃ −0.47 (Fig. 3b), which is slightly different from the value γ ≃ −0.33, predicted by the scaling approach for granular systems with constant restitution coefficient.This resembles the mimicry effect found in the binary mixtures of granular particles [60,61].For negative values of γ with larger absolute values the equipartition breaks again, but the mass-dependence of the granular temperatures becomes inverse: The granular temperature of larger particles becomes smaller [58].This is indeed observed for the force controlled driving with γ = −1 (Fig. 3b).All these findings are confirmed by the DSMC results.
As we have demonstrated above, the viscoelastic collision model yields the same temperature distribution as the simplified collision model of a constant ε.The same is true for all dissipative mechanisms and may be formulated a general theorem: Theorem: Distribution of partial granular temperatures in a driven granular mixture does not depend on the dissipation mechanism of inelastic collisions provided the size distribution of particles is steep enough.
The proof of the theorem and exact formulation of the applicability conditions are given in the section Methods.To illustrate the application of the general theorem we consider the temperature distribution in granular mixtures with other mechanisms of the dissipative collisions -the elasto-plastic model, described by Eq. ( 5) [47] and the exponential model, described by Eq. ( 6) [48].The results of DSMC simulation for granular mixtures with three different models of the restitution coefficient are shown in Fig. 4. One can see from the figure that the distribution of granular temperatures demonstrates the same slope for all mechanisms of the dissipative collisions.Moreover, for the case of the first-principle restitution coefficients -for viscoelastic and elasto-plastic dissipative mechanisms, not only the slopes of the distributions, but the distributions themselves coincide.The discrepancy for small k between the temperature distribution of the latter two models and the phenomenological model from Ref. [48] stems possibly from the unphysically steep decay of ε with the impact velocity.

CONCLUSION
We have studied kinetic properties of size-polydisperse granular mixtures where granular particles suffer pairwise inelastic collisions.These are characterized by the restitution coefficients ε, quantifying the energy losses.We consider two main mechanisms of the collisional dissipation, associated with the viscoelastic and elasto-plastic behavior of the particles material.We used the first-principle expressions for restitution coefficients for these two mechanisms, which implies the dependence of ε on the impact velocity v imp , that is, ε = ε(v imp ).We also considered a phenomenological exponential model for the impact-velocity dependent restitution coefficient ε(v imp ) and a simplified model of a constant ε that does not depend on the impact velocity.We analyzed both cases of force free and driven systems.We derived a system of equations that describe the evolution of granular temperatures of different species in the mixture and the according cooling coefficients ξ ij for the case of viscoelastic particles.We have also introduced the notion of the effective impact-velocity independent restitution coefficient ε ij , which significantly simplifies the analysis and demonstrated its adequacy and efficiency.We solved numerically the equations for the granular temperatures T k , where k specifies the granular species and found the temperature distribution T k .The theoretical studies have been accompanied by the numerical modeling of the system with the direct simulation Monte Carlo (DSMC).We observed an excellent agreement between the theoretical predictions and the numerical results.Two main conclusions follow from our work: (i) The simplified model of a constant, velocity-independent restitution coefficient fails qualitatively to describe the evolution of a granular mixture and the granular temperature distribution in a homogeneous cooling state (HCS).(ii) Temperature distribution in driven granular mixtures settles to a universal power-law distribution T k ∼ k α with the exponent α that does not depend on the dissipation mechanism of inelastic collisions.Moreover, the distribution of temperatures, obtained for the two main dissipation mechanisms -of viscoelastic and plastic energy losses, demonstrate not only coincidence of the slopes but the overall coincidence.This kinetic law is applicable for granular mixtures with a steep size distribution and may be formulated as a general theorem with a precise formulation of the conditions.
The reported results imply serious consequences for the overall kinetic behavior of granular materials -the agitated granular mixtures of very different nature may behave similarly.This result is important for fundamental science, as it helps to understand the kinetic properties of granular mixtures and perform an adequate modelling, as well as for numerous practical applications.As an immediate practical consequence of our findings, is the possibility of experimental modelling (imitation) of natural processes in a lab: One needs to care only about the size distribution of particles and the driving law, but not about particles material.Moreover, even if in the course of system evolution, the mechanism of the dissipative collisions changes (e.g. from elasto-plastic to viscoelastic) it will not affect the temperature distribution.

Derivation detail for the cooling coefficients, their sum and effective restitution coefficient
Here we present the derivation detail for obtaining ξ ki , i ξ ki and ε ik .
Cooling coefficients Due to rather small deviation of the velocity distribution function from the Maxwellian, we assume for this function the Maxwellian form, that is, In order to compute the cooling rate we multiply the Boltzmann equation (Eq.7) by m k v 2 k /2 and integrate it with Using Eq. ( 8) for the collisional integral we write, Noticing that [14] χ we recast the first integral in the r.h.s. of ( 27) into the form Since the pre-collision velocities v ′′ k and v ′′ i are related to v k and v i in the same way as v k and v i to post-collision velocities v ′ k and v ′ i , this integral may be rewritten as which finally yields, where is the difference of energy of a particle of mass m k after and before a collision.Let us introduce the center of mass velocity Using the collision rules, Eqs.
(2), one can derive ∆E k : where the restitution coefficient ε ki is given by Eq. ( 42) with the following elastic constant κ, which is a function of the Young's modulus Y , Poisson ratio ν, monomer mass m 1 and the material density of particles ρ [12,14,39].The dissipative constant A quantifies the viscous properties of the particles' material [39,46]: where η 1 and η 2 are the viscosity coefficients.Let us introduce the variable The change of energy attains the form: and the whole integral reads: with (39) Integration that refers to the first term of ∆E k [see Eq. (37] yields zero.After the integration over b we get: And Eq. ( 39) can be now presented as a sum over the following type of integrals: . Collecting all terms together, one gets Eq. ( 14).
The sum of cooling rates To compute the sum i ξ ki in Eq. ( 11) we analyze the structure of the r.h.s. of Eq. ( 20) for i ≪ k, starting from the first term.For α ≥ 1 the leading term depends on i and k as k α/2−5/6 i and for α < 1 as k −1/3 i (α+1)/2 .The next terms with 2 ≤ n ≤ 20 scale as A n/2 i −n/6 k (α−1)n/20 for α ≥ 1 and as A n/2 i −(1−α)n/20−1/6 for α < 1.For small A [which is the condition of the validity of ( 42)] these terms may be neglected for n ≥ 2 as compared to the first term, which yields Eq. ( 22).

The rigorous prove of the universality of the temperature distribution for all dissipative mechanisms
With mild assumptions we can prove that temperature distribution tends to the scaling form k α , where α does not depend on a particular model of the effective restitution coefficient.Note that effective restitution coefficient exists for every non-negative, continuous and bounded ε(v) and can be found using Eq. ( 16), although one can possibly get different values for the effective restitution coefficient and the effective squared restitution coefficient, which we denote respectively as ε ki and ε 2 ki (Here we do not consider the case of very soft particles or nano-particles, which may possess a negative restitution coefficient [64,65]).We present a proof for the Maxwell distribution, but the same techniques can be used in other cases.For simplicity we assume that m 1 = 1 and g 2 (σ ki ) = 1 (the generalization is straightforward), and rewrite (13) as First we check the convergence of Therefore, if T i increases with i slower than i, we have the convergence if In the opposite case, when Let us impose even more restrictions by the condition which holds true, if Consider now a partial sum from i = i 0 to k, which satisfies, and impose the condition which is true if ∞ i=1 n i i < ∞, and T k = O(k α ) for α from (24).Combining conditions (47) and (49) we observe that which essentially means that the full series may be replaced by its first i 0 terms with any desired accuracy ǫ (i 0 certainly depends on ǫ).Hence for k ≫ i 0 we have the same asymptotics for T k , obtained for the incomplete sum, as the one obtained for the whole series.These differ only by the factor 1 + ǫ, converging to 1 as i 0 and k increase (hereinafter it is implied that ǫ may be taken arbitrarily small).Now we illustrate that for k ≫ i 0 the terms of ξ ki (44) do converge.We have The last condition means that the terms with the negative sign in (44) disappear.Then we are left with If ε ki is continuous and has a limit for k → ∞, then 1 + ε ki is also a a value (that depends only on i) between 1 and 2 up to the factor (1 + ǫ).Solving then the stationary equation ( 12) for T k , with ξ ki from (51) leads to with α from (24), which proves the asymptotics.Note that the special case of α = 1 is to be treated with an additional care.In this case the limit of ε ik for k → ∞ may depend on the limit of T k /k α ; Eq. ( 52) for T k , with ξ ki from (51), turns into an equation for an implicit function f (x, ǫ) with x = T k /k α .In this particular case one also needs to check that ξ ki , given by Eq. ( 51) is a monotonous function of x = T k /k α .Otherwise, multiple solutions of Eq. ( 52) for lim k→∞ T k k α may exist.In all the addressed above cases, the quantity ε ik decreases slower than 1−c| v|; hence the monotony is guaranteed, so that lim k→∞ T k k α = const > 0 even for α = 1.
Substituting T k ∼ k α into ( 46) and (48) we see that the following conditions are sufficient for the convergence for the obtained asymptotics: The time between collisions and the sizes of colliding particles are chosen based on matrix Λ ij .Namely, the time between collisions is taken from the exponential distribution ∆t = − ln(rand(0, 1])/ mmax i,j=1 Λ ij and the sizes i and j of the colliding particles are chosen with probability Λ ij .
Note that matrix Λ ij has rank 6 and can be rapidly recalculated using the low-rank technique (see e.g.[66]).
2. When the sizes i and j of the colliding particles are known, we can estimate the maximum value of the relative velocity of colliding particles according to Then we use the standard technique to pick random particles and accept the collision if where e is a random unit vector.The values of |v i | max can be updated in O(log N i ) steps by keeping all velocities in a binary heap structure for each size.
To reduce the statistical noise of the simulation data for temperatures we apply the running averaging as follows.Let at time t conv the temperature of monomers ceases to change monotonously.This indicates that the system has achieved its steady state and the stochastic noise becomes the primary source of errors.We then calculate the temperatures every N collisions in the time interval (t conv , 2t conv ) and take an average.

FIG. 2 .
FIG.2.Evolution of granular temperatures T k .a) In a homogeneous cooling state.Solid lines illustrate the direct solution of the system of equations (11) with ξ ki in the form of(14). Dashed lines show the evolution of T k /T1 of particles, colliding with an effective restitution coefficient(18) with Aκ 2/5 = 0.063.b) In a uniformly heated granular gas with Γ k = Γ1 = 1 (all species are supplied with the same energy) for Aκ 2/5 = 0.063.The notations are the same as in the panel (a).

FIG. 3 .
FIG.3.Dependence of the granular temperatures T k on the reduced mass of the particle k = m k /m1 for a heated granular gas with Γ k = k γ .a) For γ = 1 (velocity controlled driving), γ = 2 and γ = 0.5.The results of the solution of the system with the effective restitution coefficient are shown with thin lines (cooling rates ξ ki , Eq. (13) with the restitution coefficient ε ki , Eq. (18)), while the full solution of the system is given by the thick lines (cooling rates ξ ki , Eq. (14)).Two solutions are practically indistinguishable.Symbols show the results of the DSMC simulations.The viscoelastic parameter is Aκ 2/5 = 0.063.b) For γ = 0 (equal distribution for the energy input for all species) and negative values of the exponent γ: γ = −0.33,−0.47, −1.The notations are the same as for the panel (a).The viscoelastic parameter is Aκ 2/5 = 0.441.

9 FIG. 4 .
FIG.4.Dependence of the granular temperatures T k on the reduced mass of particles k = m k /m1 for a heated granular gas for different models of restitution coefficients: the viscoelastic model, elasto-plastic model with V yield = 1[47] and the exponential model with δ ρ Y = 1[48] obtained in the DSMC simulations.The dashed lines indicate the slopes.