Randomly Driven Granular Fluids: large scale structure

The nonequilibrium steady state of a granular fluid, driven by a random external force, is demonstrated to exhibit long range correlations, which behave as $\sim 1/r$ in three and $\sim \ln(L/r)$ in two dimensions. We calculate the corresponding structure factors over the whole range of wave numbers, and find good agreement with two-dimensional molecular dynamics simulations. It is also shown by means of a mode coupling calculation, how the mean field values for the steady state temperature and collision frequency, as obtained from the Enskog-Boltzmann equation, are renormalized by long wavelength hydrodynamic fluctuations.


I. INTRODUCTION
Systems of granular particles, like grains of sand or more ideally glass, plastic or metal beads, exhibit different flow regimes [1], depending on the external forcing. A systematic experimental study of the rapid or collisional flow regime as compared to the quasi-static, slow or frictional regime was first performed by Bagnold [2] using an annular shear cell. Later a similar but more refined characterization was made in Ref. [3]. The possibility of coexistence of different flow regimes was observed in an experimental study of flows down an inclined chute [4].
In several more recent experimental studies of the rapid granular flow regime, more microscopic properties have been measured. In Ref. [5] the fluidization behavior of a vertically vibrated two-dimensional model granular material has been investigated using high speed photography. Patterns at the surface of a vertically vibrated granular layer, analogous to Faraday waves in molecular fluids, have been observed in Ref. [6] and stimulated the interest of many theorists [7]. An understanding of these patterns through a derivation of e.g. an amplitude equation [8] from the hydrodynamic description of the system, is still lacking, however. In Ref. [9] the effect of inelastic collisions on the formation of clusters is investigated in a system of particles rolling on a smooth surface and driven by a moving wall. Finally, Ref. [10] studies the steady state of a vertically shaken granular monolayer, and discusses clustering, inelastic collapse and long range order.
Even rapid flows of model granular materials are poorly understood in general, since complicating effects, such as gravity and interactions with boundaries, have to be taken into account. If the model granular material consists of spherical grains with a smooth surface, collisions between particles can be characterized only by their coefficient of restitution α or their inelasticity ǫ = 1 − α 2 . We assume this coefficient of restitution to be a constant, independent of the relative velocity between the colliding particles, and refer to the model as the inelastic hard sphere (IHS) model. Dissipative collisions complicate the dynamics in a nontrivial way: they may cause the system to become unstable, and give rise, for instance, to clustering; they create several new intrinsic length scales that might interfere for small inelasticity with the system size L, and for large inelasticity with the mean free path l 0 .
By driving an IHS fluid by boundaries or external fields it can reach a steady or oscillatory state. Due to the existence of these new "cooling" lengths, this state is frequently inhomogeneous, where the spatial gradients become larger at higher inelasticity. Only for small inelasticity the mean free path is well separated from the scale on which the macroscopic fields vary, and a hydrodynamic description [11] through Navier-Stokes or Burnett equations is expected to hold.
In the present paper we investigate the properties of an IHS fluid that is heated uniformly so that it reaches a spatially homogeneous steady state. This way of forcing, where a random external force accelerates a particle, was proposed by Williams and MacKintosh [12] for inelastic particles moving on a line. Peng and Ohta [13] performed simulations on a 2D version of this model. In two dimensions the model may be considered to describe the dynamics of light disks moving on an air table, a system which has been investigated experimentally in Ref. [14]. In three dimensions it can be extended to include gravitational and drag forces, making it to some extent relevant for gas-fluidized beds [15] when hydrodynamic interactions are unimportant. A similar IHS model with random external accelerations has been used by Bizon and Swinney [16] in their computer simulations to test continuum theories for vertically vibrated layers of granular material.
In the present paper we will describe the randomly driven IHS fluid in two and three dimensions, and characterize its nonequilibrium steady state (NESS). The single particle velocity distribution function in the NESS has been calculated in Ref. [17] from the Enskog-Boltzmann equation and was shown to be well approximated by a Maxwellian, except for an overpopulated tail ∼ exp(−Ac 3/2 ), where c is the velocity scaled by the thermal velocity and A ∼ 1/ √ ǫ. Computer simulations of the one-dimensional system of Ref. [12] showed the existence of long range spatial correlations in the steady state, which were addressed theoretically in in Ref. [18]. Here we will give quantitative predictions for long range correlations [19] in the two-and three-dimensional NESS. Moreover, we extend the mode coupling theory of Brito and Ernst [20] to analyze how long wavelength fluctuations in the NESS renormalize the mean field predictions of kinetic theory, and use this theory to calculate the renormalized temperature and collision frequency in the NESS.
To obtain an adequate description of the structures in steady granular flows one does not only need the equations of fluid dynamics for the average macroscopic behavior, but also the spatial correlation functions G ab (r), and their Fourier transforms, the structure factors S ab (k). Let δa(r, t) = a(r, t) − a(r, t) with (a = n, T, u α ) be the fluctuations of the slowly varying fields a(r, t), i.e. the local density, n(r, t), local temperature, T (r, t), and local flow velocity u α (r, t) (α = x, y, . . .), around their average values a(r, t) . Then the objects of interest are the correlation functions in the NESS, which are given by the limit, Here . . . is an average over some initial distribution, δa(k, t) is the spatial Fourier transform of δa(r, t), and S ab (k) that of G ab (r). Moreover we consider the unequal-time correlation functions in the NESS, defined as where F nn (k, t) is the intermediate scattering function [21]. The dynamic structure factor is then The paper is organized as follows. In section II we show how the macroscopic equations for granular flow are modified to account for the external driving/heating by the random accelerations. Section III characterizes the noise of external and internal fluctuations, and the structure factors and spatial correlation functions are calculated in sections IV and V. The latter section also presents the mode coupling calculations for the temperature and the collision frequency in the NESS. Computational details of our molecular dynamics (MD) simulations are described in section VI, and section VII compares our predictions with simulations. Some general comments and conclusions are presented in section VIII.

