Self-diffusion in freely evolving granular gases

Self-diffusion in freely evolving granular gases J. Javier Brey, M. J. Ruiz-Montero, D. Cubero, and R. Garcı ´a-Rojo Fı ´sica Teo ´rica, Facultad de Fı ´sica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain Received 18 March 1999; accepted 15 November 1999 A self-diffusion equation for a freely evolving gas of inelastic hard disks or spheres is derived starting from the Boltzmann–Lorentz equation, by means of a Chapman–Enskog expansion in the density gradient of the tagged particles. The self-diffusion coefficient depends on the restitution coefficient explicitly, and also implicitly through the temperature of the system. This latter introduces also a time dependence of the coefficient. As in the elastic case, the results are trivially extended to the Enskog equation. The theoretical predictions are compared with numerical solutions of the kinetic equation obtained by the direct simulation Monte Carlo method, and also with molecular dynamics simulations. An excellent agreement is found, providing mutual support to the different approaches. © 2000 American Institute of Physics. S1070-6631 00 00603-6 I. INTRODUCTION The properties of freely evolving granular gases have been extensively investigated in the last years. At a particle level of description, a granular medium is usually modeled as a system of inelastic hard spheres, which in the simplest case are considered to be smooth. The kinetic equations used for molecular gases, e.g., Boltzmann and Enskog equations, have been extended to account for dissipation in collisions.1,2 Once a kinetic theory description is available, it is, at least in principle, possible to proceed to a hydrodynamic description by means of the Chapman–Enskog procedure of constructing normal solutions to the kinetic equations. In this way, hydrodynamic equations, similar to the Navier–Stokes equations for molecular fluids, have been derived and explicit expressions for the transport coefficients appearing in them have been given. Although most of the derivations are restricted to the small inelasticity limit,3very recently the hydrodynamic equations have been obtained without that restriction.4Of course, there is still the problem of delimitating the range of values of the parameter characterizing dissipation for which a hydrodynamic description is valid. The simplest transport process one can think of is selfdiffusion, i.e., the diffusion of tagged particles in a system of mechanically equivalent particles, that is otherwise homogeneous. In molecular systems, self-diffusion in equilibrium has served in many ways as a prototype for more complicated transport phenomena. Here complicated refers to both from a conceptual and from a computational point of view. Self-diffusion has also been investigated in granular flows. Experimental studies include both systems with macroscopic flows5and vertically vibrated systems.6Campbell7has calculated the diffusion coefficient in a sheared cell using molecular dynamics simulations. Granular systems do not have a time-independent homogeneous state similar to the equilibrium one for molecular systems. The simplest state of a granular fluid is the so-called homogeneous cooling state HCS , in which the system is homogeneous but the temperature uniformly decreases in time as a consequence of dissipation in collisions. This will be the state we will consider. In this context, it is important to notice that the HCS is known to be unstable with respect to long-wavelength perturbations, tending to form particle clusters. This has been observed in dense systems8and also in low-density gases described by the inelastic Boltzmann equation.9The analysis of self-diffusion we will present here will be restricted to short time intervals, before the system can develop density clusters. The extension of the concept of self-diffusion for timedependent states like the HCS is not trivial. For equilibrium states of molecular fluids, there are several equivalent definitions of the self-diffusion coefficient, e.g., through a diffusion equation, as the time integral of the velocity correlation function, and from the slope of the mean-square displacement. In time-dependent states, all these quantities depend on time and the relationships between them may be rather involved. The starting point of our study will be the Boltzmann– Lorentz equation for a dilute gas of inelastic smooth hard spheres or disks. Dissipation in collisions is characterized as usual by means of a constant velocity-independent coefficient of normal restitution. The kinetic equation is particularized for a gas in the HCS and solved by means of the Chapman–Enskog procedure in Sec. II. Some details of the calculations are given in the Appendix. The analytical method we follow parallels in many points a previous work on the derivation of the Navier–Stokes transport coefficient for a simple fluid.4The main result we obtain is a diffusion equation with a transport coefficient depending on dissipation explicitly through the restitution coefficient and implicitly through the time-dependent temperature of the system. The Chapman–Enskog procedure is based on the existence of time and length scales over which there is a normal solution to the Boltzmann or Boltzmann–Lorentz equations. This is closely related with the validity of a hydrodynamic description of the system, once the accuracy of the kinetic equation has been assumed. In order to check this in the PHYSICS OF FLUIDS VOLUME 12, NUMBER 4 APRIL 2000 8761070-6631/2000/12(4)/876/8/$17.00 © 2000 American Institute of Physics context of self-diffusion in an inelastic dilute gas, in Sec. III we present numerical solutions of the Boltzmann–Lorentz equation obtained by means of the direct simulation Monte Carlo method.10 The agreement between the numerical solution and the results derived by means of the Chapman– Enskog procedure strongly supports the existence of a normal solution and, therefore, of a hydrodynamic description. The final section contains some concluding remarks and possible generalizations. II. CHAPMAN–ENSKOG SOLUTION OF THE BOLTZMANN–LORENTZ EQUATION A dilute gas of smooth hard spheres (d3) or disks (d2) of mass mand diameter is considered. Collisions between particles are inelastic and characterized by a constant coefficient of normal restitution . Some of the particles are labeled, but they are mechanically identical to the other ones. The one-particle distribution function of the the tagged particles will be denoted by fs(r,v,t), while f(r,v,t) will be used for the complete system. The time evolution of fs(r,v,t) is given by the extension to the case of smooth inelastic collisions of the Boltzmann– Lorentz equation,11 tvfsr,v,tr,v,t fsr,v,t, 1 where is the Boltzmann–Lorentz collision operator, r,v,t fsr,v,t J r,v,t fs,f d1dv1dˆ g ˆ g ˆ 2fsr,v,t f r,v1,t fsr,v,tfr,v1,t. 2 Here ˆis a unit vector along the line joining the center of particle 1 to the center of the tagged particle at contact, away from the former, g v v1is the relative velocity, is the Heaviside step function, and v,v1are the precollisional velocities leading after collision to v,v1. They are given by v v 1 2g ˆ ˆ ,v1v1 1 2g ˆ ˆ . 3 Since the distribution function f(r,v,t) obeys an independent equation, namely the inelastic nonlinear Boltzmann equation, the Boltzmann–Lorentz operator, and Eq. 1 are linear in fs(r,v,t). Tagged particles may freely exchange momentum and energy with the rest of the fluid and, therefore, these are not invariants of the collision operator J fs,f. Only the number of tagged particles is conserved. More concretely, the density nsr,t dvfsr,v,t4 obeys the conservation law nsr,t tJsr,t0, 5 with the flux of tagged particles Js, defined by Jsr,t dvvfsr,v,t. 6 The above conservation law becomes a closed hydrodynamic equation for the density of tagged particles once Jsis expressed as a functional of it, and also of the macroscopic fields defining the state of the gas as a whole. The aim here will be to obtain the lowest order in the density gradient of Js, by applying the Chapman–Enskog procedure to solve the Boltzmann–Lorentz equation. It is assumed that in the long time limit there is a normal solution to the kinetic equation 1 of the form fsr,v,t fsvnsr,t,Bi, 7 where the set Bi(r,t) denotes the hydrodynamic fields defining the macroscopic state associated to f(r,v1,t), which is assumed to also have a normal form. Although the formalism can be developed, in principle, for any state of the gas, we will restrict ourselves in the following to the so-called homogeneous cooling state HCS . This state is described by a homogeneous solution of the nonlinear Boltzmann equation, fH(v1,t), in which all the time dependence takes place through the temperature of the gas TH(t).1Therefore, it obeys the equation TH t TH fHv,t J v1,t fH,fH. 8 An evolution equation for the temperature is easily obtained from the own Boltzmann equation, THt tHTHTHt, 9 where the cooling rate H(TH) is given by H md1 /2 d112 4dd3 2nkBTH dv1,dv2v1v23 fHv1,t fHv2,t, 10 kBbeing the Boltzmann constant and nthe uniform density of the gas. The explicit form of fH(v1,t) is only known in the first Sonine approximation,1,12 fHv1,t fMv1,t1c* 4 m 2kBTH 2 1 4 2d m 2kBTH1 22d d 4, 11 c*32 1 1 2 2 9 24d8d41 30 21. 12 Here fMis the Maxwellian, fMv1nm 2kBTHt d/2 em1 2/2kBTHt. 13 In the same approximation, 877Phys. Fluids, Vol. 12, No. 4, April 2000 Self-diffusion in freely evolving granular gases H122d1 /2 d d/2 nd1kBTH m 1/2 13c* 32 . 14 Now we proceed with the application of the Chapman– Enskog algorithm to Eq. 1 . The distribution function is written as a series expansion in a formal parameter measuring the ‘‘uniformity’’ of the system, fsfs0fs12fs2. 15 Since in the present case the only possible hydrodynamical inhomogeneity is in the density ns, each factor of implies an implicit nsfactor. Consistently, the balance equation for the density, Eq. 5 , is expanded as nsr,t tj0 j1tjnsr,t, 16 with tjnsr,tJsjr,t, 17 Jsjr,t dvvfsjr,v,t. Then, to zeroth order in the kinetic equation 1 becomes TH t TH fs0fs0J fs0,fH. 18 We have taken into account that, as indicated in Eq. 7 , there is a time dependence through the temperature of the system TH(t). This dependence is not associated to any gradient in the system and, therefore, it is considered as being of zeroth order in . A comparison of Eqs. 8 and 18 shows that fs (0)(r,v,t) must be proportional to fH(v,t). Then, because of normalization, fs0r,v,tnsr,t nfHv,t. 19 It follows that the flux of tagged particles vanishes to zeroth order, i.e., Js (0) 0. To first order in the density gradient, Eq. 1 leads to an equation for fs (1) , t0nsns fs0TH t TH fs1vnsns fs0fs1. 20 The zeroth-order results imply t0ns0, ns fs01 nfH, 21 and, therefore, Eq. 20 becomes Ht THtTHtfs11 nvnsfH. 22 The solution to this equation is proportional to ns, i.e., it has the form fs1r,v,tB v nsr,t. 23 The vector B(v) will also depend on time through TH(t), but this will not be indicated explicitly. Substitution of this into Eq. 22 yields Ht THtTHtB v 1 nvfH. 24 The contribution to the particle flux of first order in the density gradient is given by Js1D nsr,t, 25 with D1 ddvv B v . 26 Upon deriving the above expression we have made use of symmetry considerations. An equation for the self-diffusion coefficient Dis easily derived from Eq. 22 , Ht THtTHDDkBTHt m, 27 where Dis a functional of B(v), D dvv B v dvv B v . 28 Since HTH 1/2 , dimensional analysis requires that also D T1/2 and D TH D 2TH , 29 so that Eq. 27 leads to DkBTHt mD H 2 1 . 30 In order to get an explicit expression for Dwe have considered a first Sonine approximation, in which only the leading term in the expansion of B(v) in Sonine polynomials is kept. Then we approximate B(v)vfM(v), with fMbeing the Maxwellian defined in Eq. 13 . Details of the analysis are given in the Appendix, where also the explicit form of fs (1) is derived. The result is D*D D0 4 12c* 32 4 3 2 . 31 Here D0is the elastic limit ( 1) of the self-diffusion coefficient at equilibrium, but with a time-dependent temperature TH(t), D0 d d/2 4d1 /2 1 nd1 kBTH m 1/2 . 32 In summary, we have derived the self-diffusion equation for a low-density granular gas in the HCS, tnsr,t D t 2nsr,t, 33 878 Phys. Fluids, Vol. 12, No. 4, April 2000 Brey where D(t) is given by Eq. 31 . As it is the case for the other Navier–Stokes transport coefficients,4Dhas a double dependence on the coefficient of restitution . First, there is the explicit dependence given by the factor D*( ) in Eq. 31 , and, second, there is an intrinsic dependence through the time evolution of the temperature TH, which is also governed by the coefficient . An expression for the self-diffusion coefficient in a system of smooth inelastic hard spheres has been derived before independently by Hsiau and Hunt13 and by Savage and Dai.14 In our notation, their result reads as D*2 1, 34 and differs from the one obtained here, except in the limit 1. Figure 1 shows the reduced diffusion coefficient D*as a function of the coefficient of restitution in the interval 0.5 1 for d3. It is seen that the deviation from the elastic functional form is quite relevant, even for values of close to unity. On the other hand, let us point out that if the distribution function of the gas as a whole is approximated by a Gaussian, i.e., we formally put c*( ) 0 in Eq. 31 , the discrepancy with the curve plotted in the figure cannot be observed on the scale used. III. DIRECT SIMULATION MONTE CARLO METHOD The direct simulation Monte Carlo DSMC method provides an algorithm to obtain numerical solutions of the Boltzmann equation. The general idea is to mimic the dynamical processes involved in the kinetic equation.10 Its extension to deal with inelastic collisions is straightforward. Since the method has been discussed extensively many times in the literature, only some specific details will be given here. An important peculiarity to be stressed is that we are considering the Boltzmann–Lorentz equation and not the general nonlinear Boltzmann equation. That means that the distribution function of the system, fH(v,t), is given, i.e., it is a needed input to solve the kinetic equation. We do not investigate whether or not a given system would actually remain in the homogeneous cooling state for the time interval considered; it is just assumed that this is the case. Of course, this avoids by construction the possibility that the system develops velocity and density instabilities that are known to occur in a freely evolving granular system. The time dependence of the diffusion coefficient in Eq. 33 can be eliminated through a change in the time and space variables. Let us introduce the Boltzmann collision frequency, 0(t), corresponding to a dilute gas with a Maxwellian distribution function at temperature TH(t), 0t2d1 /2 d/2 nd1kBTHt m 1/2 . 35 We define dimensionless time and space variables by 1 20 t dt 0t,l0t 2 m kBTHt 1/2 r. 36 The time scale provides an accurate estimation of the average number of collisions per particle in the time interval between 0 and t, while the unit of length introduced above is proportional to the mean-free path of the particles. In terms of the above variables, the self-diffusion equation reads as sl,d 4D*l 2sl, , 37 where l 2is the Laplace operator in lspace, and we have scaled the density of tagged particles with the average density of the system, i.e., sns/n. Equation 37 is a diffusion equation with a constant diffusion coefficient D*( )d/4. It follows that the mean square deviation of the ‘‘position’’ lof the tagged particles after a ‘‘time’’ interval is11 l2;d2 2D*. 38 This extends to the case of dissipative collisions the wellknown relationship between the mean square displacement of the particles and the self-diffusion coefficient in elastic systems. Equation 38 is equivalent to tr2;t2dD ,t, 39 and also to r2;d2 2 d/2 2d1n2 2 d1D*. 40 The first method we have used to measure the diffusion coefficient numerically is directly based in Eq. 38 . In the simulations we have considered N2.5 105tagged hard spheres. These particles did not see each other, so that we could say equivalently that we generated N2.5 105indeFIG. 1. The reduced self-diffusion coefficient D*as a function of the coefficient of restitution . The solid line is the theoretical prediction given by Eq. 31 , while the symbols have been obtained from the direct Monte Carlo simulation of the Boltzmann equation for a system of hard spheres by using the mean square displacement of the particles circles , and the time evolution of a sine perturbation in density of tagged particles squares . 879Phys. Fluids, Vol. 12, No. 4, April 2000 Self-diffusion in freely evolving granular gases pendent trajectories of a single tagged particle. Only the position of the tagged particles was recorded along the simulation, since, as already discussed, the distribution of the complete system was not generated by the simulation, but introduced as an input. Every time there was a collision of the tagged particle, the velocity of the other particle involved in the collision was generated from the distribution function corresponding to the HCS in the first Sonine approximation, Eq. 11 . As a consequence, there is no need either to divide the system into cells, since the ‘‘bath’’ seen by the tagged particles does not depend on their position. For the same reason, it is not necessary to specify the size of the system nor any kind of boundary conditions. In other words, the system is considered as infinite and homogeneous in all directions. Finally, let us mention that the time interval over which it is assumed that there are not correlations between free motion and collisions is t5 10 2 2kBT(0)/m1/2, where ( n2)1is the mean-free path of the particles and T(0) is the initial temperature of the system. This time interval must be chosen such that it remains small as compared with the average time between collisions. As the temperature of our system decreases monotonically in time, the above choice guarantees that this condition is verified all along the simulation. In Fig. 2 we have plotted ( l)2as a function of the cumulative number of collisions per particle, measured by the time scale defined in Eq. 36 . Several values of the coefficient of restitution in the interval 0.6 0.99 have been considered. It it observed that, after a short transient period in which each tagged particle has collided about five times, the scaled mean square displacement becomes a linear function of , as predicted by Eq. 38 . From each value of the slope in the linear regions, we have computed D*by dividing the slope by 9 2. The comparison with the Chapman– Enskog value, Eq. 31 is presented in Fig. 1. The agreement is very good in all the range of values of considered, although the discrepancy increases as the value of decreases. It must be realized that the value of the scaled diffusion coefficient D*( ) changes by an amount of the order of 50%. A second way of measuring the diffusion coefficient numerically follows as a consequence of the linear character of the self-diffusion equation. If we introduce the Fourier transformed of the density of tagged particles ns,k, the diffusion equation 37 leads to ns,kns,k0 exp sD, 41 where sD d 4D*k2. 42 We have considered a system of size L50 in the xdirection, which is the direction along which the density gradient will be studied. Consequently, periodic boundary conditions are introduced in this direction. On the other hand, no density gradient of tagged particles is allowed in the yand z directions. Suppose an initial distribution of tagged particles, given by nsx,0 n01 sin qmx, 43 with n0being a constant and qm2 /L. This density profile is compatible with the periodic boundary conditions and, therefore, Eq. 41 predicts that nsx,t nsx,t n0nsin qmx, 44 n n0exp sD,m. Here sD,mdD*km 2/4, where kmis the dimensionless wave number corresponding to qm, i.e., km2H(t)qm/H(t). To analyze the density distribution in the xdirection, the system was split into 100 layers of the same width perpendicular to the xaxis. Of course, this was to the only effect of coarse graining the density of tagged particles, since the distribution function of all particles was again given by fHin Eq. 11 . In the simulations we have verified that the profile is well fitted by a sine function with wavelength Lin all the time intervals considered. As an example, the values obtained for the amplitude of the density perturbation for 0.8 are plotted in Fig. 3 as a function of time. They decay exponentially, in agreement with the theory, and from the slope of the logarithm, the reduced self-diffusion coefficient D*( ) can be computed. The values obtained in this way for different values of the coefficient of restitution have been included in Fig. 1. As expected, they are practically the same as those obtained from the mean square displacement, then providing a self-consistent test for both the theory and the numerical method. Let us mention that the use of a sine perturbation of the hydrodynamic fields has been recently employed to measure the Navier–Stokes transport coefficients of a dilute granular gas in the linear approximation.15 IV. MOLECULAR DYNAMICS RESULTS The results discussed in the previous sections rely on the validity of the Boltzmann equation to describe a low-density granular gas. Also, the possibility of observing a system in FIG. 2. Mean square displacement of the particles as a function of time for several values of the restitution coefficient. Both quantities are measured in the dimensionless units defined in the main text. 880 Phys. Fluids, Vol. 12, No. 4, April 2000 Brey the HCS for a time period long enough as to reach the selfdiffusion regime, i.e., the existence of a time window in which the diffusion equation applies, was taken for granted. This is not at all evident due to the cluster instability exhibited by systems composed by inelastic hard particles.8,9 For these reasons, here we present and discuss some results obtained by means of molecular dynamics simulation. Of course, when using molecular dynamics, the density of the system is not asymptotically small, as required by the Boltzmann equation, and density effects can, in general, be expected. The Enskog equation provides a useful generalization of the Boltzmann equation to higher densities for a gas of hard spheres or disks.11,2 Applied to the self-diffusion problem in an inelastic gas in the HCS it has the form tvfsr,v,t gen J r,v,t fs,fH, 45 where the collision operator Jis still defined by Eq. 2 and ge(n) is the equilibrium pair correlation function at the density nof the system and a distance . It follows directly that the self-diffusion coefficient from the Enskog equation is still given by Eq. 31 , except for an obvious factor ge1(n) in the expression of D0, i.e., D0 Ed d/2 4d1 /2 1 nd1gen kBTH m 1/2 . 46 Also, it is evident that, when properly scaled, the solutions of the Boltzmann–Lorentz equation and Eq. 45 coincide. Therefore, the numerical results in the previous section can be directly translated to the Enskog equation. Let us stress that this is due to the homogeneity of the state we are considering We have simulated a system of hard disks in a square domain with doubly periodic boundary conditions. The initial condition was a uniform distribution of particles and an isotropic and homogeneous velocity distribution. Simulations were run with elastic collisions ( 1) for a period of time large enough for the system to reach a Maxwellian velocity distribution. In this way the initial conditions for the simulation of the inelastic gas ( 1) was generated. The simulations typically involved N6400 particles and the results we will present were averaged over 30 trajectories. The simulation technique was based on the ‘‘event driven’’ algorithm.16 A discussion of the application of the method to granular fluids is given in Ref. 17. In Fig. 4 we present results obtained in a very dilute system with a solid fraction 5 10 4, that is equivalent to a number density n26.25 10 4. For this density it is ge(n) 1.0008.18 The figure shows mD(4D0 E)1 t(r)2as a function of t0(0) for several values of the coefficient of restitution in the interval 0.70 1. Here 0 is the initial value of the Boltzmann collision frequency defined in Eq. 35 . After a short transient time, of the order of a few collision times, the curves reach a time-independent plateau. In the simulations, two quantities were recorded to verify whether the system stayed in the HCS, namely the root mean square of the local density and the ratio of the total macroscopic kinetic energy to the total thermal energy. In all the reported simulations both quantities remained with very low values, indicating that no significant density or momentum organization was present.19 It was also checked that the time evolution of the temperature was accurately described by the solution of Eq. 9 . According with the theoretical prediction, Eq. 39 , the steady values in Fig. 4 should correspond to D*( ), as given by Eq. 31 . The comparison between the latter and the simulation values is presented in Fig. 5. Again, a very good agreement is observed. Moreover, it is interesting to notice that the slight deviations here go in the same way as for the numerical solution of the Boltzamnn equation. For small dissipation very close to unity the theoretical curve lies below the simulation data and there is a crossover as the value of decreases. This seems to indicate that the small discrepancy is in both cases due to the approximate character of the Chapmann–Enskog solution we have obtained, FIG. 3. Time evolution of the amplitude of a sine density perturbation of tagged particles for 0.8. The density is reduced by the initial amplitude, and the time is measured in the reduced units defined in the text. The dashed line is the best fit to a straight line exponential behavior . FIG. 4. Time evolution of the reduced slope mD(4D0 E)1t(r)2for 0.7,0.8,0.9,1, from top to bottom. The steady values determine the reduced self-diffusion coefficient. The data have been obtained by molecular dynamic simulations of a system of hard disks. 881Phys. Fluids, Vol. 12, No. 4, April 2000 Self-diffusion in freely evolving granular gases namely to the use of the first Sonine approximation. In fact, if the values for the self-diffusion coefficient obtained from the numerical solution of the Boltzmann–Lorentz equation are compared directly with those following from the molecular dynamics simulation, an excellent agreement is found. This happens in spite of the fact that while in Fig. 1 we considered a system of hard spheres (d3) now we are dealing with hard disks (d2). The reason is that all the dependence of D*( ) on dcomes through c*( ), and this quantity gives a negligible contribution. V. DISCUSSION AND CONCLUSION In this paper we have discussed the simplest example of transport in a granular gas, namely self-diffusion in the homogeneous cooling state. We have combined several different approaches: numerical solutions of the kinetic equations obtained by the direct simulation Monte Carlo method, Chapman-Enskog solution of the same equations, and molecular dynamics simulations. The very good agreement found over quite a wide interval of values of the coefficient of restitution, provides a strong test of the validity and usefulness of the different approaches. It shows that kinetic theory and a hydrodynamic description are accurate to describe some states of granular fluids. In particular, the molecular dynamic results we have presented indicate that there is a relevant time interval in which the hydrodynamic limit is attained while the system stays in the HCS or, at least, very close to it, in the sense that the possible instability effects are small. The self-diffusion process in a granular fluid in the HCS can be described at a macroscopic level by a hydrodynamic equation of the same form as in the case of ordinary fluids. The main difference is in the time dependence, through the temperature of the system, of the self-diffusion coefficient. Also, there is an explicit dependence on the restitution coefficient characterizing dissipation. The results presented here are restricted to a low-density, small gradient in the density of tagged particles and to a system in the HCS. Moreover, the time must be large enough to reach the hydrodynamic regime. The extension to higher densities, beyond the range of validity of the Enskog equation, requires the use of the pseudo-Liouville equation along the lines discussed in Ref. 2. An analysis of short times and large gradients can be carried out by means of model kinetic equations, such as those discussed in Ref. 4. ACKNOWLEDGMENTS This research was partially supported by the Direccio ´n General de Investigacio ´n Cientı ´fica y Te ´cnica Spain through Grant No. PB98-1124. APPENDIX: EVALUATION OF AND 1 In this appendix some of the calculations leading to Eq. 31 are detailed. In order to evaluate velocity integrals involving the linear collision operator , it is convenient to use the relationship dvYvXvfMv d1dvdv1dˆ g ˆ g ˆ Yv 2XvfMvfHv1XvfMvfHv1 d1dvdv1dˆ g ˆ g ˆ XvfMv fHv1Yv*Yv, A1 where v*is the after-collision velocity of the tagged particle, v*v1 2ˆ g ˆ . A2 Use of Eq. A1 into Eq. 28 after approximating B(v) by the first term in the Sonine polynomials expansion, B v vfMv, A3 leads to D md1 ndkBTH 1 2dvdv1fMvfHv1 dˆ g ˆ g ˆ 2ˆ v. A4 The angular integration over ˆcan be easily carried out with the result D md1 ndkBTH 1 2 d1 /2 d3 2 dvdv1fMvfHv1gg v . A5 Using Eq. 11 for fH, the velocity integrations can be transformed into the product of Gaussian integrals by changing variables to gand G(v v1)/2. The integrals are straightforward although lengthy, and the final result is FIG. 5. A comparison of the reduced self-diffusion coefficient for a system of hard disks obtained from molecular dynamics simulations symbols with the theoretical prediction given by Eq. 31 solid line . 882 Phys. Fluids, Vol. 12, No. 4, April 2000 Brey D 2 1 d1 /2 d d/2 kBTHt m 1/2 nd11c* 64 . A6 Substitution of this expression and Eq. 14 into Eq. 30 leads to Eq. 31 directly. The first order in the gradient contribution to the distribution function is given by Eq. 23 , that in the first Sonine polynomial approximation takes the form fs1r,v,t cDfMv v nsr,t. A7 The coefficient cDcan be expressed in terms of the selfdiffusion coefficient by use of this expression into Eq. 26 , DcD ddv2fMvcDnkBTHt m. A8 The distribution function is, therefore, fs1r,v,tDm nkBTHtfMv v nsr,t. A9 1A. Goldshtein and M. Shapiro, ‘‘Mechanics of collisional motion of granular materials. 1. General hydrodynamic equations,’’ J. Fluid Mech. 282, 75 1995 . 2J. J. Brey, J. W. Dufty, and A. Santos, ‘‘Dissipative dynamics for hard spheres,’’ J. Stat. Phys. 87, 1051 1997 . 3J. T. Jenkins and M. W. Richman, ‘‘Grad’s 13-moment system for a dense gas of inelastic hard spheres,’’ Arch. Ration. Mech. Anal. 87, 355 1985 ; ‘‘Kinetic theory for plane flows of a dense gas of identical, rough, inelastic, circular disks,’’ Phys. Fluids 28, 3485 1985 ; C. K. W. Lun, S. B. Savage, D. J. Jeffrey, and N. Chepurniy, ‘‘Kinetic theories for granular flow: inelastic particles in Couette flow and slightly inelastic particles in a general flowfield,’’ J. Fluid Mech. 140, 223 1984 ; N. Sela and I. Goldhirsch, ‘‘Hydrodynamic equations for rapid flows of smooth inelastic spheres, to Burnett order,’’ ibid. 361, 41 1998 . 4J. J. Brey, J. W. Dufty, C. S. Kim, and A. Santos, ‘‘Hydrodynamics for granular flows at low density,’’ Phys. Rev. E 58, 4638 1998 . 5V. V. R. Natarajan, M. L. Hunt, and E. D. Taylor, ‘‘Local measurements of velocity fluctuations and diffusion coefficients for a granular material flow,’’ J. Fluid Mech. 304, 1 1995 , and references therein. 6O. Zik and J. Stavans, ‘‘Self-diffusion in granular flows,’’ Europhys. Lett. 16, 255 1991 . 7C. S. Campbell, ‘‘Self-diffusion in granular shear flows,’’ J. Fluid Mech. 348, 85 1997 . 8I. Goldhirsch and G. Zanetti, ‘‘Clustering instability in dissipative gases,’’ Phys. Rev. Lett. 70, 1619 1993 ; S. McNamara and W. R. Young, ‘‘Inelastic collapse in two dimensions,’’ Phys. Rev. E 50, R28 1994 ; P. Deltour and J. L. Barrat, ‘‘Quantitative study of a freely cooling granular medium,’’ J. Phys. I 7, 137 1997 . 9J. J. Brey, M. J. Ruiz-Montero, and D. Cubero, ‘‘Homogeneous cooling state of a low-density granular flow,’’ Phys. Rev. E 54, 3664 1996 . 10G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows Clarendon, Oxford, 1994 . 11P. Re ´sibois and M. de Leener, Classical Kinetic Theory of Fluids Wiley, New York, 1977 ; J. A. McLenann, Introduction to Non-Equilibrium Statistical Mechanics Prentice–Hall, Englewood Cliffs, NJ, 1989 . 12T. P. C. van Noije and M. H. Ernst, ‘‘Velocity distributions in homogeneously cooling and heated granular fluids,’’ Granular Matter 1, 57 1998 . 13S. S. Hsiau and M. L. Hunt, ‘‘Kinetic theory analysis of flow-induced particle diffusion and thermal conduction in granular material flows,’’ Trans. ASMEC: J. Heat Transfer 115, 541 1993 . 14S. B. Savage and R. Dai, ‘‘Studies of shear flows. Wall slip velocities, ‘layering’ and self-diffusion,’’ Mech. Mater. 16, 225 1993 . 15J. J. Brey, M. J. Ruiz-Montero, and D. Cubero, ‘‘On the validity of linear hydrodynamics for low density granular flows described by the Boltzmann equation,’’ Europhys. Lett. 48, 359 1999 . 16M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids Oxford University Press, Oxford, 1987 . 17S. Luding, E. Clement, A. Blumen, J. Rajchenbach, and J. Duran, ‘‘Anomalous energy dissipation in molecular dynamics simulations of grains: The ‘detachment’ effect,’’ Phys. Rev. E 50, 4113 1994 . 18N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, ‘‘Equation of state calculations by fast computing machines,’’ J. Chem. Phys. 21, 1087 1953 . 19S. McNamara and W. R. Young, ‘‘Dynamics of a freely evolving, twodimensional granular medium,’’ Phys. Rev. E 53, 5089 1996 . 883Phys. Fluids, Vol. 12, No. 4, April 2000 Self-diffusion in freely evolving granular gases