Randomly Driven Granular Fluids: collisional statistics and short scale structure

We present a molecular dynamics and kinetic theory study of granular material, modeled by inelastic hard disks, fluidized by a random driving force. The focus is on collisional averages and short distance correlations in the non-equilibrium steady state, in order to analyze in a quantitative manner the breakdown of molecular chaos, i.e. factorization of the two-particle distribution function, $f^{(2)}(x_1,x_2) \simeq \chi f^(1)(x_1) f^{(1)}(x_2)$ in a product of single particle ones, where $x_i = \{{\bf r}_i, {\bf v}_i \}$ with $i=1,2$ and $\chi$ represents the position correlation. We have found that molecular chaos is only violated in a small region of the two-particle phase space $\{x_1,x_2\}$, where there is a predominance of grazing collisions. The size of this singular region grows with increasing inelasticity. The existence of particle- and noise-induced recollisions magnifies the departure from mean field behavior. The implications of this breakdown in several physical quantities are explored.


I. INTRODUCTION
The interesting phenomena observed in recent experiments with mono-and multi-layers of granular material on vibrating plates [1][2][3][4] show the need to develop kinetic theories for rapid granular flows with mechanisms for energy input, different from those in shear flows or flows through vertical pipes. In the present article the fluidization is driven by a random external force, which gives frequent kicks to each particle in between collisions. Such a driving mechanism has recently been studied by many authors [5][6][7][8][9][10][11]. The basic physical interest is the understanding the non-equilibrium stationary states (NESS) that exist in the presence of this random force. The advantage of this fluidization mechanism, besides its potential physical realizations, lies in the fact that the NESS is linearly stable against spatial inhomogeneities.
In Ref. [10], to which we refer as paper I, we have studied the large scale structure and presented a hydrodynamic description of randomly driven granular fluids, modeled as systems of smooth inelastic hard spheres (IHS). The IHS model accounts for two essential features of granular matter: hard core exclusion and dissipative collisions [12]. The dynamics is described by a constant coefficient α of normal restitution. In collisions a fraction of the relative kinetic energy is lost, which is proportional to the inelasticity ǫ = 1 − α 2 . The stochastic external force compensates this energy loss, and drives the IHS fluid into a NESS. This stationary state, though homogeneous and stable against spatial fluctuations on large space and time scales (at least for weakly inelastic spheres), was shown to exhibit long range spatial correlations in density, velocity and granular temperature fields that extend much beyond the mean free path. In fact, the corresponding structure functions S(k) diverge as 1/k 2 as the wave number k → 0, a behavior caused by the random external force, which does not conserve momentum whereas the collisions between particles do. These long range correlations are of algebraic form, ∼ 1/r d−2 , which corresponds to ln r in two dimensions (d = 2). The existence of such extremely long range spatial correlations is one example of the many nontrivial properties of non-equilibrium stationary states in general [13,14].
Differences in the stationary states between fluids with dissipative and conservative interactions also manifest themselves in the kinetic properties of the fluid, such as the velocity distribution function, which deviates from a Maxwellian in particular in the high energy tail of the distribution. In Ref. [9] the existence of an overpopulated high energy tailf ∼ exp[−Cv 3/2 ], where C is a constant that depends on the inelasticity, has been obtained from kinetic theory. A similar behavior has been observed experimentally at high vibrational accelerations [3,4]. This observation indicates that certain features of the experiment might be reproduced by modeling the input of energy into the horizontal motion of the beads by a random external force, although other energy injection mechanisms that could be relevant to recover the large velocity tail have been put forward [15]. In similar experiments [2] with a vertically vibrating plate covered with a mono-layer of steel balls with a packing fraction around 50% the velocity distribution of the horizontal velocities has been measured, and again overpopulated non-Gaussian high energy tails have been observed. In the present paper we will investigate the kinetic properties and short scale correlations that characterize the stationary state. More specifically, we will compare molecular dynamics (MD) simulations of inelastically colliding disks with analytic predictions based on the assumption of molecular chaos.
The Boltzmann equation for dilute gases of particles that interact via short-ranged repulsive interactions is based on the assumption of molecular chaos, also called the Stosszahlansatz or mean field approximation. It assumes that the velocities of colliding particles just before collisions are uncorrelated, i.e. their pair distribution function factorizes,f (2) (x 1 , x 2 , t) =f (x 1 , t)f (x 2 , t), where x i = {r i , v i } denotes the position and velocity of particle i. Enskog's extension of the Boltzmann equation to a dense system of hard spheres [16], referred to as Enskog-Boltzmann equation, is also based on the fundamental assumption of the absence of velocity correlations. Here the assumption of molecular chaos postulates thatf (2) (x 1 , x 2 , t) = χf (x 1 , t)f (x 2 , t) for approaching particles (v 12 · r 12 < 0) just before collision (r 12 = σ + 0), where χ is assumed to be the radial pair distribution function at contact g(r 12 = σ + 0) in local equilibrium. It implies the additional assumption that spatial correlations between colliding particles just before collision are independent of their velocities, i.e. absence of position-velocity correlations. The Enskog χ factor enhances the collision frequency at higher densities. For dilute gases the assumption of molecular chaos seems to be justified. Recently Lutsko [17] and Soto and Mareschal [18] derived for a freely evolving inelastic hard disk fluid a relation between pre-and postcollision radial distribution function at contact, as a function of the angle, θ = cos −1 (v 12 ·r 12 ), between the relative velocity v 12 of the colliding particles, and their relative position at contact r 12 , and they confirmed their results by MD simulations. Their observations made it clear that further arguments are needed to clarify the meaning of the χ−factor in Enskog's formulation of the molecular chaos assumption. This will be done in section II A.
The breakdown of molecular chaos at higher densities in classical fluids with conservative forces has been extensively investigated in the sixties and seventies [19]. This breakdown is caused by sequences of correlated binary collisions, the so called ring collisions [20]. They lead to long time tails in velocity and stress autocorrelation functions [21,22], and to long range spatial correlations in NESS [13]. The quantitative effects of velocity correlations on transport coefficients at liquid densities are also significant. For instance, molecular dynamics simulations on elastic hard sphere systems at liquid densities [23] have shown that the long time tails increase the measured self-diffusion coefficient D typically by 15% to 20% with respect to the prediction of the Enskog theory, D E = D B /χ, where D B is the Boltzmann value of the self-diffusion coefficient.
A well-known example of short scale structure in granular fluids are the position-velocity correlations leading to the phenomenon of inelastic collapse [24,25], which is a divergence of the collision frequency ω in a finite time. The collapse singularity implies that an infinite number of collisions occurs within a finite time interval in a subset of (nearly) touching particles, ordered in linear arrays. The phenomenon is, however, an artifact of the assumption that the coefficient of restitution α is independent of the impact velocities, whereas on physical grounds α(v 12 ) → 1 (elastic limit), as the relative velocity v 12 vanishes. Molecular dynamics simulations have shown that the assumption of molecular chaos is also violated in undriven granular fluids in their late stages of evolution, the so called nonlinear clustering regime. For instance, the measured distribution of impact parameters is not uniform, as expected on the basis of molecular chaos, but biased toward grazing collisions [26][27][28].
As shown below, in the driven IHS fluid there is an important additional reason for the breakdown of molecular chaos, namely the strong increase in collision frequency at small relative velocities between two isolated particles, caused by the so-called noise induced re-collisions. This correction to the collision frequency, that is important at all densities, is also neglected in the molecular chaos assumption.
The main goal of this paper is to quantify, analyze and interpret the effects of the breakdown of molecular chaos in the NESS of inelastic hard spheres that are subject to a random external force between collisions. We will focus in particular on velocity-velocity correlations and position-velocity correlations between particles almost in contact, i.e. the short scale structure.
Section II presents the analytic results, based on the Enskog-Boltzmann equation, which has been modified to account for the external energy input. In section III we present molecular dynamics results for several quantities that characterize the collision processes and related short scale structure of the NESS, and make a comparison with predictions based on molecular chaos.