II. MACROSCOPIC EQUATIONS
Consider a system of inelastic disks or spheres (IHS) (d = 2, 3), driven by a heat source, which is described as a random acceleration,ξ i , Here F i is the systematic force on particle i = (1, 2, . . . , N) due to inelastic collisions. If the time constant of the heat source is much smaller than the mean free time t 0 between collisions, thenξ i (t) can be considered as Gaussian white noise with zero mean and correlation, where α, β = {x, y, . . .} denote Cartesian components of vectors or tensors. The over-line indicates an average over the noise source. It is understood that the ensemble average in Eqs.
(1) and (2), denoted by the angular brackets, also includes this noise average. To guarantee conservation of total momentum the random force has to obey the constraint iξi (t) = 0.
In thermodynamically large systems this constraint gives a correction to (5) of O(1/N), which can be neglected. The uniformly heated fluid is described by the standard macroscopic equations of fluid dynamics, where the temperature equation is supplemented with an additional source term mξ 2 0 , and a sink term Γ to account respectively for the heating and the energy loss through inelastic collisions: where ρ = mn, u the flow velocity, and 1 2 dnT the kinetic energy density in the local rest frame of the IHS fluid. The pressure tensor Π αβ = pδ αβ + δΠ αβ contains the local pressure p and the dissipative momentum flux δΠ αβ , which is proportional to ∇ α u β and contains the kinematic and longitudinal viscosities ν and ν l , defined below Eq. (A1) of appendix A. The constitutive relation for the heat flux, J = −κ∇T , defines the heat conductivity κ. For small inelasticity the transport coefficients ν, ν l , and κ are assumed to be given by the Enskog theory for a dense gas of elastic hard spheres (EHS) [22].
To lowest order in the spatial inhomogeneities the sink term, representing the energy loss through inelastic collisions, is given by [23] Γ = 2γ 0 ωT.
It is proportional to the granular temperature, to the average Enskog collision frequency (ω ∼ √ T , explicitly given in Eq. (A2) of the appendix), and to the coefficient of inelasticity ǫ = 1 − α 2 ≡ 2dγ 0 , where α is the coefficient of normal restitution. To explain the term mξ 2 0 in (6) we calculate the energy gain of a single particle due to the random force in a small time δt. This is done by formally integrating (4) and averaging over the noise source, i.e., where (5) has been used. Note that we have defined the granular temperature as twice the average random kinetic energy per translational degree of freedom, so that the Boltzmann constant k B does not appear in its definition. The above equations provide a consistent description for the heated IHS fluid at small inelasticity. The energy is not conserved in inelastic collisions, and consequently the temperature is not a hydrodynamic mode, but a kinetic mode with a relaxation rate ∝ γ 0 ω. Nevertheless, at small inelasticities (γ 0 ≪ 1), it is consistent to include temperature among the slowly changing macroscopic variables, which describe the dynamics of the system on time scales t large compared to the mean free time t 0 = 1/ω, and on spatial scales λ large compared to the mean free path, l 0 = v 0 t 0 , where v 0 = 2T /m is the thermal velocity.
At large inelasticities, where ǫ ∼ O(1), the temperature is a fast kinetic mode, that decays on the time scale t 0 , and cannot be included among the slow macroscopic variables.
In that case the IHS fluid becomes athermal, and the slow macroscopic fields only involve the density and flow field, as is the case in lattice gas cellular automata without energy conservation [24,25].
Let us consider the decay of temperature in more detail. For a homogeneous state, the fluid dynamic equations (6) will have as a solution n(r, t) = n, u(r, t) = 0, and T (r, t) = T (t), the latter satisfying For long times the system approaches a steady state with a constant temperature, determined by mξ 2 0 = 2γ 0 ωT . As ω ∼ √ T , we obtain the mean field prediction (A2), as deduced from the Enskog theory, Further symbols are defined below (A1) in the appendix. To obtain the final approach to the NESS, we linearize Eq. (9) around T E in Eq. (10). This yields an exponential approach, i.e.
In fact 3γ 0 ω can be identified as the decay rate z H (0) of the long wavelength components of the temperature fluctuations, as derived in section IV below (25). The exact time dependent solution of Eq. (9) can be obtained implicitly (t as a function of temperature), and reads where and T 0 is the temperature at t = 0.

III. NOISE CHARACTERISTICS IN THE NESS
The goal of this paper is to analyze the effects of spatial fluctuations δa(r, t) with (a = n, T, u α ) around the NESS on hydrodynamic space and time scales. As we are dealing with fluctuations, we linearize the nonlinear equations (6) around the NESS, with the result (A1) of appendix A. Moreover, to extend the average equations to fluctuating equations, valid on mesoscopic spatial and temporal scales, we need to calculate the external noise termŝ ξ ex (r, t) andθ ex (r, t) that contribute to ∂ t u and ∂ t T in (A1). These terms originate from the random accelerationξ i (t), which enters in the microscopic equations of motion (4). By starting from the microscopic expressions for the momentum and energy density one finds that the noise sources are given by the long wavelength components of These fields are again Gaussian white noise with zero mean and correlationŝ as follows from (5). Next, we argue on the basis of the hydrodynamic equations (6) that there exists, close to the NESS, a range of hydrodynamic wave numbers k ≫ k * , the so called elastic regime, where the dynamics of the fluctuations is the same as in a fluid of elastic hard spheres or disks, and is driven by internal noise that will be studied next. The validity of the hydrodynamic equations (6) and (A2) is restricted to wave numbers k ≪ 2π/l 0 (to guarantee separation of kinetic and hydrodynamic scales), and to k ≪ 2π/σ, where σ is the disk or sphere diameter (to guarantee that the Euler equations involve strictly local hydrodynamics). So, for the existence of an elastic regime in the IHS hydrodynamics, the following constraints must be satisfied Moreover, following McNamara [26] we can distinguish a dissipative regime, kl 0 ≪ ǫ [typically kl 0 < ∼ O(ǫ 2 )], and a standard regime, kl 0 ≫ ǫ [typically kl 0 > ∼ O( √ ǫ)], separated by a crossover regime, around kl 0 ∼ O(ǫ). In the dissipative regime, dissipation dominates compression effects and sound propagation, which are O(kl 0 ), as well as heat conduction, which is O(k 2 l 2 0 ). In the standard regime, dissipation effects are of the same order as heat conduction. As a consequence, the hydrodynamic modes and their propagation velocities are those of a fluid of elastic particles, while the corresponding damping rates of heat and sound modes still depend on the inelasticity. Only in the elastic regime, kl 0 ≫ k * l 0 ∼ O( √ ǫ), also these damping coefficients attain their elastic values. The above argument applies for small enough ǫ = 1 − α 2 = 2dγ 0 , where the inequalities (16) are obeyed, and an elastic regime exists and is well separated from the dissipative regime.
In the elastic regime the equations for the macroscopic deviations from the NESS are the same as those for a fluid of elastic hard spheres, deviating from thermal equilibrium. To describe fluctuating mesoscopic hydrodynamics on these length scales, one can add internal noiseξ in (r, t) andθ in (r, t), describing the rapid microscopic degrees of freedom. The noise strength of the internal fluctuations can be obtained from the fluctuation-dissipation theorem [27,23] for the EHS fluid, and is given in Fourier representation by where ν, ν l and κ are the transport coefficients for the EHS fluid, andk α is a component of the unit vectork = k/k. The effective noise in the heated IHS fluid may therefore be described by the sum of external and internal noise,ξ(k, t) =ξ ex (k, t) +ξ in (k, t), with a similar expression forθ(k, t). The noise characteristics ofξ(k, t) andθ(k, t) interpolate between two limiting behaviors and the corresponding noise strengths are given by the sum of (15) and (17).
Having specified the characteristics of the noise sources in the macroscopic equations, we conclude this section by summarizing the Langevin-type equations that describe the dynamics of the slow fluctuations. To do so, it is convenient to introduce the Fourier modes δa(k, t) exp[ık·r] of the linearized hydrodynamic equations (A1). The mesoscopic equations, valid on hydrodynamic space and time scales, i.e. t ≫ t 0 and kl 0 ≪ 1, then take the form where the components of the vector a are labeled with a = {n, T, l, ⊥}. Here a = l refers to the longitudinal velocity component u l (k, t) =k · u(k, t), and a =⊥ refers to (d − 1) transverse components of u(k, t). The matrix M ab with a, b = {n, T, l, ⊥} is given explicitly in (A3), andf a (k, t) is Gaussian white noise with nonvanishing components for a = T, l, ⊥ and correlation function The noise strength C ab (k) = δ ab C ab (k) is obtained by taking the Fourier transform of (15) together with (17), and is only nonvanishing for the following diagonal elements where the NESS condition of Eq. (9), i.e. Γ = 2γ 0 ωT = mξ 2 0 , has been used.

