Validity of the Boltzmann equation to describe low-density granular systems

Validity of the Boltzmann equation to describe low-density granular systems J. Javier Brey and M. J. Ruiz-Montero Fı ´sica Teo ´rica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain 共Received 15 July 2003; published 29 January 2004兲 The departure of a granular gas in the instable region of parameters from the initial homogeneous cooling state is studied. Results from molecular dynamics and from direct Monte Carlo simulation of the Boltzmann equation are compared. The results indicate that the Boltzmann equation accurately predicts the low-density limit of the system. The relevant role played by the parallelization of the velocities as time proceeds and the dependence of this effect on the density is analyzed in detail. DOI: 10.1103/PhysRevE.69.011305 PACS number共s兲: 45.70.⫺n, 51.10.⫹y, 05.20.Dd, 47.20.⫺k I. INTRODUCTION The study of granular systems has attracted much interest in the last years. The practical importance of these systems, which are present in many situations in daily life, together with the richness and complexity of their behavior 关1兴, has motivated many researches trying to describe and understand the physical mechanisms governing granular flows. Of course, the variety of what we call ‘‘granular media’’ makes it necessary to use different theoretical descriptions depending on the problem we wish to address. In the context of rapid, low-density, granular flows, the application of the methods of the kinetic theory for molecular 共elastic兲systems has proven to be a very useful tool. The starting point for this description has been in many cases the extension to inelastic collisions of the Boltzmann equation, which is derived under the same hypothesis as in elastic systems. From this kinetic equation, closed hydrodynamic equations with explicit expression for the fluxes and the associated transport coefficients 共up to second order in the gradients兲for inelastic hard particles have been derived 关2,3兴. This hydrodynamic picture has been found to provide an accurate description of granular systems in very different situations, driven and not driven. The Boltzmann equation has also been used to study the velocity distribution of a granular gas, modeled in most of the cases as a system of inelastic hard particles. Of course, the shape of the distribution depends on the state of the system. The simplest possibility is the so-called homogeneous cooling state 共HCS兲, the state of a homogeneous, freely evolving granular gas 关4兴. In that case, the Boltzmann equation is shown to admit a solution whose time dependence can be scaled out through its second moment. Deviations from Gaussianity in the HCS have been quantified by computing the kurtosis of the distribution 关5–7兴and by establishing the existence of exponential velocity tails 关6,8,9兴. The possible solutions of the Boltzmann equation in the case of a vibrated system in the absence of gravity have also been investigated, and the conditions for the existence of a steady solution whose space dependence is scaled out also through the second moment have been established 关10兴. In spite of its extended use and clear success in many cases, the validity of the Boltzmann equation to describe granular flows has been questioned since the early developments of the kinetic theory of granular systems. One of the first objections raised against its use was based on the tendency of these systems to form density clusters. Since it is only valid in the low-density limit, the Boltzmann equation is not suitable to describe the cluster evolution. This is not a fundamental objection, and in fact it can also be raised against the applicability of the Boltzmann equation to molecular systems in inhomogeneous states. The value of the density limits indeed the applicability of the Boltzmann description, but both in the elastic and the inelastic case. On the other hand, if we start from a homogeneous, low-density, initial configuration of a granular gas, under conditions that clusters will develop eventually, it can be expected that the Boltzmann equation will be valid to describe the first stages of the cluster formation, as long as regions with a too large density do not show up. A different objection concerns the validity of the molecular chaos hypothesis, which is on the basis of the Boltzmann description. As collisions in granular systems tend to make the particle velocities more parallel, velocity correlations may develop from the early stages of the evolution of a granular gas, making the molecular chaos hypothesis invalid. This problem was directly addressed by Soto et al. 关11,12兴, who studied the short-range velocity correlations in the HCS of a granular fluid, concluding that they were not relevant in the low-density case. Pagonabarraga et al. 关13兴also studied the validity of the molecular chaos hypothesis in a homogeneous granular system, but in this case driven by a random force. Again, it was found that, for dilute systems, deviations from molecular chaos were not significant. To put this work in a proper context, it should be taken into account that the driving force introduced by the authors induced some velocity correlations, as pointed out by them. The previous mentioned studies are restricted to homogeneous states of a granular gas. They correspond to very special situations, so the conclusions cannot be extrapolated to more general cases. A different situation was investigated recently by Nakanishi 关14兴, who studied the time evolution of the velocity distribution of a freely evolving granular fluid in conditions such that the HCS is unstable. He found that the kurtosis of the distribution did not remain stationary, but evolved towards the Gaussian value from the early stages of the evolution of the system, even before the clustering instability shows up. He argues that this is in contradiction with the predictions of kinetic theories based on the Boltzmann equation, being a clear indication of the growth of velocity correlations, which invalidates the molecular chaos hypothPHYSICAL REVIEW E 69, 011305 共2004兲 1063-651X/2004/69共1兲/011305共9兲/$22.50 ©2004 The American Physical Society69 011305-1 esis. He also found that the return to the Gaussian behavior became slower the more dilute the system, but as changing the density implies changing the ‘‘clustering time,’’ it was not clear to the author how the Boltzmann behavior could be recovered. The object of this work is to investigate the possible differences in the behavior of the velocity distribution function of an initially dilute, homogeneous, granular gas and the predictions of the Boltzmann equation when the system is under conditions such that the HCS is unstable. We will compare the results from molecular dynamics 共MD兲simulations at different average densities with those from the direct simulation Monte Carlo 共DSMC兲method, which is a method to solve numerically the Boltzmann equation 关15兴. The plan of the paper is as follows. In Sec. II the kinetic theory predictions for the HCS are briefly discussed, while the details of the simulations and the properties to be studied are given in Sec. III. The results for the DSMC and MD simulations are presented in Secs. IV and V. Some final comments and discussion are made in Sec. VI. II. KINETIC THEORY FOR THE HOMOGENEOUS COOLING STATE Let us consider a system of Nsmooth inelastic hard particles, spheres (d⫽3) or disks (d⫽2), of mass mand diameter ␴ . The loss of energy in collisions is characterized by a constant coefficient of normal restitution ␣ , which implies that the velocities of two colliding particles i,jbefore and after the collision are related by vi ⬘⫽vi⫺1⫹ ␣ 2共g• ␴ ˆ兲 ␴ ˆ, vj ⬘⫽vj⫹1⫹ ␣ 2共g• ␴ ˆ兲 ␴ ˆ,共1兲 where the primes denote velocities after the collision, g⫽vi ⫺vjis the relative velocity, and ␴ ˆa unit vector joining the centers of particles iand jat contact. Let us notice that 0 ⭐ ␣ ⭐1, and that the value ␣ ⫽1 corresponds to elastic collisions. The above rule implies that, in each collision, the component of the relative velocity in the direction of ␴ ˆis reduced in a factor ␣ , g⬘• ␴ ˆ⫽⫺ ␣ 共g• ␴ ˆ兲.共2兲 The Boltzmann equation describing the evolution of the velocity distribution f(r,v,t) of a system of freely evolving hard particles has the form 冉 ⳵ ⳵ t⫹v•“ 冊 f共r,v,t兲⫽JB关r,v 兩 f共r,v,t兲兴,共3兲 where JBis the 共inelastic兲Boltzmann collision operator, JB关r,v 兩 f共r,v,t兲兴⫽ ␴ d⫺1 冕 dv1 冕 d ␴ ˆ ␪ 共g• ␴ ˆ兲共g• ␴ ˆ兲 ⫻共 ␣ ⫺2b⫺1⫺1兲f共r,v,t兲f共r,v1,t兲. 共4兲 Here ␪ is the Heaviside step function and b⫺1an operator transforming velocities vand v1to its right into their precollisional values. As is the case with elastic particles, the derivation of this equation is based on the molecular chaos hypothesis, i.e., the factorization of the two-particle distribution function in the precollisional sphere: f(2)(x1,x2,t) ⫽f(x1,t)f(x2,t), where xi⬅ 兵 ri,vi 其 . If spatial and/or velocity correlations are present between colliding particles, this factorization does not hold and the assumption of molecular chaos fails. When a system of inelastic particles such as the one described evolves freely, its simplest possible state is the HCS. It is a homogeneous state with no fluxes, whose temperature T(t), defined as proportional to the average kinetic energy, evolves according to Haff’s law 关16兴, T共t兲⫽T共0兲 冉 1⫹t t0 冊 2,共5兲 where t0is the time characterizing the energy decay. The Boltzmann equation admits a solution f(v,t) describing the HCS, which obeys the scaling law 关5–7兴 f共v,t兲⫽nH v0共t兲d ␾ 冉 v v0共t兲 冊 ,共6兲 where nHis the homogeneous density, and v0the thermal velocity of the system, defined as v0⫽ 冑 2kBT/m, with kB the Boltzmann constant. Therefore, in the HCS all the time dependence of the velocity distribution can be scaled out through its second moment. The function ␾ is determined from the Boltzmann equation. Although its exact expression is not known, it was found that it does not deviate much from a Gaussian 关7兴. Then, it is sensible to expand it using the Sonine polynomials S(j), whose explicit expression can be found in Ref. 关17兴, ␾ 共c兲⫽e⫺c2 ␲ ⫺d/2 兺 j⭓0ajS(j)共c2兲,c⫽v v0.共7兲 Normalization and scaling imply a0⫽1 and a1⫽0. The coefficient a2is related to the kurtosis of the distribution a2⫽4 d共d⫹2兲 冋 具 c4 典 ⫺d共d⫹2兲 4 册 ,共8兲 and its value has been estimated from the Boltzmann equation up to linear order in a2关5,6兴. In the above expression, the angular brackets denote averages with the velocity distribution of the HCS. Both theoretical calculations and numerical simulations of the Boltzmann equation 关7兴show that, for J. J. BREY AND M. J. RUIZ-MONTERO PHYSICAL REVIEW E 69, 011305 共2004兲 011305-2 not too inelastic systems, a2is very small. Deviations from Gaussianity are important if we consider the tails of the distribution, which are exponential. Again, this has been verified from theoretical arguments based on the Boltzmann equation 关8,6兴and from DSMC simulations 关9兴. The time t0characterizing the kinetic energy decay in the HCS has also been computed by using the Boltzmann equation 关2兴, and its expression is t0 ⫺1⫽ ␨ *共 ␣ 兲4 ␲ (d⫺1)/2 共2⫹d兲⌫共d/2兲 冉 kBT共0兲 m 冊 1/2 nH ␴ (d⫺1),共9兲 where ␨ *( ␣ ) is a function that depends only on the coefficient of restitution ␣ , and that, for not too inelastic systems reads 关2,3兴 ␨ *共 ␣ 兲⯝2⫹d 4d共1⫺ ␣ 2兲.共10兲 In terms of the average number of collisions per particle, ␶ , Haff’s law takes the form T共 ␶ 兲⫽T共0兲e⫺2 ␨ * ␶ ,共11兲 i.e., the kinetic energy decays exponentially with the number of collisions. It is a well known result 关18,19兴that the homogeneous cooling state of a granular gas is unstable against long wavelength perturbations. In practice, this implies that if the system size is larger than a critical value Lc, it will spontaneously develop spatial inhomogeneities. The value of the critical size has also been determined from the Boltzmann equation 关3,20兴, and one obtains Lc⫽Cd共2⫹d兲⌫共d/2兲 2 ␲ (d⫺3)/2 冉 ␩ * 2 ␨ * 冊 1/2 ␭H.共12兲 Here, ␭H⫽1/(CdnH ␴ (d⫺1)) is the homogeneous mean free path, Cd⫽2 冑 2 for d⫽2 and Cd⫽ ␲ 冑 2 for d⫽3, and ␩ *is a function of ␣ related to the viscosity of a granular gas, whose expression can be found in Ref. 关3兴. Both MD 关18,19兴 and DSMC 关21兴simulations of freely evolving granular gases show that, in the development of inhomogeneities, first velocity vortices appear in the system and then the density becomes also very inhomogeneous. Of course, once the instability has set in, the velocity distribution of the system will no longer be the one of the HCS. Therefore, one expects the scaling to fail and a2to become time dependent. The question is whether the Boltzmann equation is valid to describe how a low-density granular system departs from the HCS. It might happen that, as pointed out by Nakanishi 关14兴, the development of velocity correlations are very important from the early stages of this departure, and then the Boltzmann equation fails. This is the main point we want to clarify in this work. III. COMPUTER SIMULATIONS Given that we want to check the validity of the Boltzmann equation to describe granular gases, we will compare MD simulations of initially homogeneous, low-density, granular systems with the predictions of the Boltzmann equation. Since the time-dependent analytical solution of the latter is not known, we will use the DSMC method developed by Bird 关15兴to construct numerical solutions of the Boltzmann equation under the desired conditions. It must be remembered that this method mimics the dynamics behind the Boltzmann equation by uncoupling the free flow of particles and collisions during a small enough time interval. Besides, collisions are always treated as if there were no correlations between colliding particles. In practice, this implies that space is divided in cells of size smaller than the mean free path, and particles within a cell collide with a probability proportional to their relative velocity. Technical details about the DSMC method have been extensively discussed in the literature 关15,22兴. The only point we want to stress is that, when using this method, the simulated system is by definition in the low-density limit, where the Boltzmann equation is supposed to apply, no matter how many simulated particles there are in a given space region. The simulations we will present in the following correspond to a two-dimensional system of freely evolving hard disks in a square box of side L, larger than the critical size Lc given by Eq. 共12兲. Periodic boundary conditions will be considered in all the cases. In the MD simulations, the event driven algorithm 关23兴will be used. For our purposes, it is important to study the effect of lowering the density on the evolution of the system. For that reason, two different, low values, of the density will be studied in the MD simulations, namely, nH⫽0.05 ␴ ⫺2and nH⫽0.0125 ␴ ⫺2. Nevertheless, it must be noticed that, for a given ␣ , changing the density implies changing the critical size of the system, as follows from Eq. 共12兲. As a consequence, if the systems have the same size L, the characteristic time governing the growth of instabilities will be also changed, and the departure from the HCS will be faster in the denser system. Therefore, if we want to have the same ‘‘distance to stability’’in systems with different densities, their sizes must be changed correspondingly, so L/␭H, and, therefore, L/Lc, remain unchanged. The same value of this scaled size should be considered in the DSMC simulations, again to expect the same characteristic time governing the departure from the HCS. Moreover, we have used a fixed value of the restitution coefficient, ␣ ⫽0.9, which is not too inelastic, but lies outside of what can be considered the quasielastic region. For this value of ␣ , Eq. 共12兲leads to Lc⯝23.36␭H. As we want the system to be well inside the unstable region, we have chosen the system size to be L⫽80␭H. Then, for nH ⫽0.05 ␴ ⫺2the number of particles in the system in the MD simulations was N⫽16000, while for nH⫽0.0125 ␴ ⫺2,N ⫽64000. In the DSMC simulations, N⫽2.048⫻106particles were considered, but it must be reminded that this number has only a statistical meaning. In all the cases, we started with the particles homogeneously distributed in the system and with a Gaussian velocity distribution. This situVALIDITY OF THE BOLTZMANN EQUATION TO . . . PHYSICAL REVIEW E 69, 011305 共2004兲 011305-3 ation was let to evolve with elastic collisions during several collisions per particle, until the system was equilibrated. Then, the inelasticity was switched on, and this time is taken as the origin, t⫽0, of our simulations. The initial elastic period ensures that the structure of the fluid at t⫽0 is the equilibrium one. Studied properties As the aim of the paper is to investigate the departure of the velocity distribution of the system from the HCS form, this will be one of the properties to be studied in the simulations. To be more precise, we will compute the second 具 v2 典 and fourth 具 v4 典 velocity moments of the total velocity distribution, and from them we will compute also a2, given by Eq. 共8兲. As far as the system stays in the HCS, a2will be constant. It is important to notice that, as in Ref. 关14兴,we will study the global velocity distribution of the system, even if inhomogeneities are present, causing local averages of the physical quantities to be quite different from their global values. The growth of inhomogeneities will be controlled by following the evolution of the density n(r,t), and velocity u(r,t), fields. To compute them, the system is divided into 30⫻30 square cells of size lc⬃2.7␭H, and properties are averaged in each cell. Again, while the system remains in the HCS, n(r,t) should be constant and u(r,t)⫽0共aside from statistical noise兲. Once the instability develops, vortices and also density inhomogeneities will show up. The question remains of how the departure of the velocity distribution from the HCS one is related to the development of these inhomogeneities. The possible failure of the molecular chaos hypothesis will be investigated by the study of several collisional averages, i.e., averages over colliding pairs of particles. These averages contain information about the two-particle distribution of colliding particles, which is the function whose factorization is assumed in the molecular chaos hypothesis. The first of the collisional averages we will consider is the pair distribution function at contact, g( ␴ ⫹). In an elastic system, if there are no correlations between colliding particles, as it is in fact assumed by the Boltzmann equation, g( ␴ ⫹)⫽1. In the case of a system of inelastic disks, the pair distribution function at contact can be easily computed by using 关24兴 共notice that there is a missprint in the expression for g( ␴ ⫹) provided in the cited reference兲 g共 ␴ ⫹,t兲⫽1⫹ ␣ ␣ 1 NnH ␲␴ 1 ⌬t兺 ␥ 苸⌬t 1 兩 ␴ ˆ•g ␥ 兩 ␪ 共⫺ ␴ ˆ•g ␥ 兲, 共13兲 where we are summing over all collisions ␥ taking place in the interval ⌬t, and the ␪ function implies that we are using the precollisional values of the quantities in the sum. When the system is in the HCS, it has been found 关24,11兴that gHCS( ␴ ⫹) does not depend on time, and takes the value gHCS共 ␴ ⫹兲⫽1⫹ ␣ 2 ␣ g0共 ␴ ⫹兲,共14兲 with g0the equilibrium elastic pair distribution function at contact, which in the two dimensional case is accurately given by 关25兴 g0共 ␴ ⫹兲⫽1⫹ ␲ 共25⫺4n ␲ 兲n 4共4⫺n ␲ 兲2.共15兲 In our simulations, time will be discretized so the average over collisions in Eq. 共13兲will be done not over the whole simulation time, but over collisions occurring in a given interval. In this way, we can study the evolution of the collisional averages and, in particular, of the pair distribution function at contact. Needless to say, g( ␴ ⫹,t) will only be computed in the MD simulations, as in the DSMC method its value is given. Another test of the validity of the molecular chaos hypothesis will be provided by the probability distribution of the impact parameter, b⫽( ␴ /2)sin ␹ , where ␹ is the angle formed by the impact relative velocity gand the unit vector ␴ ˆ. In a system of hard disks, and if there are no correlations, the impact parameter is uniformly distributed. In the HCS this distribution has also been measured 关11,26兴by computer simulations, and no deviations from uniformity have been found in the low-density limit. Here, we will consider the same discretization discussed above to construct the distribution of impact parameters in each time interval, and study its possible time dependence. Again, the impact parameter distribution will be studied only in the MD simulations, as the DSMC method assumes its uniformity. The possible velocity correlations will also be studied. A quantity that has been used to study them is the collisional average ⌫(t) defined as ⌫共t兲⫽1 N ␥ 兺 ␥ 苸⌬t ci ␥ •cj ␥ 兩 ci ␥ ⫺cj ␥ 兩 ␪ 共⫺ ␴ ˆ•g ␥ 兲,共16兲 with N ␥ the number of collisions in the interval ⌬t, and i,j the particles involved in collision ␥ . If there is neither macroscopic velocity field nor velocity correlations in the system, it is ⌫⫽0. Therefore, in the HCS, nonzero values of ⌫ are a clear signal of the presence of velocity correlations in the system. When the system leaves the HCS, there is another possible reason for these nonzero values. The buildup of velocity vortices implies that u(r,t) is different from zero, and then, even if there are no velocity correlations, ⌫will be different from zero. Let us notice that ⌫(t) can be measured both in molecular dynamics and DSMC simulations. Finally, the distribution of the angle ␾ formed by the velocities vi,vj, of colliding particles will also be computed. We are not aware of any analytic expression for this quantity in the elastic case, and, therefore, we will use the initial, elastic part, of the simulations to obtain the elastic distribution of the ␾ angle. The parallelization mechanism inherent to inelastic collisions may cause deviations from this behavior, even in the HCS. In any case, this distribution shows if there is a predominance of collisions between parJ. J. BREY AND M. J. RUIZ-MONTERO PHYSICAL REVIEW E 69, 011305 共2004兲 011305-4 ticles moving more or less parallel in the system. Again, this distribution will be measured both in MD and DSMC simulations. It is important to realize the different physical information behind the distributions P(b,t) and P( ␾ ,t). The former is related to the spatial distribution of the incident flux over colliding particles, while the latter contains only information about the relative direction of the velocities of colliding particles. In particular, it is easy to see that all the values of ␾ are compatible with any given value of b. IV. RESULTS The results we will present correspond, as we have already stated, to MD and DSMC simulations of a freely evolving system of hard disks with ␣ ⫽0.9 and L⫽80␭H, i.e., in conditions such that the HCS is well inside the unstable region. In the following, the mass mof the particles will be used as the unit of mass, and the initial kinetic energy per particle as the unit of energy. The results we will present have been averaged over several trajectories in all the cases in order to improve the statistics. Besides, the time evolution of the system will be expressed in terms of the number of collisions per particle, ␶ . As the system is prepared in an initially homogeneous situation, we expect it to stay in the HCS for a transient period, until the instability sets in. From that moment, the different physical properties will depart from their HCS values. In Fig. 1 we have plotted the evolution of the average kinetic energy per particle, 具 v2 典 /2, in the system. This quantity is proportional to the granular temperature as far as there is no macroscopic velocity field, i.e., while u⫽0. Also plotted is Haff’s law, Eq. 共11兲, describing the evolution of this property in the HCS. In all the cases there is an initial period in which the energy follows Haff’s law and, after that, the average kinetic energy of the system decays slower than in the HCS. The departure form Haff’s law occurs sooner in the MD simulations than in the DSMC one, but it is important to notice that, if we compare the two MD simulations, it occurs sooner in the case of larger density. Therefore, it seems that, although there are quantitative differences between the behavior of a finite density granular gas and the predictions of the Boltzmann equation, they are smaller the smaller the density, being similar the shape of the curves. Then, it seems sensible to expect that the Boltzmann equation predictions provide indeed the correct picture in the analytic low-density limit, n→0. The evolution of a2is plotted in Fig. 2. In Fig. 2共a兲, the initial evolution of this quantity is shown: there is an initial, very fast, decay from the initial Gaussian, a2⫽0, value to the HCS one, followed by a steady period for which the velocity distribution of the system seems to remain with the HCS form. The complete time evolution of the system is given in Fig. 2共b兲. Again, in all the cases we observe a similar behavior. After the initial steady period, a2grows, reaches a maximum, and afterwards it decays in time. The approximate duration of the steady period, ␶ st ,is ␶ st⬃15 for nH ⫽0.05 ␴ ⫺2, ␶ st⬃30 for nH⫽0.0125 ␴ ⫺2, and ␶ st⬃50 in the DSMC simulation. The position ␶ *of the maximum of a2is ␶ *⬃40 for nH⫽0.05 ␴ ⫺2, ␶ *⬃50 for nH⫽0.0125 ␴ ⫺2, and FIG. 1. Evolution of the average kinetic energy for the simulations discussed in the text. The units are chosen such that the initial value is equal to 1. Also shown is the theoretical prediction for the HCS, i.e., Haff’s law. Time is measured in average number of collisions per particle. FIG. 2. Time evolution of the dimensionless coefficient a2for the simulations discussed in the main text. The short time behavior 共a兲and the complete evolution 共b兲are displayed in different plots for the sake of clarity. VALIDITY OF THE BOLTZMANN EQUATION TO . . . PHYSICAL REVIEW E 69, 011305 共2004兲 011305-5 ␶ *⬃70 in the DSMC simulation. The height of the maximum is also larger the larger the density. A first conclusion that can be extracted from Fig. 2 is that the qualitative behavior of a2that follows from the Boltzmann equation is the same found in the MD simulations. Besides, the differences between the MD simulations and the DSMC one are smaller the smaller the density in the former, showing again a tendency to the Boltzmann behavior in the limit of very low density. Comparison of Figs. 1 and 2 is interesting in order to determine the sensitivity of Haff’s law to the exact shape of the velocity distribution function as measured by its fourth moment. It is found that, for ␶ ⫽ ␶ st , when a2begins to depart from the steady value, the temperature takes the value predicted by Haff’s law in all the cases, as is expected. What might be surprising is that, at the maximum ␶ *, the deviations from the Haff value are relatively small, the temperature being at most 1.3 times the Haff’s value. This is in agreement with simulations of homogeneous systems of inelastic hard disks, which show that the evolution of the temperature follows quite approximately Haff’s law, unless the fourth moment of the velocity distribution is very large as compared with the second one 关27兴. Once the validity of the Boltzmann description as a limit to the behavior of a granular gas for vanishing density has been shown to be consistent with the simulation results 共at least with regards to the behavior of the second and fourth moments of the velocity distribution兲, a natural emerging question is which is the mechanism taking the velocity distribution out of the HCS shape. With that purpose, the evolutions of the density and velocity fields have been studied. It must be noticed that, as we expect them to be inhomogeneous, and their spatial distribution changes from one realization to another, these fields cannot be averaged over different runs. Let us consider first the evolution of the density field. In all the cases we observe that the system remains homogeneous during what we have called the steady period ␶ st . Between ␶ st and ␶ *small density inhomogeneities begin to show up, but they are not very significant. The situation changes from ␶ *on: there is a very fast growth of density inhomogeneities, and in all the cases, elongated clusters of particles are formed, similar to those observed in previous studies of freely evolving granular systems. This behavior of the density field is quantified in Figs. 3 and 4. In the first of them, the evolution of the dispersion of the density fluctuations, ␦ ⫽ 冑 具 关n(r,t)⫺nH兴2 典 is plotted. Let us notice that this quantity is nonzero even in a homogeneous state, due to statistical noise, and its value in the homogeneous state depends on the number of particles per cell, which is different in each of the simulations. For that reason ␦ has been scaled with its value in the elastic part of the simulation, ␦ el .Itis observed in the figures that, in all the cases, ␦ has not increased much over its elastic value at ␶ ⫽ ␶ *, but from then on there is a very sharp growth of the density fluctuations as a consequence of the formation of clusters in the system. A similar behavior is exhibited by the maximum value of the density, nM, which is plotted in Fig. 4, scaled with the homogeneous density. While the system is in a homogeneous state, nMis a measure of the statistical noise. Let us remember that the systems were divided in 30⫻30 cells to compute the hydrodynamics fields: as the number of particles is smaller in the denser system, the noise in the fields will be larger, and that is the reason why, in the homogeneous part of the evolution, nM/nHis larger the larger the density. It must be also pointed out that, once the clustering begins, there is a time interval in which the growth of nMcan be fitted to an exponential, nM⬃nH⫹Cesn( ␶ ⫺ ␶ 1),共17兲 where Cand ␶ 1are constants whose value is not relevant for the present discussion, and sn⬃0.13 for n⫽0.05 ␴ ⫺2,sn ⬃0.163 for n⫽0.0125 ␴ ⫺2, and sn⬃0.167 in the DSMC simulation. Again, the discrepancies between the Boltzmann behavior and the one at finite density are smaller as we lower the density, the growth rate being almost the same for the lower density case and the DSMC simulation. In a freely evolving granular system, it has been shown that, when using ␶ as the time variable, the growth rate of the scaled transversal velocity mode is FIG. 3. Evolution of the average value of the density fluctuations, scaled with their value in the elastic case, for the simulations discussed in the main text. FIG. 4. Evolution of the maximum value of the scaled density for the simulations discussed in the main text. The symbols are from the simulations, and the dotted line is just a guide to the eye. J. J. BREY AND M. J. RUIZ-MONTERO PHYSICAL REVIEW E 69, 011305 共2004兲 011305-6 s⬜⫽ ␨ *⫺ ␩ * 2k2,共18兲 where kis a nondimensional wave vector defined as k ⫽2 ␲ /( 冑 ␲ nH ␴ l), with lthe wavelength of the perturbation. In a simulation, the maximum allowed wavelength, which is the one that grows faster, is l⫽L, due to the boundary conditions. For the values of the parameters used in this work, the theoretical prediction is s⬜⬃8.64⫻10⫺2. It has been argued 关18,28兴that the growth of density inhomogeneities in a freely evolving granular system is a consequence of the nonlinear coupling between the transversal velocity mode and the other hydrodynamic fields, in particular, the density. If that is indeed the case, the growth rate of the density, at least in its first stages, should be 2s⬜. Here, 2s⬜⬃0.173, which is very close to the value of snfound in the DSMC and in the MD lowest density simulation, confirming once again the picture described above. This agreement indicates the hydrodynamical character of the density fluctuations in the clustering regime. The growth of the scaled velocity field has also been followed in our simulations, and we found that vortices begin to develop quite soon in the system. To quantify this, we introduce ␦ u, the scaled average value of u2,as ␦ u ⫽ 具 n(r,t)u2(r,t) 典 . Again, due to statistical noise, ␦ uwill be different from zero when computed from a simulation, even if there are no fluxes. For that reason, in Fig. 5 we have plotted ␦ uscaled with its value in the elastic part of the simulation, ␦ u el . All the simulations show that the velocity field grows from the very early stages of the evolution. In fact, there is a first part of the growth of ␦ u/ ␦ u el , which can be quite well approximated by an exponential, which is common to all the simulations. After that, there is a saturation effect that translates into deviations from the exponential growth for larger times, occurring sooner the larger the density. In fact, the saturation occurs for times of the order of ␶ *, which is the time when the clustering is triggered. V. COLLISIONAL AVERAGES The evolution of the collisional averages in the system has been computed by discretizing the time in intervals ⌬ ␶ ⫽5. Properties are averaged over collisions taking place in each interval. Also, all the results have been averaged over several runs. In Fig. 6 we have plotted the evolution of the pair correlation function at contact, g( ␴ ⫹), from the MD simulations. Also included is the theoretical prediction in the HCS for the two values of the density displayed in the figure, calculated from Eq. 共14兲. In both cases it is found that, when the velocity distribution begins to depart from the HCS one, i.e., at ␶ ⫽ ␶ st , the pair correlation function takes still the HCS value and, in fact, deviations from it are quite small even at ␶ ⫽ ␶ *. After that, there is a very fast increase of g( ␴ ⫹) due to the formation of clusters in the system. Therefore, the increase in a2shown in Fig. 2 is not due to the development of positional correlations of colliding particles. These correlations do appear, but at rather later times. We have also studied the evolution of the impact parameter distribution P(b,t) and found that in none of the MD simulations discussed in this work it deviated significantly from uniformity for the times shown in this paper. Velocity correlations were investigated through the behavior of ⌫defined in Eq. 共16兲, and of the distribution function of the angle formed by the velocities of colliding particles, ␾ . In Fig. 7, the evolution of ⌫for the three simulations is shown. It must be noticed that, in a finite system, even if there are no correlations, ⌫takes a finite value that depends on the number of particles, which is different in the three cases. So, even in the elastic part of the simulation, ⌫is different in each case, being larger in the denser system, as it has less particles. Besides, the value of ⌫is different in an elastic system and in an inelastic one in the same conditions. Therefore, when the inelasticity is switched on at t⫽0, there is a very fast increase of ⌫to its HCS value. Then, there is a period over which ⌫does not change too much, and that lasts longer the lower the density, followed by an almost exponential increase of ⌫. It is found that the growth rate in this period depends on the density, becoming larger as n decreases. Finally, there is a slowing down in the increase of ⌫, and it seems that it saturates in the end to a value which FIG. 5. Evolution of the fluctuations of the velocity field, scaled with its value in the elastic case, for the simulations discussed in the main text. FIG. 6. Evolution of the pair distribution function at contact, g( ␴ ⫹), obtained in the two MD simulations discussed in the main text. The horizontal lines are the theoretical prediction for this function in the HCS at n⫽0.0125 ␴ ⫺2and n⫽0.05 ␴ ⫺2, from bottom to top. VALIDITY OF THE BOLTZMANN EQUATION TO . . . PHYSICAL REVIEW E 69, 011305 共2004兲 011305-7 is larger the larger the density. Also the slowing down begins sooner the larger the density. In any case, we find again that by lowering the density in the MD simulations, the behavior of the system approaches the one predicted by the Boltzmann equation. One could be tempted to conclude that the increase of ⌫is due to the development of correlations between velocities of colliding particles, this happening from the early stages of the evolution of the system. But one must be cautious when interpreting this result. Let us remember that in our system a velocity field is also being built up from the beginning of the evolution, as was shown in Fig. 5, and this leads to an increase in the value of ⌫that has nothing to do with the failure of the factorization of the two-particle distribution function of colliding particles, i.e., with a violation of the molecular chaos hypothesis. In fact, the behavior of ⌫ displayed in Fig. 7 is quite reminiscent of the one of the fluctuations of the velocity field. Also, in the DSMC simulation, it must be taken into account that the molecular chaos hypothesis is assumed in the very basis of the algorithm, so the increase of ⌫in that case cannot be due to the presence of precollisional velocity correlations. We conclude then that the growth of ⌫displayed in Fig. 7 is due to the instability of the scaled transversal velocity mode, which leads to the formation of vortices and, as a consequence, the velocities of colliding particles become more and more parallel. The existence of the parallelization mechanism is investigated also by studying the evolution of the distribution of the angle formed by the velocities of colliding particles, P( ␾ ). In Fig. 8 this distribution is plotted at different times for the MD simulation with n⫽0.0125 ␴ ⫺2. We have included the distribution at ␶ ⫽0, which is constructed from the elastic part of the simulation. While the system remains in the HCS, the distribution of the ␾ angle is indistinguishable from the elastic one, so there are no apparent angle correlations in this part of the evolution of the system. The situation changes at ␶ ⬃ ␶ st ( ␶ st⬃30 in this case兲, when deviations from the elastic distribution begin to show up. The relative number of collisions with larger relative angles of the velocity begin to decrease with respect to the elastic case. As time proceeds, this tendency becomes stronger, and for the final times considered in the simulation, most of the collisions correspond to particles that are moving almost parallel. The distortion of the angle distribution begins precisely at the time when the kurtosis of the total velocity distribution begins to depart from the HCS value. It seems sensible to conclude then that the behavior of the velocity distribution of the system is a consequence of the parallelization mechanism, which is inherent to inelastic collisions, and which, when the system is unstable, induces a collective behavior of the velocities of particles, which translates into a departure of the velocity distribution from the HCS one. The qualitative behavior of P( ␾ ) discussed for n ⫽0.0125 ␴ ⫺2also holds for the larger density MD simulation and for the DSMC one. Of course, the deviation from the elastic distribution occurs sooner the larger the density, but it is also present in the Boltzmann description. To have a clear picture of the parallelization of the velocities of colliding particles in the three cases, in Fig. 9 we have plotted the evolution of 具 ␾ 典 in the three simulations. The value of this quantity in an equilibrium, elastic fluid, is 具 ␾ 典 ⯝0.59 ␲ (⬃107°). In the three cases, the system remains with the elastic value of 具 ␾ 典 during what we have called the steady period, and afterwards it decreases with time, indicatFIG. 7. Evolution of ⌫for the simulations discussed in the main text. FIG. 8. Probability distribution of the angles of the velocities between colliding particles P( ␾ ) from the MD simulation with n ⫽0.0125 ␴ ⫺2. FIG. 9. Evolution of the average value of the angle between the velocities of colliding particles, 具 ␾ 典 , for the simulations discussed in the main text. J. J. BREY AND M. J. RUIZ-MONTERO PHYSICAL REVIEW E 69, 011305 共2004兲 011305-8 ing that there is a tendency to have more collisions with velocities that form a small angle. It must be also remembered that, in the MD simulations, we have observed that the impact parameter distribution did not show deviations from uniformity for the times considered in this work. VI. DISCUSSION The simulation of a system on freely evolving hard disks in conditions such that the HCS is unstable shows that the Boltzmann description seems to provide a valid picture for the behavior of a low-density granular gas, contrary to the statement made in Ref. 关14兴. This has been established by showing that there are no qualitative differences between the results obtained in low-density MD simulations and those that follow by using the DSMC method. It has been shown that, when the density is lowered in the suitable way in the MD simulations, i.e., leaving the characteristic time for the development of instabilities unchanged, the behavior of the system tends to the Boltzmann one. Nevertheless, it is also true that the deviations from the Boltzmann behavior are larger in a system in these conditions of instability than in a stable one, for the same values of the density and restitution coefficient. In other words, it could be said that the range of densities for which the Boltzmann equation holds depends rather strongly on the state of the system, and it cannot be given a general rule. In the situation considered here, the reason for this seems to be that, when the HCS in unstable, there is a parallelization of velocities of colliding particles that is more efficient the higher the density, although it is also present in the Boltzmann description. The role of the spatial correlations between colliding particles has been investigated in this work by studying the pair distribution at contact and the probability distribution of impact parameters. The former only deviates from the HCS value when density inhomogeneities are already developed in the system, while the latter remains always uniform. This implies that, in a low-density granular gas, the development of spatial correlations does not play a significant role in the early departure from the HCS. The above results seem to indicate that, when trying to extend the inelastic Boltzmann equation to finite higher density, the effect of velocity correlations between colliding particles must be incorporated. At least, in the physical situation considered in this paper, velocity correlations become important quite before the system develops significant spatial correlations, as measured by the pair distribution function at contact. It is sensible to expect that this effect increases as the inelasticity of the system increases. If this picture were right, kinetic equations for inelastic dense gases should not be based on Enskog-like equations, taking into account only spatial correlations, but new approximations incorporating the effect of velocity correlations in the precollisional sphere are needed. This implies a rather strong departure from the traditional methods of kinetic theory for elastic systems. ACKNOWLEDGMENTS We acknowledge partial support from the Ministerio de Ciencia y Tecnologı ´a共Spain兲through Grant No. BFM200200303 共partially financed by FEDER funds兲. 关1兴H.M. Jaeger, S.R. Nagel, and R.P. Behringer, Rev. Mod. Phys. 68, 1259 共1996兲. 关2兴J.J. Brey, J.W. Dufty, C.S. Kim, and A. Santos, Phys. Rev. E 58, 4638 共1998兲. 关3兴J. J. Brey and D. Cubero, in Granular Gases, edited by S. Luding and T. Po ¨schel, Lecures Notes in Physics 共Springer, Berlin, 2000兲. 关4兴C.S. Campbell, Annu. Rev. Fluid Mech. 22,57共1990兲. 关5兴A. Goldshtein and M. Shapiro, J. Fluid Mech. 282,75共1995兲. 关6兴T.P.C. van Noije and M.H. Ernst, Granular Matter 1,57 共1998兲. 关7兴J.J. Brey, M.J. Ruiz-Montero, and D. Cubero, Phys. Rev. E 54, 3664 共1996兲. 关8兴S.E. Esipov and T. Po ¨schel, J. Stat. Phys. 86, 1385 共1997兲. 关9兴J.J. Brey, D. Cubero, and M.J. Ruiz-Montero, Phys. Rev. E 59, 1256 共1999兲. 关10兴J.J. Brey, D. Cubero, F. Moreno, and M.J. Ruiz-Montero, Europhys. Lett. 53, 432 共2001兲. 关11兴R. Soto and M. Mareschal, Phys. Rev. E 63, 041303 共2001兲. 关12兴R. Soto, J. Piasecki, and M. Mareschal, Phys. Rev. E 64, 031306 共2001兲. 关13兴I. Pagonabarraga, E. Trizac, T.P.C. van Noije, and M.H. Ernst, Phys. Rev. E 65, 011303 共2001兲. 关14兴H. Nakanishi, Phys. Rev. E 67, 010301共R兲共2003兲. 关15兴G. Bird, Molecular Gas Dynamics and the Direct Simulations of Gas Flows 共Clarendon Press, Oxford, 1994兲. 关16兴P.K. Haff, J. Fluid Mech. 134, 401 共1983兲. 关17兴P. Resibois and M. de Leener, Classical Kinetic Theory of Fluids 共Wiley, New York, 1977兲. 关18兴I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. 70, 1619 共1993兲. 关19兴S. McNamara and W.R. Young, Phys. Rev. E 53, 5089 共1996兲. 关20兴J.J. Brey, M.J. Ruiz-Montero, and F. Moreno, Phys. Fluids 10, 2976 共1998兲. 关21兴J.J. Brey, M.J. Ruiz-Montero, and D. Cubero, Phys. Rev. E 60, 3150 共1999兲. 关22兴A. L. Garcı ´a, Numerical Methods for Physics, 2nd ed. 共Prentice-Hall, Englewood Cliffs, NJ, 2000兲. 关23兴M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids 共Oxford Science Publications, Oxford, 1987兲. 关24兴J.F. Lutsko, Phys. Rev. E 63, 061211 共2001兲. 关25兴D. Henderson, Mol. Phys. 30, 971 共1975兲. 关26兴S. Luding, ZAMM 80,9共2000兲. 关27兴J. J. Brey and M. J. Ruiz-Montero 共unpublished兲. 关28兴J.J. Brey, M.J. Ruiz-Montero, and D. Cubero, Phys. Rev. E 60, 3150 共1999兲. VALIDITY OF THE BOLTZMANN EQUATION TO . . . PHYSICAL REVIEW E 69, 011305 共2004兲 011305-9