A. Molecular chaos and Enskog approximation
The Enskog-Boltzmann equation for the single-particle distributionf (v 1 , t) in a spatially homogeneous randomly driven fluid of inelastic hard spheres of diameter σ reads in d = 2 or 3 dimensions [9]: where v 12 = v 1 − v 2 and n the number density. The Heaviside function Θ(x) restricts theσ integration to the hemisphere v 12 ·σ > 0, whereσ is the unit vector along the line of centers of the colliding spheres at contact, pointing from particle 2 to 1. In the sequelâ = a/|a| denotes a unit vector. The gain term of the collision integral describes the restituting collisions that convert the precollision velocities (v * * 1 , v * * 2 ) into (v 1 , v 2 ), while the loss term describes the direct collisions, and contains the precollision velocities (v 1 , v 2 ) leading to postcollision velocities (v * 1 , v * 2 ). The postcollision and restituting velocities have been defined in [29]. The χ factor will be discussed below.
As derived in [9], the Fokker-Planck term accounts for the external energy input, and describes diffusion in velocity space with a diffusivity proportional to the rate of energy input ξ 2 0 per unit mass. Here ξ 0 is the strength of the random driving force, which is assumed to be Gaussian white noise [9,10].
Before studying the short scale structure, we consider the single particle distribution functionf (v) in the NESS. The stationary solution of (1) is determined by the balance, mξ 2 0 = Γ, of external heating, mξ 2 0 , and internal loss of energy due to collisions, Γ. It is characterized by a time independent temperature, T = mv 2 /d , defined as the average kinetic energy per particle, and discussed in paper I. As mentioned in the introduction, this stationary solution exhibits an overpopulated high energy tailf ∼ exp[−Cv 3/2 ]. The structure of the tail distribution is determined by collisions of very energetic particles with 'thermal' particles, and can be obtained by neglecting the gain term in the Boltzmann equation [9].
In Ref.
[9]f (v) has been calculated by solving the Enskog-Boltzmann equation (1) by an expansion in Sonine polynomials. To formulate this result it is convenient to introduce a rescaled distribution function where v 0 ≡ 2T /m is the thermal velocity and d the dimensionality. This gives where the Maxwellian ϕ(c) ≡ π −d/2 exp(−c 2 ). Note that a 2 is proportional to the fourth cumulant of the scaling form f (c), i.e.
and vanishes in the elastic limit. An explicit calculation to linear order in a 2 gives [9] In the next section this prediction will be tested against molecular dynamics simulations.
Consider first the exact expression for the mean collision frequency in the homogeneous NESS, defined as wheref (2) (v 1 , v 2 , σ) is the dynamic or constrained pair distribution function with velocities aiming to collide, just before contact with r 12 = σ. Molecular chaos, also referred to as mean field theory, requires the complete factorization of the dynamic precollisional pair function, What is the meaning of the χ−factor, used in formulating the molecular chaos hypothesis? This hypothesis for dilute gases concerns the absence of correlations in precollision velocities, and in precollision positions (χ = 1).
In dense fluids on the other hand, the precollision position correlation χ is different from 1, but the precollision velocity-velocity and position-velocity correlations are still assumed to be absent.
In the literature it is common to take χ equal to the radial distribution function at contact in local equilibrium, i.e. χ = χ E ≡ g eq (r → σ + 0), which mainly accounts for precollision hard core exclusion effects. For hard disks and hard spheres the latter function is approximately given by the Verlet-Levesque (2D) and Carnahan-Starling (3D) approximations [31], where the packing fraction in d dimensions is defined as φ = n(σ/2) d Ω d /d, and Ω d = 2π d/2 /Γ(d/2) is the surface area of a d-dimensional unit sphere. In this paper we refer to the molecular chaos approximation with χ = χ E , as the Enskog approximation.
In principle, different options are open for the χ−factor. Asf (2) is the dynamic precollision pair distribution function at contact, an alternative choice for the χ in the factorized form could be the dynamic precollision radial distribution function at contact, defined as an average over the precollision hemisphere, Another option could be the unconstrained radial distribution, g(r), in the NESS, extrapolated to contact (r → σ). This function is further discussed in subsection II D.
For the randomly heated fluid under study here, the dynamics is not purely hard-sphere like. The random force acting on the particles may be expected to contribute to the value of the pair distribution function at contact. This effect will be addressed in the subsequent sections.
Equation (5) withf (v) replaced by the Maxwellian, yields for the collision frequency in the molecular chaos approximation ω mc (T ) = χω 0 (T ) , and more specifically in the Enskog approximation, Here the Boltzmann collision frequency for dilute gases is given by, and the small correction of O(a 2 ) appearing in (2) has been neglected. Spatial correlation functions in non-equilibrium stationary states are quite different from local equilibrium ones, and show long range correlations due to correlated sequences of ring collisions, also referred to as mode coupling effects [13]. In paper I we have shown the existence of very long range correlations ∼ 1/r d−2 in the randomly driven IHS fluid. The short range correlations in the NESS can in principle be obtained from the ring kinetic equation for IHS (see Ref. [29]). However, systematic methods to evaluate collision integrals and pair correlation functions at short distances using this ring kinetic theory have not yet been developed. In the section on simulation results, we return to the effects of ring collisions, and present arguments why their contributions are expected to be more important here than for elastic hard spheres.