IV. STRUCTURE FACTORS
The equal-time structure factors, introduced in (1), obey the equations of motion which follows by formally integrating (18) and using (19). The left hand side of (21) vanishes since the structure factors do not depend on time in the NESS. The resulting equation can be solved by spectral analysis, or numerically. The spectral analysis is summarized in appendix A, where w λa and v λa are respectively the a-th component of the right and left eigenvectors of the hydrodynamic matrix M, and z λ (k) is the corresponding eigenvalue.
Taking then the scalar product of (21) on both sides with left eigenvectors (A5) of the appendix yields Using the completeness relation (A6) and the fact that off-diagonal elements of C ab in (20) vanish, we obtain which is the final result for the static structure factors. The time correlation function (2) in the NESS reduces in a similar manner to The above result corresponds to the Landau-Placzek theory [21] for hydrodynamic correlations in the NESS.
To obtain more explicit results we need the explicit forms of eigenvalues and eigenvectors, which have only been calculated for small k. The eigenvalue equation can be solved numerically for any given wave number, and the results are illustrated in Fig. 1 for the two-dimensional case. Transport coefficients, equation of state p(n, T ), and pair correlation function at contact are obtained from Enskog's theory for elastic hard spheres or disks. The generic features of the spectrum in Fig. 1 are the same as McNamara's case "α = β = 0", illustrated in Fig. 6(a) of his study [26] on hydrodynamics of granular materials, which corresponds to a temperature and density independent heat source. However, neither the equation of state, nor the transport coefficients, used in Ref. [26], correspond to the heated fluid of inelastic hard spheres, used in the present simulations. In the hydrodynamic regime (kl 0 < ∼ 1) all eigenvalues are found to be negative for nonvanishing wave numbers (see appendix A). So, all modes are linearly stable.
With the help of Mathematica the structure factors in the steady state have been calculated numerically from Eq. (21) with ∂ t S ab (k) = 0 for a given wave number. The resulting structure factors are shown by solid lines in Fig. 2, and will be tested against MD simulations in section VII.
Next, we present analytic results for the dissipative regime (kl 0 ≪ γ 0 ). The eigenvalues on the largest spatial scales can be determined as an expansion in powers of k at a fixed value of γ 0 , with the results In later applications the explicit form of the eigenvectors of M(k) is needed to lowest order in k. They read The coefficients g(n), v D , and D λ with λ = {⊥, H, ±} are calculated in the appendix. In the dissipative regime there are two propagating sound modes (λ = ±) with a propagation speed v D and a damping constant D S . There is a kinetic heat mode (λ = H), with a long wavelength relaxation rate z H (0) = −3γ 0 ω, in agreement with Eq. (11). Therefore, on the largest spatial scales the temperature deviations have decayed to zero, and temperature gradients do not exist; there is no heat conduction. In addition, there are (d − 1) transverse velocity or shear modes (λ =⊥), which are purely diffusive. The corresponding diffusivity, D ⊥ = ν, has the same form as for EHS. In (A7)-(A10) the coefficients are expressed explicitly in terms of thermodynamic quantities and transport coefficients.
In the standard regime, kl 0 ≫ γ 0 (and γ 0 small) the eigenvalues for shear and sound modes are to leading nonvanishing order the same as for the EHS fluid, where the sound waves propagate with the adiabatic sound speed v S of the elastic fluid, which is larger than the propagation speed v D in (A7). The damping of the sound and heat modes, on the other hand, are larger than in the elastic fluid due to the inelastic collisions. In the elastic regime, defined in section III, where kl 0 ≫ √ γ 0 , all transport coefficients are equal to their EHS values.
In summary, the eigenvalue spectrum z λ (k) for the uniformly heated IHS fluid is quite different from that of the freely evolving IHS fluid, linearized around the homogeneous cooling state [26]. In the free case, all shear modes (λ =⊥) and the heat mode (λ = H) are unstable in the dissipative regime (small k), and propagating modes do not exist for kl 0 ≪ γ 0 . Moreover, there exist in this regime a stable diffusive density mode and a kinetic temperature mode, which combine into two propagating modes for kl 0 ∼ O(γ 0 ), where crossover occurs from the dissipative to the standard regime. In the heated case, however, all modes are linearly stable, and the sound modes remain propagating in the dissipative regime down to k = 0.

