Energy fluctuations in the homogeneous cooling state of granular gases

Energy fluctuations in the homogeneous cooling state of granular gases J. Javier Brey, M. I. García de Soria, P. Maynar, and M. J. Ruiz-Montero Física Teórica, Universidad de Sevilla, Apdo. de Correos 1065, E-41080 Sevilla, Spain (Received 25 March 2004; published 29 July 2004) Starting from the hierarchy of equations for microscopic densities in phase space, a general theory for fluctuations and correlations in a dilute granular gas of hard particles is developed. Then, the particular case of the homogeneous cooling state is addressed. Explicit expressions for some distributions describing the presence of velocity correlations and their dynamics are obtained. These correlations are inherent to the dissipative dynamics of the collisions. The implications for the behavior of the total energy of the system are analyzed and the results are expressed in terms of a fluctuation-dissipation theorem. The theoretical predictions are shown to be in agreement with results obtained by molecular dynamics simulations, which also indicate that energy fluctuations are well described by a Gaussian distribution. DOI: 10.1103/PhysRevE.70.011302 PACS number(s): 45.70.⫺n, 51.10.⫹y, 05.20.Dd I. INTRODUCTION The understanding of the microscopic and macroscopic physical mechanisms involved in rapid granular flows has experimented a great advance in the last years [1,2]. Kinetic theory and nonequilibrium statistical mechanics methods have been extended to include inelasticity of collisions, which is the main feature of grain interactions and is associated with their macroscopic character. In these studies, the most widely employed idealized model for granular fluids is a system of smooth hard spheres (or disks in two dimensions)whose collisions are characterized by a constant coefficient of normal restitution. This coefficient measures the decrease in magnitude of the normal component of the relative velocity of the two colliding particles. Starting from the dynamics of the particles, it is possible to derive the corresponding Liouville equation providing the basis for the application of many-body methods [3,4].In particular, the inelastic Boltzmann equation describing the time evolution of the one-body distribution function is obtained in the dilute limit. This kinetic equation has been extensively used to address many fundamental issues and to compute averages of dynamical properties. Primary among these are the hydrodynamic equations, with explicit expressions for the transport coefficients, which have been derived by using a generalization of the Chapman-Enskog algorithm [5,6]. More recently, the hydrodynamic description has been identified in the context of the eigenfunctions and eigenvalues of the linearized inelastic Boltzmann collision operator [7,8]. To put these results in a proper context, it must be pointed out that there is no proof of time scale separation in rapid granular flows. At the level of the Boltzmann equation, this would require that the rest of the spectrum of the linearized operator is bounded away from the hydrodynamic modes, having a larger real part. This property has only been established for some simplified model Boltzmann collision operators [8,9]. On the other hand, granular hydrodynamics has been applied with success to many different situations, and the Chapmann-Enskog expressions for the transport coefficients have been confirmed by direct Monte Carlo simulations [10]. Most of the work carried out up to date has focused on kinetic equations for the one-body distribution and the information following from them, while the knowledge about fluctuations and correlations in granular gases is much more limited. First studies were based on mesoscopic fluctuating hydrodynamics [11]and also on model equations [12], and focused on the initial buildup of spatial correlations in the development of the clustering instability starting from an homogeneous system. Precollisional velocity correlations have been analyzed in detail at a fundamental level starting from the microscopic equations of motion and considering correlated sequences of collisions [13]. The system was initially in a state of thermal equilibrium and the results are valid at the short time and length scales. Then, although evidently relevant, they do not provide too much insight into the structure of the velocity correlations in any particular macroscopic state of the granular fluid. The main goal of this paper is to formulate a theory of fluctuations for dilute granular gases and to apply it to one of the simplest possible situations, exploiting the scaling properties of the particular state under consideration. This will allow us to obtain explicit expressions for some distributions characterizing velocity fluctuations in the system. Special attention will be paid to the possibility of expressing the results in the form of a macroscopic theory of fluctuations, including fluctuation-dissipation relations. In order to develop a theory for equal and different time fluctuations, the hierarchy method will be used [13,14]. The general idea of the method is to derive from the Liouville equation hierarchies of coupled equations for the distribution functions describing the fluctuations. Then the hierarchies are closed by using the same kind of approximations as needed to derive the corresponding kinetic equation, in our case the Boltzmann equation. In this way, a unified formalism provides the usual kinetic equation as well as evolution equations for the one-time and two-time correlations. An isolated fluid with elastic collisions tends to the equilibrium Gibbs state. The single particle distribution function is Maxwellian and there are no correlations between the velocities of the particles. Equilibrium fluctuations are well understood since long ago [15,16]. For inelastic collisions, there is no stationary state for an isolated system. The simplest state is the so-called “homogeneous cooling state” (HCS)with a monotonically decreasing temperature. Of PHYSICAL REVIEW E 70, 011302 (2004) 1539-3755/2004/70(1)/011302(12)/$22.50 ©2004 The American Physical Society70 011302-1 course, this is a consequence of the dissipation of energy in collisions. For the same reason, the single particle distribution function deviates from the Maxwellian, both at small and large velocities [17–19]. Here, we will investigate a third implication of the dissipative character of collisions: the presence of velocity correlations, even in the low density limit. It must be realized that the existence of velocity correlations does not invalidate by itself the Boltzmann equation, as long as they are small enough, since they are obtained within a self-consistent theoretical frame, expected to be valid in the low density limit. The study of fluctuations in the HCS is of primary importance for the development of a general theory of fluctuations in granular systems. The HCS defines the reference state from which macroscopic hydrodynamic equations can be derived [5], thus playing in this context a similar role to the equilibrium state in molecular systems. In the same way, a general theory of hydrodynamic fluctuations requires, as a necessary first step, the knowledge of the fluctuations and correlations occurring in the HCS. It should be stressed that no conceptual difficulty arises from using the HCS, which is a time-dependent state, as a reference state. In fact, the reference state is determined by the kinetic equation itself and cannot be arbitrarily chosen [20]. In the next section, hierarchies of equations for the onetime and two-time reduced distribution and correlation functions are derived. The closure relations obtained in the low density limit, leading in particular to the Boltzmann equation, are also discussed. The specific case of the HCS is considered in Sec. III, where the scaling property of this state is used to simplify the equations and eliminate trivial time dependencies. With an appropriate change of variables, the equations looks like corresponding to a stationary, but non- Gibbs, state. In Secs. IV and V, marginal distributions for the one-time and two-time velocity correlations are introduced and their hydrodynamic parts are evaluated. The derived expressions are then used to obtain the second moment of the total energy fluctuations and the two-time energy correlation function. The theoretical predictions are compared with molecular dynamics simulation results and a quite good agreement is found over a wide range of values of the restitution coefficient. Moreover, the simulation data show that the scaled energy fluctuations are Gaussian within the numerical errors. Finally, Sec. VI contains a summary of the main results, and also a discussion of their connection with a fluctuation-dissipation relation. II. EVOLUTION EQUATIONS FOR THE REDUCED DISTRIBUTION FUNCTIONS We consider a dilute gas of smooth inelastic hard spheres 共d=3兲or disks 共d=2兲of mass mand diameter ␴ . Let Xj共t兲 ⬅兵Rj共t兲,Vj共t兲其 denote the position and velocity of particle j at time t. Both Rj共t兲and Vj共t兲are parametric functions of the initial positions and velocities of all particles. For t⬎0 the time evolution of a phase function A关兵Xj共t兲其兴 is defined by [3,4] A关⌫共t兲兴 =etL共⌫兲A关⌫共0兲兴,共1兲 where ⌫⬅兵Xj;j=1,...,N其denotes the N-particle phase point. The generator of the dynamics Lis L共⌫兲=L共0兲共⌫兲+1 2兺 j N 兺 i⫽j N T共Xi,Xj兲,共2兲 with the first term generating free streaming, L共0兲=兺 j=1 N Vj· ⳵ ⳵ Rj,共3兲 and the second terms describing velocity changes in collisions, T共Xi,Xj兲= ␴ d−1 冕 d ␴ ˆ⌰共−Vij · ␴ ˆ兲兩Vij · ␴ ˆ兩 ␦ 共Rij − ␴ 兲 ⫻关b ␴ 共Vi,Vj兲−1兴.共4兲 Here ␴ = ␴␴ ˆ,d ␴ ˆis the solid angle element for ␴ ˆ,Rij⬅Ri −Rj,Vij⬅Vi−Vj,⌰is the Heaviside step function, and the operator b ␴ 共Vi,Vj兲changes functions of Vi,Vjto the same functions of the scattered velocities given by Vi→Vi ⬘⬅b ␴ Vi=Vi−1+ ␣ 2共 ␴ ˆ·Vij兲 ␴ ˆ, Vj→Vj ⬘⬅b ␴ Vj=Vj+1+ ␣ 2共 ␴ ˆ·Vij兲 ␴ ˆ.共5兲 Inelasticity of collisions is characterized by means of the coefficient of normal restitution ␣ taking values in the interval 0⬍ ␣ 艋1. Microscopic densities in phase space Fs共x1,x2,...,xs,t兲 are defined by F1共x1,t兲=兺 j=1 N ␦ 共x1−Xj共t兲兲,共6兲 F2共x1,x2,t兲=兺 i N 兺 j⫽i N ␦ 共x1−Xi共t兲兲 ␦ „x2−Xj共t兲…,共7兲 etc. The lower-case variables xi⬅兵ri,vi其are field variables referring to a particular point in ␮ space. A hierarchy of equations for the above densities is obtained from Eq. (1). The first two equations, valid for t⬎0, are 冋 ⳵ ⳵ t+L共0兲共x1兲 册 F1共x1,t兲= 冕 dx2T ¯ 共x1,x2兲F2共x1,x2,t兲,共8兲 冋 ⳵ ⳵ t+L共0兲共x1,x2兲 册 F2共x1,x2,t兲=T ¯ 共x1,x2兲F2共x1,x2,t兲 + 冕 dx3关T ¯ 共x1,x3兲+T ¯ 共x2,x3兲兴F3共x1,x2,x3,t兲,共9兲 where L共0兲共x1,x2兲=L共0兲共x1兲+L共0兲共x2兲,L共0兲共xi兲=vi· ⳵ ⳵ ri, 共10兲 and T ¯ 共xi,xj兲is the binary collision operator BREY et al. PHYSICAL REVIEW E 70, 011302 (2004) 011302-2 T ¯ 共xi,xj兲= ␴ d−1 冕 d ␴ ˆ⌰共vij · ␴ ˆ兲兩vij · ␴ ˆ兩 ⫻关 ␣ −2 ␦ 共rij − ␴ 兲b ␴ −1共vi,vj兲− ␦ 共rij + ␴ 兲兴. 共11兲 The operator b ␴ −1共vi,vj兲replaces all the velocities viand vj appearing to its right by the precollisional values v1 *and v2 *, vi→vi *⬅b ␴ −1vi=vi−1+ ␣ 2 ␣ 共 ␴ ˆ·vij兲 ␴ ˆ, vj→vj *⬅b ␴ −1vj=vj+1+ ␣ 2 ␣ 共 ␴ ˆ·vij兲 ␴ ˆ.共12兲 The averages of the microscopic densities Fs共x1,x2,¯,xs,t兲over the probability distribution function ␳ 共⌫,0兲characterizing the initial state of the system are the usual one-time reduced distribution functions, fs共x1,...,xs,t兲⬅具Fs共x1,¯,xs,t兲典,共13兲 where we have introduced the notation 具G典⬅ 冕 d⌫G共⌫兲 ␳ 共⌫,0兲.共14兲 The one-time reduced distributions obey a hierarchy of equations that follow directly by averaging the hierarchy for the microscopic densities. This is the generalization for inelastic collisions of the BBGKY hierarchy. In particular, from Eqs. (8)and (9)it is trivially obtained [3]: 冋 ⳵ ⳵ t+L共0兲共x1兲 册 f1共x1,t兲= 冕 dx2T ¯ 共x1,x2兲f2共x1,x2,t兲, 共15兲 冋 ⳵ ⳵ t+L共0兲共x1,x2兲 册 f2共x1,x2,t兲=T ¯ 共x1,x2兲f2共x1,x2,t兲 + 冕 dx3关T ¯ 共x1,x3兲+T ¯ 共x2,x3兲兴f3共x1,x2,x3,t兲.共16兲 Two-time reduced distribution functions can also be defined in terms of the microscopic densities as fr,s共x1,...,xr,t;x1 ⬘,...,xs ⬘,t⬘兲⬅具Fr共x1,...,xr,t兲 ⫻Fs共x1 ⬘,...,xs ⬘,t⬘兲典. 共17兲 It will be assumed that t⬎t⬘⬎0 for concreteness. Again, evolution equations for these quantities are derived from the hierarchy for the densities Fr. The lowest order distribution, f1,1, is seen to obey the equation 冋 ⳵ ⳵ t+L共0兲共x1兲 册 f1,1共x1,t;x1 ⬘,t⬘兲 = 冕 dx2T ¯ 共x1,x2兲f2,1共x1,x2,t;x1 ⬘,t⬘兲,共18兲 which has to be solved with the initial condition f1,1共x1,t⬘;x1 ⬘,t⬘兲= ␦ 共x1−x1 ⬘兲f共x1,t⬘兲+f2共x1,x2,t⬘兲,共19兲 that follows directly from the definition in Eq. (17). With regard to the above equations for the reduced distribution functions, it must be reminded that although the streaming operators also generate unphysical trajectories for unphysical initial conditions, the latter have a vanishing weight in the initial probability distribution and, therefore, the dynamics in all the equations for the reduced distribution functions is in fact well behaved. It is convenient to introduce also correlation functions through the usual cluster expansions. From the one-time reduced distributions, one-time correlations gs共x1,...,xs,t兲are defined by f2共x1,x2,t兲=f1共x1,t兲f1共x2,t兲+g2共x1,x2,t兲,共20兲 f3共x1,x2,x3,t兲=f1共x1,t兲f1共x2,t兲f1共x3,t兲+f1共x1,t兲g2共x2,x3,t兲 +f1共x2,t兲g2共x1,x3,t兲+f1共x3,t兲g2共x1,x2,t兲 +g3共x1,x2,x3,t兲,共21兲 etc. Similarly, two-time correlation functions hr,scan be defined. In particular, h1,1 and h2,1 are introduced through f1,1共x1,t;x1 ⬘,t⬘兲=f1共x1,t兲f1共x1 ⬘,t⬘兲+h1,1共x1,t;x1 ⬘,t⬘兲, 共22兲 f2,1共x1,x2,t;x1 ⬘,t⬘兲=f1共x1,t兲f1共x2,t兲f1共x1 ⬘,t⬘兲 +g2共x1,x2,t兲f1共x1 ⬘,t⬘兲 +h1,1共x1,t;x1 ⬘,t⬘兲f1共x2,t兲 +h1,1共x2,t;x1 ⬘,t⬘兲f1共x1,t兲 +h2,1共x1,x2,t;x1 ⬘,t⬘兲.共23兲 Evolution equations for the correlation functions are obtained from the equations for the corresponding reduced distributions. Thus from Eqs. (15)and (16)one gets 冋 ⳵ ⳵ t+L共0兲共x1兲 册 f1共x1,t兲= 冕 dx2T ¯ 共x1,x2兲关f1共x1,t兲f1共x2,t兲 +g2共x1,x2,t兲兴,共24兲 冋 ⳵ ⳵ t+L共0兲共x1,x2兲 册 g2共x1,x2,t兲 =T ¯ 共x1,x2兲关f1共x1,t兲f1共x2,t兲+g2共x1,x2,t兲兴 + 冕 dx3共1+P12兲T ¯ 共x1,x3兲关g2共x1,x2,t兲f1共x3,t兲 +g2共x2,x3,t兲f1共x1,t兲+g3共x1,x2,x3,t兲兴,共25兲 where the permutation operator P12 interchanges the labels of particles 1 and 2 in the quantities to its right. In an analogous way, from Eq. (18)it results ENERGY FLUCTUATIONS IN THE HOMOGENEOUS…PHYSICAL REVIEW E 70, 011302 (2004) 011302-3 冋 ⳵ ⳵ t+L共0兲共x1兲 册 h1,1共x1,t;x1 ⬘,t⬘兲 = 冕 dx2T ¯ 共x1,x2兲关h1,1共x1,t;x1 ⬘,t⬘兲f1共x2,t兲 +h1,1共x2,t;x1 ⬘,t⬘兲f1共x1,t兲+h2,1共x1,x2,t;x1 ⬘,t⬘兲兴. 共26兲 The last equation holds for t⬎t⬘and the initial condition is h1,1共x1,t⬘;x1 ⬘,t⬘兲=g2共x1,x1 ⬘,t⬘兲+ ␦ 共x1−x1 ⬘兲f1共x1,t⬘兲. 共27兲 All the equations derived in this section are exact consequences of the dynamics of the particles as defined by the generator given in Eq. (2). In fact, they are the extension to inelastic hard particles of the equations discussed in Ref. [14]. In the next sections, the low density limit and the particular case of the HCS will be addressed. III. LOW DENSITY FLUCTUATIONS AROUND THE HOMOGENEOUS COOLING STATE In order to obtain closed kinetic equations for the reduced distribution functions, some kind of closure relations for the hierarchies of equations derived in the previous section are required. Here we are interested in the limit of a dilute inelastic gas and the method we will apply is based on the general assumption that the order in the density of correlation functions increases as the number of involved particles increases. More precisely, what it is assumed is that the above is true for the precollisional part of the correlation functions when two of the particles are at contact [21]. This method was discussed in detail by Ernst and Cohen for the hard sphere gas with elastic collisions [14]and it has been already applied to the study of one-time reduced distributions for inelastic hard particles [4,22]. Closely related results follow from a formal small parameter expansion of the reduced distribution functions [20,23]. From a physical point of view, the underlying idea of these kinetic theories is that no hypothesis are necessary to derive equations for fluctuations and time-correlations beyond those required to derive the kinetic equation for the one-particle distribution function, i.e., the Boltzmann equation in the dilute gas. Let us stress that this seems to be a quite general property, beyond the range of validity of the Boltzmann equation, as emphasized in Ref. [24]. Then, in the low density limit, the term involving T ¯ 共x1,x2兲g2共x1,x2,t兲in Eq. (24)is assumed to be negligible as compared with the one involving T ¯ 共x1,x2兲f1共x1,t兲f1共x2,t兲. Furthermore, in the same limit the spatial separation between the centers of colliding particles can be neglected, so that the above mentioned equation reduces to 冋 ⳵ ⳵ t+L共0兲共x1兲 册 f1共x1,t兲=J关x1,t兩f1兴,共28兲 where J关x1,t兩f1兴= 冕 dx2 ␦ 共r12兲T ¯ 0共v1,v2兲f1共x1,t兲f1共x2,t兲,共29兲 T ¯ 0共v1,v2兲= ␴ d−1 冕 d ␴ ˆ⌰共v12 · ␴ ˆ兲兩v12 · ␴ ˆ兩 ⫻关 ␣ −2b ␴ −1共v1,v2兲−1兴.共30兲 This is the Boltzmann equation for smooth inelastic hard spheres [3,17]. Similarly, Eq. (25)reduces to 冋 ⳵ ⳵ t+L共0兲共x1,x2兲−K关x1,t兩f1兴−K关x2,t兩f1兴 册 g2共x1,x2,t兲 = ␦ 共r12兲T ¯ 0共v1,v2兲f1共x1,t兲f1共x2,t兲,共31兲 where we have introduced the operator K关xi,t兩f1兴⬅ 冕 dx3 ␦ 共ri3兲T ¯ 0共vi,v3兲共1+Pi3兲f1共x3,t兲. 共32兲 Therefore, two particle correlations are generated from the uncorrelated product of solutions of the Boltzmann equation through inelastic binary collisions. Of course, the solutions of this equation must be consistent with the approximation leading to the Boltzmann equation. The low density limit of the equation for the two-time correlations, Eq. (26), follows by neglecting the term proportional to h2,1 on the right-hand side and also the distance between colliding particles, so that it is given by 冋 ⳵ ⳵ t+L共0兲共x1兲−K关x1,t兩f1兴 册 h1,1共x1,t,x1 ⬘,t⬘兲=0. 共33兲 Finally, an equation for h1,1共x1,t;x1 ⬘,t兲is obtained by taking time derivative in Eq. (27)and using Eqs. (28)and (31), 冋 ⳵ ⳵ t+L共0兲共x1,x1 ⬘兲−K关x1,t兩f1兴−K关x1 ⬘,t兩f1兴 册 h1,1共x1,t;x1 ⬘,t兲 =⌫共x1,x1 ⬘,t兲,共34兲 with ⌫共x1,x1 ⬘,t兲= ␦ 共x1−x1 ⬘兲J关x1,t兩f1兴−共K关x1,t兩f1兴 +K关x1 ⬘,t兩f1兴兲 ␦ 共x1−x1 ⬘兲f1共x1,t兲 + ␦ 共r1−r1 ⬘兲T ¯ 0共v1,v1 ⬘兲f1共x1,t兲f1共x1 ⬘,t兲.共35兲 It is worth pointing out that the above derivation makes no reference to whether collisions are elastic or inelastic. In fact, when it is compared with the analysis carried out in Ref. [14], the only difference is in the replacement of the elastic binary collision operators by the inelastic ones and this does not affect the density dependence of each of the contributions to the equations. This is also the case when the kinetic equations are obtained from a formal expansion in a dimensionless density parameter [20,23]. Therefore, at least at a superficial level, it seems that the theoretical basis for low density kinetic theory is the same both for elastic and inelastic systems. Of course, the actual range of validity of the kinetic equations, i.e., how small the density of the system must be so that the above kinetic equations provide an accurate description, might depend on the inelasticity of the system, and also on the particular state being considered. BREY et al. PHYSICAL REVIEW E 70, 011302 (2004) 011302-4 The inelastic Boltzmann equation (28)admits a special solution fHCS共v,t兲describing the HCS for which all the time dependence of the distribution function occurs through the granular temperature (defined from the average kinetic energy). This solution has the scaled form [17] fHCS共v,t兲=nHv0 −d共t兲 ␹ HCS共c兲.共36兲 Here nHis the uniform number of particles density, v0共t兲⬅ 冋 2kBTHCS共t兲 m 册 1/2,共37兲 kBbeing the Boltzmann constant, and ␹ HCS共c兲an isotropic function of scaled velocity c=v/v0共t兲. The above scaling property implies that the temperature obeys the equation ⳵ THCS共t兲 ⳵ t+ ␨ HCS共t兲THCS共t兲=0. 共38兲 The cooling rate ␨ HCS共t兲is a given bilinear functional of fHCS, and it is proportional to THCS共t兲1/2. Its explicit expression will not be given here [7,17]. Substitution of Eq. (36) into the Boltzmann equation (28)and use of Eq. (38)leads to a closed equation for ␹ HCS共c兲, ␨ 0 2 ⳵ ⳵ c·共c ␹ HCS兲=Jc关c兩 ␹ HCS兴,共39兲 where ␨ 0=ᐉ ␨ HCS共t兲 v0共t兲 =共1− ␣ 2兲 ␲ 共d−1兲/2 2⌫ 冉 d+3 2 冊 d 冕 dc1 冕 dc2c12 3 ␹ HCS共c1兲 ␹ HCS共c2兲 共40兲 is the time-independent dimensionless cooling rate of the HCS and Jc关c兩 ␹ HCS兴= 冕 dc1T ¯ 0共c,c1兲 ␹ HCS共c兲 ␹ HCS共c1兲,共41兲 T ¯ 0共c,c1兲= 冕 d ␴ ˆ⌰关共c−c1兲· ␴ ˆ兴共c−c1兲· ␴ ˆ关 ␣ −2b ␴ −1共c,c1兲−1兴. 共42兲 The operator b ␴ −1共c,c1兲is again defined by Eqs. (12), but replacing the velocities viand vjby cand c1, respectively. It is useful to introduce dimensionless time and length scales by s= 冕 0 tdt1 v0共t1兲 ᐉ,l=ᐉ−1r,共43兲 respectively, where ᐉ⬅共nH ␴ d−1兲−1 is proportional to the mean free path of the gas. Thus, s共t兲is a measure of the average number of collisions per particle in the interval 共0,t兲. In terms of s, Eq. (38)becomes ⳵ THCS ⳵ s=− ␨ 0THCS共s兲,共44兲 showing that the temperature of the HCS decreases exponentially on this time scale. The equivalent expression in the original time scale is given by Haff’s law [25], THCS共t兲=THCS共0兲 关 1+1 2 ␨ HCS共0兲t 兴 −2.共45兲 In order to particularize the equations for the two-particle correlations, the concept of the HCS as introduced above will be extended. It will be assumed that in the HCS the N-particle distribution function of the system, ␳ HCS共⌫,t兲,is given by a scaling function for which all the time dependence occurs through the scaling of the velocities with the thermal velocity v0共t兲[3], ␳ HCS共⌫,t兲=关v0共t兲兴−dN ␳ HCS *共兵rij,ci其兲.共46兲 The scaled distribution ␳ HCS *is invariant under space translations. Therefore, ␳ HCS共⌫,t兲represents a spatially homogeneous fluid. The above scaling implies a similar scaling for all the reduced distribution functions and, in particular, the one given by Eq. (36)for the one-particle distribution function. Then, Eq. (31)when applied to the HCS becomes 冋 c12 · ⳵ ⳵ l12 −⌳共c1兲−⌳共c2兲 册 g ˜ HCS共l12,c1,c2兲 = ␦ 共l12兲T ¯ 0共c1,c2兲 ␹ HCS共c1兲 ␹ HCS共c2兲,共47兲 where ⌳共ci兲is the linearized Boltzmann collision operator in the scaled representation, ⌳共ci兲= 冕 dc3T ¯ 0共ci,c3兲共1+Pi3兲 ␹ HCS共c3兲− ␨ 0 2 ⳵ ⳵ ci·ci, 共48兲 and the scaled two-particle one-time correlation function is defined by g ˜ HCS共l12,c1,c2兲=nH −1ᐉdv0 2d共t兲g2,HCS共r12,v1,v2,t兲.共49兲 This quantity does not depend on sas a consequence of the scaling property. In a similar way, Eq. (33)yields 冋 ⳵ ⳵ s+c1· ⳵ ⳵ l1−⌳共c1兲 册 h ˜ HCS共l1−l1 ⬘,c1,s−s⬘;c1 ⬘兲=0, 共50兲 with h ˜ HCS共l1−l1 ⬘,c1,s−s⬘;c1 ⬘兲=nH −1ᐉdv0 d共t兲v0 d共t⬘兲 ⫻h1,1,HCS共x1,t;x1 ⬘,t⬘兲.共51兲 The initial condition for this equation is h ˜ HCS共l1−l1 ⬘,c1,0;c1 ⬘兲⬅h ˜ HCS共l1−l1 ⬘,c1;c1 ⬘兲 =g ˜ HCS共l1−l1 ⬘,c1,c1 ⬘兲 + ␦ 共c1−c1 ⬘兲 ␦ 共l1−l1 ⬘兲 ␹ HCS共c1兲. 共52兲 ENERGY FLUCTUATIONS IN THE HOMOGENEOUS…PHYSICAL REVIEW E 70, 011302 (2004) 011302-5 In the above expressions, we have made explicit the translational invariance of the HCS. Finally, Eq. (34)becomes 冋 c1· ⳵ ⳵ l1+c1 ⬘· ⳵ ⳵ l1 ⬘−⌳共c1兲−⌳共c1 ⬘兲 册 h ˜ HCS共l1−l1 ⬘,c1;c1 ⬘兲 = ␦ 共l1−l1 ⬘兲⌫ ˜ 共c1,c1 ⬘兲,共53兲 ⌫ ˜ 共c1,c1 ⬘兲=−关⌳共c1兲+⌳共c1 ⬘兲兴 ␦ 共c1−c1 ⬘兲 ␹ HCS共c1兲 +T ¯ 0共c1,c1 ⬘兲 ␹ HCS共c1兲 ␹ HCS共c1 ⬘兲.共54兲 Equations (47)and (53)describe the correlations between equal time fluctuations. Their solutions are related through Eq. (52). In both equations, the distribution function of the HCS is needed and has to be found from Eq. (39). The same happens with Eq. (50)for the unequal time fluctuations. Nevertheless, the analytical form of ␹ HCS共c兲is only partially known. In particular, an expansion about the Gaussian in Sonine polynomials has been used. To first order, ␹ HCS共c兲is approximated by ␹ HCS共c兲=e−c2 ␲ d/2关1+a2共 ␣ 兲S共2兲共c2兲兴,共55兲 where S共2兲共c2兲=c4 2−d+2 2c2+d共d+2兲 8.共56兲 The coefficient a2is related with the fourth moment of ␹ HCS共c兲by a2共 ␣ 兲=4具c4典 d共d+2兲−1; 具c4典⬅ 冕 dcc4 ␹ HCS共c兲.共57兲 Substitution of Eq. (55)into Eq. (39)leads to a closed equation for the coefficient a2, and by neglecting nonlinear terms in a2one gets [17,18] a2共 ␣ 兲=16共1− ␣ 兲共1−2 ␣ 2兲 9+24d+共8d−41兲 ␣ +30 ␣ 2−30 ␣ 3.共58兲 Using Eq. (55), the following approximate expression for the dimensionless cooling rate is found: ␨ 0=冑2 ␲ 共d−1兲/2共1− ␣ 2兲 ⌫ 冉 d 2 冊 d 冋 1+ 3 16a2共 ␣ 兲 册 .共59兲 The above expressions have been shown to agree with numerical results obtained by means of the direct simulation Monte Carlo method in the thermal velocity region, i.e. for velocities cof the order of a few units [26]. Let us mention that the shape of the scaling distribution has also been analyzed for large velocities and found to have an exponential form [18,19]. IV. ENERGY FLUCTUATIONS IN THE HCS To study the effect of velocity fluctuations on the behavior of the total energy of the system in the HCS, we introduce a marginal velocity correlation function by ␸ HCS共c1,c2兲⬅ 冕 dl12h ˜ HCS共l12,c1;c2兲.共60兲 An equation for this distribution is obtained by integrating Eq. (53), 关⌳共c1兲+⌳共c2兲兴 ␸ HCS共c1,c2兲=−⌫ ˜ 共c1,c2兲.共61兲 From the definition of g2in Eq. (20)and the scaling of the several distributions defined in the previous section it follows that 冕 dc1 ␸ HCS共c1,c2兲= 冕 dc2 ␸ HCS共c1,c2兲=0. 共62兲 Instead of trying to find the complete solution of Eq. (61), we will restrict ourselves to compute what we will call the hydrodynamic part of ␸ HCS共c1,c2兲, defined as follows. We consider the eigenvalue problem associated with the linearized homogeneous Boltzmann operator [7,8] ⌳共c兲 ␰ i共c兲=␭i ␰ i共c兲.共63兲 The hydrodynamic part of the spectrum of ⌳is defined by those eigenvalues that coincide with the eigenvalues of the balance equations for the number density, momentum, and temperature following from the homogeneous linearized Boltzmann equation. Such eigenvalues are ␭1=0, ␭2= ␨ 0 2,␭3=− ␨ 0 2.共64兲 The corresponding eigenfunctions of ⌳for each of the above eigenvalues are found to be [7,8] ␰ 1共c兲= ␹ HCS共c兲+ ⳵ ⳵ c·关c ␹ HCS共c兲兴, ␰ 2共c兲=− ⳵␹ HCS共c兲 ⳵ c, ␰ 3共c兲=− ⳵ ⳵ c·关c ␹ HCS共c兲兴,共65兲 respectively. The eigenvalue ␭2is, therefore, d-fold degenerated. The scalar product of two functions f共c兲and g共c兲is defined as 具f兩g典⬅ 冕 dc ␹ HCS −1 共c兲f*共c兲g共c兲,共66兲 with f*being the complex conjugate of f. The eigenfunctions ␰ ␤ given in Eq. (65)are not orthogonal, as a consequence of the operator ⌳being non-Hermitian in the associate Hilbert space. On the other hand, it is easily verified that the set of functions ␰ ¯ 1共c兲= ␹ HCS共c兲, ␰ ¯ 2共c兲=c ␹ HCS共c兲, ␰ ¯ 3共c兲= 冉 c2 d+1 2 冊 ␹ HCS共c兲共67兲 verify the biorthogonality condition BREY et al. PHYSICAL REVIEW E 70, 011302 (2004) 011302-6 具 ␰ ¯ ␤ 兩 ␰ ␤ ⬘典= ␦ ␤ , ␤ ⬘,共68兲 ␤ , ␤ ⬘=1,2,3. The definition of the scalar product given above is trivially generalized to two-velocity functions by 具f共c1,c2兲兩g共c1,c2兲典 ⬅ 冕 dc1 冕 dc2 ␹ HCS −1 共c1兲 ␹ HCS −1 共c2兲 ⫻f*共c1,c2兲g共c1,c2兲.共69兲 Then, a projector operator Pis defined as Pf共c1,c2兲=兺 ␤ =1 3 兺 ␤ ⬘=1 3 ␰ ␤ 共c1兲 ␰ ␤ ⬘共c2兲具 ␰ ¯ ␤ 共c1兲 ␰ ¯ ␤ ⬘共c2兲兩f共c1,c2兲典 共70兲 and the “hydrodynamic part” of ␸ HCS共c1,c2兲is by definition ␸ HCS 共h兲共c1,c2兲⬅P ␸ HCS共c1,c2兲.共71兲 Next we apply Pon both sides of Eq. (61)and, in order to get a closed equation for ␸ HCS 共h兲, we make the approximation P⌳共ci兲=P⌳共ci兲P.共72兲 A theoretical estimation a priori of the accuracy of this approximation would require to know more about the spectrum of ⌳共c兲and its adjoint than it is available at present. Therefore, it will be taken here as a working hypothesis to be evaluated later on by comparing the predictions it leads to with the results from particle simulations of the system. In this way, we obtain from Eq. (61) 关⌳共c1兲+⌳共c2兲兴 ␸ HCS 共h兲共c1,c2兲=−PT ¯ 0共c1,c2兲 ␹ HCS共c1兲 ␹ HCS共c2兲 +关⌳共c1兲+⌳共c2兲兴P ␦ 共c1−c2兲 ␹ HCS共c2兲.共73兲 By writing ␸ HCS 共h兲共c1,c2兲=兺 ␤ , ␤ ⬘=1 3 a ␤␤ ⬘ ␰ ␤ 共c1兲 ␰ ␤ ⬘共c2兲,共74兲 the above equation can be easily solved for the coefficients a ␤␤ ⬘. Some details of the calculations are given in the Appendix. The result is: ␸ HCS 共h兲共c1,c2兲=a33 ␰ 3共c1兲 ␰ 3共c2兲,共75兲 where a33共 ␣ 兲=d+1 2d+d+2 4da2共 ␣ 兲+b共 ␣ 兲.共76兲 In the first Sonine approximation, a2共 ␣ 兲is given by Eq. (58) and b共 ␣ 兲=2+d−6d2−共10 − 15d+2d2兲 ␣ −2共2+7d兲 ␣ 2+2共10 − d兲 ␣ 3 6d共2d+1兲−2d共11 − 2d兲 ␣ +12d ␣ 2−12d ␣ 3.共77兲 The structure of the result for ␸ HCS in Eq. (75)is a consequence of the strict conservation of the number of particles and total momentum. Moreover, in the elastic limit ␣ →1, a2 vanishes and b=−共1+d兲/2d, so that a33 also vanishes, reflecting that the energy is also strictly conserved in this limit, and there are no velocity correlations. On the other hand, for ␣ ⬍1, velocity correlations are induced by the inelasticity of collisions and, as a consequence, fluctuations of the total energy E=兺j NmVj 2/2 show up. These fluctuations have an intrinsic dissipative character. To compute their second moment, we start by realizing that the energy can be expressed as E关⌫共t兲兴 = 冕 dxmv2 2F1共x,t兲,共78兲 with F1共x,t兲being the microscopic phase space density defined in Eq. (6). The instantaneous deviations of E关⌫共t兲兴 from its average value in the HCS are given by ␦ E共⌫,t兲⬅E关⌫共t兲兴 −具E共t兲典HCS,共79兲 具E共t兲典HCS ⬅ 冕 d⌫E共⌫兲 ␳ HCS共⌫,t兲.共80兲 Then, using Eq. (78)it is easily seen that 具关 ␦ E共t兲兴2典=m2 4 冕 dx1 冕 dx2v1 2v2 2h1,1,HCS共x1,t;x2,t兲 =NkB 2THCS 2共t兲 冕 dc1 冕 dc2c1 2c2 2 ␸ HCS共c1,c2兲 =NkB 2THCS 2共t兲d2具 ␰ ¯ 3共c1兲 ␰ ¯ 3共c2兲兩 ␸ HCS 共h兲共c1,c2兲典 =NkB 2THCS 2共t兲e共 ␣ 兲,共81兲 where e共 ␣ 兲=d2a33共 ␣ 兲.共82兲 As expected, this function monotonically increases as the value of ␣ decreases, vanishing in the limit of elastic collisions. To verify the above theoretical predictions, we have performed Molecular Dynamics (MD)simulations of a freely evolving system of inelastic hard disks, i.e. d=2, in a square box of size Lwith periodic boundary conditions. The event ENERGY FLUCTUATIONS IN THE HOMOGENEOUS…PHYSICAL REVIEW E 70, 011302 (2004) 011302-7 driven algorithm [27]has been used. Three different, low values, of the density have been considered, namely nH =0.02 ␴ −2, 0.01 ␴ −2, and 5⫻10−3 ␴ −2. The number of particles used for each of the above densities is N=1000,2000, and 4000, respectively. These values guarantee that, in all the simulations to be reported in the following, the size of the system is smaller than the critical size, Lc, beyond which the HCS becomes unstable and velocity vortices and spatial inhomogeneities are developed [28]. An estimation of Lchas been made by using the hydrodynamic equations obtained from the Boltzmann equation [5,29]. Using this expression, one gets that for the most unfavorable situation to be reported ( ␣ =0.7),itisL/Lc⬇0.86. Therefore the system is inside the stable homogeneous region. Moreover, we have checked in all the simulations that the system actually stays in the HCS, by monitoring the local velocity and density fluctuations. For each value of the density and the restitution coefficient a number of trajectories have been generated. By averaging over these trajectories, the time evolution of the average energy and, therefore, of the granular temperature THCS共t兲, has been obtained. In all cases it was very well fitted by Haff’s law, Eq. (45). The average value of E2共t兲has been obtained in the same way. Afterwards, we have verified that, when the time-dependent energy deviation is scaled with THCS 2共t兲, it fluctuates around a stationary average value, as predicted by Eq. (81). An example, corresponding to ␣ =0.85 and nH=0.02 is shown in Fig. 1. In this case, the results have been averaged over 800 trajectories. In Fig. 2, the results for the steady values of ␴ E 2 ⬅具共 ␦ E兲2典HCS/NkB 2THCS 2are plotted as a function of the coefficient of restitution ␣ for 0.7艋 ␣ 艋1 for nH=5⫻10−3 ␴ −2. Also plotted is the theoretical prediction, i.e., the function e共 ␣ 兲given by Eq. (82). A quite good agreement is observed between the theoretical and numerical results. As expected, the amplitude of the fluctuations monotonically increases as ␣ decreases, vanishing in the limit of elastic collisions. The error bars in the figures have been determined from the deviations of the instantaneous values of the scaled energy dispersion from its stationary average value. Similar results have been found for the other two larger densities considered in this paper, although discrepancies between theory and simulations of the order of 10%, probably due to density effects, are observed for ␣ =0.7. We have also investigated the shape of the probability distribution for the energy fluctuations. For each trajectory of the system, the trace of 共E共t兲−具E共t兲典HCS兲/具E共t兲典HCS has been partitioned into nonoverlapping bins of value 0.001 and the frequency distribution has been built using the data from all the trajectories corresponding to the same values of Nand ␣ . The resulting normalized distribution is shown in Fig. 3 for the case nH=0.01 ␴ −2 and ␣ =0.9, where it is seen that it is very well fitted by a Gaussian distribution (solid line). Simi- FIG. 1. Scaled dispersion of the energy fluctuations ␴ E 2 ⬅N具共 ␦ E共s兲兲2典/具E共s兲典2, as a function of the dimensionless time s defined in the main text, for a system of inelastic hard disks with ␣ =0.85 and nH=0.02 ␴ −2. The angular brackets average has been taken over 800 trajectories. FIG. 2. Average steady values of the scaled second moment of the energy fluctuations ␴ E,st 2, as a function of the restitution coefficient. The symbols are from MD simulations of a system of hard disks of density nH=5⫻10−3 ␴ −2 while the solid line is the theoretical prediction given in Eq. (82). FIG. 3. Normalized distribution of the relative energy fluctuations for a system of inelastic hard disks with nH=0.01 ␴ −2. the symbols are from MD simulation and the solid line the fit to a Gaussian. BREY et al. PHYSICAL REVIEW E 70, 011302 (2004) 011302-8 lar results are found for the other values of ␣ and nHconsidered along this work. This suggests the existence of a physical potential governing the amplitude of the energy fluctuations in the HCS, similar to the equilibrium state of ordinary fluids, but we have not been able to identify it. V. ENERGY TIME CORRELATION FUNCTION IN THE HCS In this section, the time correlation function of the total energy of the HCS, will be studied. To begin with, a global two-time velocity correlation function is introduced by ␺ HCS共c1,s−s⬘;c1 ⬘兲⬅ 冕 dl1h ˜ HCS共l1−l1 ⬘,c1,s−s⬘;c1 ⬘兲, 共83兲 where s艌s⬘. Then, from Eq. (50)it follows that 冋 ⳵ ⳵ s−⌳共c1兲 册 ␺ HCS共c1,s−s⬘;c1 ⬘兲=0, 共84兲 whose formal solution is ␺ HCS共c1,s−s⬘;c1 ⬘兲=e共s−s⬘兲⌳共c1兲 ␸ HCS共c1,c1 ⬘兲.共85兲 As in the previous section, we are interested in the hydrodynamic part of ␺ HCS, ␺ HCS 共h兲共c1,c2兲⬅P ␺ HCS共c1,c2兲.共86兲 Using again the approximation given in Eq. (72), Eq. (85) can be rewritten as ␺ HCS 共h兲共c1,s−s⬘;c1 ⬘兲=e共s−s⬘兲⌳共c1兲 ␸ HCS 共h兲共c1,c1 ⬘兲共87兲 and substitution of Eq. (75)yields ␺ HCS 共h兲共c1,s−s⬘;c1 ⬘兲=a33e−关 ␨ 0共s−s⬘兲/2兴 ␰ 3共c1兲 ␰ 3共c1 ⬘兲.共88兲 This is our primary result for the spatially integrated hydrodynamic part of the two-time velocity correlation function in the HCS. It shows that the correlation decays exponentially with a rate determined by the cooling rate. In the following, we are going to employ the above expression to compute the HCS time correlation function CEE共t,t⬘兲for the total energy, defined by CEE共t,t⬘兲⬅具E共t兲E共t⬘兲典HCS −具E共t兲典HCS具E共t⬘兲典HCS,共89兲 with t艌t⬘艌0. Here, 具E共t兲E共t⬘兲典HCS = 冕 d⌫E关⌫共t兲兴E关⌫共t⬘兲兴 ␳ HCS共⌫,0兲.共90兲 Taking into account Eqs. (83)and (51), it is easy to show that the above definition can be transformed into CEE共t,t⬘兲=NkB 2THCS共t兲THCS共t⬘兲 ⫻ 冕 dc1 冕 dc1 ⬘c1 2c1 ⬘2 ␺ HCS共c1,s−s⬘;c1 ⬘兲 =NkB 2THCS共t兲THCS共t⬘兲d2 ⫻具 ␰ ¯ 3共c1兲 ␰ ¯ 3共c1 ⬘兲兩 ␺ HCS 共h兲共c1,s−s⬘;c1 ⬘兲典 =NkB 2THCS共t兲THCS共t⬘兲e共 ␣ 兲e−关 ␨ 0共s−s⬘兲/2兴.共91兲 This result is consistent with Eq. (81), since the former reduces to the latter for t=t⬘and, therefore, s=s⬘. In Fig. 4 results from MD simulations for a system with ␣ =0.85 and n=0.02 ␴ −2 are shown. They have been averaged over 800 trajectories. In the figure, CEE *共s,s⬘兲 ⬅CEE共t,t⬘兲/NkB 2e共 ␣ 兲THCS共t兲THCS共t⬘兲is plotted as a function of s−s⬘for several values of s⬘, namely s⬘=15.16,31.12, and 47.07. It is seen that the results do not depend on the value of s⬘as implied by Eq. (91). This is a consequence of the scaling property of the HCS, as discussed in Sec. III. Moreover, there is a quite good qualitative and quantitative agreement with the exponential decay predicted by the equation, given by the solid line in the figure and that has been evaluated by using Eq. (59)for the cooling rate ␨ 0. This is more clearly observed in Fig. 5, where a logarithmic scale has been used. For large times the simulation data become very noisy, but there seems to be some evidence of a systematic deviation from the exponential behavior, becoming stronger as the value of ␣ decreases. Whether this is the case and its possible relationship with hydrodynamic mode coupling effects FIG. 4. Normalized energy time autocorrelation function CEE *共s,s⬘兲as a function of s−s⬘for three different values of s⬘. The dimensionless time scale shas been defined in the main text. The symbols are from MD simulations of a system of inelastic hard disks with ␣ =0.85 and nH=0.02 ␴ −2. The solid line is the theoretical prediction as given by Eq. (91). ENERGY FLUCTUATIONS IN THE HOMOGENEOUS…PHYSICAL REVIEW E 70, 011302 (2004) 011302-9