B. Collisional averages
In hard sphere systems there are many properties that involve the pair distribution function of particles just before collision. To study these, it is convenient to introduce the collisional average ... coll for a quantity A in the NESS, defined as In the sequel it is more convenient to work with a rescaled pair distribution functionf (2) To express the collisional averages (11) in rescaled variables one replacesf (2) by f (2) These objects can be conveniently computed in event driven molecular dynamics algorithms for hard sphere systems [30]. Collisional averages are defined for particles that are about to collide (i.e. |r 12 | = σ + 0), and can be calculated from kinetic theory using the molecular chaos assumption, possibly supplemented with the Enskog approximation at higher densities.
Collisional averages of great importance are the collisional energy loss per unit time, d 2 nΓ, and the excess hydrostatic pressure, p − nT , resulting from collisional transfer of momentum. With a minor generalization to d dimensions we obtain from Ref. [32] the exact expression for the pressure in the NESS: The second equality is obtained by introducing the collisional average (11) and expressing its denominator in terms of the collision frequency given by (5). In fact, inserting (6) into the first line of (12) allows one to carry out theσ−integration, and the right hand side becomes proportional to the rescaled velocity average dc 1 dc 2 g 2 f (c 1 )f (c 2 ) = 2 without any further assumption about neglecting the term proportional to a 2 in (2). This argument is special for the pressure, as other moments involve the knowledge of the complete f (c). Indeed, the generic collisional average becomes |g ·σ| m coll = 2 m/2 Γ( 1 2 m + 1), independent of dimensionality, assuming molecular chaos (6) and replacing f (c) by the Maxwellian ϕ(c) (the contributions coming from a 2 are quite small and can be neglected; they have been computed in [9]). Finally, the pressure can be expressed as, Different choices for χ yield different approximations. For instance, with χ = χ E we obtain the Enskog approximation p E (T ) for the pressure of IHS.
In the elastic limit p E (T ) at α = 1 gives the standard equation of state for elastic hard disks or spheres. Notice that the pressure for IHS is only defined kinetically as the momentum flux, which leads to (12). A statistical mechanical derivation of the equation of state from the partition function or free energy for the IHS fluid does not exist.
In a similar manner we obtain the exact expression for the collisional damping rate, where γ 0 = (1 − α 2 )/2d is the dimensionless damping constant introduced in Refs. [9,10]. The last equality (14) expresses the balance between the energy input due to the white noise, and the collisional loss of energy in the NESS, and determines the temperature T in the NESS. By specializing this equation to the Enskog mean field approximation, f (2) = χ E f f , where |g ·σ| 2 coll = 2, we obtain the approximate result, and similar relations for different choices of χ. It is convenient to define a reference temperature T E through the relation, or more explicitly Moreover the definition of T E combined with the NESS condition Γ(T ) = mξ 2 0 implies the relation Γ(T ) = Γ E (T E ), and consequently, In the sequel we will also use a reference frequency ω E without any argument to denote, Although we will discuss the simulation results in detail in the next section, it is of interest already at this point to note that for these systems, as shown in Fig. 1, the ratio of the kinetic temperature and the reference temperature, T /T E , is only somewhat larger than 1 for all α, that it approaches 1 in the elastic limit (α → 1), and that it monotonically increases with decreasing α (see Fig. 1). The same figure shows that the ratio, ω/ω E also approaches 1 for α → 1, with a steep increase to a value 5.6 as α → 0. Further discussion of these points is postponed till section III. If the Enskog factorization,f (2) = χ Eff , would be exact, then T = T E and ω = ω E . A third quantity of interest, the precollisional χ (−) factor, defined in (8), can also be expressed as a collisional average using (11), Before closing this section a caveat about internal consistency is appropriate. To obtain consistent theoretical predictions for the pressure p or dissipation rate Γ, it is paramount that both factors, ω and |g · σ| m coll , be calculated using identical approximations for f (2) . For instance, the mean field or molecular chaos approximation for the dissipation rate, Γ E (T ) = 2γ 0 ωT , -an expression commonly used in granular hydrodynamic equations -should necessarily be combined with ω E (T ) in (9). Any improved theoretical calculation for ω without a concomitant correction to the mean field result for |g · σ| m coll is inconsistent .

C. Velocity distributions
We study a variety of collisional averages and corresponding probability distributions. By choosing A(c 1 , c 2 , σ) = δ(|g| − g) we obtain the probability P r (g) that two colliding particles have a relative speed |c 12 | = g. From here on we only quote results for two dimensions. Analytic calculations are based on the molecular chaos assumption (6) in combination with (2). Inspection of (11) shows that under this assumption the collisional averages are independent of the χ−factor. Straightforward algebra gives for the constrained g−distribution, Similarly we obtain the probability distribution for the center of mass velocity, G ≡ 1 2 (c 1 + c 2 ), It equals the unconstrained equilibrium distribution function apart from a small term of O(a 2 ). Furthermore, the probability that the precollision speed |c i | of one of the colliding particles (i = 1, 2) has a value v is whereas the unconstrained distribution is ∼ v exp(−v 2 ). In evaluating this collisional average we neglect the a 2 contribution, and carry out the constrainedσ integration. To calculate the remaining integral dc 2 c 12 ϕ(c 2 ), we change integration variables to g expressed in polar coordinates {g, φ}, and use the relation π 0 dφ exp(−2c 1 g cos φ) = πI 0 (2c 1 g). The subsequent g−integration follows from (6.618.4) in Ref. [33]. Using the asymptotic expressions I 0 (x) ∼ I 1 (x) ∼ exp(x)/ √ 2πx for the modified Bessel functions of the zeroth and first order, we obtain the high energy behavior, In a similar manner we obtain the following velocity moments and correlations, using the molecular chaos assumption, Here c * i are the postcollision velocities, as defined in paper I. The sum of the third and fourth equality depends only on the center of mass velocity, i.e. G 2 coll . In the elastic limit α → 1, the average energy of a particle that is about to collide, c 2 coll = (5/4) c 2 , is above the mean energy per particle, c 2 , which equals unity.
In the molecular chaos approximation an average like (c 1 · c 2 ) m g n coll with {m, n} integers, is in general non-vanishing, except in the special case n = −1. Then (c 1 · c 2 ) m /g coll reduces to an unconstrained average, proportional (c 1 · c 2 ) m , which vanishes for odd values of m. Additional information about the relative orientation of the incoming velocities can be obtained from the distribution of the angle ψ 12 , defined by c 1 · c 2 = c 1 c 2 cos ψ 12 . A numerical calculation (again neglecting a 2 corrections) gives cos ψ 12 coll ≃ −0.233, which is close to the value −0.2, estimated from c 1 · c 2 coll ≃ c 2 1 coll cos ψ 12 coll using the above results. A very sensitive probe for studying the violation of molecular chaos is the probability distribution P (b) of impact parameters, b = |ĝ ×σ| = sin θ, where θ = cos −1 (ĝ ·σ) is the angle of incidence. It is defined as the collisional average and P (b) can be easily computed in a molecular dynamics experiment. As long as molecular chaos holds, the distribution of b is independent of the functional form of f and we obtain straightforwardly which reduces in two dimensions to the uniform distribution, In order to analyze molecular chaos breakdown in more detail, we have introduced a collection of moments M nm and their dimensionless counterparts B nm for n, m = {0, 1, 2 . . .}, (see definition below), to analyze in detail the possible breakdown of the molecular chaos factorization (6). These moments M nm (T ) of the pair distribution at contact are defined as, which are averages over the precollision hemisphere, where θ = cos −1 (ĝ ·σ). Let M E nm (T E ) denote the same quantity evaluated in Enskog's formulation of the molecular chaos approximation, and evaluated at the reference temperature T E , i.e. evaluated withf (2) where We prefer to normalize the reduced moments by M E nm (T E ), because its analytic form is given explicitly. One could also normalize by is that the computation requires the simulated values of the kinetic temperature T . The collisional averages v n 12 | cos θ| m coll expressed in terms of these new moments give, We first observe that the average collision frequency ω, defined in (5), is proportional to M 11 (T ), so that with ω E defined in (19). This implies that the reduced moments B nm (T ) can all be expressed in collisional averages, i.e.
The average · · · E coll is defined through (30) with M nm (T ) replaced by M E nm (T E ), and calculated in (A.2). It represents the collisional average, evaluated with the Enskog factorization f (2) = χ E f f and also taken at the reference temperature T E . Note that the equality (32) consists of two factors, ω and · · · coll , which are measured separately in event driven MD simulations.
We also observe that the equality Γ(T ) = Γ E (T E ), explained above (18), implies that Furthermore we have for the excess pressure, p ex (T ) ≡ p(T ) − nT , and for the dynamic pair correlation at contact χ (−) , In the appendix we present a more complete set of relations for the B nm .
In the next section MD simulations will show that the predicted deviation from a Maxwellian (see (2)) in the (unconstrained) velocity distribution of a single particle can be observed for small inelasticities. However larger deviations are found between the observed constrained probability distributions and averages, and the corresponding kinetic theory predictions given by (24), based on molecular chaos. Consequently the small corrections resulting from a 2 in (2) can be neglected in most cases.