V. EFFECTS OF LONG RANGE CORRELATIONS
In this section we study static and dynamic structure factors and corresponding correlation functions at the largest spatial scales. Moreover, we show by means of a mode coupling calculation, how average properties, which were calculated in section II on the basis of a mean field theory (i.e. the Enskog-Boltzmann equation) are renormalized by spatial fluctuations.
The static structure factors have been calculated in (23). In the dissipative regime (kl 0 < ∼ γ 0 ), the relevant eigenvalues (25) and eigenmodes (26) are discussed below (26). The dominant singularity at small wave number of the structure factors S ab (k) ∼ O(1/k 2 ) in (23) originates from pairs of transverse modes, where 2z ⊥ (k) = −2νk 2 , and from anti-parallel sound modes, where z + (k)+z − (k) = −2D S k 2 . We start with the transverse structure factor, where only the shear modes in (26) contribute, and deduce from the equations above, where the relation (20) has been used for kl 0 ≪ γ 0 . The structure factors S ab (k) for a, b =⊥ derive their dominant small-k behavior from two anti-parallel sound modes and we obtain with the help of (23), (25), (26), and (20) at small k, where the sound damping constant in the dissipative regime D S has been calculated in (A9). It depends on the inelasticity. The contribution of the internal noise is subdominant in this regime. In a similar manner we find where all nonvanishing coefficients labeled (ab) = (ll, nn, T T, nT ) are listed in Table I. The remaining structure factors are of O(1) as k → 0. Next, we consider the spatial correlation functions G ab (r), which are the inverse Fourier transforms of S ab (k). The small-k behavior of S ab (k), obtained above, enables us to calculate the large r behavior of the spatial correlation functions. The calculations are given in appendix B. One finds for the leading large r behavior in three-dimensional systems, and in two-dimensional systems, valid for r ≪ L, where L is the linear dimension of the system. The subleading large r corrections to (31) are constant terms, independent of r. In the calculations given in appendix B, these constants depend on a cutoff wave vector k min = 2π/L, used to evaluate the divergent k integrals occurring in the Fourier inversion of S ab (k). To calculate their precise values the subleading small-k corrections to (27) and (28) are required. The long wavelength behavior of the time dependent correlation function F ab (k, t) in (24) can be evaluated in a similar manner. We quote the result in terms of the Laplace transform,F ab (k, z), from which the dynamic structure factor (3) follows. In the dissipative regime (kl 0 ≪ γ 0 ) we find to leading order for small wave numbers, In a similar manner we obtain where all nonvanishing coefficients B ab are listed in Table I. The dynamic structure factor (3) then becomes It contains only Brillouin peaks, coming from the sound modes. There is no central Rayleigh peak, because the heat mode is not a slow, but a fast kinetic mode in this regime. In the elastic regime S nn (k, Ω) has the standard Rayleigh and Brillouin lines of the EHS fluid.
The existence of long range spatial correlations shows that the NESS is quite different from a thermal equilibrium state [19]. In fact, the spatial fluctuations also modify ("renormalize") the mean field predictions for the averages and the particle distribution functions. In appendix C a mode coupling calculation is presented to estimate the renormalization effects on the average energy per particle E = (1/N) i 1 2 mv 2 i and average collision frequency ω in the NESS, and we recall their mean field values, i.e. E E = 1 2 dT E and ω E given by (A2), i.e. ω E ∝ nχ(n) √ T E , where T E is given in (10). As it turns out, the fluctuation contributions, δE and δω, are finite and well-behaved in three dimensions, but logarithmically divergent in the system size L in two dimensions, so that d = 2 is the upper critical dimension. The mode coupling calculations of appendix C yield then in two dimensions, for large L, where C E and C ω are calculated in appendix C. The argument of the logarithm γ 0 L/l 0 is an estimate for the ratio of the values for the right (k ∼ γ 0 /l 0 ) and left (k min = 2π/L) boundaries of the dissipative k range, where the small-k behavior in (27) to (29) is valid. The logarithmic correction becomes only appreciable for large systems with a size L, much larger that the so called homogeneous cooling length l T = l 0 /γ 0 , which diverges in the elastic limit. Then the renormalization corrections for small inelasticity vanish as δT ∼ ǫ ln(ǫL/l 0 ) and δω ∼ ǫ 2 ln(ǫL/l 0 ). Here we have used the relations C E ∼ ǫ and C ω ∼ ǫ 2 for ǫ → 0, as can be deduced from the results in appendix C. A similar mode coupling theory has been recently used in Ref. [20] for freely evolving granular fluids to calculate the long time decay of the energy, which deviates from Haff's cooling law due to inhomogeneities in the hydrodynamic fields, and good agreement between theory and simulations was found. Before concluding this section we compare the theoretical predictions for the structure factors, with and without internal noise in Eqs. (20) and (21), as shown in Fig. 2. Inspection of (27) and (28) shows that only the external noise determines their dominant small-k behavior. The question then arises what are the effects of internal noise, and is it meaningful to include it in the theoretical description? The answer is affirmative, as we will show below. The steady state solution of (21) without internal noise clearly behaves at small wavenumbers as k −2 , but the k independent plateau values, shown by the full theory (solid lines in Fig. 2), are missing. These plateau values represent a very short distance correlation ∼ δ(r). Calculation of the plateau value for S αβ (k) yields (1/n 2 V ) i v iα v iβ = (T /ρ)δ αβ , i.e. the self-correlation term (i = j) in the definition (1) of S αβ (k), or more explicitly in Eq. (40) below. Then, addition of this plateau value to the numerical solution of (21) without internal noise yields the structure factors S without αβ (k), shown as dotted lines in Fig. 2. For the transverse velocity structure factor the dotted line coincides with the solid line in Fig. 2, as can be shown analytically from (20) and (21). Only the structure factor for the longitudinal flow field, S without (k) differs appreciably from S (k) in the relevant intermediate regime, 0.2 < ∼ kσ < ∼ 0.5. In section VII the simulation results will be compared with the theoretical predictions with and without internal noise. As it turns out, the comparison shows convincingly that the theory with/without internal noise agrees/disagrees with the simulations.

