Universal reference state in a driven homogeneous granular gas

# Universal reference state in a driven homogeneous granular gas

## Abstract

We study the dynamics of a homogeneous granular gas heated by a stochastic thermostat, in the low density limit. It is found that, before reaching the stationary regime, the system quickly “forgets” the initial condition and then evolves through a universal state that does not only depend on the dimensionless velocity, but also on the instantaneous temperature, suitably renormalized by its steady state value. We find excellent agreement between the theoretical predictions at Boltzmann equation level for the one-particle distribution function, and Direct Monte Carlo simulations. We conclude that at variance with the homogeneous cooling phenomenology, the velocity statistics should not be envisioned as a single-parameter, but as a two-parameter scaling form, keeping track of the distance to stationarity.

## I Introduction

A granular gas may be viewed as a collection of macroscopic particles undergoing dissipative collisions. This very ingredient –inelasticity– gives rise to a rich phenomenology (1); (2); (3); (4), the understanding of which requires statistical mechanics tools: Kinetic theory has proven to be powerful at a microscopic/mesoscopic level of description, while at a macroscopic scale, hydrodynamic equations have been derived (2) and put to the test. Yet, the relevance and consistency of a hydrodynamic framework, which is a central question, is still elusive (5); (6); (7).

Due to collisional dissipation, the granular temperature, defined as the variance of velocity fluctuations, decays monotonically in time in an isolated granular system (8). 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 is encoded in the temperature. This is the so-called homogeneous cooling state (HCS) which has been widely studied in the literature (9); (10). In such a homogeneous situation, the dynamics involves two time-scales: the kinetic or fast scale –a few collisions per particle– in which the scaling regime has not been reached yet and where the “microscopic” excitations relax, and the following hydrodynamic or slow scale in which the memory of initial conditions has been lost, and the velocity distribution evolves through the granular temperature (11). Considering non homogeneous states, this separation of time scales opens the possibility of a hydrodynamic or coarse-grained description in terms of the density, velocity and temperature field, and the HCS then plays the role of the reference state when the Chapmann-Enskog method is applied (12).

On the other hand, several studies and experiments in granular matter deal with stationary states, which are reached under the action of some energy driving, often realized by a moving boundary, or by an interstitial medium that acts as a thermostat, see e.g. (13); (14); (15). In all these cases, the energy injected compensates for the energy lost in collisions. From a theoretical point of view, a minimalistic approach is to consider the system as driven by some random energy source, which can be implemented in different ways (16). 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 (17); (10); (18); (19); (20); (21); (22); (23); (24); (25); (26); (27). In the low density limit, the distribution function of the homogeneous state has been characterized (10). Hydrodynamic equations have been derived via the Chapmann-Enskog expansion (21) and fluctuating hydrodynamics have been put forward in order to understand the large scale structure found in the stationary state (18), or fluctuations of global quantities (25). We stress that in all the studies pertaining to the hydrodynamics of a system heated by the stochastic thermostat (18); (21), the stationary state played the role of “reference” state, as the HCS happens to be in the undriven case. It was therefore assumed that a one parameter scaling holds for the velocity probability distribution. The objective in this work is to analyze this point critically in the homogeneous case at the level of the Boltzmann equation. We will study, for arbitrary initial conditions, the type of state the system evolves into in a kinetic time scale. Surprisingly, 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 the instantaneous temperature and its stationary value.

The outline is as follows. Section II opens with a definition of the model, and a summary of relevant previously known results. The key question addressed lies in the scaling form of the velocity distribution close to the steady state. Does it depend only on the suitably reduced velocity variable, as is the case in the HCS, or is another parameter relevant, that would encode the distance to stationarity ? We will argue in section II.1 that a single parameter scaling form is inconsistent. We shall then show in sections II.2 and II.3 that a consistent two parameter scaling form can be identified. Its properties will be characterized by complementary numerical and analytical tools. Conclusions and perspectives will finally be discussed in section III.

## Ii Scaling form of the velocity distribution function

The system of interest is a dilute gas of smooth inelastic hard particles of mass and diameter , which collide inelastically with a coefficient of normal restitution independent of the relative velocity (3). If at time there is a binary encounter between particles and , having velocities and respectively, the post-collisional velocities and are

 V′i = Vi−1+α2(^σ⋅Vij)^σ, V′j = Vj+1+α2(^σ⋅Vij)^σ, (1)