D. Radial distributions
The static or unconstrained radial distribution function in the spatially homogeneous IHS fluid is defined as, It may be averaged over all directions of r because of statistical isotropy. The unconstrained radial distribution function at contact is defined as the extrapolation, Y = g(r → σ + 0). By splitting theσ−integration into a precollision (ĝ ·σ < 0) and postcollision hemisphere (ĝ ·σ > 0), we obtain Y as sum of two terms, The definitions of Y (−) , Y (+) follow from (36) by adding respectively factors Θ(−ĝ ·σ) and Θ(ĝ ·σ) under the integral sign in (36). The dummy integration variables in Y (+) represent the postcollision velocities, (c * 1 , c * 2 ), corresponding to the precollision ones, (c 1 , c 2 ).
On the other hand, we have the dynamic precollision correlation χ (−) , defined in (8), and a similar postcollision one, χ (+) , defined by replacing Θ(−ĝ ·σ) in (8) by Θ(ĝ ·σ). They are related by continuity of flux. Because the incident flux of (c 1 c 2 )−pairs just before collision is equal to the scattered flux of (c * 1 c * 2 )−pairs just after collision, we have inside dynamic averages the equality, where g n = g ·σ = g cos θ. The reflection law, g * n = α|g n |, for inelastic collisions in combination with the continuity of the flux and (8) yields at once, In principle, equations (37) and (20) provide two alternative routes to compute the precollisional pair correlation at contact: the first one, denoted by Y (−) , can be implemented numerically as a static or unconstrained average, namely by extrapolation to contact of the pair correlation function for pairs aiming to collide. The second one, denoted by χ (−) , can be computed as a dynamic collisional average, calculated from f (2) (c 1 , c 2 , σ) at contact. It is important to stress that the dynamic χ (−) is calculated as a time average over the subset of colliding pairs at contact, and the static Y (−) as a time average over all pairs, satisfying the relation, g ·σ < 0 and extrapolated to r → σ + 0, i.e. calculated from f (2) (c 1 , c 2 , r), where the limit is taken after all integrations have been performed. This may lead to different results, because the integrand contains the the function f (2) which turns out to be singular near r = σ and v 12 small (see discussion in subsection III C). Physically, it is also clear why the averages in the NESS need not be the same. For instance, the relation (38) may not hold for the limiting (r → σ) values, Y (−) and Y (+) , because non mean-field effects (in particular, the 'rotation-induced' recollisions discussed at the start of section III, or noise induced recollisions, see below) may result in differences between the two methods to evaluate χ (−) and Y (−) . The reason is that the validity of (38), expressing flux continuity for the limiting values (r → σ + 0), is questionable in the presence of the external random force. When the kicking frequency is much larger than the collision frequency (situation considered here), a pair of particles may indeed be put in contact under the action of the random force only. We will investigate possible numerical differences between χ (−) and Y (−) in the next section on MD simulations.
In paper I we have studied the long range spatial correlation functions G ab (r) of the density field n(r, t) and the flow field u(r, t) in the NESS. These functions are closely related to (36), i.e.
and, in the notation of paper I, where · · · is an average over the N -particle non-equilibrium steady state and the static velocity correlation, c 1 · c 2 (r), is defined as, The correlation functions G ab (r) above are very long ranged, decaying like r 2−d for large distances. In the first part of this subsection we have introduced the static correlations, Y , Y (±) , and the dynamic ones, χ (±) . In (A.8) and (A.9) we have done the same for the dynamic counter parts c 1 · c 2 dyn and c 1 · c 2 (−) dyn of the static correlation c 1 · c 2 (r → σ), introduced in (41).
In the next section the short range behavior of these functions will be studied by MD simulations.

III. SIMULATION RESULTS FOR THE NESS
To investigate the short scale structure characterizing the NESS and the validity of molecular chaos we will present in this section MD simulation results, and compare these with our theoretical predictions whenever possible. The details of the simulations of the randomly driven inelastic hard disk system have been reported elsewhere [10]. We will work in the limit in which the kicking frequency of the external random force is much larger than the collision frequency. This is the limit in which the Fokker-Planck term in (1) models the random energy input through the random kicks. The external random force will, in principle, have a quantitative influence on the short range structure of the fluid. There is only one important difference with respect to the simulations carried out in [10]. There, the random rotation proposed in [36] was implemented to avoid inelastic collapse at high inelasticities (α < 0.5). This procedure amounts to rotating the relative velocity g by a small random angle after each collision. Consider the completely inelastic situation α = 0 for the sake of the argument. After each collision, the vector g lies exactly at the border of the precollisional hemisphere (g·σ = 0), so that if the aforementioned random angle has equi-probable positive and negative values, the rotation procedure will lead to a recollision with probability 1/2. This leads to a spurious increase of the number of collisions by a factor ∞ n=0 1/2 n = 2 (the recollision can itself induce a recollision with probability 1/2 etc. . . so that the frequency of collision effectively doubles !). When α is small but non vanishing, this effect is still present but weaker. This is clearly an artificial violation of molecular chaos that has been discarded in the present work: for α < 0.5, we have also implemented the rotation procedure, but if a small rotation leads to a recollision, a new angle is drawn until the pair separates. In this way, we reduce an important source of correlations (the effect is dramatic on all the low order moments B nm , not only on the collision frequency; in particular, the moments with n ≤ 1 that correspond to collisional averages of negative powers of g, are strongly biased toward bigger values if the "rotation-induced recollisions" are present). After applying this new rule, we are then left only with correlations induced by the hard sphere dynamics plus the ones induced by the noise itself (see below).

A. Cumulants
First we focus on the single particle velocity distribution function f averaged over all particles, which deviates from a Maxwellian distribution due to the inelasticity of the collisions. In the previous section we presented predictions for these deviations, assuming molecular chaos. The resulting expression given by (4) for the fourth cumulant a 2 of the distribution as well as the prediction for its overpopulated tail are in perfect agreement with 3D Direct Simulation Monte Carlo (DSMC) results over the whole region of inelasticities [34]. As DSMC itself invokes molecular chaos, this observation merely justifies the approximations made in the analytic calculation. Information about the validity of molecular chaos can only be obtained from a comparison with molecular dynamics simulations, and Fig. 2 shows this comparison for the fourth cumulant a 2 as a function of the coefficient of restitution α in 2D. The simulation results are in agreement with (4) for small inelasticity, but start to deviate significantly from the theoretical prediction at α = 0.6. These deviations, together with the perfect agreement between the theoretical prediction and DSMC results, provide direct evidence for the breakdown of molecular chaos for α < ∼ 0.6. The theoretical result (4) is independent of the density. As a 2 represents only a small correction in (2), one needs a large number of collisions and a large number of particles to reach sufficient statistical accuracy. So, high densities, for which one can use linked lists [30], are well suited. The data in Fig. 2 are typically obtained at high packing fractions (φ = 0.63 at α = 0.92 and 0.7, and φ = 0.55 at α = 0.6) and N = 10201 particles. At low densities and weaker inelasticities (α > ∼ 0.6) we are unable to collect enough statistics to measure the small correction, represented by a 2 . Simulations at higher inelasticities did not show any density dependence of a 2 in the range 0.2 < ∼ φ < ∼ 0.6, suggesting that the cumulant expression (3), obtained from the Enskog-Boltzmann equation also applies at liquid densities. Similar results as those displayed in Fig. 2 have been observed for the 3-dimensional version of the present model [11].