VI. COMPUTATIONAL DETAILS
In the two subsequent sections we describe molecular dynamics simulations performed to verify our theoretical predictions. In this section we present the computational details, before testing our theoretical results against simulations in the next section.
The model studied here has been extensively used in computer simulations for the freely evolving case (no forcing) [28][29][30]23], as well as for the randomly accelerated case in one [12] and two dimensions [13]. We consider a system of N inelastic hard disks having diameter σ in a two-dimensional square cell of length L, with periodic boundary conditions. The disks interact via inelastic collisions with coefficient of normal restitution α. For a colliding pair (i, j) of particles having equal masses the postcollision velocities are: where v ij = v i − v j , the asterisk denotes velocities after collision, andσ is a unit vector along the line connecting the centers of particle j and particle i. The energy loss in consecutive collisions, which is proportional to ǫ = 1 − α 2 , is compensated by a periodic (in time) and instantaneous perturbation of all velocities by a random amount. After every time step ∆t, the velocity of each particle is modified according to where the components of the vectors ϕ i are taken from a random distribution of zero mean and variance ϕ 0 (in practice a Gaussian or a flat function of finite support). The time step ∆t of this "heating" or "kicking" is chosen much smaller than the mean time between successive collisions of a tagged particle (typically a factor 10 4 smaller), so that the system under scrutiny reduces to that described in section II. The opposite limit, where ∆t is much bigger or comparable to the mean free path, was considered in Ref. [31] (in the presence of an additional external damping or drag force), and in that limit clustering was observed.
The relation between ϕ 0 and ξ 0 , the variance of the noise term in Eq. (4), can be deduced from the energy fed into the system by the kicks. This yields straightforwardly, Between the heating events, the motion of the disks is free, which enables us to implement an event-driven molecular dynamics scheme with a linked-list method [32]. The CPU time however scales like N 2 because the lists in all cells need to be updated after each heating event.
Four parameters determine the state of the system: the inelasticity ǫ = 1 − α 2 , the packing fraction φ = πNσ 2 /(4L 2 ), the reduced lower wave number cutoff k min σ = 2πσ/L, and the heating rate ξ 2 0 . The values of N investigated in this article vary between 10 3 and 10 4 , and we shall restrict our attention to high packing fractions for which the use of linked lists implies the most significant reduction of computer time. We consider the cases of moderate inelasticities (0.6 < α < 1) and complete inelasticity (α = 0). For the latter case, inelastic collapse occurs, i.e. the collision frequency involving only a small number of correlated particles diverges, as observed first by McNamara and Young [29] in freely evolving fluids of inelastic hard disks. Also in our system, for high enough inelasticity the heating seems never sufficient to prevent the inelastic collapse. For α < 0.5, the inelastic collapse has been avoided by introducing a slight modification of collision rule (36), as proposed in [30]: in each collision, the velocities are first computed according to the standard procedure ; the relative velocity v * 12 is then rotated by a random angle smaller than a maximum value Θ (typically less than a few degrees), keeping the center of mass velocity fixed. Note that this modified collision rule does not change the total energy loss of the colliding pair, and does not introduce any spurious drag or forcing on the particles.
The structure factors in the NESS have been computed for wave vectors compatible with the periodic boundary conditions, i.e. of the form (2π/L)(n x , n y ). We have obtained the density-density structure factors and the velocity-velocity structure factor, defined as, Here the averages are taken in the spatially uniform NESS. The fluctuation δg α in the momentum density, δg α (k) = i mv iα exp(−ık · r i ), and those in the flow field are related as δg α = ρδu α , where ρ is the average mass density in the steady state. The second rank tensor S αβ (k) is isotropic and can be split into a longitudinal and transverse part, wherek = k/k. From Eq. (39), it appears that the knowledge of S nn requires the computation of an O(N) quantity. We can rewrite the velocity-velocity structure factor so that its computation also increases linearly with the number of particles. For example, for the longitudinal part of S αβ (k) in (41) we have, In practice, the different S aa (k) have been computed for every k lying in the disk |k| < k max = 6π/σ, then averaged over shells of thickness k min = 2π/L, to achieve better accuracy. Moreover, the statistics in the NESS has been increased by averaging over time. Note that the above procedure, which gives insight into the microscopic to large scale structure of the system, does not require the knowledge of the hydrodynamic (coarse-grained) density and velocity fields.

A. Approach and characterization of the NESS
Before addressing the question of the large scale structure of the inelastic fluid, we investigate the validity of the macroscopic description given in section II and V. The former section gives the mean field results for the steady state temperature T E in (10) and collision frequency ω E in (A2), based on the Enskog-Boltzmann equation. The latter section and appendix C show how the long range spatial fluctuations renormalize these mean field values and lead to estimates in (35) for the corrections δT = T − T E and δω = ω − ω E , using a mode coupling calculation.
When the system is initially prepared in a configuration having a temperature T 0 different from the steady state temperature T E , the time dependence predicted by the mean field result (12) is in good agreement with the numerical data. This is shown in Fig. 3 for a system with small inelasticity. When the initial temperature T 0 is much larger than T E , the heating at short times is dominated by the inelastic dissipation, and Eq. (12) becomes Haff's homogeneous cooling law for a freely evolving system where t 0 = 1/ω E (T 0 ) is the mean free time in the initial state. This can be seen from the asymptotic expansion of the function f defined in (13), These analytic results for short times are confirmed by MD simulations, as shown in Fig. 4. Moreover, for initial temperatures T 0 ≪ T E , Eq. (12) predicts at short times a linear increase of T , as in a heated fluid of elastic hard spheres. Our simulations confirm this behavior. Fig. 5 shows that the measured kinetic energy per particle T is larger than the temperature T E predicted on the basis of mean field theory. This effect is noticeable albeit small in the results reported in Figs. 3 and 4, which correspond to the nearly elastic limit. For the densities studied here, we observe (see Fig. 5) that the correction δT is positive and decreases with decreasing inelasticity. The positive excess in temperature is already present for small inelasticity (see e.g. Fig. 4), and vanishes as ǫ → 0. Note that the above results are at variance with those reported by Peng and Ohta [13], who find that T /T E does not depend on α. The measured collision frequency ω is also larger than the Enskog estimate ω E , as shown in Fig. 5. The excess δω increases with increasing inelasticity.
We first observe that the simulation data in Fig. 5, where both ω > ω E and T > T E , cannot be explained consistently as a possible over-or underestimation of the IHS pair correlation function χ IHS at contact by its value for elastic disks. On the basis of (10) and (A2) we note that an overestimation of χ IHS would increase the collision frequency (A2), and decrease the temperature (10). These trends are at variance with the observations. The discrepancy between the measured and predicted temperatures and collision frequencies is more likely due to large scale fluctuations, at least for not too large inelasticities (α > ∼ 0.6). These long range spatial fluctuations renormalize the mean field Enskog values for the temperature T E and collision frequency ω E by amounts δT and δω. The theoretical estimates (35), based on mode coupling arguments, show the correct trends for the dependence of these corrections on the inelasticity, i.e. δT ∼ ǫ and δω ∼ ǫ 2 as ǫ → 0, and give a rough estimate of the magnitude of these terms as illustrated in the inset of Fig. 5. We also point out that the system size considered in Fig. 5 is much too small for the asymptotic theory (35) to be applicable. For instance at α = 0.8, we deduce from the data (l 0 = 1.1σ, L = 60σ) in the caption of For large inelasticity (α < 0.5), the temperature T and collision frequency ω depend on the maximum angle Θ of the random rotations, used to avoid the collapse singularity, as explained in section VI. In fact, ω diverges at Θ = 0 (inelastic collapse), in agreement with the observations of Peng and Ohta [13]. Intuitively, one might expect that a larger randomization of the postcollision velocities (larger Θ), would more effectively destroy the correlations leading to inelastic collapse, and consequently would decrease the deviations in T and ω from the mean field and the mode coupling predictions. This happens indeed for the temperature, but not for the collision frequency. At packing fraction φ = 0.07, the collision frequency is indeed monotonically decreasing with increasing Θ (see Table II). However, at φ = 0.2 it diverges at Θ = 0, decreases to a finite value at Θ = 5 o , and increases again to its value at Θ = 10 o (see Fig. 5). We have no explanation for this behavior.

