Institute of Physics, University of Potsdam, Am Neuen Palais, D-14469 Potsdam, Germany Institute for Plasma Research, University of Maryland, College Park, Maryland 20742

**Olaf Petzschmann, Udo Schwarz, Frank Spahn, Celso
Grebogi, and Jürgen Kurths**

**March 6, 2001
**

Clustering of granular gases is caused by the dissipation occurring in inter-particle collisions. We use the generalized dimension to characterize the structural evolution and the index-of-dispersion in order to quantify the typical length scale of the density clusters of force-free granular assemblies obtained by 3-D numerical computer simulations. We derive then an expression for the expected length-scales of cluster formation from stability analysis in the hydrodynamical approximation. We find a good agreement between our theoretical prediction and numerical experiments.

74.80.-g, 47.20.-k, 05.20.Dd, 02.50.Fz

In recent studies much attention has been devoted to the dynamics of granular gases [1, 2, 3, 4]. The difference of such systems to conventional gases is their dissipative nature, i.e., the energy loss occurring in inter-particle collisions. This dissipation causes the granules to form clusters, which results from an increase in the collision frequency in regions of fluctuative density enhancements. As more and more mechanical energy is dissipated, there is a corresponding increase in the effective cooling and a decrease in pressure. The consequence is that the granular matter migrates into regions of smaller pressure and higher densities, resulting in the cluster formation [4, 5].

In this Letter, we argue and demonstrate that the decisive quantity for the rate of granular clustering is the coefficient of restitution , which is the ratio between the normal components of the rebound velocity and the impact velocity of the colliding particles. In general, is a function of the impact velocity . Intuitively, a lower , corresponding to larger dissipation per collision, leads to a more intense clustering of the granular gas. We validate this notion by numerical and analytical investigations of the cooling process in three-dimensional configuration space. We characterize and quantify the clustering of granular assemblies through a statistical analysis of the phase-space data obtained from event-driven particle simulations. We introduce two statistical measures to accomplish that:

(i) To characterize the increase of inhomogeneity in the particle density with time, we define a quantity in terms of the generalized dimensions. This quantity reflects the -dependence of clustering formation. We find a monotonous dependence of the the rate of clustering on the dissipation.

(ii) To quantify the characteristic size of inhomogeneities, we define another quantity in terms of the index-of-dispersion. We find that the dependence of cluster size on , as obtained from our 3-D simulations, agrees well with our theoretical predictions derived from a stability analysis using a hydrodynamical approximation [e.g. Eq. 6]. =

**Figure:**
Two-dimensional projection of the 3-D configuration space slice
obtained after about 15 collisions per particle, where , the particle diameter is 2 cm and cm.

The reason for cluster formation is that the fluctuations of the *
homogeneous cooling state* (*HCS*)
are unstable [6]. In this case, the molecular
chaos assumption breaks down, the assumption which serves as the
foundation for the formulation of a kinetic theory of Boltzmann-type.
Numerical simulations have shown [6]
that the beginning of the cooling process is well described by the
*HCS*,
while, after a certain critical time, depending on the number of
collisions per particle, the cooling process is slowed down
by velocity fluctuations.
After this critical moment, the cooling is almost no longer determined
by particle-particle but by
cluster-cluster interactions. This would imply that the clusters
have already been formed during this *HCS*-phase.

In our simulations, we use a 3-D event-driven code with
identical spherical particles of diameter
.
Using periodic boundary conditions,
the particle ensemble is simulated in a cubic box of
side length , chosen to correspond
to a filling factor of about 0.02.
Initially, the particles are randomly placed in the
configuration space having initial velocities
with randomly chosen directions according to an initial
distribution function .
We ensure that the initial conditions are homogeneous in
density and isotropic in the velocities.
For convenience, we restrict ourselves to the case
of constant coefficients of restitution fixing its values at 0.1,
0.2, ..., 0.9.
In our simulations only translational degrees of freedom
are considered and, thus, only the normal coefficient of
restitution is taken into account. According to this, the
relative velocity between the collision
partners is projected onto
the direction of the particle centers,
represented by the unit vector . The relative velocity after the collision is then
given by .
We simulate the assemblies of different for the same absolute
time, taking care that the system is still in the *HCS*-phase.