B. Radial distribution function
Next, we present results for the static or unconstrained radial distribution function g(r), and in particular its extrapolated values at contact Y , Y (−) , Y (+) . Here g(r) is essentially the density-density correlation function, whose long range behavior was studied in Ref. [10]. Figure 3 shows the measured values of g(r) for short distances and packing fraction φ = 0.2, at different inelasticities. At small inelasticities (α ≃ 0.9), g(r) resembles the radial distribution function for elastic hard disks (EHD). At higher inelasticities, deviations start to appear: the first and second maximum in the measured g(r) are enhanced with respect to their EHD values at the same density. Moreover, the functional shape also deviates from the corresponding pair distribution of EHD at an appropriately chosen higher density; e.g. if this density is chosen such that the value of the second maximum of the pair distribution of EHD coincides with the simulation result for IHD, the observed value at contact would still be underestimated by the EHD pair distribution.
It seems worthwhile to compare these results with existing experiments on granular fluids in which the pair distribution g(r) has been measured. In the experiment of Ref. [2] on a vertically vibrated thin granular layer, g(r) has been measured at φ = 0.46. In the fluidized ('gas-like') phase, it follows the equilibrium result for elastic hard disks almost identically. This result may be compared to our simulations for a randomly driven fluid of inelastic disks at α = 0.9, corresponding to the value for stainless steel balls used in the experiment. It would be of interest to measure experimentally how g(r) in the fluidized phase depends on the inelasticity, and see if a behavior similar to that of Fig. 3 is observed. It is also interesting to note that the pair correlation function g(r) in a non-Brownian suspension of spherical particles, fluidized between two vertical parallel plates, shows an enhanced value at contact as well [35].
In Fig. 4a we show the value at contact, Y , obtained by extrapolation from g(r) at φ = 0.05, together with the extrapolated values for approaching and receding pairs, Y (−) and Y (+) respectively. For α > ∼ 0.8 no significant deviations are found from the Verlet-Levesque value χ E = 1.084 for elastic hard disks at the same density. More surprising is the value of g(r = σ) at large inelasticities, reaching a value around 40 for α → 0. This property, combined with the observation that the first and second maximum in g(r) are shifted to smaller r-values, and are larger (up to 20% at small α) than the corresponding hard disks values, may be interpreted as a tendency to cluster, i.e. to stay in continuously rearranging configurations with large density inhomogeneities. We return to this point in subsection III G. Figure 4 also shows the dynamic correlation, χ (−) = χ E B 00 , measured as a collisional average. Fig. 4b compares the static ratio Y (−) /Y (+) with the dynamic one, χ (−) /χ (+) = α, and also shows the ratio Y (−) /χ (−) . The plots clearly show that the dynamic and static correlations, χ (±) and Y (±) , are different. For α < ∼ 0.5 the differences are large, and for α > ∼ 0.6 both functions are about equal. For the case of a freely evolving IHS fluid, Soto and Mareschal [18] have recently observed a similar behavior, and explained it in terms of the effect of the increase of grazing collisions on the effective χ (−) . In the randomly driven IHD fluid the same effects are present. All correlations, χ (−) , Y (±) are large, especially at small α. This is caused by the divergence of f (2) (c 1 , c 1 , σ) at small g and small cos θ, which corresponds to grazing collisions and will be further discussed in the next subsection. As a result of noise-induced recollisions, collisions with small g and small cos θ are oversampled; consequently, a dynamical average involving negative powers of g cos θ such as χ (−) is expected to be larger than its static counterpart Y (−) . This feature can be observed in Fig. 4b. Moreover, in the absence of recollisions, we would expect Y (−) = αY (+) as a result of plain hard sphere dynamics. However, in the heated system, the flux continuity, as expressed in Eq. (38), is no longer satisfied for the extrapolated Y 's; some pairs are put in pre-collision configuration under the action of the random force, which leads to Y (−) > αY (+) . The breakdown of Y (−) = αY (+) signals the inelasticity beyond which noise-induced correlations become relevant. It is furthermore possible to consider situations where the recollisions dominate the dynamics, e.g. at small α by allowing rotation-induced recollisions. In this extreme limit, we expect Y (−) ≃ Y (+) , the pairs being put in pre-or post-collision configuration essentially at random. In the same limit, the population of colliding pairs with small g and small cos θ is enhanced, leading to a more pronounced discrepancy between dynamic and static averages (i.e. a much smaller ratio Y (−) /χ (−) than observed in Fig. 4b).