where is the relative velocity and is the unit vector pointing from the center of particle to the center of particle at contact. Between collisions, the system is heated uniformly by a white noise acting independently on each grain (10); (18); (20); (21); (22); (24); (25) so that the one-particle velocity distribution, , then obeys the Boltzmann-Fokker-Planck equation (10); (28). For a homogeneous system this equation reads

 ∂∂tf(v1,t)=σd−1∫dv2¯T0(v1,v2)f(v1,t)f(v2,t)+ξ202∂2∂v21f(v1,t), (2)

where is the dimension of space, measures the noise strength, and is the binary collision operator

 ¯T0(v1,v2)=∫d^σΘ(v12⋅^σ)(v12⋅^σ)(α−2b−1σ−1). (3)

Here we have introduced the operator which replaces the velocities and by the precollisional ones and given by

 v∗1 = v1−1+α2α(^σ⋅v12)^σ, v∗2 = v2+1+α2α(^σ⋅v12)^σ. (4)

### ii.1 One-parameter scaling or beyond ?

It is an observation from numerical simulations, that for a wide class of initial conditions the system reaches a stationary state (10); (18); (20). Assuming that total momentum is zero, i.e. , the state is characterized by an isotropic stationary distribution, . Let us define the scaled distribution function, , by

 fs(v)=nvdsχs(c),c=vvs, (5)

where is the density, is the thermal velocity and is the stationary temperature, . As is rather close to a Maxwellian distribution, a reasonable strategy is to perform an expansion in Sonine polynomials (29). In the so-called first Sonine approximation, the steady state function then reads (10)

 χs(c)≈χM(c)[1+as2S2(c2)], (6)

where is the Maxwellian distribution with unit temperature, , is the second Sonine polynomial, and is the kurtosis of the distribution. Within this approximation, the distribution function can be calculated, with the result (10)

 as2(α)=16(1−α)(1−2α2)73+56d−24dα−105α+30(1−α)α2, (7)

and a stationary temperature

 Ts=m⎡⎢⎣dΓ(d/2)ξ2o2πd−12(1−α2)nσd−1⎤⎥⎦2/3. (8)

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, its precise value being immaterial). It is clear that the system will reach the stationary state in a hydrodynamic scale. The ensuing question is two-pronged. First of all, is the dynamics compatible with a universal scaling form –once the memory of initial condition is lost– that would provide a consistent solution to the Boltzmann equation, or is memory only washed out strictly speaking at the steady state point? Second, assuming such a scaling regime exists in some vicinity of the steady state, what is the minimal number of parameters required for its description? By analogy with unforced (HCS) phenomenology, a single parameter scaling might be anticipated:

 Extra open brace or missing close brace (9)

is defined from the instantaneous temperature . We note that this scaling property holds for the Gaussian thermostat as well (19), where the particles are accelerated between collisions by a force proportional to its own velocity (30). Moreover, in the stochastic thermostat case, the one parameter scaling (9) was implicitly assumed, and it seemed to be yield reasonable predictions at least close to the stationary state, see (18); (21). Nevertheless, when the form (9) is inserted in the Boltzmann equation, Eq. (2), we obtain

 ξ202v30(t)∂2∂c21χs(c1)+1v20(t)dv0(t)dt∂∂c1⋅[c1χs(c1)]=−nσd−1∫dc2¯T0(c1,c2)χs(c1)χs(c2), (10)

which is inconsistent: The left hand side depends on time while the right hand side does not. We conclude that such a solution should be ruled out, except when stationarity is reached and . The hope is to capture the post-kinetic time dependence of , through a more involved functional form, that would be free of the above inconsistency.

As is customary, we shall seek for a normal solution (29), which in the present homogeneous case means that should only depend on time via the instantaneous granular temperature . In conjunction with dimensional analysis, this leads to a function that should only depend on and , which we write as

 f(v,t)=nvd0(t)χ(c,β),withβ≡vsv0(t) (11)

and again . Note that we have assumed isotropy, depending on and not on the full vector , but this assumption can be easily relaxed. Note also that equivalent expressions can of course be chosen, such as .

If a state such as (11) holds, it represent a strong constraint on the form of the velocity distribution. The corresponding dynamics can be partitioned in a first rapid stage –that we do not attempt to describe– where initial conditions matter, and a subsequent universal relaxation towards stationarity, where only the distance to the steady state is relevant, through the dimensionless inverse typical velocity .