B. Fluctuations in the NESS
In this section we analyze the effects of inelasticity on the large distance behavior of the fluid. We have computed the structure factors in the spatially homogeneous NESS of the inelastic hard disk fluid as explained in the previous section, and have focused either on values of α close to the elastic limit, where the theoretical description is supposed to apply, or on values close to 0, in order to test how large deviations from the theory might be. Local mean field values, like T E and ω E , reach their steady state values rapidly. However the time scale needed for the structure factors S ab (k) and the contributions of spatial fluctuations, δT and δω, to reach their steady state values are diffusive, and increase as k −2 min ∼ L 2 with system size. We have checked in the simulations that the large scale behavior of the structure factors was properly equilibrated before accumulating the data used to compute the averages.
First of all, we have tested the isotropy of tensor (40) by checking that the average vanishes for k values compatible with the periodic boundary conditions. In Fig. 6 we show the density-density and the relevant components of the velocity-velocity structure factors. For elastic hard disks the plateau values of S aa (k) around kσ ≃ 2 extend all the way down to k = 0. The excess correlations in the dissipative regime (k < ∼ γ 0 /l 0 ), which for S ⊥ (k) extend up to k ≃ √ γ 0 /l 0 , are characteristic of the randomly driven inelastic fluid. Fig. 6 shows that for small inelasticities the agreement between simulations and theory is quite reasonable. The structure factors diverge at small scales like k −2 , in agreement with the theoretical predictions (27)- (29). The packing fraction has been chosen fairly high (φ = 0.63) but lower than the two-dimensional random close-packing of monodisperse disks φ RCP ≃ 0.82 [33], and inside the liquid region of the phase diagram for elastic hard disks.
In addition to the gain in computer time, such a packing fraction leads to a mean free path l 0 smaller than the particle diameter σ (e.g. l 0 ≃ 0.095σ for φ = 0.63). Therefore, the hydrodynamic regime will hold up to typical particle diameters, enlarging the range of wave vectors where comparison between simulations and theoretical predictions is feasible. Moreover, for dense systems, a marked density structure is to be expected at the molecular scale, especially at wavelengths close to σ (kσ ≃ 2π). Fig. 7 shows that for α close to 1, this structure is indistinguishable from the structure for elastic hard disks. Note that S ll (k) and S ⊥ (k), although quite structureless for kσ > 1, show a weak and broad peak correlated with the maximum of S nn (k). This feature is more pronounced as the inelasticity increases, as shown in Figs. 10 and 13. As the molecular structure of the fluid has not been taken into account in the long wavelength hydrodynamic approach of section IV, the present theory cannot explain the structure in Figs. 10 and 13, and the structure factors predicted by the theory reach a plateau in the elastic regime (kl 0 > ∼ √ γ 0 ), given by where χ T is the isothermal compressibility of the elastic hard disk fluid. For α close to 1, the above limiting behavior is observed numerically for S ll or S ⊥ , and for S nn only in the limit of small packing fraction, where the molecular structure disappears. When the inelasticity is increased, the structure factors exhibit the same k −2 behavior at large scale, but the theoretical expressions are less accurate (see Fig. 8). However, the theoretical curves are based on the Enskog estimate T E for the granular temperature T , whereas for the system corresponding to Fig. 8, T ≃ 1.45T E . Our mode coupling theory predicts here T NESS ≃ 1.41T E . Fig. 9 displays the comparison between theory and simulation when the measured granular temperature is taken as an input for the hydrodynamic description. It appears that the large scale correlations (for which the present theory has been constructed) are well described by the theory, as long as the temperature is corrected from the mean field Enskog prediction (10) to the measured value T . In the case of S ⊥ , the amplitude only depends on the shear viscosity. The good agreement of the amplitude when the temperature is rescaled, while keeping for the shear viscosity the elastic hard disk value, suggests that the dependence of the shear viscosity on the inelasticity could be attributed only to the change in temperature. At the molecular scale, S ll (k) and S ⊥ (k) appear to be correlated to the density-density structure factor (see Fig. 10), in marked contrast to the elastic situation, where a plateau value would be reached. Such an effect is beyond the scope of our hydrodynamic approach, and is currently under investigation.
Surprisingly, in the case of complete inelasticity (α = 0), the theoretical structure factors give a reasonable picture of the large wavelengths correlations in the fluid (especially for the longitudinal velocity component S ll (k), as shown in Fig. 11). When the theoretical structure factors are deduced from the measured granular temperature T ≃ 1.5T E and not from T E (our mode coupling theory predicts T NESS ≃ 1.77T E ), the agreement for S ⊥ (k) improves, but the mismatch for S ll (k) increases (see Fig. 12). At small scales, the density correlations differ significantly from the elastic ones (Fig. 13) and the velocity structure factors exhibit the oscillatory behavior already present in Fig. 10, with peak positions locked in on the peaks in S nn (k). As can be expected from Figs. 11 and 12, the large scale correlations are compatible with the expected k −2 law (see Fig. 14) unlike the results of Peng and Ohta who report a k −1.4 asymptotic behavior for α = 0. However, a log-log plot such as Fig. 14 does not allow an accurate evaluation of the scaling exponents. The k −2 law is better inferred from the direct comparison with theory (Figs. 6,9,12).
Before concluding this section we compare the structure factors for the longitudinal flow fields, obtained from MD simulations with two different theoretical predictions in Fig. 2, obtained by including or excluding internal noise. First observe that all parameters in Fig. 2 and Fig. 6 are identical, as well as units on both axes. The simulation results for S (k) in Fig. 6 are in excellent agreement with the theoretical prediction (dashed line), which corresponds to the solid line in Fig. 2 (internal plus external noise). The dotted line (without internal noise) for S without (k) in Fig. 2 disagrees with the simulations in the relevant interval 0.2 < ∼ kσ < ∼ 0.5. Hence, inclusion of internal noise extends the validity of the asymptotic theory to intermediate wave numbers.

