Collision statistics of driven granular materials

Collision statistics of driven granular materials Daniel L. Blair and A. Kudrolli arXiv:cond-mat/0212336v1 [cond-mat.soft] 13 Dec 2002 Department of Physics, Clark University, Worcester, MA 01610, USA (Dated: February 1, 2008) We present an experimental investigation of the statistical properties of spherical granular particles on an inclined plane that are excited by an oscillating side-wall. The data is obtained by high-speed imaging and particle tracking techniques. We identify all particles in the system and link their positions to form trajectories over long times. Thus, we identify particle collisions to measure the effective coefficient of restitution and find a broad distribution of values for the same impact angles. We find that the energy inelasticity can take on values greater than one, which implies that the rotational degrees play an important role in energy transfer. We also measure the distance and the time between collision events in order to directly determine the distribution of path lengths and the free times. These distributions are shown to deviate from expected theoretical forms for elastic spheres, demonstrating the inherent clustering in this system. We describe the data with a two-parameter fitting function and use it to calculated the mean free path and collision time. We find that the ratio of these values is consistent with the average velocity. The velocity distribution are observed to be strongly non-Gaussian and do not demonstrate any apparent universal behavior. We report the scaling of the second moment, which corresponds to the granular temperature, and higher order moments as a function of distance from the driving wall. Additionally, we measure long time correlation functions in both space and in the velocities to probe diffusion in a dissipative gas. PACS numbers: 81.05.Rm, 05.20.Dd, 45.05.+x 45.70.Mg, 45.70.-n, 51.10.+y I. INTRODUCTION Granular material represent a type of matter not well defined by conventional means. Although each granular particle is obviously solid, an assemblage of these particles show distinctly non-solid behavior when subjected to external forces [1]. In the rapid flow regime, the interaction between the grains is collisional and the system resembles a dense granular gas. Indeed, the kinetic theory for dense gases formulated by Chapman and Enskog [2] have been modified to include the dissipative nature of the collisions [3, 4]. However, a number of approximations have to be made in any calculation that can be only validated by experiments. Furthermore, even if key assumptions such as equipartition breakdown [5, 6, 7], it is important to have a measure of the failure to guide further development. Energy has to be constantly supplied from an external source to observe a steady state in granular gas systems. Therefore, model experiments consist of granular particles inside a container where energy is continuously injected at a side-wall [8, 9, 10]. Thus gradients are present in experimental granular systems, which implies that care must be taken when comparing results to non-equilibrium kinetic theory [11, 12, 13]. With advances in high speed image acquisition, it is now possible to obtain positions of particles several times between collisions. However, particle positions and velocities can be obtained accurately only in two dimensions by direct imaging thus forcing certain constraints on the geometry of the system. One of the first such experiments to investigate velocity distribution functions (VDFs) utilized an apparatus in which particles are vibrated vertically inside a narrow transparent box [8, 14, 15]. Maxwellian statistics were reported for the vertical and horizontal velocity components of the particles parallel to the plane of the transparent side-walls. Additional interactions in these systems arise due to collisions between particles and the side-walls [14]. Following this work, Wildman et al. [16] were able to do long time particle tracking to measure diffusion constants by interpreting mean square displacement data over a very broad range of density. More recently, in a similar apparatus, Rouyer and Menon [17] report that their VDFs have a universal form that can be parameterized by a single variable, the granular temperature. A different method of energy injection utilizes large flat container that is vibrated vertically to excite a sub mono-layer of particles [10, 18, 19]. The velocity of the particles in the horizontal plane are measured and are found to follow a non-Gaussian distribution. However, the impact of the velocity gradient in the vertical direction on the observed distributions are not taken in to account because these components cannot be measured. Our experiment is a variation of the vertically vibrated apparatus. Spherical particles are constrained to roll on an inclined two dimensional surface. This geometry allows for a direct investigation of the interplay between energy injected at the side-wall and the dissipation through inelastic collisions. In addition, the inclination reduces the effects of gravity, therefore minimizing shock waves. This system has been used to demonstrate clustering and collapse when the inter-particle collision frequency is much greater than particle-driving wall collision frequency [9]. Recent works have explored a full range VDFs, from very near Gaussian behavior to highly non-Gaussian distribution functions, as well as velocity correlations [20, 21]. 2 In addition to analytical techniques and experiments, several groups have utilized computer simulations of inelastic hard spheres with both Molecular Dynamics [12, 22, 23, 24, 25, 26, 27] and Direct Simulation MonteCarlo [28, 29, 30, 31, 32] techniques to investigate the statistical properties of granular gases. Using DSMC simulations, Baldassarri et al. [28] have found velocity and density distributions that are qualitatively similar to our previous experimental results [21]. Recent work by Brey and Ruiz-Montero [27] investigate how the second and fourth moments of the VDFs scale as a function of distance from the driving wall, which until now, have not been experimentally tested. In this paper, we report on the statistical properties of a gas of inelastic particles constrained to two-dimensions. An inclined geometry reduces the gravitational acceleration acting on each particle which results in lower mean velocities. The combination of slow dynamics and high speed imaging allows us to accurately identify the particle trajectories and collision events. By using velocities before and after a collision event, we measure the normal coefficient of restitution. We find that these quantities are found to be broadly distributed for the same impact parameters. By calculating the distance and time between collision events we measure the distributions of free paths and times. We find that these distributions do not follow the result found from kinetic theory. The path and time distributions have an overpopulation of short distance and time bins, demonstrating the inherent clustering present in granular gases. We propose an empirical form that captures the distributions, which is then used to calculate the mean free path and free time as a function of density. The particle trajectories are also used to measure the mean square displacement, velocity auto-correlation, and diffusion rates. The distribution of particle velocities are measured with a variation in density of an order of magnitude and show distinctly nonGaussian behavior with no apparent universal form. We compare our results to recent experiments, as well as theoretical and simulation treatments of equivalent systems. The paper has the following structure. In Section II we present the experimental apparatus and imaging methods. Section III provides the overall system characteristics such as the density distributions and coefficients of restitution and inelasticities. We then present our analysis of the trajectories of the particles in Sec. IV. Finally, in Sec. V we summarize our results in the context of granular kinetic theory and simulations. II. EXPERIMENTAL METHODS The experimental configuration [Fig. 1], consists of a 100 σ x 60 σ glass plane that is inclined at an angle β with the horizontal. The particles are stainless steel with diameter σ = 3.175 mm and a high degree of sphericity (δσ/σ = 10−4 ). The number of particles, measured in number of mono-layers Nl across the driving wall, is var- Camera (a) Soleniod β x y Fcn. Gen. Amplifier Computer (b) 19 cm Driving Wall 31 cm FIG. 1: (a) Schematic diagram of the experimental setup. The inclined plane is a smooth glass surface, the side-walls and driving wall are stainless steel so that the particle-boundary collisions approximate those between particles. The driving is produced by a solenoid connected to the lowest side-wall. The angle of inclination β, can be varied from β = 0 − 8 degrees, the values of β we have chosen are 2◦ ± 0.1◦ and 4◦ ± 0.1◦ . (b) An image of the system taken from above. The bottom right corner is considered the origin of our coordinate system (0, 0). The white bars allow us to track the position of the driving wall. ied between Nl = 1–5 in steps of one layer, (viz. from Np = 100 − 500 in steps of 100, where Np is the number of particles). The energy source is an oscillating sidewall, driven by a solenoid, that is located as shown in Fig. 1(a). The driving signal is a 10 Hz pulse with a velocity during each pulse of ∼ 40 cm s−1 . The driving frequency and amplitude were chosen to ensure that no phase dependence on the center of mass is observed (at frequencies below 2 Hz the particle positions are phaselocked with the driving). The signal is produced with a computer interfaced Aglient Technologies 33120-A waveform generator that is and subsequently amplified by an 3 7 respectively [see Fig. 1]. A typical particle trajectory is shown in Fig. 2(a). Multiple collision events can be distinguished with nearly straight paths between each event. A particle that freely rolls on the inclined plane will follow a parabolic trajectory [see Fig. 2(b)]. The particle trajectory is given by y (cm) (a) 3.5 y(x) = 0 0 3.5 7 5 x2 g sin(β), 7 2vx2 (1) where, g is the acceleration due to gravity, vx is measured from the width of the parabola, and the 75 factor is due to the moment of inertia for a solid sphere. x (cm) (b) III. 11.5 SYSTEM CHARACTERISTICS y (cm) A. 9.5 7.5 particle position / 2 2 y(x) = g (x /v x) 8 9 10 11 x (cm) 12 13 FIG. 2: (a) Linked particle positions over 1365 time steps (Nl = 3). We can determine particle collision events with a high degree of accuracy from trajectories such as this. (b) The parabolic path of a particle. We use fixed values for vx and g = 980 cm s−2 to measure g ′ = 75 g sin β. The fit gives β = 2.2◦ , the deviation from the measured value of β is 9%. HP 6824A Amplifier. The inclination of the plane can be varied between β = 0◦ −8◦ , for our experiments the angle was fixed at β = 2◦ ± 0.1◦ or 4◦ ± 0.1◦ . In the extreme case of β ≪ 1◦ , the particles essentially cease to interact with the energy source and cluster at the side opposite of the driving. The particles are imaged using a Kodak MotionCorder SR1000 high speed digital camera. We measure the positions of all particles contained in the apparatus for 1365 frames at 250 frames per second at full spatial resolution of 512 × 480 pixels. These digital images are then transfered to a computer and analyzed using a centroid method that allows us to resolve each particle to subpixel accuracy. After each particle is located the positions are then connected in time to form continuous trajectories for 5.46 s. Our coordinate system is such that the x, y axes are parallel and perpendicular to the driving Density Distributions The results presented will be given in terms of the number of single layers across the cell, Nl and the angle of inclination β, which the determine the area fraction φ [see Table I]. We measure φ by defining a region of interest (ROI) that is centered about the peak in ρ(y), [Fig 3(b)] whose extent in the y-direction is limited to ±10% of ρ(y). The ROI scheme excludes all particles that are within 3σ of the side walls to ensure that clustering due to the side-walls does not affect our results [see Fig. 3(a)]. The over-plotted box in Fig. 1(b) demonstrates the ROI definition for φ = 0.13. A more stringent division of the system in the y-direction will be used when the behavior of the temperature, pressure and kurtosis are discussed in Sec IV. The form of the density in the ydirection is similar to that found in Refs. [14, 36] however, we find that the form of the tails of ρ(y) at higher values of Nl β deviate from Boltzmann distribution. This implies that the law of isothermal atmospheres breaks down for granular systems as we shall also see when we discuss the scaling of the granular temperature in Sec. IV E. Nl β φ Np 1 2 3 4 4 5 2.0 2.0 2.0 2.0 4.0 4.0 0.022 0.068 0.138 0.191 0.302 0.581 100 200 300 400 400 500 TABLE I: Experimental values of the number of layers Nl , the angle of inclination β, and the resulting measured value of φ. Np the number of particles in the system is given for clarity. 4 0.4 (a) 0.3 ρ (x) where v̂ = v/|v| the unit vector of the calculated velocity. If 20◦ ≤ ψ ≤ 180o the proximity of all particles at the same time instant is checked. If a particle is found within a radius σ + ∆σ, whose velocity also satisfy Eq. (2) it is considered as a candidate for a collision. To assure that re-collisions are not occurring, we maintain a record of the identity of the previous collision partner. We then ensure that those particles can re-collide if and only if the partner particle has undergone a collision with yet a third particle. If particles pass these requirements then a collision has occurred. To extend the algorithm to include collisions with the boundary walls we first check if Eq. (2) is satisfied. We then check if the particle’s center is within σ + ∆σ of a boundary and it’s velocity component perpendicular to the wall is reversed. Nl = 1, β = 2 Nl = 2, β = 2 Nl = 3, β = 2 Nl = 4, β = 2 Nl = 4, β = 4 Nl = 5, β = 4 0.2 0.1 0 0 10 20 30 x (cm) 0 10 C. Coefficient of Restitution ρ (y) (b) −1 10 −2 10 0 10 20 y (cm) α=− FIG. 3: (a) The density ρ(x) versus x for all Nl . The obvious clustering due to inelastic collisions at the side-walls is demonstrated here. Also, as Nl is increased the system becomes more inhomogeneous in across the cell. This effect is most likely due to the onset of clustering instabilities that have been recently discussed [33, 34, 35]. (b) The aerial density plots ρ(y) for each Nl and β on a log-linear plot. The density φ is measured in a particular area by integrating ρ(y) over that region of interest. The total are under each curve corresponds to the average area fraction for that particular Nl . The solid line shows an exponential fit over the tail of the distribution of Nl = 200. However, we will demonstrate that the isothermal atmosphere is not obeyed for any density. B. Particle Collisions We identify collision events from the trajectories by using the following algorithm. Velocities are constructed as finite differences vj = ∆x ∆t , where ∆x = x(tj ) − x(ti ) and the subscripts i, j represent positions separated by the time difference ∆t = 4 ms. All velocity vectors are compared sequentially to find direction changes given by ψ = cos −1 (v̂i · v̂j ), The loss of energy in a collision is given by the coefficient of restitution. Two particles that undergo an inelastic collision with a relative velocity between parti∗ · σ̂ = cles v12 = v1 − v2 , will obey the reflection law v12 −α v12 · σ̂, where α is the normal component of the restitution coefficient and σ̂ is the unit vector connecting the centers of the particles. Having an efficient method for collision identification, we are able to measure the relative velocities of two particles before and after collision events. The coefficient of normal restitution during a binary collision is given by (2) ∗ (v̄12 · σ̂) , (v̄12 · σ̂) (3) where the over-bar denotes average over three pre/postcollisional velocities measured in the ROI described above [see Fig. 1(b)]. The angle between the relative velocities of two colliding particles is given by ∗ θ = cos−1 (v̄12 · v̄12 ). (4) Thus we can characterize the coefficient of restitution as a function of θ. The probability distributions P (α) for 60◦ ≤ θ ≤ 180◦ for each Nl β, are shown in Fig. 4(a– f). Data for θ < 60◦ suffers from a lack of statistics and therefore is not included. Each graph represents the probability of the inelasticity having a value α for a range of θ + ∆θ, where ∆θ = 2◦ . P (α) follows a very broad distribution of values over all θ, and have a decreasing mean value as function of φ [see Fig. 5(a)]. Thus we find that the coefficient of restitution can have a broad distribution of values for the same impact angle. We also measured the energy loss due to a collision as a function of Nl β. The ratio of the magnitudes of the relative velocities before and after a collision, η= ∗ |v̄12 | , |v̄12 | (5) 5 1.0 (a) α 0.8 0.6 0.4 0.2 0.0 0.1 0.2 0.3 φ 0.4 0.5 0.6 0 10 (b) -1 P(η) 10 -2 10 -3 10 -4 10 -5 10 FIG. 4: The distribution of the normal component of restitution α versus 60◦ ≤ θ ≤ 180◦ , the relative angle of incidence between particle velocities. (a) Nl = 1 , β = 2.0, (b) Nl = 2 , β = 2.0, (c) Nl = 3 , β = 2.0, (d) Nl = 4 , β = 2.0, (e) Nl = 4 , β = 4.0, (f) Nl = 5 β = 4.0. The value of the z-axis for each graph is the probability of a collision giving a value of α in a range of θ + ∆θ, where ∆θ = 2◦ determines the energy restitution coefficient, (η 2 = α2 if all θ are averaged). Figure 5(b) shows the distributions of measured values of η shifted for clarity. We find that a peak exists at a value that is consistent with α2 . Furthermore there exists a power law tail for values of η > 1, which has been interpreted as a random inelasticity [37]. The appearance of a tail at high η implies that the rotational degrees of freedom are actively transferring energy to translational motion during a collision. IV. A. RESULTS Distributions of Paths and Times We measure the distribution of paths lengths from the the geometric distance between collision events defined in our ROI at each Nl β [Fig. 6(a–f)]. By basic kinetic theory arguments [38], the distribution of path lengths -1 10 0 10 η 1 10 FIG. 5: (a) The mean value of the distributions of α shown in Fig. 4 averaged over 60◦ ≥ θ ≥ 180◦ , as a function of the average covering fraction φ. The bars indicate the spread in the distribution. (b) The distribution of energy inelasticities given by Eq. (5) for (⊲)Nl = 1 , β = 2.0, (⊳)Nl = 2 , β = 2.0, (△)Nl = 3 , β = 2.0, (⋄)Nl = 4 , β = 2.0, ()Nl = 4 , β = 4.0, (◦)Nl = 5 β = 4.0. Each distribution is shifted vertically for clarity. for an elastic hard-sphere gas (and by a similar treatment the distribution of free times) is given by, √ √ P (l) = (2 2 φ) e −2 2 φ l . (6) The distribution therefore should follow a simple exponential form depending only on the density. However it is clear from the dashed lines in Fig. 6(a–f) that the simple form given by Eq. (6) does not describe the behavior over all l. The distributions of times between collisions P (τ ) [Fig. 7(a–f)] is also measured and shows similar behavior to that of the path length distributions, that is an overpopulation of the short time bins. This should be expected from the simple relationship between the displacement and the time. However, it is noteworthy to mention that the ratio of l/τ versus the path length l, is not a constant over all values of l, implying that the average speed of the system depends on the distance or time between collisions. 6 -5 10 -2 10 -2 10 10 P(l) 10 P(l) l (cm) 10 1 -5 10 -2 10 -3 10 l (cm) 1 10 -3 10 -3 10 -4 -4 10 l (cm) 10 0 20 0 10 10 0 -1 -3 10 l (cm) 1 10 P(l) -5 10 -2 10 0.6 τ (s) 1.2 0 0 10 (d) 10 0 -1 10 -2 10 10 0 15 -5 10 -2 10 -3 10 (c) l (cm) 10 0 10 -1 -4 -3 10 -5 l (cm) -5 10 0 9 6 0 10 (e) 10 10 P(l) -3 (f) 10 -5 10 -2 10 l (cm) 10 0 10 10 10 0 0.3 τ (s) -1 -3 10 (e) -1 10 -2 10 l (cm) 10 0 10 1 l (cm) 2 -3 10 10 -5 10 1 l (cm) 2 FIG. 6: The probability distributions of path lengths P (l) versus l, on a log-linear scale, and inset log-log scale (a) Nl = 1 , β = 2.0, (b) Nl = 2 , β = 2.0, (c) Nl = 3 , β = 2.0, (d) Nl = 4 , β = 2.0, (e) Nl = 4 , β = 4.0, (f) Nl = 5 , β = 4.0. The dashed line shows the theoretical form given by Eq. (6) derived for elastic particles, and the solid line is an empirical fit given by Eq. (7a). Table II shows the fit parameters. Elastic hard spheres will have a mean free path that is simply l̄ = v̄τ̄ , where v̄ and τ̄ are the average speed and collision time respectively. Also, the mean free path can be derived directly from the distribution of path lengths, R∞ l̄ = 0 l P (l) dl, where P (l) is given by Eq. 6. Grossman et al. [6] have interpolated how the mean free path for a granular system should be modified to account for higher collision rates due to increased density. Although the interpolation gives a qualitatively accurate correction for passing between the high and low density limits, the actual distribution of path lengths has not been measured or calculated for a granular gas. We have found an empirical form that well describes the measured distributions of path lengths and free times, P (l) = a (l)−b e−c l P (τ ) = a (τ )−b e−c τ -3 10 -4 -4 10 0 0.26 -2 10 10 10 -5 10 -5 10 0 (f) 10 -2 -5 10 -5 0 0 -4 -4 10 10 0.6 10 -2 10 0.13 τ (s) -5 0 4 0 -1 -2 10 l (cm) 3 P(τ) -1 2 0 0 P(l) 10 1 P(l) 3 10 P(τ) -5 0.4 -3 10 10 10 10 0 0.2 τ (s) 10 -4 10 -4 (d) -2 10 -4 10 1.2 10 -2 -2 10 0.6 τ (s) 0 10 P(τ) 10 (c) l (cm) 10 P(l) -1 5 0 P(l) 10 10 -4 -5 0 P(l) 10 10 -3 10 P(l) (b) -2 -2 10 -2 10 -1 10 (a) P(τ) P(l) 10 10 0 (b) -1 P(l) -1 -1 0 10 0 P(τ) 10 (a) P(τ) 0 10 (7a) (7b) where a, b, c for the path lengths and free times are shown in Table II for all Nl β. This form appears to capture both the short l and τ power-law behavior. In the dilute case the form returns to the theoretical prediction for larger path lengths. 0 FIG. 7: τ , on a 2.0, (c) 4.0, (f) (7b) 0.14 τ (s) 0.28 -6 10 0 The probability distributions of free times P (τ ) versus log-linear scale (a) Nl = 1 , β = 2.0, (b) Nl = 2 , β = Nl = 3 , β = 2.0, (d) Nl = 4 , β = 2.0, (e) Nl = 4 , β = Nl = 5 , β = 4.0. The solid line is a fit given by Eq. From the distribution of path lengths and free times, we calculate the mean free path and time by utilizing the fitting form and it’s parameters. The ratio of the mean free path to the mean collision time should determine the average speed v̄ in the ROI where the distributions are measured. We have taken the ratios of the integrated distributions, R∞ ¯l l P (l) dl 0 , v̄ = = R ∞ τ̄ 0 τ P (τ ) dτ (8) and compared that to the average of the speed distribution hvx,y i in the same ROI. Figure 8 the both measurements for all Nl β. The agreement is within 10% over the entire range of Nl β indicating that the proposed forms in Eqs. (7a,b) quantitatively capture the behavior of the distributions. B. Velocity Autocorrelation The velocity autocorrelation function (VAF) is computed for the x components of the velocities within an 7 1 Nl = 1, β = 2 Nl = 2, β = 2 Nl = 3, β = 2 Nl = 4, β = 2 Nl = 4, β = 4 Nl = 5, β = 4 (a) Cv (t) −1 <v>, v (cm s ) 24 12 0.5 0 0 0 0.2 0.4 0 0.6 0.5 1 φ 1.5 t (s) 80 60 Cx2 (t) FIG. 8: The average speed measured for each Nl β. (◦) v̄ obtained from Eq. (8), and (⋄) hvi measured from the mean of the speed distribution. The two independent measures give similar values over the entire density range. 40 (b) 20 ROI by using the following [39]: Np ,Ns tmax X X 1 Cv (t) = vij (to ) · vij (to + ∆t) (9) Np Ns tmax i,j=0 0 0 0.5 1 t (s) 1.5 2 2.5 ∆t=1 Cv (t) = hv(0)2 ie −t/τE , 24 (c) 2 −2 Dx (cm s ) where, Ns is the total number of data sets, Np is the number of particles and tmax is the total number of time origins. Figure 9(a), shows the measured values of the VAF normalized by hv(0)2 i in our system. In simulations of hard sphere fluids, Alder and Wainwright [40] first found that the form for the VAF was strongly depended on the density of the system. For very low densities the the characteristic form of the correlation function was given simply by 12 0 0 0.2 (10) where τE is the Enskog collision time. If the density of the system is increased however, the form of Eq. (10) breaks down and Cv (t) can become negative with long range tails due to the caging of particles by their neighbors. We find, that the lowest density case becomes, and remains negatively correlated after the decay from hv(0)2 i [Fig 9(a)]. This appears to be in contradiction with Ref. [40] but is due to a finite size effect. That is, the particles are interacting frequently with the side-walls Nl β a(a) b(b) c(c) 1 2 3 4 4 5 2 2 2 2 4 4 0.031(0.0027) 0.025(0.0025) 0.008(0.0003) 0.003(0.0008) 0.002(0.0002) 0.002(0.0005) 0.428(0.511) 0.603(0.665) 0.932(1.203) 1.258(1.209) 0.970(1.043) 1.096(1.384) 0.154(2.969) 0.393(6.024) 0.742(8.306) 1.489(16.91) 3.171(26.45) 2.887(28.33) TABLE II: Fitting parameters for Eqs. (7a,b). The values are arranged a(a) for P (l)(P (τ )), respectively. The (· · · ) correspond to the values for P (τ ) φ 0.4 0.6 FIG. 9: (a) The velocity auto correlation function Cv (t). (b) The mean square displacement in the x-direction for each Nl β for t = 0 − 2.5 s. Inset: The short time behavior indicated by the box in the main figure. (c) The diffusion constants calculated for each φ, where (◦) corresponds to the numerical integration of Cv from (a), and (⋄) corresponds to the least squares fit of Cx2 from (b). The dashed line shows the kinetic theory result for a fixed temperature, which is given by the average measured value over this range of φ. at low densities, which reverse the sign of velocity vectors, thus leading to the observed anti-correlation. The predominance of the sidewall interactions are screened for the intermediate densities due to the increased number of particle-particle collisions, therefore no anti-correlations are observed. C. Mean Square Displacement To determine the mean square displacement of the xcomponent of the particle positions [Fig 9(b)], by using the sequential time zeros for each trajectory [39] given by 8 Np ,Ns tmax X X 1 Cx2 (t) = |xij (to ) − xij (to + ∆t)|2 . Np Ns tmax i,j=0 solid line in Fig. 9(c) shows the form of Eq. (14) with the temperature T given by the granular temperature. The granular temperature is defined by ∆t=1 (11) Where Ns is the total number of data sets, Np is the number of particles and tmax is the total number of time origins. For these measurements, the ROI is allowed to increases in size along the y-direction as Nl decreases. We have chosen to make this increase to ensure that particles at low Nl have had the opportunity to undergo a collision while under consideration. The long time behavior of Cx2 for each Nl β [Fig. 9(b)], displays linear dependence on time, indicating diffusive behavior. However for Nl = 1, Cx2 clearly shows a crossing from one linear regime to another, which may be a possible indication of finite system size for low density. For short times, [Fig. 9(b)Inset] the behavior is ballistic as indicated by the quadratic increase of Cx2 in time. As Nl is increased, the range of the ballistic regime dramatically decreases indicating a decrease in the Enskog collision time τE . The ballistic and diffusive regimes are consistent with what is expected for kinetic theory of elastic, finite-sized particles. D. Self-Diffusion The self diffusion constant D, can be determined for a system of particles by either evaluating the time integral of the velocity autocorrelation function, Z ∞ D= Cv (t)dt, (12) 0 or using the relationship between the mean square displacement of the particles and time over long times, D = lim t→∞ 1 Cx2 (t), 2dt (13) where d is spatial dimension. From kinetic theory [38], the diffusion constant of a two dimensional gas is calculated as:  1/2 σ πT D= (14) 8φg(σ) m where, g(σ) is the radial correlation function at contact [41] given by, g(r = σ) = 16 − 7φ . 16(1 − φ)2 (15) By numerically integrating the curves in Fig. 9(a), and performing a least squares fit to the data in Fig. 9(b) after the ballistic regime, we obtain the self-diffusion constant [see Fig. 9(c)]. We find that the values for the self-diffusion from Eqs. (12,13) are self consistent. The Tx,y = 1 m[hvx2 i + hvy2 i], 2 (16) where, m is the mass of the particles and h · · · i denote averages over the component distributions [see Sec. IV E]. The other constants in Eq. 14 are determined from system parameters. The theory for the diffusion of elastic particles, given by Eq. (14) closely matches our results all φ. Thus we show that the effects of inelasticity on the self-diffusion are small. E. Velocity Distributions The distribution of the x- and y-components of the particle velocities are plotted in Figs. 10–12(a–f). The distributions correspond to velocities that are measured within a region of interest. The ROI is defined by making a narrow slice across the y-direction that is centered upon the peak in ρ(y) while excluding particles lying within a distance of 3σ from the side walls. We utilize this ROI to ensure that large gradients in ρ(y), and the clustering produced by the side-walls, do not affect the measured VDFs. Each distribution correspond to ∼ 2 × 106 unique velocities that are found within our ROI. The velocities of elastic particles follow a distribution given by the Maxwell-Boltzmann form, P (v) = (2πkB T )−d/2 e−v 2 /2kB T (17) where, d is the dimensionality of the system and T is the temperature of the heat bath that the system is in contact with. Hence, if a system of particles is at equilibrium, its temperature determined by the width of the distribution of particle velocities. Equation (17) is fit to the data for the x-components of the velocities, and is shown on both linear and logarithmic scales [Fig. 10,11]. We observe that the form given by Eq. (17) displays deviations both at low and high velocities. The distributions of velocities are normally displayed in a log-linear fashion to accentuate the tails of the VDF, however this suppresses the deviations at low velocities. By plotting the distributions on a linear scale we display the more statistically significant deviations from Eq. (17). In a recent experimental work [17], a two dimensional collection of particles is driven into a steady state. Using analysis techniques that are similar to ours, the authors proposed a governing form for the VDF given by R(v) = A e−B |vx /Tx | −1.5 , (18) where A and B are constants and Tx is the x-component of the granular temperature defined in Eq. (16). They claim to have seen a universal VDF which they parameterized by a single value, regardless of the system density 9 -1 -1 0.03 10 10 0.03 (a) (b) (a) -2 (b) -2 10 10 0.01 P (vx) 0.02 P (vx) P (vx) P (vx) 0.02 -3 10 -4 0.01 -4 10 10 -5 0 -20 0 0 -20 20 -1 -5 10 0 -60 20 10 -40 20 40 60 -60 (d) (c) -2 10 P (vx) P (vx) P (vx) P (vx) 20 40 60 -1 10 0.02 0 10 0.04 0.04 -20 -1 10 0.06 -40 vx (cm s ) -1 0.08 (c) 0 vx (cm s ) vx (cm s ) 0.06 -20 -1 -1 vx (cm s ) -3 10 -3 10 (d) -2 -3 10 10 -4 -4 10 0.02 -5 10 -5 0 0 -20 20 -1 10 0 -6 20 -60 -1 vx (cm s ) -40 0.06 40 (f) P (vx) 0.04 20 -1 0 -20 (e) 10 -3 10 -4 10 20 -1 FIG. 10: The velocity distribution functions P (vx ) versus vx on a linear-linear scale. (a) Nl = 1 , β = 2.0, (b) Nl = 2 , β = 2.0, (c) Nl = 3 , β = 2.0, (d) Nl = 4 , β = 2.0, (e) Nl = 4 , β = 4.0, (f) Nl = 5 , β = 4.0. The solid curves are a least squares fit to a Gaussian form given by Eq. (17). Note that the deviation from a Gaussian distribution extends all the way to the lowest velocity bins. Each distribution correspond to ∼ 2 × 106 unique velocities that are found within the ROI defined in the text. or the value of the inelasticity of the particles. From our VDFs, whose corresponding densities range over an order of magnitude and where the average inelasticity varies by nearly a factor of two, we cannot find any single parameter fit that describes the overall form. The VDFs for the y components, P (vy ) versus vy for each Nl β in our ROI are also measured [see Fig. 12(a– f)]. The VDFs are highly skewed by the asymmetry in the driving against the direction of gravity. To identify the effects that the asymmetry in P (vy ) has upon P (vx ), we have separated the vx distributions by the sign of vy , i.e. P (vx | + vy ; −vy ). We have found that the form for these conditional distributions are not affected by the sign of vy , however we do note that their widths differ, with hvx i+vy < hvx i−vy . We also measure the x-component of the granular temperature Tx [see Eq. (16)], to probe the scaling behavior of the velocity distributions. Figure 13(a) shows the measured granular temperature as a function of distance from the driving wall. At low densities (Nl ≤ 2), Tx (y) initially increases and then decays. In contrast, for (Nl ≥ 3), 10 -60 40 60 (f) -2 -3 10 10 -4 -5 10 -6 0 20 -1 -5 vx (cm s ) 0 vx (cm s ) 10 0 -20 -1 -2 0.02 vx (cm s ) -40 10 10 0.02 0 -20 10 -60 60 -1 0.06 P (vx) 20 10 0.04 P (vx) 0 vx (cm s ) 0.08 (e) -20 -1 vx (cm s ) P (vx) 0 -20 -6 -40 -20 0 20 40 60 10 -60 -1 -40 -20 0 20 40 60 -1 vx (cm s ) vx (cm s ) FIG. 11: The velocity distribution functions P (vx ) versus vx on a log-linear scale.(a) Nl = 1 , β = 2.0, (b) Nl = 2 , β = 2.0, (c) Nl = 3 , β = 2.0, (d) Nl = 4 , β = 2.0, (e) Nl = 4 , β = 4.0, (f) Nl = 5 , β = 4.0. The solid curves are a least squares fit to a Gaussian form given by Eq. (17). Here the apparent deviation in the tails of the distribution functions are present. Tx (y) has a distinct minimum. We note that Tx (y) never reaches a constant value and the minimum (maximum) does not correspond to the peak in ρ(y). To further show the non-universality of the VDFs, we plot the kurtosis as a function of distance from the driving wall. The kurtosis is obtained by the following: γ= hvx4 i hvx2 i2 . (19) If the velocity distribution is a Gaussian then γ = 3, shown by the solid line in Figure 13(b), and if the distribution is given by Eq. (18) then γ = 3.576. We find that the the measured values for γ exceed the value for a Gaussian and also vary as a function of distance from the driving. This analysis is consistent with our previous results [21] and recent MD simulations of Brey and Ruiz-Montero [27] that closely mimic our experiment. 10 -1 -1 P (vy) -3 10 -3 10 -4 -4 10 10 -5 -5 -90 -60 -30 0 30 10 -120 60 -1 -90 -60 -30 0 30 60 -1 -1 2 0 10 P (vy) 10 -3 10 10 -3 8 -4 10 -5 10 Nl = 1, β= 2 Nl = 2, β= 2 Nl = 3, β= 2 Nl = 4, β= 2 Nl = 4, β= 4 Nl = 5, β= 4 (b) 10 -4 10 (d) -2 γx (c) 10 4 vy (cm s ) 10 -2 6 -1 vy (cm s ) P (vy) (a) −2 10 10 P (vy) (b) -2 Tx (g cm s ) (a) -2 10 -120 8 10 10 6 10 -90 -60 -30 0 30 -5 -90 60 -60 0 30 60 vy (cm s ) -1 -1 10 (e) P (vy) 10 10 -3 10 10 (f) -2 10 -3 10 -5 -5 -6 -60 -30 0 30 60 10 -90 -60 -1 -30 0 30 60 -1 vy (cm s ) vy (cm s ) FIG. 12: The velocity distribution functions P (vy ) versus vy on a log-linear scale.(a) Nl = 1 , β = 2.0, (b) Nl = 2 , β = 2.0, (c) Nl = 3 , β = 2.0, (d) Nl = 4 , β = 2.0, (e) Nl = 4 , β = 4.0, (f) Nl = 5 , β = 4.0. The large skewness in the distributions for the negative values of vy is due primarily to the driving from the bottom wall. Particles that are moving in the −y direction are leaving the moving wall. F. 6 12 -4 -4 -90 0 y (cm) 10 10 4 2 10 -2 P (vy) -30 -1 -1 vy (cm s ) Equation of State The equation of state for ideal gases relates the pressure to the temperature and the density, P = nkB T, (20) where n is the number density and kB is Boltzmann’s constant. If we assume that kinetic theory is valid for a granular gas, we can immediately relate the average squared speed, hv 2 i of the particles to the temperature, mhv 2 i = kB T, (21) for each degree of freedom. Figure 14(a) shows the pres2 sure, P = nm 2 hvy i as a function of distance from the driving wall. Due to the effects of gravity on the particles, the density should follow the well known atmospheric law, ρ(y) = ρ0 e−mgy/ Ty , (22) which assumes a constant temperature. We find that the temperature is not constant for any Nl β [see Fig. 13]. This is also consistent with our observations of the density distributions in Sec. III A where ρ(y) deviates from the form of Eq. (22). FIG. 13: (a) The granular temperature 12 mhvx2 i as a function of distance from the driving wall for each Nl β. If the isothermal atmosphere condition was satisfied these would be constant values for all y above the peak in ρ(y). For values of Nl > 2 the temperatures follow a non-monotonic form that has a distinct minimum. (b) The kurtosis, γx measured from P (vx ) as a function of the distance from the driving. The values given by a Gaussian (solid line) and the form proposed in Eq. (18) (dashed) are only attained very far from the energy source. Momentum balance implies that the gradient of the pressure is related to n by the following equation: dP = −nmg dy (23) where, m is the mass of a particle and g is the acceleration of gravity. Due to the non-vanishing gradient of temperature the general form of the pressure gradient must be taken into account dP dn dTy = Ty +n . dy dy dy (24) We find that Eq. (24) is indistinguishable from the numerical derivative of d(nTy )/dy. We have measured the pressure gradient acting on a particle held at a particular y by evaluating,   Ty d(nTy ) = g′, (25) − nm dy where g ′ = 57 g sin(β). Figure 14(b) shows the LHS of Eq. (25) for β = 2, and Fig. 14(c) for β = 4. The solid lines correspond to the values of g ′ with β = 2, 4 found by utilizing Eq. (1) and the data in Fig. 2(b). We find that the measured values systematically overestimate the actual values for g ′ in the region where the density reaches 11 40 Nl = 1, β = 2 Nl = 2, β = 2 Nl = 3, β = 2 Nl = 4, β = 2 Nl = 4, β = 4 Nl = 5, β = 4 −1 −2 P (g cm s ) (a) 30 20 10 (1/mρ) dP/dy 0 (b) −15 −30 (1/mρ) dP/dy −45 (c) −15 −30 −45 −60 0 6 12 y (cm) FIG. 14: (a) The pressure P (y) = ρT as a function of y, the distance from the driving wall. (b) The ratio of the mass 1 dP density to the granular pressure force mρ as a function of dy distance from the driving wall. The solid line corresponds to g ′ for β = 2◦ (c) The same as (b) for β = 4◦ . The obvious deviation for regions of high density show a breakdown of the simple treatment of the granular equation of state in regions of high density. it’s maximum. Our interpretation assumes a dilute gas, therefore the deviations near the peak in ρ(y) are not surprising. To incorporate the effects due to increased density, we have obtained the pressure from the interpolated equation of state derived by Grossman et al. [6]: P = nTy nc + n nc − n (26) where, nc is the close packing number density. Utilizing Eq. 26 to calculate the gradient of the pressure [Eq. (23)]. However, we find that the disagreement persists between the LHS of Eq. (25) and the measured value of g ′ near the peak in ρ(y). V. SUMMARY AND CONCLUSION In this paper we have presented a statistical analysis of an inelastic gas that is constrained to two dimensions. Utilizing high speed digital image processing we perform long time tracking over a broad range of densities. Not surprisingly, we observe that the statistical properties of inelastic gases deviate from expectations of the kinetic theory for smooth elastic particles. The most apparent discrepancies are found in the distribution of free paths and times and the distribution of particle velocities. To characterize our system we measure the effective coefficient of restitution from the relative pre- and postcollision velocities of particles undergoing binary collisions. We find that the normal component of restitution and the energy inelasticity are not single-valued, but have a distribution of values even for the same impact parameters. The mean value of the normal components of restitution systematically decreases with the system density. We also find that the energy inelasticity can take on values greater than unity, demonstrating a transference of energy from the rotational to linear degrees of freedom. In a recent numerical work Barrat and Trizac have measured the projected one-dimensional coefficient of restitution Ref. [37, 42] and the energy inelasticity. Their interpretation is that the coefficient of restitution and the energy restitution are random variables that characterize collisions, consistent with our findings. The distribution of path lengths and free times are shown to have an overpopulation of the short distance and time bins. We have proposed an empirical form in Equations (7 a,b) which capture the overall behavior of the observed distributions of paths and times. By integrating these distribution functions, we are able to measure the mean free path and mean time. The average speed obtained from the speed distribution and from the mean free path and time are in close agreement. Inspired by these finding, Paolotti et al. [26] have reported similar results for the mean free time in a simulation that mimics our system. Particle diffusion constants are measured from two independent long time averaged correlation functions. The mean square displacement and the velocity autocorrelation function are calculated. By then performing least squares fitting and numerical integration to these quantities respectively, the self diffusion over a broad range in density is calculated. We find that the diffusion constants are similar to that of a two dimensional gas over this density regime. Therefore, long-time averaged correlation functions seem to accurately capture the diffusive properties of granular gases. We find that the distribution of particle velocities perpendicular to the direction of driving does not have a universal form, but depends on both the density and the inelasticity. In addition, we find a distinct asymmetry in the VDFs parallel to the driving direction. We measure the granular temperature as a function of distance from the driving source and find non-monotonic behavior. For low densities, the granular temperature has a distinct maximum and for high densities there exists a distinct minimum. The temperature inversion at higher densities has recently been described via granular hydrodynamics 12 by Ramírez and Soto [43]. However, the crossover from T (y) having a maximum (for low densities), to T (y) with a minimum at high densities has not been discussed in any kinetic or hydrodynamic models. By using kinetic theory and simple hydrodynamics we have tested the force balance between the gradient of the pressure exerted by a granular gas on a particle and the force due to gravity [Eq. (23)]. We find strong deviations in the regions of high density. A simple hydrodynamic form, that describes the behavior over all densities, is not yet available. [1] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Rev. Mod. Phys. 68, 1259 (1996). [2] S. Chapman and T. Cowling, The Mathematical Theory of Non-Uniform Gases (Cambridge, Cambridge University Press, 1970). [3] J. T. Jenkins and S. B. Savage, J. Fluid Mech. 130, 187 (1983). [4] P. K. Haff, J. Fluid Mech. 134, 401 (1983). [5] S. McNamara and W. R. Young, Phys. Fluids A 4, 496 (1992). [6] E. L. Grossman, T. Zhou, and E. Ben-Naim, Phys. Rev. E 55, 4200 (1997). [7] L. P. Kadanoff, Rev. Mod. Phys. 71, 435 (1999). [8] S. Warr, G. T. H. Jacques, and J. M. Huntley, Powder Tech. 81, 41 (1994). [9] A. Kudrolli, M. Wolpert, and J. P. Gollub, Phys. Rev. Lett. 78, 1383 (1997). [10] J. S. Olafsen and J. S. Urbach, Phys. Rev. Lett. 81, 4369 (1998). [11] T. P. C. van Noije and M. H. Ernst, Granular Matter 1, 57 (1998). [12] I. Goldhirsch, M.-L. Tan, and G. Zanetti, J. Scient. Comp. 8, 1 (1993). [13] A. Puglisi, V. Loreto, U. M. Marconi, A. Petri, and A. Vulpiani, Phys. Rev. Lett. 81, 3848 (1998). [14] S. Warr, J. M. Huntley, and G. T. H. Jacques, Phys. Rev. E 52, 5583 (1995). [15] S. Warr and J. M. Huntley, Phys. Rev. E 52, 5596 (1995). [16] R. D. Wildman, J. M. Huntly, and J.-P. Hansen, Phys. Rev. E 60, 7066 (1999). [17] F. Rouyer and N. Menon, Phys. Rev. Lett. 85, 3676 (2000). [18] J. S. Olafsen and J. S. Urbach, Phys. Rev. E 60, R2468 (1999). [19] W. Losert, D. Copper, J. Delour, A. Kudrolli, and J. P. Gollub, Chaos 9, 682 (1999). [20] A. Kudrolli and J. Henry, Phys. Rev. E 62, R1489 (2000). [21] D. L. Blair and A. Kudrolli, Phys. Rev. E 64, 050301 (2001). [22] T. P. C. van Noije, M. H. Ernst, E. Trizac, and I. Pagonabarraga, Phys. Rev. E 59, 4326 (1999). [23] S. J. Moon, M. D. Shattuck, and J. B. Swift, Phys. Rev. E 64, 031303 (2001). [24] I. Pagonabarraga, E. Trizac, T. P. C. van Noije, and M. H. Ernst, Phys. Rev. E 65, 011303 (2002). [25] J. S. van Zon and F. C. MacKintosh, cond-mat/0205512 (2002). [26] D. Paolotti, C. Cattuto, U. Marconi, and A. Puglisi, cond-mat/0207601 (2002). [27] J. J. Brey and M. J. Ruiz-Montero, cond-mat/0210158 (2002). [28] A. Puglisi, V. Loreto, U. MariniBettoloMarconi, and A. Vulpiani, Phys. Rev. E 59, 5582 (1999). [29] J. J. Brey, M. J. Ruiz-Montero, R. García-Rojo, and J. W. Dufty, Phys. Rev. E 60, 7174 (1999). [30] A. Baldassarri, U. M. Marconi, A. Puglisi, and A. Vulpiani, Phys. Rev. E 64, 011301 (2001). [31] W. A. M. Morgado and E. R. Mucciolo, condmat/0204084, Physica A to appear. [32] A. Baldassarri, U. M. Marconi, A. Puglisi, and A. Vulpiani, Phys. Rev. E 64, 011301 (2001). [33] E. Livne, B. Meerson, and P. V. Sasorov, Phys. Rev. E 65, 021302 (2002). [34] E. Khain and B. Meerson, Phys. Rev. E 66, 021306 (2002). [35] M. Argentina, M. G. Clerc, and R. Soto, Phys. Rev. Lett. 89, 044301 (2002). [36] S. Luding, E. Clément, A. Blumen, J. Rajchenbach, and J. Duran, Phys. Rev. E 49, 1635 (1994). [37] A. Barrat and E. Trizac, cond-mat/0205546 (2002). [38] C. E. Hecht, Statistical Thermodynamics and Kinetic Theory (Dover, Dover Publications New York, 1998). [39] M. P. Allen and D. Tildesley, Computer Simulation of Liquids (Oxford, Oxford University Press, 1996). [40] B. A. T. Wainwright, Phys. Rev. Lett. 18, 988 (1967). [41] D. Henderson, Molec. Phys. 30, 971 (1975). [42] A. Barrat and E. Trizac, Phys. Rev. E 66, 051303 (2002). [43] R. Ramírez and R. Soto, cond-mat/0210471 (2002). We acknowledge stimulating discussions with H. Gould, J. Tobochnik, and A. Puglisi. This work was supported by the NSF under Grant No. DMR-9983659, and by the donors of the Petroleum Research Fund.