To put the above scenario to the test, we have performed Direct Monte Carlo Simulations (DSMC) (31) of hard disks () of unit mass and unit diameter, that collide inelastically with the collision rule given by equation (II). The thermostat is implemented following previous investigations (18) and the results have been averaged over trajectories. For a given value of the inelasticity, we thus solve the time-dependent Boltzmann equation for different initial conditions and analyze whether, after some kinetic transient, 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 worked with the cumulants of the scaled distribution, . In terms of the velocity moments, , we have measured the kurtosis of the distribution

 a2=dd+2⟨v4⟩⟨v2⟩2−1, (12)

which is proportional to the fourth cumulant of and the quantity

 a3=−d2(d+2)(d+4)⟨v6⟩⟨v2⟩3+3dd+2⟨v4⟩⟨v2⟩2−2, (13)

which can be viewed as the reduced sixth cumulant. If our scaling is correct, we expect that the cumulants quickly collapse for different initial conditions, as a function of . In Fig. 1, we have plotted and versus for . The initial conditions are either Maxwellian distributions with three different temperatures , significantly above the steady state value , or asymmetric distributions made up of three possible velocities with different probabilities

 f(vx,vy,t=0) = 36δ(vx+8D/3)δ(vy+8D/3)+26δ(vx−4D/3)δ(vy−4D/3) (14) + 16δ(vx−16D/3)δ(vy−16D/3).

Here, the parameter is chosen to match the initial desired temperature (chosen the same as in the Gaussian initial condition). All the quantities are measured every collisions, so that each four consecutive points in figure 1 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 () is reached following a universal route. In figure 1, those data points associated to the Gaussian initial distribution approach the scaling curve from above (circles) while those for the initial asymmetric case (14) approach the scaling curve from below (squares). A very similar behavior can be seen in Fig 2 for and four different initial temperatures, again such that , which ensures that (we have also probed the regime obtained with , where similar conclusions hold, see e.g. Fig. 4 below). We have started either with a Maxwellian distribution, as above, or with a distribution in which all the velocities have the same probability density in a square centered on (referred to as the “flat” case). We emphasize that the initial transient is fast: memory of the initial condition is lost after at most or collisions per particle, a phenomenon that cannot be appreciated from the figures.

Borrowing ideas from the extended self-similarity technique (34), we put to the test the possibility of an enhanced universality, by plotting as a function of , see Fig. 3. In doing so, it appears that the universal part of the versus curve seen in Fig. 1 is not enhanced by the reparametrization : different initial conditions do not lead to a data collapse, beyond the interval that was already evidenced in Fig. 1. However, a given functional form (say Maxwellian) leads to a unique path in the - plane, which is already a non trivial point, and furthermore, a plot like Fig. 3 leads to a significantly reduced scatter of points that Fig. 1. It is therefore more amenable to chart out the universal regime sought for. In these figures, the first measure reported after the dynamics has acted on the initial conditions is for a time of collisions per particle for the Maxwellian distribution and of around collisions per particle for the flat and asymmetric distributions. The present results establish numerically the existence of a universal non-trivial scaling regime, for which we now seek analytical characterization.

### ii.3 Analytical approach

For the sake of analytical progress, it is convenient to change variables in the Boltzmann equation (2), from the set , to . In these variables, the scaled distribution function fulfills

 [μ(β)−μ(1)β3]{∂∂c1⋅[c1χ(c1,β)]+β∂∂βχ(c1,β)} =∫dc2¯T0(c1,c2)χ(c1,β)χ(c2,β)+12μ(1)β3∂2∂c21χ(c1,β), (15)

where

 μ(β)=−12d∫dc1∫dc2(c21+c22)¯T0(c1,c2)χ(c1,β)χ(c2,β). (16)

We note that here and in contrast with Eq. (10), the equation is fully consistent as it appears as a change of variables where simply plays the role of time. Nevertheless, proving that for any “reasonable” initial 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 distribution function. As in the stationary state, the distribution will be worked out in the first Sonine approximation

 χ(c,β)≈χM(c)[1+a2(β)S2(c2)], (17)

where the kurtosis, , has been defined in (12) and, by definition, we have

 ∫dcχ(c,β)=1,∫dccχ(c,β)=0,∫dcc2χ(c,β)=d2. (18)