Compared to the 2-D case (see e.g. Fig. 1 in[5] or Fig. 3
in[4]), a 3-D granular assembly shows a much
weaker clustering, when we project, for visualization purpose the 3-D
particle locations onto a
two-dimensional plane, as illustrated in Fig. 1.
Thus, the inhomogeneities in Fig. 1 are
difficult to distinguish from those ones of an ensemble with purely
elastic particle interactions.
Therefore, we use measures from nonlinear dynamics and
statistics to characterize and quantify the spatial inhomogeneities
in our simulations.
First, we split the box in *m* = 2, 4, 8, ..., 4096 sub-boxes.
The equal length *L* of the diagonal of each sub-box is given by
in units of the smallest grid
size , with .
Second, we characterize the cluster formation quantitatively, using
the generalized dimension, based on the Rényi entropies
[7] defined as
with ,
the probability for
particles to be in the box labeled by
.
Recently, we have shown that the generalized dimension , defined
in terms of , is an appropriate tool to
characterize the time evolution of the degree of clustering
[5, 8]. The
generalized dimension is defined by
for each time step *t* and restitution coefficient .
As *q* increases from zero, accentuates more and more the degree
of clumping (high density regions) while, as *q* decreases from zero,
accentuates better the degree of depletion (low density regions).
To characterize the increase of
inhomogeneity with time, we calculate from simulation data
how changes with
*q*, i.e., .
Then, to get a representative measure of the change of inhomogeneity
with time, we introduce the quantity (for )

where and are the times at the beginning and at the end
of the simulation, respectively.
The value , based on the change of the generalized dimension
with *q* reflects the
-dependent structure formation.
Figure 2 shows
the dependence of on indicating
that an initially homogeneous distribution of
particles changes the greatest the smaller the coefficient of
restitution is.
For elastically colliding granules (), one gets corresponding to an unchanged homogeneous particle distribution
during the time evolution.
=

**Figure:**
vs. [Eq. (1)].
The diamonds are the simulation data points, the
solid line is the result of a quadratic regression fit.

This behavior agrees with predictions [4, 5] that,
as the dissipation rate increases,
the clustering becomes more intense.
However, it is not possible to quantify, using this method, the
length-scales of the inhomogeneities in the
particle density of the simulations.
Therefore, we
introduce a new measure defined in terms of the index-of-dispersion
[9].
The test statistic of the index-of-dispersion is given by
with the particle density in the box,
the standard deviation
,
and the number of degrees of freedom *m* (total number of boxes).
This test statistic is that of a goodness-of-fit test
for the hypothesis that the particles are independently and uniformly
distributed in the whole simulation box.
=

**Figure:**
The evolution of the index-of-dispersion [Eq. (2)] in the
critical region vs. the sub-box diagonal
*L* and the simulation time *t* for different restitution coefficients
. The value of is color coded as indicated in
the color bar. Black indicates a
Poissonian distribution (), values higher than 0
(colored) indicate a significant lumpy distribution of the particles
in the configuration space. White corresponds to the highest
deviation from a Poissonian particle distribution.

For a homogeneous particle distribution, the index-of-dispersion satisfies
,
which is characteristic for a Poissonian distribution
of particles.
A higher index-of-dispersion points towards
inhomogeneities, while is typical for a
regular distribution of particles, as it is found, e.g., in lattices.
Therefore, we define the normalized quantity

where we relate
the index-of-dispersion to the critical value
. The significance level of the
-statistics is chosen to be .
Hence, is a measure of the inhomogeneity,
i.e., of the degree of clustering.
The Heaviside function ensures that purely regular
features are excluded from the statistical treatment.
The time evolution of in the critical region () for different values of ,
obtained from the 3-D simulations,
is shown
in Fig. 3.
Here the quantity is color coded and
plotted vs. time *t* and length scale *L*.
The clumping starts at the smallest scales *L*
after a characteristic time.
At later times, inhomogeneities at larger scales evolve. Figure
3 shows the tendency
of the values as a function of *L*, and time.
Due to the finite number of particles used in the simulations,
one observes, in addition to the density inhomogeneities caused by
the dissipation,
fluctuations .

The most prominent length-scale of the
inhomogeneities, under these conditions, is obtained as follows.
From the index-of-dispersion, which is in
the critical region of the corresponding -test statistics
(Fig. 3), we find the following expression for the critical length
scale by averaging the value *L* of space and time, as

Note that the size of the inhomogeneities detected by this procedure
[Eqs. (2)-(3)] are half of the size of
the typical wave length.
Using Eq. (3), we derive the typical wave number
based on the numerical simulations.
Figure 4 shows , Eqs. (3), as a function
of for a finite simulation time.

To interpret and better understand the results of our computer
simulations,
we use a stability analysis in the hydrodynamical
approximation, which describes well the physics of a granular gas
during the *HCS*-phase [1, 2, 3, 4, 5, 6].
The equations of hydrodynamics for a force-free granular gas read
[10, 11] as