C. Equation of state. Molecular chaos breakdown
To what extent does the extrapolated static radial distribution function, Y = g(r → σ), describe the nontrivial dependence in the NESS of collision frequency in (5), collisional damping in (14) and pressure in (12) on the inelasticity? If molecular chaos holds, the latter quantities depend, according to Eqs. (6), (13), (9) and (15), on the precollisional pair function at contact, χ (−) , where the particles are aiming to collide. This function differs from the extrapolated static Y (−) at high inelasticities (see Fig. 4). Consider first the collision frequency in the molecular chaos approximation, ω mc = χω 0 (T ) above (9), with χ = χ (−) the dynamic correlation, i.e.
where we have used Eqs. (9) and (35). This is an extremely poor approximation, as can be seen from Fig. 1, which shows that the measured value ω/ω E approaches 5.6 as α → 0, whereas B 00 is essentially divergent. Next we replace χ (−) in (43) by its static counterpart Y (−) , shown in Fig. 4. This yields Its limiting value for α → 0 is about a factor 3 too large when compared to ω/ω E . We conclude that all mean field approximations for the collision frequency, including the Enskog approximation ω E (T )/ω E = T /T E , break down for α < 0.6. Fig. 5 shows the pressure of the IHD fluid, compared with the molecular chaos prediction given by (13), taking for χ either the Enskog approximation χ E in (7), or the simulated Y (−) , or χ (−) . The Enskog approximation, accounting for the short range geometric exclusion effects in the precollision state, gives a reasonable description of p(T ) for all α, while both the static Y (−) and the dynamic χ (−) give an extremely poor description except for α > 0. 8. Consistent with this conclusion is the good estimate for the temperature T E in the NESS, obtained by balancing the energy dissipation rate Γ E (T E ) in (16) with the energy input from the random force, as shown in Fig. 1. Moreover the collisional energy loss, Γ(T )/Γ E (T ) = (T E /T ) 3/2 in (18), is in agreement with MD simulations over the whole α-interval within 30%. All other mean field approximations with ω E replaced by ω mc (T ) or ω stat (T ) give very poor results for Γ(T ) = mξ 2 0 . How can these paradoxical results be reconciled? Let us compare the individual definitions of χ (−) , ω, p and Γ, which all contain factors |g · σ| n f (2) (c 1 , c 2 , σ) with n = 0, 1, 2, 3. To find a possible explanation of these paradoxical results, we test the following scenario: the molecular chaos assumption (6) only breaks down at very small relative velocities g, and more precisely, at very small g n = g ·σ = g cos θ, which is the component of g, parallel to the line of centers of the colliding particles (physical arguments for this scenario will be offered in subsection III F where we discuss the noise induced recollisions). On the basis of this scenario the singularity in f (2) at small g makes the dynamic correlation χ (−) = B 00 χ E ( shown in Fig. 6) very much larger than χ E , essentially divergent as α → 0. In calculating the collisional frequency from (5) the extra factor g n in the integrand makes the small g n −singularity integrable, giving a finite correction to the Enskog collision frequency, also for α → 0 (see Fig. 1). The contributions of the small g n −singularity in p(T ) and Γ(T ) are essentially suppressed by extra factors of |g ·σ| n . This possibility has been analyzed systematically by measuring the behavior of the moments B nm (T ), which are useful tools to investigate the breakdown of the molecular chaos assumption. We have made n and/or m small in order to analyze the nature of the singularities in f (2) near small relative velocities and near grazing collisions, as displayed in Fig. 7. All deviations of these quantities from unity give a quantitative measure for the violation of the molecular chaos assumption. In the elastic limit we have carefully checked for a large number of cases that the reduced moments B nm tend to 1. Fig. 6 shows the values of different moments B nn (T ), and one can clearly see how the deviations from the elastic limit rapidly decrease as n increases to n = 3, after which they start to increase slowly. For larger n-values, the moments are reasonably close to unity, but statistical inaccuracy precludes any definite conclusion about the large n behaviour.
Further evidence for the above scenario is shown in Fig. 7, where we display two sequences of moments B nm . To draw some further conclusions from Figs. 6 and 7 we note that the integrands in B nn , B 0n and B n0 , as defined in (28), contain apart from f (2) respectively the factors g n | cos θ| n , | cos θ| n , g n . The reduced moments B 01 and B 10 contain again very large contributions from the divergence of f (2) near vanishing g n = g cos θ. Fig. 7 suggests that the presence of equal powers of g and cos θ in B nn simultaneously suppresses the large contributions from the singularities at g = 0 and cos θ = 0.
We conclude that the numerical results, displayed in Figs. 6 and 7, give support to the previous scenario, showing that molecular chaos breaks down only in a very small portion of the phase space, around g n = g ·σ = g cos θ = 0. The size of this 'pocket' in phase space increases as α decreases. Therefore, only those collisional quantities that contain low powers of g and cos θ (such as χ (−) and ω) will be very sensitive to this breakdown as the inelasticity increases, while physical quantities involving higher powers of g and cos θ, such as the temperature, pressure or energy dissipation will be well approximated by their molecular chaos counterparts.

D. Velocity correlations at contact
In the previous subsection we have considered the pair distribution function f (2) (c 1 , c 2 , σ) in the precollision state, and have examined how molecular chaos is broken down, and which physical quantities are most sensitive to it. Now we will analyze the effect of the breakdown of molecular chaos on collisional statistics.
We show in Fig. 8 different velocity collisional averages at φ = 0.05. In the simulations, these quantities are obtained by averaging over successive collision events in the steady state. We first observe that the simulation results in Fig. 8 approach for α → 1 the analytic results for elastic spheres, calculated in (24). At small inelasticities, the simulation data follow the trends of the theoretical prediction with systematic deviations depending on the quantity considered. For instance, the behavior of the center of mass motion G 2 coll is close to the analytical prediction of (24) in the whole range of α values. This indicates that the center of mass velocity G is not correlated with the relative velocity g. Consequently, f (2) (c 1 , c 2 , σ) in the collisional average (11) factorizes, and we can expect the contributions in numerator and denominator in (11) coming from G−integrations to cancel. Consistent with this behavior, we observe that the two curves in Fig. 8 The reduced moments have been measured independently (see Figs. 6 and 7), and used to calculate the expressions (44) and (45). The results have been plotted in Fig. 8 as dashed and dashed-dotted lines, which agree very well with the direct measurements of these quantities as collisional averages, shown in Fig. 8 respectively as squares and circles. In deriving (44) and (45) we have again used that the velocity variables G and g are statistically uncorrelated. The present results strongly support this assumption. The correlation ĉ 1 ·ĉ 2 coll = cos ψ 12 coll , also plotted in Fig. 8, cannot be expressed in B-moments. However, the approximate relation already employed to show that cos ψ 12 coll ≃ c 1 · c 2 coll / c 2 coll in section II C, holds for the simulation data over the whole range of inelasticities. As the system becomes more inelastic, the typical "temperature" of colliding particles (defined as the collisional average c 2 coll ) decreases and even becomes lower that the unconstrained average c 2 that defines the temperature. On the other hand, as already noted below (24), c 2 coll = 5/4 > 1 in the elastic limit. This decrease of c 2 coll is directly related to the increase of the small g-portion of phase space where molecular chaos is violated. At small α, most of the collisions occur between particles with small and even vanishing relative velocities. An extreme example is the inelastic collapse, mentioned in the introduction.
The correlation function c 1 · c 2 /g coll for the freely evolving IHD fluid has been simulated by Soto and Mareschal, and was shown to be small, but non-vanishing [18].
We have also investigated r−v−correlations by measuring the expectation value of c 1 ·c 2 (r) for two particles separated by a distance r, as defined in (42). The results are shown in Fig. 9 and Fig. 10. The plot shows an intermediate range of r values with an exponentially decaying correlation. It is again of interest to compare the extrapolation of the static correlation c 1 · c 2 (r → σ) with its dynamic counter part, c 1 · c 2 dyn , calculated at collision. The results, derived in (A.8) and (A.9) of the appendix, read for hard disks, The first term on the RHS represents the precollision part, Figure 10 compares the extrapolation c 1 · c 2 (r → σ) (circles) of the static correlation with its dynamic analogs (46) and (47). The numerical data for both correlations agree well for α > ∼ 0.8, but for α < ∼ 0.5 the dynamic correlation (solid line) is substantially larger than the static one. This is consistent with the difference between χ (−) and Y (−) observed in Fig. 4. For comparison the dynamic precollision correlation (dashed line) is also shown . It should be noted that the divergence of f (2) at small g and small cos θ implies in particular that B 00 ≫ B 20 > B 22 , so that (46) predicts that the dynamic correlation at contact c 1 · c 2 dyn should increase at α → 0 and saturate close to 1/2. By the same arguments its precollision part in (47) approaches the same limit. This can be observed in Fig. 10.
The velocity correlation c 1 · c 2 coll in (44) involve the reduced moments b 31 and b 11 . Consistent with the scenario, developed in subsection III C, the divergence of f (2) (c 1 , c 2 , σ) near g = 0 and cos θ = 0 is largely suppressed in these higher moments, which remain finite for α → 0, where b 11 ≃ 4b 31 . Consequently c 1 · c 2 coll does not approach the value 1/2 as α → 0, but a value close to 0.3, as can be deduced from Fig. 8.