In expansion (17), we neglect contributions in and higher order. This is justified as long as the inelasticity is not too strong, and is backed up here by the fact that , as can be seen in Figs. 1 and 2. Of course, inclusion of higher order terms in (17) would improve the accuracy of the subsequent calculation.

Inserting (17) into the Boltzmann equation (II.3), taking the fourth velocity moment while neglecting nonlinear terms in , we obtain the following evolution equation for the cumulant

 14β(1−β3)ddβa2(β)=(1−B−β3)a2(β)+Bas2, (19)

where the parameter depends on dissipation and space dimension according to

 B=73+8d(7−3α)+15α[2α(1−α)−7]16(1−α)(3+2d+2α2)+as2[85+d(30α−62)+3α(10α(1−α)−39)]. (20)

Eq. (19) is an inhomogeneous linear differential equation that can be integrated. It exhibits two singular points at , , and we start with the interval . In this case the general solution of the associated homogeneous equation reads

 aH2(β)=K(1−β3)43Bβ4(B−1) (21)

and a particular solution can be obtained by variations of parameters. The general solution will then be the sum of these two contributions. The ensuing depends on the initial conditions through , and since our purpose here is to extract the universal behaviour of as a function of , we note that the contribution (21) fades rapidly as approaches unity (as where can be large ; note that it diverges in the elastic limit ). The universal behaviour is consequently encoded is the particular solution, and we finally have

 a2(β)=as2[1+1−β3B−12F1(−13,1;4B−13;β3)],0<β<1. (22)

where is the hyper-geometric function (32). This expression is well behaved in all the interval . An analogous analysis can be performed for . Following similar lines, we identify the universal solution to be

 a2(β)=−4Bas27β3(1−1/β3)4B32F1(73,1+4B3;103;1β3). (23)

Clearly, the same technical procedure can be applied to the higher order cumulants. For the sake of simplicity, we restrict to the function , that we wish to test against simulation results.

In order to compare the above theoretical predictions to the simulation data, attention should be payed to the fact that analytical computation of velocity moments or cumulants is plagued by non-linear effects that have been discussed in the literature (19); (33). This results in some error in the calculation of the steady state value , and we can also expect to suffer from a similar inaccuracy, that may be of the order of 10 or 20%. To circumvent this (somewhat minor) drawback, we take appearing in Eqs. (22) and (23) from the Monte Carlo simulations, and we adjust to match the measured function . This procedure provides us with Fig. 4, where the agreement between the functional forms (22) and (23) with Monte Carlo is excellent. It is of course important to check a posteriori that and thereby obtained are close enough to our predictions. The precise values are reported below:

 α=0.95α=0.8as2 from DSMC −0.0171−0.0150as2 from Eq. (???) −0.0157−0.0135B from DSMC 5.172.32B from Eq. (???) 6.432.61 (24)

## Iii Conclusion and perspectives

To summarize, we have studied the dynamics of a system of inelastic hard spherical grains, heated homogeneously by a stochastic thermostat. We have restricted the analysis to low density systems, amenable to a Boltzmann equation treatment. We have found that generically, after a kinetic transient, the system evolves into a scaling solution that no longer depends on initial conditions, before the steady state is finally reached. The relevant scaling form is not of the single parameter family as is the case in the homogeneous cooling state, but involves 2 parameters. The velocity distribution function was calculated in the so-called first Sonine approximation, which provides a very good agreement with the Monte Carlo simulations.

At this point several questions arise: What is the counterpart of the universal state brought to the fore at two-particle level, or even -particle ? What is the effect of density, and of a change in the driving mechanism ? Do the hydrodynamic relations (worked out say at Chapman Enskog level) depend on the structure of this state? The complete answer to all these interrogations requires further studies, but we expect that similar scaling forms should occur at higher densities, and with different thermostats as long as a steady state can be reached. In this respect, the Gaussian thermostat is presumably singular, since it can be mapped onto the free cooling case. The questions pertaining to hydrodynamics are more subtle; How the universal -scaling behaviour discussed in the present work impinges, as a reference state, on transport properties, should be explored.

## Iv Acknowledgments

This research was supported by the Ministerio de Educación y Ciencia (Spain) through Grant No. FIS2008-01339 (partially financed by FEDER funds).

### References

