Which reference state for a granular gas heated by the stochastic thermostat?

AIP Conference Proceedings 1501, 1009 (2012); https://doi.org/10.1063/1.4769652 1501, 1009 © 2012 American Institute of Physics. Which reference state for a granular gas heated by the stochastic thermostat? Cite as: AIP Conference Proceedings 1501, 1009 (2012); https://doi.org/10.1063/1.4769652 Published Online: 27 November 2012 P. Maynar, M. I. García de Soria, and E. Trizac Which Reference State for a Granular Gas Heated by the Stochastic Thermostat? P. Maynar∗,M.I.GarcíadeSoria ∗and E. Trizac† ∗Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, Spain †Laboratoire de Physique Théorique et Modèles Statistiques, UMR CNRS 8626, Université Paris-Sud, 91405 Orsay, France Abstract. We study the time evolution of the one-particle distribution function of a granular gas heated by the stochastic thermostat. It is found that, for a wide class of homogeneous initial conditions, memory of the initial conditions is rapidly erased and the system approaches the stationary state through a hydrodynamic state in which all the time dependence of the distribution function goes through the temperature. As a consequence, the simple one parameter scaling velocity distributions reported earlier in this context do not suffice, and a two parameter scaling is required. We find a very good agreement between the theoretical prediction and DSMC simulation results. Keywords: Granular gases, Boltzmann equation, driven systems PACS: 45.70.-n, 05.20.Dd, 05.60.-k, 51.10.+y INTRODUCTION A granular system may be defined as an ensemble of macroscopic particles which collide inelastically, so that part of the kinetic energy of the grains is dissipated in a collision. This simple ingredient gives rise to a rich phenomenology, which is of interest not only for practical reasons but also because of the theoretical challenges thereby brought to the fore [1, 2, 3]. Statistical mechanics and kinetic theory have proved to be powerful tools to study granular matter at a microscopic/mesoscopic level of description, while at a macroscopic level hydrodynamic equations have been widely used [2]. Due to the non conservativity of the interactions, which has shed doubts on the consistency of a hydrodynamicdescription [4, 5, 6], the granular temperature, defined as the second velocity moment of the one-particle distribution function, decays monotonically in time in an isolated granular system [7]. In the fast-flow regime, it has been shown numerically that for a wide class of initial conditions the system reaches a homogeneous state in which all the time dependence of the one-particle distribution function goes through the temperature. This is the so-called homogeneous cooling state (HCS) which has been widely studied in the literature [8, 9]. Actually, in a homogeneous situation, this state is reached very quickly (some collisions per particle). We have then two time scales: the kinetic or fast scale in which the scaling regime has not been reached yet and where the “microscopic” excitations evolve, and the hydrodynamic or slow scale in which the temperature changes [10]. Considering non homogeneous states, this separation of time scales renders possible a hydrodynamic description of the system in terms of the density, velocity and temperature field, and the HCS plays the role of the reference state when the Chapman-Enskog method is applied [11]. In this work we will discuss if this kind of behavior is also present in a driven granular system. As energy is injected, a stationary state is reached in which the energy supplied by the thermostat is compensated by the energy lost in collisions. For the hard particle model, one of the most used homogeneous heating is the so-called stochastic thermostat, which consists in a white noise force acting on each grain [12, 9, 13]. In the low density limit, the distribution function of the homogeneous state has been characterized [9]. Hydrodynamic equations have been derived via the Chapman-Enskog expansion [14] and fluctuating hydrodynamics have been successfully applied in order to understand the large scale structure found in the stationary state [15]. Let us remark that in all the studies concerning hydrodynamics for a system heated by the stochastic thermostat, the stationary state played the role of “reference” state as the HCS in the undriven case although, as stressed by Lutsko [16], this is far from being trivial (the zeroth order in the gradients distribution is not simply related to the stationary state distribution). The objective in this work is to analyze this point critically at the level of the Boltzmann equation. We will study, for arbitrary homogeneous initial conditions, the type of state the system reaches in a kinetic scale. We find that a universal state is reached in a kinetic scale (universal in the sense that it is independent of the initial conditions) but that depends on the quotient between 28th International Symposium on Rarefied Gas Dynamics 2012 AIP Conf. Proc. 1501, 1009-1016 (2012); doi: 10.1063/1.4769652 © 2012 American Institute of Physics 978-0-7354-1115-9/$30.00 1009 the instantaneous temperature and the stationary temperature. Although the majority of the results has been published elsewhere [17], here a more exhaustive simulation study is provided, and the differences with the special case of the Gaussian thermostat are discussed in more detail. The plan of the paper is as follows. In the next section, the model is introduced and the main results considering the stationary state are summarized. We then argue that the most general homogeneous “hydrodynamic” distribution may depend on an additional dimensionless parameter. The distribution is calculated in the first Sonine approximation finding a very good agreement with DSMC simulation results. Finally, we discuss in the last section the possible implications of the existence of this state, focusing particularly on the hydrodynamic description. THE MODEL The system considered is a dilute gas of Nsmooth inelastic hard particles of mass mand diameter σ which collide inelastically with a coefficient of normal restitution α independent of the relative velocity. If at time tthere is a binary encounter between particles iand j, with velocities Vi(t)and Vj(t)respectively, the postcollisional velocities V i(t) and Vj(t)are V i=Vi−1+ α 2( σ ·Vij) σ , Vj=Vj+1+ α 2( σ ·Vij) σ ,(1) where Vij ≡Vi−Vjis the relative velocity and  σ is the unit vector pointing from the center of particle jto the center of particle iat contact. Between collisions, the system is heated uniformly by adding a random velocity to the velocity of each particle independently with certain frequency and with a given probability distribution. Let us define the jump distribution, PΔt(Δv), as the probability density that a particle experiments a jump Δvin the time interval Δt, that will be considered symmetric in the velocity space. If the variance of this distribution is small compared to the velocity scale in which the one-particle distribution, f(r,v,t), varies then the evolution equation for this distribution is the Boltzmann-Fokker-Planck equation [9, 18]. For a homogeneous system, this equation reads ∂ ∂ tf(v1,t)= σ d−1dv2¯ T0(v1,v2)f(v1,t)f(v2,t)+ ξ 2 0 2 ∂ 2 ∂ v2 1 f(v1,t),(2) where dis the dimension, ξ 2 0=limΔt→01 dΔtdxx2PΔt(x),and ¯ T0is the binary collision operator ¯ T0(v1,v2)=d σ Θ(v12 · σ )(v12 · σ )( α −2b−1 σ −1).(3) Here we have introduced the operator b−1 σ which replaces the velocities v1and v2by the precollisional ones v∗ 1and v∗ 2 given by v∗ 1=v1−1+ α 2 α ( σ ·v12) σ , v∗ 2=v2+1+ α 2 α ( σ ·v12) σ .(4) Let us remark that, under the conditions noted above, the evolution equation does not depend on the details of the distribution PΔt, but only on its second moment, through the coefficient ξ 2 0. It is known numerically that for a wide class of initial conditions the system reaches a stationary state. Assuming that total momentum is zero, i.e. dvv f(v,0)=0, the state is characterized by an isotropic stationary distribution, fs(v). Let us define the scaled distribution function, χ s,by fs(v)= n vd s χ s(c),c=v vs ,(5) where nis the density, vs≡2Ts/mis the thermal velocity and Tsis the stationary temperature, d 2nTs= dv1 2mv2fs(v). As this distribution is very close to a Maxwellian distribution, an expansion in terms of Sonine polynomials [19] does make sense. In the so-called first Sonine approximation the function reads [9] χ s(c)≈ χ M(c)[1+as 2S2(c2)],(6) 1010 where χ Mis the Maxwellian distribution with unit temperature, S2(c2)=d(d+2) 8−d+2 2c2+1 2c4, is the second Sonine polynomial, and as 2is the kurtosis of the distribution. Within this approximation, the distribution function can be calculated, obtaining for the kurtosis [9] as 2( α )= 16(1− α )(1−2 α 2) 73+56d−24d α −105 α +30(1− α ) α 2,(7) with a stationary temperature Ts=mdΓ(d/2) ξ 2 0 2 π d−1 2(1− α 2)n σ d−12/3 .(8) THE REFERENCE STATE Now, let us consider an initial condition with a temperature that differs appreciably from the stationary temperature (we also assume that total momentum is zero). The system will reach the stationary state in a hydrodynamic scale. As said, the question is whether this evolution will take place in two steps: a first quick evolution, whose dynamics depends on the initial conditions, followed by some isotropic universal state, independent of the initial condition, through which the final stationary state is reached. A first guess to this universal state could be an isotropic state with a scaling similar to the one found for the HCS, i.e. f(v,t)= n vd 0(t) χ s(c),wherec≡v/v0(t),v0(t)≡2T(t)/mand T(t)is the instantaneous temperature, d 2nT(t)=dv1 2mv2f(v,t). Let us remark that this scaling property, as will be seen later, is actually true in the Gaussian thermostat, where the particles between collisions are accelerated by a force proportional to its own velocity [20]. Moreover, in the stochastic thermostat case there are also some studies where it was implicitly assumed and where it seems to be valid, at least close to the stationary state, see [15] and references there in. Nevertheless, when this scaling form is substituted in the Boltzmann equation, Eq. (2), we obtain ξ 2 0 2v3 0(t) ∂ 2 ∂ c2 1 χ s(c1)+ 1 v2 0(t) dv0(t) dt ∂ ∂ c1 ·[c1 χ s(c1)] = −n σ d−1dc2¯ T0(c1,c2) χ s(c1) χ s(c2).(9) Except in the steady state where v0is time independent, such a form is not consistent (the left hand side depends on time while the right hand side does not) and we can conclude that such a solution is not possible. Another possibility is to try with the most general normal solution in which all the time dependence of the distribution function goes through the instantaneous temperature, T(t). In our case, we can define a quantity with temperature dimensions apart from the instantaneous temperature. This can be done with the strength of the noise parameter, ξ 0, and for concreteness, we take the stationary temperature given above. Then, by dimensional analysis the distribution function is of the form f(v,t)= n vd 0(t) χ (c, β ),(10) wherewehaveintroducedthe dimensionless parameter β ≡vs v0(t)and we have redefined c≡v v0(t). The scaled distribution function, χ , is a universal isotropic function of the scaled velocity cand the dimensionless parameter β . If this state really exists and plays the role of reference state, we expect to have a dynamics divided in a first quick transient that depends on the initial conditions, and a second regime where the system has forgotten the initial condition and evolves to the stationary state in such a way that it only feels the “distance” to stationarity through the parameter β ,i.e. f(v,t|f0)−→ n vd 0(t) χ (c, β )−→ fs(v).(11) Let us remark that a similar two parameter scaling occurs in the uniform shear flow of granular gases [21, 22], where theroleof β is played by the dimensionless shear rate a∗=a n σ d−1v0(t). To see whether this state actually exists, we have performed DSMC simulations [23] of N=1000 hard disks (d=2) of unit mass and unit diameter that collide inelastically with the collision rule given by Eq. (1), specified by the coefficient of normal restitution, α . The thermostat is implemented by a distribution which is constant in a square of length Δ=0.1 centered in the origin, and a given Δtfor each set of parameters. The results have been averaged over 105trajectories. For a given value of the inelasticity, the objective is to solve the Boltzmann equation 1011 0 0.2 0.4 0.6 0.8 1 β -0.05 -0.04 -0.03 -0.02 -0.01 0 a2 0 0.2 0.4 0.6 0.8 1 β -0.02 -0.015 -0.01 -0.005 0 a3 FIGURE 1. Coefficients a2and a3from Monte Carlo simulations, for α =0.80. The initial conditions are Maxwellian (circles) and flat (squares) distributions. In both cases, the initial temperatures are T0/Ts=7.3,21.8,36.4,51. The dotted lines are the theoretical predictions for the stationary values (a2from Eq. (7) and a3from [24]) 0 0.2 0.4 0.6 0.8 1 β -0.02 -0.01 0 0.01 0.02 0.03 0.04 a2 0 0.2 0.4 0.6 0.8 1 β -0.06 -0.04 -0.02 0 0.02 a3 FIGURE 2. Same as Fig 1, but for α =0.6, T0/Ts=15.5,38.8,77.7 and two initial conditions, Maxwellian (circles) and asymmetric (squares) distributions. for different initial conditions and see if, after some time, all the time dependence of the distribution function goes through the dimensionless parameter β . As it is difficult to measure the complete distribution function with the desired accuracy, we have measured the cumulants of the scaled distribution, χ (c, β ). In terms of the velocity moments, vl≡1 ndvvlf(v,t), we have measured the kurtosis a2=d d+2 v4 v22−1,(12) which is proportional to the fourth cumulant of χ (c, β )and the quantity a3=−d2 (d+2)(d+4) v6 v2+3d d+2 v4 v22−2,(13) which can be viewed as the reduced sixth cumulant. If the scaling is correct, we expect that, after some transient, the cumulants collapse for different initial conditions as a function of β . In Fig 1 we have plotted a2and a3as a function of β for α =0.80. The initial velocity distributions are Maxwellian distributions with four different temperatures, T0, significantly above the steady state value, Ts, and a distribution in which all the velocities have the same probability density in a square centered on v=0(referred to as the “flat” case) and the same initial temperatures as in the Maxwellian case. All the quantities are measured every 250 collisions, so each four consecutive points in the figure 1012 correspond to a time span of one collision per particle. It can be seen that, after some transient, memory of the initial condition is forgotten, so that the stationary distribution ( β =1) is reached following a universal route. The same behavior can be seen in Fig 2, for α =0.60 and three different initial temperatures, again T0Ts.We have started either with a Maxwellian distribution, as above, or with an asymmetric distribution made up of three possible velocities with different probabilities f(vx,vy,t=0)=3 6 δ (vx+8D/3) δ (vy+8D/3)+ 2 6 δ (vx−4D/3) δ (vy−4D/3) +1 6 δ (vx−16D/3) δ (vy−16D/3),(14) where the parameter Dis chosen to match the initial desired temperature, that are the same as in the Maxwellian initial condition. Note that although the figures emphasize the initial transient, the initial condition is forgotten very quickly (3 or 4 collisions per particle maximum). It can also be seen that the sign of the steady state a2,as 2, significantly affects the pathway followed in the a2- β plane: while the Maxwellian data in Fig. 1 approach the scaling form from above (as 2<0), all data in Fig. 2 approach the scaling form from below (as 2>0). In principle, a more than one parameter scaling form is expected to hold for any driving mechanism in which a stationary state is reached in the long time limit. Nevertheless, we will see now that the Gaussian thermostat is an exceptional case. We have performed the same kind of simulations but heating the system by accelerating each particle with a force proportional to its own velocity, Fi=mw0Vi. The time evolution of the cumulants, a2and a3,in the β scale is represented in Fig. 3, for a system of N=1000 and α =0.80. It is seen that the curves collapse only once the stationary value for a2is reached, in agreement with a one-parameter scaling function (no dependence on β ). The results are for a Maxwellian initial distribution (circles) and an asymmetric distribution (squares) and three initial temperatures (T0/Ts=4.8,9.6,28.9). This is obtained for any value of the parameter w0. In fact, such a behavior could have been suspected due to the mapping between the HCS-Gaussian thermostat. Since in the former there is no stationary temperature, there is no room for a two parameter scaling function in both cases (even if one could think that it is possible in the Gaussian thermostat case). To see this point more clearly, let us see that the Boltzmann equation for the Gaussian thermostat admits a solution of the form f(v,t)= n vd 0(t) χ s(c). The homogeneous Boltzmann equation in this case is [26] ∂ ∂ tf(v1,t)= σ d−1dv2¯ T0(v1,v2)f(v1,t)f(v2,t)− ω 0 ∂ ∂ v1 ·[v1f(v1,t)].(15) Inserting the scaling form into Eq. (15) and taking into account the evolution equation for the temperature, we obtain the fully consistent equation n σ d−1dc2¯ T0(c1,c2) χ s(c1) χ s(c2)= ω 0 vs ∂ ∂ c1 ·[c1 χ s(c1)],(16) that coincides with the equation for the scaled distribution in the HCS, with the proper definition of the dimensionless cooling rate [9]. Focusing again on the stochastic thermostat problem, and in order to analyze the problem theoretically, we change variables in the Boltzmann equation, Eq. (2), from {t,v},to{ β ,c}. In these variables, the scaled distribution function fulfills [ μ ( β )− μ (1) β 3] ∂ ∂ c1 ·[c1 χ (c1, β )]+ β∂ ∂βχ (c1, β ) =dc2¯ T0(c1,c2) χ (c1, β ) χ (c2, β )+1 2 μ (1) β 3 ∂ 2 ∂ c2 1 χ (c1, β ),(17) where μ ( β )=−1 2ddc1dc2(c2 1+c2 2)¯ T0(c1,c2) χ (c1, β ) χ (c2, β ).(18) Let us note that in this case, in contrast with Eq. (9), the equation is fully consistent as it appears as a change of variables where β simply plays the role of time. Of course, it only makes sense if the change of variables is well defined, i.e. the temperature evolves monotonically in time. Nevertheless, to prove that for any “reasonable” initial 1013 0 0.2 0.4 0.6 0.8 1 β -0.06 -0.04 -0.02 0 a2 0 0.2 0.4 0.6 0.8 1 β -0.05 -0.025 0 0.025 a3 FIGURE 3. Coefficients a2and a3for a system with α =0.80 heated by the Gaussian thermostat. The initial condition is either Maxwellian (circles) or asymmetric (squares), and the initial temperatures are T0/Ts=4.8,9.6,28.9. The dotted lines are the theoretical predictions, and the dashed lines are the simulations values reported in [24] for hard disks under the Gaussian thermostat. condition, the system forgets the initial condition and reaches a universal state is a formidable task. For this reason we limit ourselves to the simplified problem of deriving an approximate expression for this universal distribution function. As in the stationary state, the distribution is always close to a Maxwellian distribution and it makes sense to remain in the first Sonine approximation χ (c, β )≈ χ M(c)[1+a2( β )S2(c2)],(19) where, the kurtosis a2has been defined in (12) and, by definition, we have dc χ (c, β )=1,dcc χ (c, β )=0,dcc2 χ (c, β )=d 2.(20) If the approximate distribution (19) is substituted in the Boltzmann equation for the scaled distribution, Eq. (17), and we take the fourth velocity moment, neglecting the nonlinear terms in a2, we obtain the following evolution equation for a21 4 β (1− β 3)d d β a2( β )=(1−B− β 3)a2( β )+Bas 2,(21) where Bis the following function of α B=73+8d(7−3 α )+15 α [2 α (1− α )−7] 16(1− α )(3+2d+2 α 2)+as 2[85+d(30 α −62)+3 α (10 α (1− α )−39)].(22) Eq. (21) is an inhomogeneous linear differential equation that can be integrated. As the equation has three singular points, β =0, β =1and β →∞, let us first solve it in the interval (0,1). In this case the general solution of the associated homogeneous equation is aH 2( β )=K(1− β 3)4 3B β 4(B−1)(23) and a particular solution can be calculated by variations of parameters ap 2( β )=as 21+1− β 3 B−12F1−1 3,1; 4B−1 3; β 3,(24) where 2F1is the hypergeometric function [25]. As B−1>0 for each α , we obtain that each solution of (21) tends to as 2in the limit β →1, but for K=0 the solution diverges when β →0. This enables us to identify the universal a2( β ) as the unique solution that is finite in all the interval (0,1) a2( β )=as 21+1− β 3 B−12F1−1 3,1; 4B−1 3; β 3,0< β <1.(25) 1014 00.5 11.5 2 β 0 0.05 0.1 0.15 0.2 0.25 a2 α=0.1 α=0.3 α=0.6 00.5 11.5 22.5 β -0.03 -0.025 -0.02 -0.015 -0.01 -0.005 0 a2 α=0.80 α=0.95 FIGURE 4. Kurtosis as a function of β , for systems with α =0.1, 0.3, and 0.6(left)and α =0.80, and 0.95 (right). The points are simulation results starting from a Maxwellian distribution and the solid lines are the theoretical predictions. For the values of T0/Ts, see the main text. Moreover, if the initial condition is such that the first Sonine approximation is a good approximation for all times, we expect that Eq. (21) accurately describes the time evolution of the kurtosis and the universal solution be reached as (1− β 3)4 3B/ β 4(B−1). An analogous analysis can be done for β >1. Following similar ideas, we identify the “universal” solution to be a2( β )=−4Bas 2 7 β 3(1−1/ β 3)4B 3 2F17 3,1+4B 3;10 3;1 β 3.(26) In order to compare the above theoretical predictions to the simulation data, attention should be paid to the fact that analytical computation of velocity moments or cumulants is plagued by nonlinear effects that have been discussed in the literature [26, 27]. This results in some error in the calculation of the steady state value as 2, and we can also expect Bto suffer from a similar inaccuracy, that may be of the order of 10 or 20%. To circumvent this drawback, we take as 2 appearing in (25) and (26) from the Monte Carlo simulations, and we adjust Bto match the measured function a2( β ). This has been done for several values of the inelasticity. In Fig. 4 (left) we have plotted the results for α =0.1, 0.3, and 0.6, with initial temperatures T0/Ts=0.01,104.7, T0/Ts=0.01,98.7, and T0/Ts=0.008,77.7, respectively, and in Fig. 4 (right) for α =0.80 and 0.95, with initial temperatures T0/Ts=0.11,51 and T0/Ts=0.11,44 respectively. The solid lines are the theoretical predictions, for β <1and β >1. Two values of the initial temperature are used to generate separately the β <1 branch (associated to large T0)andthe β >1 branch (obtained for small T0). Let us note that the agreement is excellent for the whole range of β considered and all the values of the inelasticity. We have checked that as 2and Bthereby obtained are close enough to our predictions. DISCUSSION We have studied the dynamics of a homogeneous system of hard inelastic particles heated homogeneously by the stochastic thermostat. We have restricted to low densities, for which the time evolution of the one-particle distribution function is given in terms of the homogeneous Boltzmann equation. We find that for general initial conditions, the system forgets the initial condition on a kinetic time scale, and then evolves towards the stationary state in such a way that all the time dependence of the distribution function goes through the temperature. This universal state plays, for heated systems, a similar role as the homogeneous cooling stage for the free cooling problem. The distribution function is calculated in the first Sonine approximation finding a good agreement with the DSMC simulation results. Interesting questions remain open: Is this state reached for other sorts of homogeneous heating? Does it appears at the two-particle, or even at the N-particle level? What happens for arbitrary densities? Do the hydrodynamic equations depend on the structure of this state? Of course the complete answers to all these questions require further studies but we can guess some of the responses with the intuition gained here. In principle the only necessary ingredient to have this kind of scaling is to reach a homogeneous stationary state, in order to define the dimensionless time, β . Then, in principle, the scaling might be extended to other driving mechanisms. In this sense, the Gaussian thermostat is a 1015 particular case where the scaling is simpler (it does not depend on β ) due to the mapping with the free cooling case. The two-parameter scaling may also be present at higher densities. In fact, molecular dynamics simulation results indicate this possibility, at least for moderate densities [28]. From the point of view of hydrodynamics, the existence of this state should definitively be exploited. The key point to derive hydrodynamic equations from the Boltzmann equation is to have two time scales, in the first transient the system evolves to a local version of the reference state followed by a hydrodynamic regime. In this respect, the role of the reference state is played by the β -state, which is identified as the universal state reached after a general homogeneous perturbation. ACKNOWLEDGMENTS This research was supported by the Ministerio de Educación y Ciencia (Spain) through Grant No. FIS2008-01339 (partially financed by FEDER funds). It is a pleasure to thank A. Santos and A. Astillero for valuable discussions. REFERENCES 1. J. W. Dufty, J. Phys.: Condens. Matter 12, A47-A56 (2000). 2. I. Goldhirsch, Annu. Rev. Fluid Mech. 35, 267-293 (2003). 3. I. S. Aranson and L. S. Tsimring, Rev. Mod. Phys. 78, 641-692 (2006). 4. L. P. Kadanoff, Rev. Mod. Phys. 71, 435-444 (1999). 5. I. Goldhirsh, “Granular Gases: Probing the Boundaries of Hydrodynamics” in Granular Gases, edited by T. Poeschel and S. Luding (Springer, Berlin, 2001), pp. 79-99. 6. J. W. Dufty, Advances in Complex Systems 4, 397-406 (2001); cond-mat/0109215. 7. P. K. Haff,J. Fluid. Mech. 134, 401-430 (1983). 8. A. Goldshtein, and M. Shapiro, J. Fluid. Mech. 282, 75-114 (1995). 9. T. P. C. van Noije, and M. H. Ernst, Granular Matter 1, 57-64 (1998). 10. J. J. Brey, J. W. Dufty, and M. J. Ruiz-Montero, “Linearized Boltzmann Equation and Hydrodynamics for Granular Gases” in Granular Gas Dynamics, edited by T. Poeschel and N. Brilliantov (Springer, Berlin, 2003), pp. 227-250. 11. J. J. Brey, J. W. Dufty, C. S. Kim, and A. Santos, Phys.Rev.E58, 4638-4653 (1998). 12. D. R. M. Williams, and F. C. MacKintosh, Phys. Rev. E 54, R9-R12 (1996). 13. A. Puglisi, V. Loreto, U. M. B. Marconi, and A. Vulpiani, Phys. Rev E 59, 5582-5595 (1999). 14. V. Garzó, and J. M. Montanero, Physica A 313, 336-356 (2002). 15. T. P. C. van Noije, M. H. Ernst, E. Trizac, and I. Pagonabarraga, Phys. Rev. E 59, 4326-4341 (1999). 16. J. F. Lutsko, Phys.Rev.E73, 021302-1-021302-19 (2006). 17. M. I. García de Soria, P. Maynar, and E. Trizac, Phys. Rev. E 85, 051301-1-051301-16 (2012). 18. N. G. van Kampen, Stochastic Processes in Physics and Chemestry (North-Holland, Amsterdam, 1992). 19. P. Résibois and M. de Leener, Classical Kinetic Theory of Fluids (John Wiley, New York, 1977). 20. This property can be understood from the mapping there exists between this thermostat and the HCS. See for example J. Lutsko, Phys.Rev.E63, 061211-1-061211-14 (2001), or J. J. Brey, M. J. Ruiz-Montero, and F. Moreno, Phys. Rev. E 69, 051303-1-051303-13 (2004) 21. A. Astillero, and A. Santos, Europhys. Lett. 78, 24002-p1-24002-p6 (2007). 22. A. Astillero, and A. Santos, Phys. Rev. E 85, 021302-1-021302-16 (2012). 23. G. Bird, Molecular Dynamics and the Direct Simulation of Gas Flows (Clarendon, Oxford, 1994). 24. A. Santos and J. M. Montanero, Granular Matter 11, 157-168 (2009). 25. A. F. Nikiforov and V. B. Uvarov, Special Functions of Mathematical Physics (Birkhäuser Verlag, Basel, 1988). 26. J. M. Montanero and A. Santos, Granular Matter 2, 53-64 (2000). 27. F. Coppex, M. Droz, J. Piasecki, and E. Trizac, Physica A 329, 114-126 (2003). 28. M. I. García de Soria, N. Khalil, P. Maynar, and E. Trizac (in preparation). 1016