E. Grazing collisions
The data in Fig. 8 for c 1 · c 2 coll , cos ψ 12 coll , √ 1 − b 2 coll and b coll clearly illustrate that the violation of molecular chaos strongly increases with increasing inelasticity. Consider first the average, This average remains at a plateau value 1/2 for α > ∼ 0.5 , which is determined by the uniform distribution P (b) corresponding to molecular chaos in two dimensions. Recall that the value 1/2 holds regardless of the functional form of the velocity distribution function f . It is thus a good probe for molecular chaos breakdown. Moreover, from its trend we can also estimate the way in which such a breakdown takes place. Specifically, as the inelasticity increases the average value increases by about 50%, which indicates a strong bias toward grazing collisions. To illustrate this, we model the normalized distribution of impact parameters as a uniform background and a 'half' delta-peak at b = 1, i.e. P (b) = 1 − p + 2pδ (1 − b), where p is the fraction of grazing collisions. This yields the average b coll = 1 2 (1 + p), which implies, according to Fig. 8, that at α = 0.0, 0.1 and 0.3 respectively a fraction of 50%, 35% and 5% is grazing at φ = 0.05. This qualitative picture is supported in a more quantitative manner in Fig. 11, which shows the measured P (b), which is strongly peaked near grazing collisions (b = 1). At small inelasticity all impact parameters are equally probable as expected on the basis of molecular chaos, and consistent with Fig. 8. Only for α < ∼ 0.5 deviations become significant: upon decreasing the coefficient of restitution, collisions with a larger impact parameter occur more frequently, implying an increase of the frequency of grazing collisions. The behavior of P (b) is then fully consistent with the divergence of f (2) at small cos θ, discussed in subsection III C.
To avoid inelastic collapse for α ≤ 0.5 the postcollision velocities of colliding pairs are rotated over a small random angle as described in Refs. [36,10], with the important restriction mentioned at the beginning of section III. Alternative algorithms to avoid inelastic collapse are described in Ref. [37]. For α > 0.5 no such rotation was applied. To check if the deviations of the impact parameter for α < ∼ 0.5 are due to this applied rotation, we have also performed simulations where even for α > 0.5 a random rotation was applied. Regardless of the applied random rotation, we found b coll close to 1 2 for α > ∼ 0.5. Both Figs. 8 and 11 show that for α < ∼ 0.5 molecular chaos is strongly violated, and that the violation is weaker in the small inelasticity regime. The average √ 1 − b 2 coll supports the same conclusions. The data for c 1 · c 2 coll and cos ψ 12 coll in Fig.8 are consistent with the predominance of grazing collisions at large inelasticities. They show the average relative angle between the velocities of the incoming particles, which has a strong α−dependence and no plateau value near the elastic limit. Near α = 1 the particles are on average on approaching trajectories with cos ψ 12 coll ≃ −0.25 and ψ 12 coll ≃ 105 o . As α decreases, cos ψ 12 coll increases linearly to a value 0.50, while ψ 12 coll approaches 60 o , at α = 0. This corresponds to collisions of more or less parallel-moving pairs of particles, where faster particles overtake slower ones.
Figs. 12 and 13 show the distribution of relative orientations of incoming velocities. The distribution of angles between the incoming particles (ψ 12 ) shows moderate deviations from what is expected for an elastic system in the range 0.5 < ∼ α < 1. As an analytic expression for elastic disks is not available, deviations are compared with the simulation results for elastic hard disks (in the absence of a random external force). At α = 0.5 the frequency of collisions of parallel-moving particles is strongly increased, a trend which is enhanced upon increasing the inelasticity. Finally, the probability distribution P (c 1 · c 2 ) is shown in Fig. 13. When the inelasticity increases, this distribution becomes more peaked around the origin, as the colliding particles on average move more slowly relative to each other. In the mean time, the typical angle ψ 12 decreases, which causes this peak to shift to positive values.
There is strong evidence that the effects of ring collisions are considerably enhanced in fluids with dissipative interactions, such as granular flows, where relative kinetic energy is lost in binary collisions. As a result the postcollision velocities {v * 1 , v * 2 } will be on average more parallel than the precollision ones {v 1 , v 2 } [24], i.e. the trajectories are less diverging than in the elastic case, and there is a much larger {r 3 , v 3 } phase space, in which particle 3 will knock, say, particle 1 back to recollide with particle 2.
This increase of phase space is confirmed by gathering recollision statistics. We have counted the fraction of recollisions as a function of α, as shown in Table I. The column labeled R 1 (recollisions between two partners mediated by a third particle) shows that at a packing fraction φ = 0.2 in the elastic case (α = 1) only a fraction of 6.7% of all collisions is a recollision. This frequency gradually increases to about 15% at α = 0.4.
In the randomly driven IHS fluid there is the additional effect of noise induced recollisions that do not require the intervention of other particles. This type of recollision (denoted R 0 ) occurs with high probability when the relative velocity after collision is so small that it can be simply reversed by a random kick. At α = 0.6 the frequency of noise-induced recollisions is about 4%, and it increases to 52% at α = 0 (see column R 0 in table I). The effect is of importance at all densities, because it does not require the mediation of a third particle. Indeed, at a low packing fraction of 1% and in the completely inelastic case α = 0, the frequency of R 0 -like events is still 34%, while R 1 -like events have dropped to 5%. Moreover, we have verified that inclusion of rotation-induced recollisions modifies most of the collisional quantities we have analyzed, increasing their deviations with respect to the molecular chaos prediction.
At present, more quantitative theories or estimates of the effect of both types of recollisions and other ring collisions on the short range behavior of the pair distribution function f (2) (x 1 , x 2 ) are lacking. A natural way to incorporate the noise-induced recollisions into a kinetic theory description would be to include them into an effective two-particle scattering operator, which transforms an asymptotic precollision state of two independent particles into an asymptotic postcollision state, without involving intermediate two-particle scattering states, as in the present case. This may lead to an instantaneous Boltzmann collision term (without memory effects), provided the mean free time and the time between random kicks are very well separated (dilute gases). Such a description would suppress the recollions of type R 0 , and make the violation of molecular chaos less severe, say comparable to the freely evolving IHD fluid.

