Mesoscopic theory of critical fluctuations in isolated granular gases

Mesoscopic Theory of Critical Fluctuations in Isolated Granular Gases J. Javier Brey,*A. Domı ´nguez, M. I. Garcı ´a de Soria, and P. Maynar Fı ´sica Teo ´rica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain (Received 24 November 2005; published 21 April 2006) Fluctuating hydrodynamics is used to describe the total energy fluctuations of a freely evolving gas of inelastic hard spheres near the threshold of the clustering instability. They are shown to be governed only by vorticity fluctuations that also lead to a renormalization of the average total energy. The theory predicts a power-law divergent behavior of the scaled second moment of the fluctuations, and a scaling property of their probability distribution, both in agreement with simulations results. A more quantitative comparison between theory and simulation for the critical amplitudes and the form of the scaling function is also carried out. DOI: 10.1103/PhysRevLett.96.158002 PACS numbers: 45.70.n, 05.20.Dd, 51.10.+y Granular gases are assemblies of macroscopic particles evolving independently between inelastic collisions [1]. The methods of nonequilibrium statistical mechanics, kinetic theory, and hydrodynamics have been successfully extended to describe the observed macroscopic behavior and also, although in a much more limited form, the fluctuations around it [2]. The lack of energy conservation makes these systems behave quite differently from molecular fluids. A simple widely used model for them consists of smooth inelastic hard spheres (IHS’s), with momentum conserving dynamics. Inelasticity is characterized by means of a constant coefficient of normal restitution . Recently, molecular dynamics (MD) simulation results have been reported for the total energy fluctuations of a two-dimensional freely evolving IHS gas, near the threshold of the clustering instability [3]. The dimensionless second moment was found to exhibit a power-law divergent behavior with the distance to the instability. Also, the scaled cooling rate was found to tend to zero according to a power law, although in a weak way. Besides, the distribution function for the energy fluctuations, when properly scaled, turned out to be independent of the parameters defining the system. This was associated with a scaling property of the distribution. Quite remarkably, the scaling function was very well fitted by the same expression as several equilibrium and nonequilibrium molecular systems [4,5]. The main goal of this Letter is to provide an explanation for the above results on the basis of fluctuating hydrodynamics [6]. Consider an isolated system of NIHS’s of mass mand diameter . The total (kinetic) energy ~ Eof the system can be expressed in the form [7] 2~ EtZdrd~ nr;t~ Tr;tm~ nr;t~ u2r;t;(1) where dis the dimension of the system, ~ nr;tthe number density field, ~ Tr;tthe temperature field, and ~ ur;tthe flow field. The tildes indicate that all the quantities are understood as fluctuating variables. In the following, systems in the homogeneous cooling state (HCS) will be considered. At a macroscopic level, this state is characterized by a constant uniform density nH, a vanishing flow field uH0, and a uniform time-dependent temperature obeying the law [8] @tTHtHTHTHt, where H/ THt1=2is the cooling rate. Moreover, we will restrict ourselves to the region in which the amplitudes of the fluctuations of the fields around their HCS values remain small on the average [see below Eq. (12)]. Then retaining up to quadratic order in the deviations, Eq. (1) yields ~ Et~ EtEHt 1 2ZdrdnH~ Tr;td~ nr;t~ Tr;t mnHj~ ur;tj2:(2) Here, EHtdNTHt=2,~ nr;t~ nr;tnH, ~ ur;t~ ur;t, and ~ Tr;t~ Tr;tTHt.Itis now convenient to introduce dimensionless position, l, and time, s, scales by lr=l0and ds vHtdt=l0, respectively, where vH2THt=m1=2is the thermal velocity and l0nHd11is proportional to the mean free path. Moreover, dimensionless fields are defined by l;s~ nr;t=nH,!l;s~ ur;t=vHt, and l;s~ Tr;t=THt. Then, Eq. (2) takes the form s0s V1 V2X kksks2 dj!ksj2;(3) where s~ Es=EHs,VLdis the volume of the system in the new units, and the Fourier transforms of the fields have been introduced. It is assumed that after a time of the order of the mean free time, the system reaches a regime in which all its energy is stored in the hydrodynamic modes. In this regime, the hydrodynamic fields are expected to be described at a mesoscopic level by fluctuating hydrodynamic equations. Here, they will be assumed to be linear Langevin equations obtained by linearizing the Navier-Stokes equations for a granular gas around the HCS. Moreover, it is postulated that the noise terms are defined by the same properties as for molecular, elastic PRL 96, 158002 (2006) PHYSICAL REVIEW LETTERS week ending 21 APRIL 2006 0031-9007=06=96(15)=158002(4)$23.00 158002-1 ©2006 The American Physical Society gases [9]. This is not expected to be true, except in the nearly elastic limit, i.e., when is very close to unity. Consequently, the theory will be restricted in the following to this limit. Thus, the transversal flow field or vorticity field, !k?, obeys the following equation in the scaled variables [6,9]: @s=2k2!k?sk?s:(4) In this expression, HTHtl0=vHtand  HTHt=mnHl0vHt,Hbeing the shear viscosity. The noise term k?sis Gaussian, with hk?sk0?s0i  V2 Nss0k;k0k2l;(5) lbeing the unit tensor in the subspace perpendicular to k, and the angular brackets denoting average over the noise realizations. A main advantage of using the scaled variables is that the coefficients in Eq. (4) and the strength of the noise are time independent, contrary to what happens in the original variables. The equation shows that !k?grows in time for those values of ksuch that ?k=2 k2>0. Although this does not imply by itself that the HCS is linearly unstable, due to the time-dependent scaling of the velocity introduced above, simulation results and nonlinear analytical analysis of the Navier-Stokes equations have shown that this growth is the origin of the clustering instability [10,11]. The minimum value of k for a system of linear extent L, measured in the lscale, is kmin 2=L. Then, for given values of the other parameters, the system becomes unstable if L>L c, with Lc22=1=2.ForL<L c, the HCS is stable and the long time solution of Eq. (4) is !k?sZs 1 ds0ess0?kk?s0:(6) From this expression, it is easily obtained h!k?s!k0?s0i   V2k2 2N?kess0?kk;k0l;(7) for ss01. This shows that as Lapproaches Lcfrom below, the amplitudes of the fluctuations of the transversal components of the velocity increase very fast due to contributions from values of kclose to kc. For the same reason the decay of these fluctuations becomes very slow. This is not the case for the fluctuations of the other hydrodynamic fields, whose Langevin equations are decoupled from Eq. (4) [6]. Therefore, it seems possible to consider a range of values of f L LcL=Lcwhere the fluctuations of !k;?, although still small, dominate over the fluctuations of density and temperature. But, although this is true for components with k>0, some care is needed when analyzing Eq. (3), since it involves 0s. The Langevin equation for sis obtained from the linearization around the HCS of the macroscopic average equation for the total energy, @tEtd 2Zdrnr;tHn; TTr;t:(8) The result is @sss 3 2V0s:(9) Here, the dependence of the cooling rate on the temperature has been taken into account. Moreover, the noise term discussed in Ref. [12], associated with the localized character of the energy dissipation, has been omitted. Although it can be expected to be negligible far from the instability in the quasielastic limit, this may not be the case near the instability. Equation (9) shows a coupling between the fluctuations of the volume averaged temperature and those of the total energy. Use of Eq. (3) into Eq. (9) and neglecting contributions from the density and longitudinal velocity fluctuations gives @ss 2s !s; !s 6 V2dX kj!k?sj2; (10) valid in the region f L 1. The long time limit of the average value of sis, therefore, hist lim s!1h!si  3d1 Nd X k k2 ?k:(11) Since we are considering f L 1, the sum over kin the above expression is dominated by the 2dmodes with the largest wavelength, for which ?kmin’f L. Using this into Eq. (11), it follows that there is a renormalization by fluctuations of the average total energy of the HCS, Eth~ Eti, given by EtEHt13d1 nHLd cf L1:(12) Consistency of the theory we are developing requires that nHLd cf L11, a condition involving the inelasticity and the distance to the instability. Similarly, there is also a renormalization of the temperature of the HCS, Tt h~ Ttist, that can be evaluated directly from the long time limit of the average of Eq. (9), TtTHt1h0ist VTHt12d1 nHLd cf L1: (13) Alternatively, an effective temperature Teftcan be defined as Teft2Et=Nd. Of course, the form of the renormalized law for the temperature depends on the definition used for the latter. In Ref., [3], what was actually measured was  ef efTefl0=vHTef, with ef defined by @tTef efTefTef. Then, it is found PRL 96, 158002 (2006) PHYSICAL REVIEW LETTERS week ending 21 APRIL 2006 158002-2  ef "13d1 nHLd cf L1#1=2 :(14) This result predicts that near the clustering instability threshold, 2 ef 2=2Af L1with A3d 1=nHLd c, that is just the behavior observed in Ref. [3]. Define f Es~ EsEs=Esshist EHs=Esand !s !sh!ist. Note that we are considering deviations from the renormalized average values, i.e., including the fluctuations effects, and not from the macroscopic bare values. A standard calculation using Eq. (7) and exploiting the Gaussian character of the noise, gives that at the instability threshold and for ss01it is h!s!s0ist 9d1 n2 HL2d cdf L2ess0=sc;(15) where sc2f L1is a divergent ‘‘critical’’ relaxation time. Now Eq. (10) can be easily solved with the result hf Esf Es0ist h!s!s0ist;(16) valid for ss01. Thus below the instability, the scaled total energy fluctuations decay with the same rate as the fluctuations of the kinetic energy associated with the transversal modes of the velocity. For ss0, Eq. (16) yields 2 Ehf E2ist A2 f L2;(17) with A2 9d1=n2 HL2d cd. Therefore, close to the instability point, the relative dispersion of the total energy fluctuations Epresents a divergent behavior with a critical exponent 1, and an amplitude Adepending on nH and (through the value of the critical length Lc). Again, this is the same behavior as reported in Ref. [3] from MD simulations. To carry out a more detailed check of the theory presented here, we have performed MD simulations of twodimensional systems with different values of and nH(see Table I). In all cases, the dependence on f L of both the cooling rate and the dispersion of the total energy, i.e., the exponents in the power laws (14) and (17), was in agreement with the theoretical predictions. This was illustrated in Figs. 1 and 2 of Ref. [3] and no more details will be given here. The comparison between the predicted critical amplitudes and the MD results given in Table I can be considered as satisfactory, in the sense that the theory correctly predicts the order of magnitude of the amplitudes, especially taking into account the smallness of the quantities being measured. Next, let us proceed to investigate the form of the probability distribution of the energy fluctuations. Particularization of Eq. (4) for the modes with the smallest possible value of kin the limit f L 1gives @sf L!k?sk?s;(18) where it is understood that jkjkmin. Define a new time scale d f Lds, and a new transversal velocity field by ! k?!k?=Ld c1=2 E. Equation (18) becomes @1! k? k?;(19) with h k? k0?0i  d1=2 6d11=2k;k00l:(20) Equation (19) implies that the probability distribution for ! k?with jkjkmin near the clustering instability depends only on the dimension dof the system. In fact, since the noise term  k?is Gaussian, it is trivial to write the long time form of this distribution using Eq. (7) with ss0, Pst! k?22 !d1=2e!2 k?=22 !;(21) with 2 !d1=2=12d11=2. In the time scale , and keeping only the dominant modes, Eq. (10) reads f L@y1 2y6 dX jkjkminj! ?kj2;(22) where y=Eand the sum is restricted to vectors kwith jkjkmin. From the comparison of Eqs. (19) and (22) it is seen that, on the scale and in the threshold of the instability, ydecays much faster than the dominant components of ! k?. Consequently, for large the solution of Eq. (22) is y6 dPjkjkmin j! ?kj2, where the probability distribution of the modes ! ?kis given by Eq. (21). Since the latter does not depend on the parameters of the system other than the dimensionality, the same property follows for the probability distribution of both yand the variable f E Edd11=26 dX jkjkminj! k?j2:(23) TABLE I. Comparison between the predicted and MD values for the critical amplitudes of the cooling rate Aand the total energy dispersion A. All the values of the amplitudes have been multiplied by 103. nH2A theory AMD Atheory AMD  0.02 0.9 0.88 1.11 0.62 0.6 0.02 0.8 1.62 3.63 1.15 1.5 0.1 0.98 1.06 0.50 0.75 0.47 0.1 0.95 2.4 2.38 1.7 1.45 0.2 0.98 1.97 2.34 1.4 1.34 0.2 0.95 4.59 7.78 3.24 3.6 PRL 96, 158002 (2006) PHYSICAL REVIEW LETTERS week ending 21 APRIL 2006 158002-3 This is equivalent to saying that the probability distribution for f E verifies the scaling relation Pf E 1 E ff E E;(24) where fis a scaling function. This is just the property assumed in Ref. [3] and verified by MD simulations. Since the probability distribution function for ! k?is known, it is possible to numerically generate the probability distribution function Pf E. The result is shown in Fig. 1. Also plotted is the function f EKexexa;xbsf E;a=2; (25) with K2:14,b0:938, and s0:374, that fits extremely well the MD results for EPf E[3]. It is important to remark that fluctuations in a large number of equilibrium and nonequilibrium systems exhibiting selforganized criticality as well as confined turbulent flows present the same kind of behavior [4,5]. Although the agreement between both plotted curves is not so bad for positive values of f E, strong discrepancies are observed for negative values. A major source for them is easily identified from Eq. (23), that for d2implies f E=E 2 p, while smaller values are found in the MD simulations. Since Eq. (23) is a consequence of Eq. (9), it seems plausible that in order to elaborate a more accurate theory the intrinsic noise associated with the cooling rate must be taken into account. In summary, we have developed a mesoscopic theory for the fluctuations of the total energy of an isolated granular gas near the threshold of the clustering instability. The theory describes accurately the qualitative behavior obtained in MD simulations, namely, the divergent behavior of the dimensionless second moment and the decrease of the apparent cooling rate. Also, it is consistent with the observed scaling property of the probability distribution function of the fluctuations. On the other hand, it seems clear that a more refined formulation is needed in order to get a more satisfactory quantitative agreement, especially for the distribution function. Also, it should be interesting to check whether a similar behavior occurs in more realistic models of granular gases in which depends on the relative velocity and, although present, the clustering instability seems to be a transient phenomenon [13]. This research was supported by the Ministerio de Educacio ´n y Ciencia (Spain) through Grant No. FIS2005-01398 (partially financed by FEDER funds). *Electronic address: [email protected] [1] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Rev. Mod. Phys. 68, 1259 (1996). [2] I. Goldhirsch, Annu. Rev. Fluid Mech. 35, 267 (2003). [3] J. J. Brey, M. I. Garcı ´a de Soria, P. Maynar, and M. J. Ruiz- Montero, Phys. Rev. Lett. 94, 098001 (2005). [4] S. T. Bramwell, P. C. W. Holdsworth, and J.-F. Pinton, Nature (London) 396, 552 (1998). [5] S. T. Bramwell, K. Christensen, J.-Y. Fortin, P. C. W. Holdsworth, H. J. Jensen, S. Lise, J. M. Lo ´pez, M. Nicodemi, J.-F. Pinton, and M. Sellitto, Phys. Rev. Lett. 84, 3744 (2000). [6] L. Landau and E. M. Lifshitz, Fluid Mechanics (Peragamon Press, New York, 1959). [7] R. Brito and M. H. Ernst, Europhys. Lett. 43, 497 (1998). [8] P. K. Haff, J. Fluid Mech. 134, 401 (1983). [9] T. P. C. van Noije, M. H. Ernst, R. Brito, and J. A. G. Orza, Phys. Rev. Lett. 79, 411 (1997). [10] I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. 70, 1619 (1993); S. McNamara and W. R. Young, Phys. Rev. E 50, R28 (1994). [11] J. J. Brey, M. J. Ruiz-Montero, and D. Cubero, Phys. Rev. E 60, 3150 (1999). [12] J. J. Brey, M. I. Garcı ´a de Soria, P. Maynar, and M. J. Ruiz- Montero, Phys. Rev. E 70, 011302 (2004). [13] N. Brilliantov, C. Saluen ˜a, T. Schwager, and T. Po ¨schel, Phys. Rev. Lett. 93, 134301 (2004). -6 -4 -2 0 -2 0 2 4 ln [σE P(δ E)] δ E ~ ~ / σE FIG. 1. Probability density function of the relative total energy fluctuations EPf Efor a system of inelastic hard disks. The broken line is the theoretical prediction derived in this Letter and the solid line Eq. (7). PRL 96, 158002 (2006) PHYSICAL REVIEW LETTERS week ending 21 APRIL 2006 158002-4