where the material derivative is given by ; ":" and ""
denote the double inner product and the dyadic product,
respectively. The quantities ,
, and *T* are the density, the average particle velocity, and
the temperature, respectively.
The heat flux is ,
and the pressure tensor is denoted by
with the shear tensor
,
where *p* is the pressure and is the unity tensor.
The constitutive relations in the dilute limit, based on an expansion
of the distribution function around the
Maxwell-Boltzmann distribution, read for a constant coefficient of
restitution [2]
for the viscosity,
for the heat conductivity, and
for the cooling due to inelastic collisions. Here,
is the particle diameter and
the mass and the Boltzmann constant are
.

In order to analyze the stability of an initially homogeneous granular
gas, we linearize the hydrodynamic Eqs. (4) and
expand the state variables about the ground state
with . For a force-free granular gas,
is constant in time and space, , and the ground state temperature is
a function of time, due to the permanent cooling.
In a periodic or infinite system, the eigenfunctions of the linearized
Eqs. (4) are plane waves , where is the wave vector and is the growth rate. Without loss of
generality, one can
choose . The resulting eigenvalue problem
gives a system of linear equations , where is the
coefficient matrix
with the solubility condition , determining a characteristic fifth order polynomial,
which can easily be reduced to third order.
If one of its roots takes , an
instability is found.
The linear stability analysis yields a condition, for which an
initially homogeneous granular ensemble becomes monotonically
unstable [5],

The subscript 0 denotes the ground state about which the linear
stability analysis is carried out.
To quantify the derivatives in Eq. (5), we apply
the equation of state of an ideal gas with as an independent
variable .
Using this, we get , ,
, and
.
This yields an expression for the critical wave number
.
According to this equation, unstable wave-like oscillations with a
characteristic length scale only occur if
, which
quantifies the pressure instability mentioned above.
It has been shown that clusters are formed under this condition
either for a constant [4] or variable [5] .

In the dilute
limit, the temperature evolution of an initial homogeneous
granular gas with
and an initial temperature is given by
with , and therefore,
, which, together with
the constitutive relations for and
yields the critical wave number

Equation (6) indicates that clustering sets in after a certain
time , which is the order of the average collision
time for .
On the other hand, for a given finite time *t*, there is a
characteristic coefficient of restitution , where
has a pole and above which no clustering is expected.
Obviously, there is no clustering for . Especially, for () one
gets [dotted line in
Fig. 4], which is the result found by
Goldhirsch & Zanetti [4],
where the derivative has been
assumed to be a negative constant.
Note that the term has almost
no influence on the shape of .
Figure 4 shows given by Eq. (6) as
a function of for a finite time and for .

In Fig. 4, the wave number derived from numerical simulations (stars and solid line for the fit) is compared with the corresponding theoretical relation for [Eq. (6)] in dependence on (dashed line). We find a surprisingly good agreement between the theoretical approach and the numerical experiments. This confirms the validity of the theory to explain cluster-formation, including the sensitive quantitative properties such as the typical cluster-size. The rather small deviation of the theoretical curve from the data points is due to fluctuations associated with the finite number of particles, and the averaging procedure. The dotted line corresponds to the wave number of the homogeneous cooling for , indicating that no dramatic change of cluster formation can be expected under homogeneous cooling state conditions.

=

**Figure:**
Critical wave number vs. . The
stars are the data points obtained from the simulations, the solid
line shows the curve for the quadratic regression fit.
The dashed line presents the result of the stability
analysis [Eq. (6)] using the same parameters as in the
simulations, while the dotted one is the result for .

In conclusion, our theoretical investigations, based on a stability analysis of the hydrodynamical equations, Eqs. (4), yield typical cluster sizes depending on the dissipation, represented by the wave number [Eq. (6)]. Using the generalized dimension and the index-of-dispersion, we characterize and quantify the structural evolution of density inhomogeneities. Statistical analyses of 3-D numerical simulation data, based on the index-of-dispersion, yield typical length-scales of cluster formation [Eq. (3)] as a function of the restitution , which is in good agreement with our theoretical predictions [Eq. (6)]. However, the homogeneous cooling state is a transient phase for the evolution of an initially homogeneous force-free granular gas. The linearized hydrodynamic approximation, as well as the constitutive relations, might not hold for granular assemblies with very inhomogeneous particle distributions and for . Therefore, the analysis of the cooling regime beyond the homogeneous cooling state, both through theoretical investigations and numerical simulations, is needed.

This work was supported by the "Deutsche Forschungsgemeinschaft"
(grant *Sp 384/5-1*), the "Deutscher Akademischer
Austauschdienst", and the "Alexander von Humboldt Stiftung".

Tue Mar 6 08:27:03 MET 2001