VIII. CONCLUSION
We have presented a theory for the large scale dynamics of a granular fluid that is driven into a nonequilibrium steady state by a random external force. Our description combines the macroscopic equations of motion for the hydrodynamic fields, accounting for energy dissipation through inelastic collisions and uniform heating, together with the fluctuating forces. The long range character of the spatial correlation functions is determined by the small wave number divergence ∼ k −2 of the corresponding structure factors. This k −2 behavior is typical for systems which combine conserving deterministic dynamics (conservation of particle number and momentum in collisions) with nonconserving noise [34], thus violating the fluctuation-dissipation condition, and is generic for rapid granular flows that are driven by external noise. We also draw attention to the analogy of our equations of motion for the fluctuating fields to the Edwards-Wilkinson model [35] that was proposed for growth of a granular surface. In that case the dynamic variable is a scalar field, namely the height of the surface, which obeys a similar equation of motion as any of the (d − 1) components of the transverse velocity field, u ⊥α (k, t), in our case. The only difference is that in the Edwards-Wilkinson model there is only nonconserving noise, whereas in our case both nonconserving (external) and conserving (internal) noise are present.
We have tested our predictions for the structure factors against molecular dynamics simulations and have demonstrated that there is quantitative agreement over the whole wave number range, if internal fluctuations are taken into account. In our two-dimensional simulations we have found deviations in the steady state temperature and collision frequency that grow with the inelasticity and the system size. For not too large inelasticities (α > ∼ 0.6), we have explained these deviations in terms of mode coupling effects of the long range fluctuations.
The phenomenological mode coupling theory, proposed in [20] and extended here to the driven IHS fluid, starts from the same basic ingredients as in the case of elastic fluids [36], where that theory was used to calculate the long time tails of the velocity autocorrelation function and other current-current time correlation functions. For the elastic case the mode coupling theory can be derived from the ring kinetic theory in the low density limit [37], which accounts for dynamic correlations built up by sequences of correlated binary collisions, leading to nonlocal effects in space and time. Such collision sequences correct the mean fieldtype Boltzmann or Enskog kinetic equations for the errors induced by the breakdown of the molecular chaos assumption. The ring kinetic theory for rapid granular flows of IHS has been developed in Ref. [38], but has not yet been used to derive the present phenomenological mode coupling theory from the more fundamental kinetic theory, valid in the low density limit.
In detailed balance models, such as elastic hard spheres, dynamic correlations created by correlated collision sequences lead to long time tails, which imply that transport coefficients in two dimensions diverge as ln L for large systems [36,37].
Non-detailed balance models, such as IHS fluids, generically exhibit long range spatial correlations [19]. In randomly driven IHS fluids, as studied in this paper, these correlations between densities and flow fields at distant points in the fluid behave as 1/r in 3D and ln r in 2D, and already modify (renormalize) the mean field Enskog-Boltzmann values for steady state properties, such as the temperature and collision frequency. In 2D systems these renormalization corrections, δT and δω, exhibit the ln L divergence, which in the case of detailed balance models appears only in the transport coefficients [19].
At larger inelasticity (α < ∼ 0.6), molecular chaos is also violated due to the presence of short range velocity-velocity correlations (see Fig. 13), which are beyond our mode coupling theory. A detailed investigation of the small scale structure, which for large inelasticity clearly deviates from an equilibrium structure, will be reported in a subsequent publication. It is surprising that our description, which is based on the Enskog theory and neglects any dependence of transport coefficients on inelasticity, even at α = 0 predicts the long range structure reasonably well, provided that the temperature is not taken as the mean field APPENDIX A: In this appendix we derive expressions for the transport coefficients that govern the decay of fluctuations in the steady state. Linearization of the macroscopic equations (6) around the NESS, defined in (9) and (10), gives the deterministic part of the following set of equations: The noise termsξ(r, t) andθ(r, t) have been discussed in section III. The pressure p is assumed to be that of EHS, p = nT (1 + Ω d χnσ d /2d), where Ω d = 2π d/2 /Γ(d/2) is the d-dimensional solid angle, and χ(n) the equilibrium value of the pair correlation function of EHS of diameter σ and mass m at contact. The kinematic and longitudinal viscosities ν and ν l , as well as the heat conductivity κ are also assumed to be approximately equal to the corresponding quantities for EHS, as calculated from the Enskog theory [22], where ρν = η and ρν l = 2η(d − 1)/d + ζ are expressed in shear viscosity η and bulk viscosity ζ. The collisional energy loss in (7), Γ = 2γ 0 ωT , is proportional to the collision frequency as obtained from the Enskog theory. In two dimensions we use the Verlet-Levesque approximation χ = (1 − 7φ/16)/(1 − φ) 2 . By taking spatial Fourier transforms δa(k, t) = drδa(r, t) exp(−ık · r) in (A1), one obtains the mesoscopic equation (18) with the hydrodynamic matrix It contains the coefficients In the body of the paper we need the eigenvalues z λ (k) of the asymmetric matrix M ab (k), and its right and left eigenvectors, which are obtained from where M T is the transpose of M. Here λ = ± labels the sound modes, λ = H the heat mode, and λ =⊥ labels (d − 1) degenerate shear or transverse velocity modes. The eigenvectors form a complete biorthonormal basis, i.e. Consequently all modes of the heated IHS fluid are linearly stable. There is no clustering instability [23,28] and no instability in the flow field [23], as in the freely evolving IHS fluid.
In the dissipative regime (kl 0 ≪ γ 0 ) the eigenvalue equation is solved by an expansion in powers of k, and one finds to dominant orders the eigenvalues in the form (25) and eigenvectors in the form (26). The eigenmodes to dominant nonvanishing order in k are listed in (26). There are (d − 1) transverse velocity or shear modes (λ =⊥), which are purely diffusive with a diffusivity D ⊥ = ν; there are two propagating modes (λ = ±) with a speed of propagation v D , and sound damping constant D S , and a kinetic heat mode (λ = H) with a nonvanishing z H (0).
For later reference we also express these coefficients in thermodynamic quantities and transport coefficients. The speed of sound v D in the dissipative regime (kl 0 ≪ γ 0 ) is: It satisfies the inequality v D < v S , where v S is the adiabatic speed of sound in the standard regime (γ 0 ≪ kl 0 < 1), The damping constant of the sound modes is and the dispersion relation for the kinetic heat mode contains the positive constant The ratios v 2 D /T , v 2 S /T , D S / √ T , and D H / √ T are independent of temperature.

