Simulation study of the Green-Kubo relations for dilute granular gases

Simulation study of the Green-Kubo relations for dilute granular gases J. Javier Brey and M. J. Ruiz-Montero Física Teórica, Universidad de Sevilla, Apdo. de Correos 1065, E-41080 Sevilla, Spain (Received 24 June 2004; published 1 November 2004) The Green-Kubo relations for dilute granular gases are employed to compute their transport coefficients by means of the direct simulation Monte Carlo method. This requires not only to follow the dynamics of the system, but also to identify some modified fluxes appearing in the time-correlation functions. The results are compared with those obtained from the Boltzmann equation by means of the Chapman-Enskog procedure in the first Sonine approximation. A good agreement is found for the shear viscosity over a wide range of inelasticities. Nevertheless, for the two transport coefficients associated with the heat flux, significant discrepancies appear for strong inelasticity. Their origin is discussed, showing that they are partially due to the presence of velocity correlations in the homogeneous cooling state of a dilute granular fluid. DOI: 10.1103/PhysRevE.70.051301 PACS number(s): 45.70.⫺n, 51.10.⫹y, 05.20.Dd I. INTRODUCTION Hydrodynamics has been extensively used with clear success to describe the behavior of low-density, rapid granular flows [1–4]. From a theoretical point of view, the appropriate context to address fundamental issues is provided by the kinetic theory and nonequilibrium statistical mechanics methods. This includes the existence itself of a macroscopic description analogous to the one provided by the Navier-Stokes equations for molecular gases, the form of these equations, and the explicit expressions of the transport coefficients appearing in them. The prototypical idealized model for a granular gas is a system of smooth inelastic hard spheres or disks, the inelasticity being characterized by a constant coefficient of normal restitution. For this model, hydrodynamic equations have been derived starting from the Boltzmann equation and using the Chapman-Enskog procedure [5,6]. The accuracy of some of these results has been confirmed via the direct Monte Carlo simulation (DSMC)method, at least for moderate dissipation [7]. Nevertheless, and to put them in a proper context, it must be emphasized that the method used in the derivation is formal, in the sense that it does not determine the range of validity of the obtained hydrodynamic description [8]. A limitation of the explicit expressions for the transport coefficients as derived in the works mentioned above, is that it leads to rather complicated differential equations. Then, in practice, one has to resort to expansions in orthogonal polynomials, restricting the evaluation to the lowest orders, without any solid justification about the accuracy of such approximation. Recently, the transport coefficients of a granular gas following from the Boltzmann equation, have been expressed in the form of low density Green-Kubo relations [9,10]. They involve averages, with the one-particle distribution function of the homogeneous cooling state, of the product of two one-particle dynamical properties computed at different times. The time dependence is defined by means of a linear Boltzmann collision operator. As expected, they differ from those for molecular systems in many relevant ways, due to the energy dissipation in collisions. Although much more complicated than the expressions for molecular fluids, the Green-Kubo relations for a dilute granular gas can also be transformed in a form that is suitable for evaluation by means of N-particle simulation techniques. The general strategy to be followed has been discussed in detail in Ref. [12]. It is based on the property that the dynamics of a granular system in the time-dependent homogeneous cooling state (HCS)can be exactly transformed into a different dynamics around a stationary state [13]. Moreover, in order to transform the one-particle problem into an equivalent N-particle one, the same assumptions as needed to derive the Boltzmann equation are used. In particular, it is assumed that velocity correlations of colliding particles are negligible in the HCS. In this paper, we present simulation results obtained by employing the above method. This is relevant for several reasons. The accuracy of truncating the polinomial expansion, as carried out when deriving analytical expansions for the transport coefficients by the Chapman-Enskog method, is not known a priori. There is no reason to expect the same level of errors as in the case of elastic, molecular systems. A second, and more fundamental, possible source of discrepancy between the simulation results and the ChapmanEnskog predictions, can be the presence of relevant velocity correlations in the HCS, even in the very dilute limit. Although it is not easy to disentangle in practice both effects from the simulation data, it will be shown that some relevant information can be obtained. Moreover, the analysis presented here provides information about the decay of the correlation functions between the fluxes and the dynamical variables coupled to them. The latter are in fact closely related with the eigenfunctions of the linearized Boltzmann equation corresponding to the hydrodynamic modes of a dilute granular gas [8,10]. The fast enough decay of the correlation functions is a necessary condition for the existence of a hydrodynamic description. Finally, although the form of the Green-Kubo relations for dense granular fluids is not known, it can be expected, on the basis of what happens in molecular systems, that their structure will not differ too much from their dilute limit. Consequently, the present analysis may enlighten the study of denser systems. Green-Kubo relations for arbitrary densities have been derived in Ref. [14], by considering the linear response funcPHYSICAL REVIEW E 70, 051301 (2004) 1539-3755/2004/70(5)/051301(10)/$22.50 ©2004 The American Physical Society70 051301-1 tion to spacial perturbations of the HCS given by linear combinations of the local densities of mass, momentum, and energy. The low density limit of the expressions derived in this way differ from those being used in this work, which are consistent with the formal Chapman-Enskog result (before introducing any polinomial expansion). The origin of this discrepancy is discussed in Ref. [9]and is related with the form of the eigenfunctions of the linearized inelastic Boltzmann collision operator [10]. The plan of the paper is as follows. In the next section, the Green-Kubo relations for dilute granular gases are shortly reviewed, as well as the steady representation of the HCS and its implementation in N-particle simulations. In Sec. III, the DSMC method to be used in the simulations is described. Besides, the results for the velocity distribution function needed for the identification of the modified fluxes appearing in the Green-Kubo relations are reported. The evaluation of the transport coefficients is addressed in Sec. IV. All the involved time correlation functions are found to decay in an exponential way. Moreover, a fairly good agreement over a wide range of values of the inelasticity is found between the simulation results for the shear viscosity and the theoretical predictions obtained from the Boltzmann equation in the first Sonine approximation. Nevertheless, the presence of relevant velocity correlations manifests itself very clearly for small values of the restitution coefficient. For the transport coefficients associated to the heat flux, although the agreement is fairly good at low and moderate inelasticities, systematic deviations occur for strong dissipation. The accuracy of some analytical approximations for the modified fluxes present in the correlation functions is discussed as well. Also included is a comparison of the simulation results for the transport coefficient coupling heat flux and density gradient with some recent measurements directly based on the hydrodynamic description of a vibrated granular gas [11]. Finally, Sec. V contains a short summary of the results and some additional comments. II. GREEN-KUBO EXPRESSIONS FOR THE TRANSPORT COEFFICIENTS The expressions for the pressure tensor, Pij共r,t兲, and heat flux, q共r,t兲, to Navier-Stokes order for a dilute granular gas of dimension dare given by [5,6] Pij =p ␦ ij − ␩ 冉 ⳵ ui ⳵ rj+ ⳵ uj ⳵ ri−2 d ␦ ij ⵱·u 冊 ,共1兲 q=− ␬ ⵱T− ␮ ⵱n,共2兲 where pis the pressure, Tthe temperature, uthe velocity flow, and nthe number of particles density. Moreover, ␩ is the shear viscosity, ␬ the (thermal)heat conductivity, and ␮ another transport coefficient that vanishes in the elastic limit and will be referred to as the diffusive heat conductivity. Explicit expressions for the above transport coefficients have been derived from the Boltzmann equation for smooth inelastic hard spheres 共d=3兲and disks 共d=2兲of mass mand diameter ␴ , by using the Chapman-Enskog procedure, eigenfunction expansions, and also linear response theory, finding equivalent results. Moreover, it was shown [9,10]that the transport coefficients can be written in the form of GreenKubo relations. A particularly useful representation for N-particle simulations is obtained by exploiting an exact mapping of the HCS of a granular fluid onto a steady state [13]and assuming that, if there are the (one-time)velocity correlations present in the HCS, their effect can be neglected when computing the two-time correlation functions. The details of this steady representation for a dilute gas have been discussed in Ref. [12], where the particular case of the selfdiffusion coefficient was addressed. The analysis of the transport coefficients considered here proceeds in exactly the same way and, therefore, we directly quote the final expressions: ␩ 共T兲=nmv0共T兲 v ˜ 0,stN 冕 0 ⬁ dt具⌬xy共v,t兲⌽2,xy共v/v ˜ 0,st兲典ste− ␻ 0t,共3兲 ␬ 共T兲=nkBv0共T兲 N 冕 0 ⬁ dt具⌺ ˜ x共v,t兲⌽3,x共v/v ˜ 0,st兲典ste ␻ 0t,共4兲 ␮ 共T兲=mv0 3共T兲 N 冕 0 ⬁ dt具⌺ ˜ x共v,t兲⌽3,x共v/v ˜ 0,st兲典st共e ␻ 0t−1兲 +mv0 3共T兲 2v ˜ 0,stN 冕 0 ⬁ dt具⌺ ˜ x共v,t兲vx典st.共5兲 In the above expressions, ␻ 0is an arbitrary positive constant, Nis the number of particles in the system, kBis the Boltzmann constant, v0⬅共2kBT/m兲1/2 is the thermal velocity, and v ˜ 0,st⬅共2kBT ˜ st/m兲1/2, where T ˜ st will be identified below. Moreover, we have used the definitions ⌬xy共v兲=vxvy,共6兲 ⌺ ˜ x共v兲= 冉 v2 v ˜ 0,st 2−d+2 2 冊 vx,共7兲 ⌽2,xy共c兲=−cx ⳵ ln ␹ HCS共c兲 ⳵ cy,共8兲 ⌽3,x共c兲=−cx 2 冋 d+c· ⳵ ln ␹ HCS共c兲 ⳵ c 册 .共9兲 The function ␹ HCS共c兲is defined from the one-particle distribution of the HCS, fHCS共v,t兲, through fHCS共v,t兲=nv0 −d关T共t兲兴 ␹ HCS共c兲,c=v v0关T共t兲兴,共10兲 and it is an isotropic function of the vector c. The angular brackets in Eqs. (3)–(5)denote averages defined by 具a共v,t兲b共v兲典st = 冕 dr 冕 dvf ˜ st共v兲a共v,t兲b共v兲共11兲 with J. J. BREY AND M. J. RUIZ-MONTERO PHYSICAL REVIEW E 70, 051301 (2004) 051301-2 f ˜ st共v兲=nv ˜ 0,st −d ␹ HCS 冉 v v ˜ 0,st 冊 .共12兲 Finally, the time dependence of the dynamical variables is given by a共v,t兲=et⌳ ¯ sta共v兲.共13兲 Here ⌳ ¯ st is some linear operator involving both a Boltzmann collision term and also a streaming contribution proportional to ␻ 0[12]. Its explicit form will not be needed here. The relevant point for the analysis to be carried out in this paper is that, if velocity correlations in the HCS are assumed to be negligible in the low density limit, in this limit expression (11)is equivalent to the time-correlation function CAB,st共t兲=具A共t兲B典N,st −具A典N,st具B典N,st,共14兲 where 具A典N,st = 冕 d⌫ ␳ st共⌫兲A共⌫兲,具B典N,st = 冕 d⌫ ␳ st共⌫兲B共⌫兲, 共15兲 具A共t兲B典N,st = 冕 d⌫ ␳ st共⌫兲A共⌫,t兲B共⌫兲,共16兲 with ⌫denoting a point in the phase space of the system, ⌫⬅兵Ri,Vi;i=1,...,N其, and A共⌫兲=兺 i=1 N a共Vi兲,B共⌫兲=兺 i=1 N b共Vi兲.共17兲 Moreover, A共⌫,t兲⬅A关⌫共t兲兴 is generated from A共⌫兲by a modified particle dynamics consisting of an accelerating streaming between collisions, ⳵ ⳵ tRi共t兲=Vi共t兲,共18兲 ⳵ ⳵ tVi共t兲= ␻ 0Vi共t兲,共19兲 while the effect of a collision between particles iand jis to instantaneously alter their velocities according to Vi→Vi ⬘=Vi−1+ ␣ 2共 ␴ ˆ·Vij兲 ␴ ˆ, Vj→Vj ⬘=Vj+1+ ␣ 2共 ␴ ˆ·Vij兲 ␴ ˆ,共20兲 where Vij⬅Vi−Vjand ␴ ˆis the unit vector pointing from the center of particle jto that of particle iat contact. The parameter ␣ is the coefficient of normal restitution characterizing the inelasticity of collisions. It is defined in the interval 0⬍ ␣ 艋1 and is considered here as a constant, independent of the relative velocities. Under this dynamics, the system is expected to reach a steady state after a short transient period [12,13,15]. In the steady state, the energy dissipated in collisions is balanced by the effect of the acceleration between them. The function ␳ st共⌫兲in Eqs. (15)and (16)is the N-particle distribution function corresponding to this steady state. Moreover, T ˜ st, introduced implicitly above through v ˜ 0,st, is the temperature parameter of the steady state, i.e., dNkBT ˜ st/2=具E共⌫兲典N,st, with Ebeing the total kinetic energy of the system. The value of T ˜ st is related with the cooling rate of the HCS, ␨ HCS共t兲,by[12,15] T ˜ st = 冉 2 ␻ 0 ␨ ¯ 冊 2, ␨ ¯ = ␨ HCS共t兲 THCS 1/2 共t兲.共21兲 Since all the time dependence of ␨ HCS共t兲occurs through the temperature and it is proportional to THCS 1/2 , it follows that ␨ ¯ does not depend on time. The dynamics defined by Eqs. (18)–(20)can be easily implemented in particle simulations. Then, molecular dynamics (MD)simulations could be used to evaluate the transport coefficients as given by Eqs. (3)–(5). Of course, in this case one should keep in mind that Eqs. (11)and (14)are expected to be equivalent only in the low density limit. Another possibility is to employ the DSMC method [16], which is specially designed to simulate the N-particle dynamics of a system in the low density limit. An important advantage of the DSMC method in the present context, as compared with MD simulations, is that it allows to particularize the dynamics of the system for the case of homogeneous situations, therefore eliminating the spontaneous development of the spacial inhomogeneities following from the long wavelength hydrodynamic instability exhibited by the HCS [17]. As already mentioned in the Introduction, the structure of Eqs. (3)–(5)differs from the standard forms for the GreenKubo expressions for molecular (elastic)systems in several ways. First, the averages are taken over the velocity distribution corresponding to the steady state reached by the system under the modified dynamics. This distribution is different from the Maxwellian for all ␣ ⬍1. Second, the time correlation functions appearing in the expressions are not constructed from the momentum and energy fluxes, ⌬xy and ⌺ ˜ x, alone. Each of them is paired with another function, a “modified flux,” which is related with the derivative of the velocity distribution of the HCS. Third, the time evolution of the dynamical variables is not defined in terms of the particle Newton equations of motion, but includes a friction term. Finally, the time integrals contain, in addition to the time correlation functions, exponential in time factors, due to the collisional cooling of the HCS. III. THE SIMULATION METHOD In order to evaluate the transport coefficients from Eqs. (3)–(5), we have used the DSMC method to simulate the N-particle dynamics of a dilute granular gas [16,18]. Since we are interested in computing averages and time correlations of position-independent properties in a homogeneous state, the positions of the particles play no role in the simulations, and it is enough to consider just one cell in configuration space. In other words, every pair of particles in the system can collide with a probability depending only on their SIMULATION STUDY OF THE GREEN-KUBO …PHYSICAL REVIEW E 70, 051301 (2004) 051301-3 relative velocity. Consequently, neither the size of the system nor boundary conditions must be specified. In the simulations to be reported here, we have considered a system of N=104hard disks 共d=2兲. The results will be expressed in the following units. The unit of mass is the mass mof a particle and the unit of length is ᐉ=共n ␴ d−1兲−1, which is proportional to the mean free path. The unit of time is ᐉ关2kBT ˜ 共0兲/m兴−1/2, where T ˜ 共0兲is the initial scaled temperature. Moreover, we set kB=1, implying that in our units it is T ˜ 共0兲=1/2. Starting from a Maxwellian velocity distribution, the system is allowed to evolve with the dynamics defined by Eqs. (18)–(20)until it reaches a steady state. Then, all the statistical averages of interest are accumulated. Moreover, the results to be presented have been averaged over a number of different trajectories of the system, typically 6000, in order to increase the statistical accuracy.Along the simulations, the behavior of the total momentum of the system must be controlled, since it is unstable due to the presence of the friction term in the scaled dynamics, and round-off numerical errors propagate exponentially in time. This difficulty is eliminated by computing the total momentum at regular time intervals and subtracting it evenly from the momentum of each particle. A practical important point is the choice of the parameter ␻ 0. In principle, its value is arbitrary and determines the value of the temperature of the steady state, T ˜ st, as established by Eq. (21), and also the rate at which this steady state is approached [12,15]. On the other hand, inspection of Eqs. (3)–(5)shows that the numerical evaluation of the correlation functions appearing in the expressions of the transport coefficients is simplified if v ˜ 0,st=1. In the units we are using, this is equivalent to T ˜ st=1/2 or ␻ 0= ␨ ¯ /2冑2. The problem is that the expression of ␨ ¯ is only partially known. In the socalled first Sonine approximation, it is given by [19,20] ␨ ¯ ⬇ ␨ ¯ 共1兲=2 ␲ 共d−1兲/2共1− ␣ 2兲 ⌫ 冉 d 2 冊 ᐉd 冉 kB m 冊 1/2 冋 1+ 3 16a2共 ␣ 兲 册 , 共22兲 with a2共 ␣ 兲=16共1− ␣ 兲共1−2 ␣ 2兲 9+24d+共8d−41兲 ␣ +30 ␣ 2−30 ␣ 3.共23兲 In fact, the analytical expression of the distribution function of the HCS, ␹ HCS, that is needed to construct the “modified fluxes” appearing in the expressions of the transport coefficients, is only known in the same approximation, in which it reads ␹ HCS共c兲⬇ ␹ HCS 共1兲共c兲=e−c2 ␲ d/2关1+a2共 ␣ 兲S共2兲共c2兲兴,共24兲 where S共2兲共c2兲=c4 2−d+2 2c2+d共d+2兲 8.共25兲 Then, what has been done is the following. For each value of the coefficient of restitution ␣ , a preliminary series of simulations has been carried out, with the parameter ␻ 0set to ␻ 0= ␨ ¯ 共1兲/2冑2. These simulations were used to determine ␹ HCS共c兲and also the actual value of ␨ ¯ , from the measured value of T ˜ st through Eq. (21). Afterwards, in the second series of simulations, the value of ␻ 0is fixed by the same expression as before, but now using for ␨ ¯ the result obtained in the previous simulations. This guarantees that T ˜ st=1/2 and v ˜ 0,st=1 within the numerical errors. Once the steady state is reached, the time correlation functions are measured. Now we describe the results from the first series of simulations. The expressions of the transport coefficients, Eqs. (3)–(5), contain velocity derivatives of ␹ HCS共c兲that, due to the isotropic property of this function, can be easily related to ⳵ ln ␹ HCS共c兲/ ⳵ c. To measure this quantity in the simulations, the range of chas been partitioned into nonoverlapping bins of value ⌬c=8⫻10−2, and the frequency distribution has been built from the simulation data, measured once the system is in the steady state. This provides ␹ HCS, and afterwards its logarithm is computed also numerically. In Figs. 1 and 2, ⳵ ln ␹ HCS/ ⳵ cis plotted as a function of cfor ␣ =0.9 and ␣ =0.6, respectively. The circles are the results from computing the numerical derivative directly from the raw simulation data, while the solid line has been obtained by carrying out an interpolation of the numerical data for ln ␹ HCS共c兲to a smaller bin value 共⌬c=5⫻10−3兲before computing its derivative. For comparison, we have also included in the figures the derivative for the Gaussian distribution (dashed line)as well as for the first Sonine approximation, i.e., from Eq. (24)(dot-dashed line). As expected, the Sonine approximation describes quite well the behavior of the distribution function for cⱗ2.5 (thermal region), while the discrepancy grows very fast for larger values of the velocity. In Fig. 2, it is observed that FIG. 1. Plot of 共ln X兲⬘⬅ ⳵ ln ␹ HCS/ ⳵ cas a function of cfor ␣ =0.9. The circles are the numerical derivative of the simulation results, the solid line the fitted function, while the dashed and dotdashed lines are the Gaussian and first Sonine approximations for this quantity, respectively. Quantities are measured in the dimensionless units defined in the main text. J. J. BREY AND M. J. RUIZ-MONTERO PHYSICAL REVIEW E 70, 051301 (2004) 051301-4 ⳵ ln ␹ HCS/ ⳵ ctends to a constant value for large c, consistently with the known exponential decay of ␹ HCS for large velocities [21,22]. This behavior is not observed in Fig. 1 because the velocity range for which the exponential decay shows up increases very fast as ␣ approaches unity. As already indicated, from these simulations we also determined the actual value of ␨ ¯ from the measured value of T ˜ st and the value to be used for ␻ 0in the second series of simulations. Let us mention that the relative discrepancy between the value obtained in this way and the prediction of the fist Sonine approximation, Eq. (22), was always smaller than 0.2%. In the next section, these values as well as the interpolated results for ⳵ ln ␹ HCS/ ⳵ c, exemplified by the solid lines in the above two figures, will be used to evaluate the transport coefficients. IV. TRANSPORT COEFFICIENTS Since in the simulations where the correlation functions are measured, ␻ 0is set to ␻ 0= ␨ ¯ /2冑2, the steady temperature is the same as the initial one. Of course, the velocity distribution changes from the initial Gaussian to its steady form. In Fig. 3, the time evolution of the temperature in the steady state is plotted for several values of the coefficient of restitution, namely ␣ =0.9, 0.6 and 0.3. Time is measured by the accumulated number of collisions per particle, and the origin has been taken once the system is in the steady state. It is seen that, as predicted, the temperature fluctuates around its initial value (note the very small vertical scale used in the figure). The physical origin of these fluctuations has been discussed in Ref. [23]. There, it was shown that they are intrinsically associated to the inelasticity of collisions and that their amplitude increases as ␣ decreases. A. The shear viscosity Let us define a reduced dimensionless shear viscosity ␩ * by ␩ *= ␩ 共T兲 ␩ 0共T兲,共26兲 where ␩ 0共T兲=共d+2兲⌫共d/2兲共mkBT兲1/2 ␴ −共d−1兲 8 ␲ 共d−1兲/2 共27兲 is the elastic shear viscosity in the first Sonine approximation. Use of Eq. (3)yields ␩ *= 8冑2 ␲ 共d−1兲/2 共d+2兲⌫共d/2兲ᐉv ˜ 0,st 冕 0 ⬁ dtJ ␩ 共t兲e− ␻ 0t,共28兲 with J ␩ 共t兲=1 N具⌬xy共v,t兲⌽2,xy共v/v ˜ 0,st兲典st.共29兲 In the simulations, we have measured the function J ␩ 共t兲, using the stationarity of the HCS in the scaled dynamics, i.e., that 具a共v,t+t0兲b共v,t0兲典st =具a共v,t兲b共v兲典st,共30兲 for arbitrary t0艌0. This allows us to average over many different samplings along each trajectory of the system. Figures 4 and 5 show, in a logarithmic representation, two typical correlation functions, corresponding to ␣ =0.95 and 0.5, obtained in this way. The symbols are the results from the simulations, while the solid line is a fit to an exponential function. It is seen in the figures that the decay of J ␩ 共t兲is very well fitted by an exponential at least until it decays two orders of magnitude from its initial value. In fact, for the times where relevant deviations from the exponential behavior are observed, the statistical noise is too large as to make a precise statement about whether the deviations are something more than just noise. In any case, the contribution to FIG. 2. The same as Fig. 1 but for ␣ =0.6. FIG. 3. Time evolution of the temperature, measured in the units defined in the text, in the steady state. Time is measured in terms of the number of accumulated collisions per particle, ␶ . The solid line corresponds to ␣ =0.9, the dotted line to ␣ =0.6, and the dashed line to ␣ =0.3. SIMULATION STUDY OF THE GREEN-KUBO …PHYSICAL REVIEW E 70, 051301 (2004) 051301-5 the time integral in Eq. (28)corresponding to times where the numerical data for J ␩ 共t兲differ significantly from the exponential is negligible. For these reasons, in order to calculate the coefficient of shear viscosity, we have fitted J ␩ 共t兲to an exponential, J ␩ 共t兲=J ␩ 共0兲e−␭ ␩ t.共31兲 Then, Eq. (28), after particularizing for d=2 and the units defined in Sec. III, leads to ␩ *=2冑2 ␲ J ␩ 共0兲 ␭ ␩ + ␻ 0.共32兲 In Fig. 6, the values of ␩ * obtained in this way are plotted as a function of the coefficient of normal restitution ␣ . They are represented by the black circles. The solid line is the theoretical prediction derived from the Boltzmann equation by using the Chapman-Enskog procedure in the first Sonine approximation [7]. Moreover, for comparison purposes, the simulation results that follows from making, in the expression of J ␩ given in Eq. (29), each of the two approximations discussed in the context of the velocity distribution, are also displayed. More precisely, the expression for the modified flux ⌽2,xy in the first Sonine approximation has been employed, i.e., ⌽2,xy共c兲⬇−cx ⳵ ln ␹ HCS 共1兲共c兲 ⳵ cy⬇2⌬xy共c兲 冋 1−a2 冉 c2−d+2 2 冊 册 , 共33兲 with a2共 ␣ 兲given by Eq. (23). In the last transformation, we have neglected nonlinear in a2contributions, consistently with the approximation leading to Eq. (23)[19,20]. The simulation results for ␩ * in this approximation are indicated by squares in the figure, while triangles are used for those corresponding to the Gaussian approximation [equivalent to formally set a2=0 in Eq. (33)]. The latter agrees with the result that is obtained by linear response methods and constructing the response function for a spatial perturbation of the HCS coupling only to the local densities of mass, momentum and energy [14]. The first conclusion following from the analysis of Fig. 6 is that the analytical expression derived in Ref. [7]fits quite well the simulation data with no approximations for the modified flux over the whole range of values of ␣ considered. In fact, for 0.65ⱗ ␣ 艋1, the results obtained in the different approximations are close, and their ␣ -dependence shows the same trend. Let us remark that for ␣ =0.7 the simulation results based on the different approximations for FIG. 4. Time evolution of the correlation function J ␩ 共t兲for ␣ =0.95. The circles are the results of the simulation, while the solid line is the best fit to an exponential. Quantities are measured in the units defined in the main text. FIG. 5. The same as in Fig. 4 but for ␣ =0.5. FIG. 6. Dimensionless reduced shear viscosity coefficient, ␩ *, as a function of ␣ . The solid line is the theoretical prediction derived in Ref. [7], while the symbols are from the simulations: the circles are obtained using the true (from the simulation)velocity distribution for the modified flux, the triangles with the Gaussian approximation, and the squares correspond to the first Sonine approximation. J. J. BREY AND M. J. RUIZ-MONTERO PHYSICAL REVIEW E 70, 051301 (2004) 051301-6 the modified fluxes are almost indistinguishable on the scale used in the figure. This is not surprising, since the fourth moment of ␹ HCS is known to coincide with that of the Maxwellian approximation for a value of ␣ very close to 0.7 [24]. On the other hand, for smaller values of the restitution coefficient, the N-particle Green-Kubo expression for the shear viscosity is overestimated if the Gaussian approximation is used for the modified flux, and underestimated when the first Sonine approximation is used for it. Quite interestingly, in the latter case ␩ * exhibits a maximum for ␣ ⯝0.4, a behavior that is qualitatively different from the numerical results with the right expression for the modified flux. There is a point deserving some additional comments. It can be wondered why the analytical results show a better agreement with the DSMC data for the exact Green-Kubo expression than the simulation results obtained by using the first Sonine approximation for the modified flux, since in the analytical derivation the expansion in Sonine polynomials is also used, and only the first order is kept. A possible reason is that, in the Chapman-Enskog procedure, the Sonine expansion is carried out not only at the level of the velocity distribution of the HCS, but also when computing the dynamics of the fluxes. The results in Fig. 6 suggest that both approximations together lead to some kind of self-consistency improving the accuracy of the results. Nevertheless, let us point out that things are in fact much more complicated since, for low values of ␣ , the simulation results clearly indicate the presence of relevant velocity correlations in the HCS, giving a nonvanishing contribution to the Green-Kubo relations. A more detailed discussion of this is delayed to the final section of the paper. B. The (thermal) heat conductivity The dimensionless reduced heat conductivity is defined as ␬ *= ␬ 共T兲 ␬ 0共T兲,共34兲 where ␬ 0共T兲=d共d+2兲2⌫共d/2兲kB共kBT兲1/2 ␴ −共d−1兲 16共d−1兲 ␲ 共d−1兲/2m1/2 共35兲 is the elastic limit. Then, from Eq. (4)it is obtained that ␬ *=16共d−1兲冑2 ␲ 共d−1兲/2 d共d+2兲2⌫共d/2兲ᐉ 冕 0 ⬁ dtJ ␬ 共t兲e ␻ 0t,共36兲 with J ␬ =1 N具⌺ ˜ x共v,t兲⌽3,x共v/v ˜ 0,st兲典st.共37兲 As it was the case for J ␩ 共t兲, the simulation results also show an exponential decay of J ␬ 共t兲. Two typical examples are given in Figs. 7 and 8 corresponding to ␣ =0.95 and ␣ =0.5, respectively. Therefore, we have fitted again the simulation data to an exponential function, J ␬ 共t兲=J ␬ 共0兲e−␭ ␬ t.共38兲 Substitution of Eq. (38)into Eq. (36), after particularizing for d=2 and the units we are using, yields ␬ *= 冉 ␲ 2 冊 1/2 J ␬ 共0兲 ␭ ␬ − ␻ 0.共39兲 Here we have assumed that ␭ ␬ ⬎ ␻ 0, otherwise the time integral would diverge and the thermal conductivity would not exist. This condition has been fulfilled in all the cases we have considered. The results obtained for ␬ * are plotted as a function of the coefficient of restitution in Fig. 9, where also the theoretical prediction obtained by the Chapman-Enskog procedure in the first Sonine approximation [7]is shown. Although in both cases the value of the heat conductivity increases as ␣ decreases, there is a relevant quantitative discrepancy, the theoretical prediction growing much faster than the simulation data for ␣ ⱗ0.65. This is probably due to the fact that FIG. 7. Decay of the time correlation function J ␬ , measured in the units defined in the main text, for ␣ =0.95. The circles are from the simulations, while the solid line is the best exponential fit. FIG. 8. The same as Fig. 7 but for ␣ =0.5. SIMULATION STUDY OF THE GREEN-KUBO …PHYSICAL REVIEW E 70, 051301 (2004) 051301-7 the fluxes appearing in the expression for the heat conductivity involve higher velocity moments, as compared with the expression for the shear viscosity. The first Sonine approximation seems not to be able to accurately describe the behavior of these moments. Of course, when interpreting the above discrepancies, it must be kept in mind that a part of them can be due to the presence of velocity correlations in the HCS, as already mentioned. We again refer to the next section for a further discussion of this. The results coming from the other approximations discussed in the context of the shear viscosity are not illuminating and will not be addressed here. C. The diffusive heat conductivity The dimensionless reduced coefficient ␮ * associated to ␮ 共T兲is defined by ␮ *= n T ␬ 0共T兲 ␮ 共T兲,共40兲 so that substitution of Eq. (5)gives ␮ *=32冑2共d−1兲 ␲ 共d−1兲/2 d共d+2兲2⌫共d/2兲ᐉ 冋 冕 0 ⬁ dtJ ␬ 共t兲共e ␻ 0t−1兲 +1 2v ˜ 0,stN 冕 dt具⌺ ˜ x共v,t兲vx典st 册 .共41兲 The simulation results show that the second term on the right hand side is negligible, as compared with the first term. Even more, it is seen to identically vanish within the numerical precision of the simulation data. Since we have already shown that J ␬ 共t兲can be accurately described by Eq. (38),itis obtained that ␮ *=共2 ␲ 兲1/2J ␬ 共0兲 ␻ 0 共␭ ␬ − ␻ 0兲␭ ␬ ,共42兲 where we have particularized for d=2 and the units used in the simulations. Moreover, it has been assumed again that ␭ ␬ ⬎ ␻ 0in all the simulations being described. The simulation results for ␮ * as a function of ␣ are displayed in Fig. 10, where they are represented by the black circles. The solid line is, as in the previous figures, the analytical expression obtained by the Chapman-Enskog procedure in the first Sonine approximation [7]. A systematic discrepancy is observed for ␣ ⱗ0.7. The theoretical prediction grows much faster than the simulation results from the Green-Kubo expression, a similar behavior to that exhibited by the (thermal)heat conductivity in Fig. 9. The value of the transport coefficient ␮ * has been measured recently by direct application of Eq. (2)at the position of the temperature minimum presented by a steady vibrated granular gas in presence of gravity [11]. From DSMC measurements of the heat flux and the temperature at the above minimum, the values of ␮ * were obtained. They are represented by triangles in Fig. 10. Although it is not completely clear, one is tempted to say that these “hydrodynamic” measurements of ␮ * show a similar behavior to the Green-Kubo expression evaluated in this paper. In particular, both grow slower than the first Enskog approximation when ␣ decreases below ␣ ⱗ0.75. Unfortunately, this hydrodynamic method to measure ␮ * cannot be extended to arbitrarily small values of ␣ . In the considered steady state, there is a coupling between inelasticity and gradients, so that for small values of ␣ the system develops strong gradients and the Navier-Stokes approximation is no longer valid. V. SUMMARY In this paper, the Green-Kubo expressions for the transport coefficients of a dilute granular gas composed of smooth FIG. 9. The dimensionless reduced coefficient of thermal conductivity, ␬ *, as a function of ␣ . The symbols are from the simulations, while the solid line is the theoretical prediction derived in Ref. [7]. FIG. 10. The dimensionless reduced coefficient of diffusive heat conductivity ␮ * as a function of ␣ . The circles are from the simulations using the Green-Kubo expression, and the solid line is the theoretical prediction derived in Ref. [7]. The triangles are from an independent study of this transport coefficient in vibrated systems [11]. J. J. BREY AND M. J. RUIZ-MONTERO PHYSICAL REVIEW E 70, 051301 (2004) 051301-8 inelastic hard disks have been evaluated by means of the DSMC method. This N-particle algorithm is designed to mimic the dynamics of a low density gas. The structure of the Green-Kubo formulas is strongly modified by the inelasticity of collisions. Of particular relevance for their evaluation is that the correlation functions involve, in addition to the usual microscopic fluxes, other dynamical variables that are expressed in terms of velocity derivatives of the distribution of the HCS. Therefore, their exact analytical expressions are not known, and their values have to be obtained from the simulations themselves. The simulation technique used here takes advantage of the existence of a mapping between the homogeneous cooling state of a dissipative hard-sphere model and the steady state reached by the system under a modified dynamics. It has been found that, in the steady state representation, the time correlation functions of the fluxes and their paired dynamical variables decay exponentially in time. Moreover, in the case of the transport coefficients associated with the heat flux, the relaxation time is small enough as to compensate the explicit exponentially growing factor coming from the energy dissipation in collisions. This nontrivial result provides further support to the validity of the hydrodynamic description for dilute granular gases even in the case of quite strong inelasticity. Proving that the formal Green-Kubo relations for inelastic gases are amenable to computer simulation is one of the aims of this paper. The results for the transport coefficients show that all of them increase as the value of the coefficient of restitution ␣ decreases, in agreement with the predictions following from the Chapman-Enskog solution of the Boltzmann equation in the first Sonine approximation. Nevertheless, significant quantitative discrepancies occur in the small ␣ region in the case of the transport coefficients defining the heat flux, i.e., the thermal and diffusive heat conductivities, ␬ and ␮ , respectively. The Sonine approximation predicts a much more rapid increase than the one observed in the simulations. To properly evaluate this discrepancy, two main features must be taken into consideration. First, it must be realized that the way in which the Sonine approximation is introduced in the Chapman-Enskog scheme affects both, the initial form of the dynamical variables and its time evolution. The second point to be considered has a deeper physical origin. The time-correlation functions appearing in the Green-Kubo like relations derived from the Boltzmann equation are single-particle correlation functions in the HCS, i.e., they describe the time correlations of dynamical properties of just one-particle in that state, as generated by the linear Boltzmann operator. In order to transform the above expressions into others involving the N-particle dynamics needed for particle simulations, Eqs. (3)–(5), the assumption was made that (one-time)velocity correlations are negligible in the HCS [12]. It is possible to partially investigate whether this property actually holds as follows. Consider the initial value of J ␩ 共t兲defined in Eq. (29). As discussed in Sec. II what is actually computed in the simulations is J ␩ ⬘共0兲=1 N兺 i N 兺 j N 具⌬xy共vi兲⌽2,xy共vj/v ˜ 0,st兲典N,st =J ␩ ⬘共1兲共0兲+J ␩ ⬘共2兲共0兲,共43兲 where J ␩ ⬘共1兲共0兲=1 N兺 i N 具⌬xy共vi兲⌽2,xy共vi/v ˜ 0,st兲典N,st,共44兲 J ␩ ⬘共2兲共0兲=1 N兺 i N 兺 j⫽i N 具⌬xy共vi兲⌽2,xy共vj/v ˜ 0,st兲典N,st.共45兲 The diagonal component, J ␩ ⬘共1兲共0兲, is identically the same as J ␩ 共0兲. In fact, it can be easily evaluated and, in the units used in this paper and with our choice for ␻ 0,itisJ ␩ ⬘共1兲共0兲=1/2. The other component, J ␩ ⬘共2兲共0兲, only differs from zero if velocity correlations are present in the system, and it has been assumed to be negligible. We have measured the values of FIG. 11. Simulation results for the correlation function J ␩ ⬘共0兲 defined in Eq. (43)(filled circles), and its diagonal part, J ␩ ⬘共1兲(empty circles), as a function of the restitution coefficient ␣ . Quantities are measured in the units defined in the text. FIG. 12. Simulation results (filled circles)for the N-particle correlation function J ␬ ⬘共0兲. The empty circles are its diagonal part, J⬘共1兲. Quantities are measured in the units defined in the main text. SIMULATION STUDY OF THE GREEN-KUBO …PHYSICAL REVIEW E 70, 051301 (2004) 051301-9