G. Cold dense inhomogeneities
In Ref. [10] we have shown by analyzing the Fourier modes of the granular hydrodynamic equations, which are valid for small inelasticities (say α > 0.7), that the NESS in a randomly driven IHS fluid is linearly stable against spatial inhomogeneities. Consequently, when observed over sufficiently long times, the NESS should be spatially homogeneous. However, it was also shown that the NESS exhibits strong fluctuations, resulting in long range spatial correlations in density, flow field and granular temperature. The observation of density inhomogeneities for large inelasticities has already been reported by Peng and Ohta [8]. These density inhomogeneities, as shown by the snapshot of the density in Fig. 14, are not quasi-static, as in the freely evolving case [38,36,39,24], but seem to behave as dynamic assemblies of particles that dissolve and re-assemble again. Also for a uniform shear flow, dynamical density inhomogeneities have been reported [40]. The existence of density inhomogeneities was already suggested by the static pair distribution functions g(r), which showed an enhancement of the first few maxima as compared to their elastic values (see Fig. 3).
In Fig. 8 we show that the mean energy c 2 1 coll , of particles aiming to collide, is above the mean, c 2 = 1, for small inelasticity. It decreases from its elastic value 5 c 2 /4 with decreasing α, then crosses the mean value value c 2 = 1 at α ≃ 0.2, and further decreases to approximately 0.7 c 2 at α = 0.
It is interesting to observe that in the strong dissipation range, the mean kinetic energy or granular temperature of particles that are about to collide is lower than the average temperature. We combine this observation with Figs. 3a,b of Peng and Ohta [8], which show that essentially all collisions occur inside "cold" regions of high densities. This last observation applies even more so to undriven IHS fluids [38,41]. We expect that, also in the randomly driven IHS fluid, the majority of collisions takes place inside cold high density regions.
If the predominance of cold particles in strongly inelastic collisions, c 2 coll < c 2 is indeed a signal for the appearance of density inhomogeneities, then Fig. 8 suggests that at a packing fraction φ = 0.05 density inhomogeneities may occur for α < ∼ 0.2. This is indeed confirmed by the snapshots in Figs. 14. In Fig. 15 we illustrate the existence of cold inhomogeneous dense regions for α = 0.2 and φ = 0.2. The particles with a less (more) than median kinetic energy are shown on the left (right). The formation of inhomogeneities is more clear for the colder particles. The temporal evolution of these regions show that they dissolve after some time, while new inhomogeneous regions appear. The formation of "living" inhomogeneous regions can be understood using the hydrodynamic picture put forward in [10], where it was shown that the structure factor behaves as S(k) ∼ k −2 , implying density correlations increasing with distance as ln(r) in two dimensions. These long range spatial correlations induce a slowing down of the dynamics, as in critical phenomena. This, in turn, implies the slow decay of density perturbations, that could lead to visible density inhomogeneities as the kicking frequency is reduced (in this respect, see Refs. [6]). We can also expect that upon decreasing the forcing frequency, the dynamics should be closer to its "free cooling" counterpart so that well defined clusters are then likely to appear.
More details about the predominance of cold particles, among those involved in collisions, can be seen in Fig.  16, which shows the constrained probability distribution P (c), defined in subsection II C and obtained from MD-simulations at different inelasticities. For α < ∼ 0.5 the distribution has significantly shifted to smaller impact velocities. For the completely inelastic case, collision events involving "immobile" particles are more than twice as frequent as for the elastic case. The second moment of the distribution displayed in Fig. 16 decreases when increasing the inelasticity. In fact, all functional forms with simulation data at different α can essentially be collapsed onto a single universal curve (the elastic one) by plotting T (α)P (c|α) as a function of c/ T (α), where T (α) = c 2 coll is the mean temperature of a particle at collision. The collapse plot is shown in Fig. 16b. This data collapse confirms the concept of cold dense regions dominating the energy dissipation. This could point to a possibly relevant two fluid picture of a "hot" dilute background gas coexisting with continuously rearranging configurations of "cold" dense regions.

IV. CONCLUSION
We have performed extensive MD simulations to study the kinetic properties and short range correlations in the non-equilibrium steady state of a randomly driven fluid of inelastic hard disks, as a model for fluidized granular material. The MD results have been compared with kinetic theory predictions derived from the Enskog-Boltzmann equation, properly modified with a Fokker-Planck diffusion term ξ 2 0 (∂/∂v) 2 to account for the applied random driving force [9].
It appears that the kinetic theory predictions, based on molecular chaos, are essentially in agreement with the MD results for small inelasticities (α > ∼ 0.5) at φ = 0.05. For larger inelasticities the deviations from the molecular chaos predictions start to become manifest: the radial distribution function at contact differs strongly from its local equilibrium form; there is a predominance of grazing collisions. When increasing φ, the effects of the inelastic collisions become relevant at smaller inelasticities; e.g. at φ = 0.2 and φ = 0.5, we observe already significant deviations for α ≤ 0.7.
To avoid inelastic collapse of the system at low α, we have implemented a modified rotation procedure (see the beginning of section III). In its original version, this procedure induces dramatic violations of molecular chaos. It could then be argued that the important deviations of low order reduced moments B nm are also spurious consequences of the above rotation procedure. However, we checked that circumventing the collapse by applying elastic collisions when the relative velocity of a pair is below a certain cutoff [37], also induces very important violations of molecular chaos (quantified by B 00 for instance), unless the cutoff is chosen unphysically high.
Sequences of ring collision processes, which lead to the breakdown of molecular chaos in classical fluids with conservative interactions, are strongly enhanced in fluids with dissipative interactions, like rapid granular flows. We have analyzed how molecular chaos is broken, i.e. essentially only through pairs of colliding particles at very small relative velocities. This means that molecular chaos is violated only in a small portion of phase space, implying that only certain physical properties will be sensitive to this violation. This explains why quantities like the collision frequency, or the pair distribution function at contact are very sensitive to the inelasticity parameters, while others like the pressure or the energy dissipation rate are well approximated by their Enskog prediction. Disentangling the effects of hard disk and noise-induced correlations remains an interesting point to explore. The studies performed in a freely evolving IHS fluid also shows the predominance of grazing collisions at long times. The fact that we have observed an analogous behavior for this homogeneous steady state indicates that the mechanism of breakdown of molecular chaos in granular fluids through grazing collisions is generic for this type of fluids.
The extra feature of noise-induced recollisions, which do not require mediation of a third particle, will further enhance the violation of the molecular chaos assumption. A natural way to develop a kinetic theory for randomly driven fluids, thereby presumably restoring the validity of the molecular chaos assumption in the dilute gas case, could be to include the noise-induced recollisions in an effective two-particle scattering operator. It would be of interest to study its properties, either analytically or by simulating a two-particle inelastic collision in the presence of external noise. An additional theoretical complication here is the validity of the Boltzmann equation (1) with Fokker-Planck diffusion term due to the fact that there are two limits involved when dealing with hard spheres in combination with external white noise. The actual properties of the effective collision operator depend on the order in which both limits are taken. In the simulations, one always takes the hard sphere limit first, while the white noise is approximated by discrete kicks which are applied to the particles at discrete times.
In Ref. [10] we have calculated the equal-time spatial correlations of the fluctuations in the hydrodynamic densities in the NESS. Here we have focused on the dynamic properties of these enhanced fluctuations, in particular of the dynamic inhomogeneities observed in the density field. The collisional velocity moments, introduced in section II and measured in MD simulations, reveal that the dense regions consist mostly of particles colder than average. This is clearly shown in the velocity distribution P (v|α) of particles which are about to collide.
The MD simulations have been carried out in the limit in which the time interval between the external random kicks is much shorter than the mean free time between collisions. In this limit, regions with density larger than average are not seen to survive for a long time. Rather, they form, dissolve and reappear elsewhere. The spatial correlations analyzed in [10] show long-range correlations, which imply also a slowdown in the temporal decay of density perturbations. Therefore, we expect than the decrease of the kicking frequency will be accompanied by the appearance of apparent clusters. This fact, together with the shape modification of the velocity distribution P (v|α) (see Fig. 16) suggests the picture of a two fluid model, in which a "hotter" more dilute background gas coexists with continuously rearranging configurations of "cold" dense clusters. This point remains open for subsequent investigation; for example, it would be interesting to analyze separately the collisional statistics in the dense and dilute regions to assess the role of density fluctuations. under the integral sign in Eq.(41) by its value at contact, f (2) (c 1 , c 2 , σ). We proceed in the same fashion as in Eqs.(36)- (39), and split the numerator in (42) in a pre-and postcollision part, as done in subsection II D. One finds after a lengthy calculation, Here the first term on the RHS is its precollision part, i.e. In section III these quantities are compared with MD simulations.