APPENDIX B:
In this appendix we calculate the tails of the spatial correlation functions G ab (r), which is done by Fourier inversion of S ab (k). Consider first the tensor fields G αβ (r) and S αβ (k) in Eq. (1) with δa α (r, t) = u α (r, t) (α, β = x, y, . . .) being the components of the flow field. Both tensor fields are isotropic and can be split into longitudinal and transverse components, i.e.
where the small k behavior of S ⊥ (k) and S (k) is given in (27) and (28). By contracting the second line above withr αrβ we obtain Contraction of the second line of (B1) with (δ αβ −k αkβ ) yields in a similar manner an expression for G ⊥ (r).
In three dimensions the k integral can be performed explicitly and yields dk (2π) 3 exp(ık · r) The resulting large r behavior of G αβ (r) is given in (30).
In two dimensions the integral in (B2) over the azimuthal angle yields for large r dk (2π) 2 exp(ık · r) In two dimensions the k integral diverges for k → 0. It should in fact be restricted to k ≥ k min = 2π/L, which is the smallest allowed wave number when periodic boundary conditions are imposed. To obtain the last equalities one needs the small z behavior of the Bessel functions J ν (z) ≃ (z/2) ν /Γ(ν + 1) [39]. Combining these results with (B2) yields for large r The remaining spatial correlation functions G ab (r) involving a, b = {n, T }, are scalar fields. Their large r behavior is given by where B ab (k) is given in Table I. The subleading corrections of O(r 0 ) depend on the cutoff k min . To evaluate these terms requires the subleading small k behavior of S ab (k) in Eqs. (27) to (29).

APPENDIX C:
In this appendix we present a mode coupling calculation to estimate the contributions of the long wavelength fluctuations in the NESS to some quantity h. Examples are the particle distribution functions, the energy per particle E, and the collision frequency ω. The fluctuations are correlated over large distances, as a consequence of sequences over dynamically correlated collisions, the so called ring collisions [40], as shown in section V by explicit calculation of their spatial tails.
To calculate their contributions to average quantities like h, one may solve the ring kinetic equations [38], or estimate these quantities from a more phenomenological mode coupling approach, as developed in Ref. [36]. The basic assumption made there is that the state of the system rapidly decays to a state of local equilibrium, described by the fluctuating hydrodynamic fields a(r, t) = {n(r, t), u α (r, t), T (r, t)}.
For the quantity under consideration this can be implemented by representing h = (1/V ) drh(r), and approximating h(r) by its value in local equilibrium, i.e.
where the average is taken over the fluctuating hydrodynamic fields a(r) in the NESS. To carry out the average over the fluctuations we expand h l (a) in powers of the fluctuations δa = a − a around the NESS, yielding where summation convention for repeated indices has been used. Here A ab is the matrix of second derivatives ∂ 2 h l (a)/∂a∂a at a = a . The asterisk indicates that k integrals are restricted to the long wavelength range, k < γ 0 /l 0 , the so called dissipative range, discussed below (16). In this range the structure factors have the form S ab (k) ≃ E ab /k 2 on account of (27) to (29). For dimensionality d ≥ 3 the fluctuation contribution, δh = h NESS −h l ( a ), is convergent at small k, and gives only small well-behaved corrections to h l ( a ). However, for d = 2, the k integral diverges logarithmically at small k (where k > ∼ 2π/L), and the excess, δh, is given by Consequently, the fluctuation contribution δh in 2D systems is a singular function of the system size, L, that diverges in the thermodynamic limit. We first apply the above results to the energy per particle, E = (1/N) dre l (a(r)), where e l (a) = 1 2 ρu 2 + d 2 nT is the energy density in local equilibrium. From e l (a) the expansion coefficients corresponding to h l ( a ) and A ab in (C2) can be calculated, to yield in d = 2: δE ≃ 1 2n * dk (2π) 2 [ρS ⊥ (k) + ρS ll (k) + 2S nT (k)] .
Inserting (27) to (29) then yields δE ≃ γ 0 ω E T E 4πn where the coefficient B nT is listed in Table I. For the collision frequency the analog of (C1) is ω = (1/N) dr n(r)ω l (a(r)) with ω l (a) ∝ nχ(n) √ T given in (A2). This gives in two dimensions: where the coefficients B ab are given in Table I, and A nn = ∂ 2 (nω)/∂n 2 , A nT = ∂ 2 (nω)/∂n∂T , and A T T = ∂ 2 (nω)/∂T 2 at a = a . The same method can be used to calculate other averages, as well as particle distribution functions. For instance, for the single particle distribution function the starting point would be and the above procedure can be applied at once.   ) as a function of Θ, for a totally inelastic system (α = 0) with N = 1600 particles, and packing fraction φ = 0.07.