1. J. W. Dufty, J. Phys.: Condens. Matter 12, A47 (2000).
2. I. Goldhirsch, Annu. Rev. Fluid Mech. 35, 267 (2003).
3. N. Brilliantov and T. Pöschel, Kinetic Theory of Granular Gases (Clarendon Press, Oxford, 2004).
4. I. S. Aranson and L. S. Tsimring, Rev. Mod. Phys. 78, 641 (2006).
5. L. P. Kadanoff, Rev. Mod. Phys. 71, 435 (1999).
6. I. Goldhirsh, in Granular Gases, edited by T. Poeschel and S. Luding (Springer, Berlin, 2001).
7. J. W. Dufty, Advances in Complex Systems 4, 397-406 (2001); cond-mat/0109215.
8. P. K. Haff, J. Fluid. Mech. 134, 401 (1983).
9. A. Goldshtein and M. Shapiro, J. Fluid. Mech. 282, 75 (1995).
10. T. P. C. van Noije, and M. H. Ernst, Granular Matter 1, 57 (1998).
11. J. J. Brey, J. W. Dufty, and M. J. Ruiz-Montero, in Granular Gas Dynamics, edited by T. Poeschel and N. Brilliantov (Springer, Berlin, 2003).
12. J. J. Brey, J. W. Dufty, C. S. Kim, and A. Santos, Phys. Rev. E 58, 4638 (1998).
13. A. Prevost, D.A. Egolf, and J.S. Urbach, Phys. Rev. Lett. 89, 084301 (2002).
14. A. E. Lobkovsky, F. Vega Reyes and J. S. Urbach, Eur. Phys. J. Spec. Top. 179, 123 (2009).
15. A. Puglisi, A. Gnoli, G. Gradenigo, A. Sarracino and D. Villamaina, J. Chem. Phys. 136, 014704 (2012).
16. R. Cafiero, S. Luding and H. J. Herrmann, Phys. Rev. Lett. 84, 6014 (2000).
17. D. R. M. Williams and F. C. MacKintosh, Phys. Rev. E 54, R9 (1996).
18. T. P. C. van Noije, M. H. Ernst, E. Trizac, and I. Pagonabarraga, Phys. Rev. E 59, 4326 (1999).
19. J. M. Montanero and A. Santos, Granular Matter 2, 53 (2000).
20. I. Pagonabarraga, E. Trizac, T. P. C. van Noije and M. H. Ernst, Phys. Rev. E 65, 011303 (2002).
21. V. Garzó, and J. M. Montanero, Physica A 313, 336 (2002).
22. P. Visco, A. Puglisi , A. Barrat, E. Trizac, F. van Wijland, J. Stat. Phys. 125, 533 (2006).
23. P. Visco, A. Puglisi, A. Barrat, F. van Wijland and E. Trizac, Eur. Phys. J. B 51, 377???387 (2006)
24. M.I. Garcia de Soria, P. Maynar, E. Trizac, Mol. Phys. 107, 383 (2009).
25. P. Maynar, M.I Garcia de Soria, and E. Trizac, Eur. Phys. J. Spec. Top. 179, 123 (2009).
26. K. Vollmayr-Lee, T. Aspelmeier and A. Zippelius, Phys. Rev. E 83, 011301 (2011).
27. for a variant of the model, see A. Puglisi, V. Loreto, U. M. B. Marconi, and A. Vulpiani, Phys. Rev E 59, 5582 (1999).
28. N. G. van Kampen, Stochastic Processes in Physics and Chemestry (North-Holland, Amsterdam, 1992).
29. P. Résibois and M. de Leener, Classical Kinetic Theory of Fluids (John Wiley, New York, 1977).
30. This property can be easily understood, since it can be seen that there is a mapping between this thermostat and the HCS. See for example J. Lutsko, Phys. Rev. E 63, 061211 (2001), or J. J. Brey, M. J. Ruiz-Montero, and F. Moreno, Phys. Rev. E 69, 051303 (2004)
31. G. Bird, Molecular Dynamics and the Direct Simulation of Gas Flows (Clarendon, Oxford, 1994).
32. A. F. Nikiforov and V. B. Uvarov, Special Functions of Mathematical Physics (Birkhäuser Verlag, Basel, 1988).
33. F. Coppex, M. Droz, J. Piasecki and E. Trizac, Physica A 329, 114 (2003).
34. R. Benzi, S. Ciliberto, R. Tripiccione, C. Baudet, F. Massaioli, and S. Succi, Phys. Rev. E 48, 29 (1993).
100023