\documentstyle[multicol,aps,psfig]{revtex}
\begin{document}
\pagestyle{myheadings}
\title{Dynamics of Freely Cooling Granular Gases}
\author{Xiaobo Nie$^{1}$, Eli Ben-Naim$^{2}$, and Shiyi Chen$^{1,2,3}$} 
\address{
${}^{1}$Department of Mechanical Engineering, The Johns Hopkins University,
Baltimore, MD 21218\\
${}^{2}$Center for Nonlinear Studies and Theoretical Division,
Los Alamos National Laboratory, Los Alamos, NM 87545\\
${}^{3}$Peking University, Beijing, China}
\maketitle

\begin{abstract}

  We study dynamics of freely cooling granular gases in two-dimensions
  using large-scale molecular dynamics simulations. We find that for
  dilute systems the typical kinetic energy decays algebraically with
  time, $E(t)\sim t^{-1}$, in the long time limit. Asymptotically,
  velocity statistics are characterized by a universal Gaussian
  distribution, in contrast with the exponential high-energy tails
  characterizing the early homogeneous regime.  We show that in the late
  clustering regime particles move coherently as typical local
  velocity fluctuations, $\Delta v$, are small compared with the
  typical velocity, \hbox{$\Delta v/v\sim t^{-1/4}$}.  Furthermore,
  locally averaged shear modes dominate over acoustic modes.  The
  small thermal velocity fluctuations suggest that the system can be
  heuristically described by Burgers-like equations.
\vspace{2pt}
\\ {PACS:} 45.70.Mg, 47.70.Nd, 05.40.-a,  02.50.-r, 81.05.Rm
\end{abstract}
\begin{multicols}{2}
  Freely evolving granular media, i.e., ensembles of hard sphere
  particles undergoing dissipative inelastic collisions, exhibit
  interesting collective phenomena including clustering,
  vortices, shocks, and multiple dynamical regimes.  This dissipative
  nonequilibrium gas system can be described by kinetic theory in
  the early homogeneous phase where density fluctuations and velocity
  correlations are relatively small. However, quantitative
  characteristics of the late clustering regime such as the typical
  length and velocity scales, particle velocity distributions, 
  as well the corresponding continuum
  theory remain open questions
  \cite{Haff,Goldhirsch1,Goldhirsch2,McNamara1,McNamara2,Brey1,Deltour,Noije,Brito,Luding1,Luding2,Ben-Naim,Chen,Trizac,Orza,Esipov,Brey2,Brilliantov,Huthmann}.

  Recent molecular dynamics simulations in one-dimension have shown
  that the asymptotic energy decay is universal and that clustering
  corresponds to the formation of shocks in the inviscid Burgers
  equation \cite{Ben-Naim}. In this paper we study long time
  asymptotic dynamics of freely evolving granular gases in two
  dimensions using Molecular Dynamics (MD) simulations.  Treating the
  particles as identical, totally undeformable hard disks, we
  developed an efficient event-driven algorithm
  \cite{Alder,Rapaport,Marin} that allows us to probe the system well
  into the clustering regime. Our main result is that the asymptotic
  dynamics of dilute granular gases is universal as it is independent
  of the degree of dissipation.
  
In the simulations, $N=10^6$ identical disks of radius $R=0.15$ and
mass $m=1$ are placed in a two-dimensional system of linear dimension
$L=10^3$ with periodic boundary conditions.  The particle
concentration, $c=1$, so the mass density (or volume fraction) is
$\alpha=c\pi R^2=0.0707$. Initially, particles are distributed
randomly in space and their velocities are drawn from a Gaussian
distribution with zero mean and unit variance.  To ensure random
initial conditions, the system is first evolved under elastic
collisions only, with each particle undergoing $10^2$ collisions on
average. Then, time is reset to zero and each particle experiences up
to an average of $4 \times 10^4$ inelastic collisions.  In such a
collision the particle velocity, ${\bf v}$, changes according to
\begin{equation}
\label{rule}
  {\bf v}\to {\bf v}-\frac{1}{2}(1+r)({\bf g}\cdot {\bf n})\,{\bf n},
\end{equation}
with ${\bf g}={\bf v}-{\bf v'}$ the relative velocity of the colliding
particles and ${\bf n}$ a unit vector connecting the centers of the
two particles. In each collision, the normal component of the relative
velocity is reduced by the restitution coefficient, $0\leq r\leq 1$,
and the energy dissipation equals ${1\over 2}(1-r^2)({\bf g}\cdot{\bf
n})^2$. In the following, time is quoted in units of
$t_0=\frac{3}{20cR}\sqrt{m/E_0}$ which is proportional to the initial
mean free time, where $E_0$ is the average particle energy initially.

Inelastic collapse, the formation of a cluster via an infinite series
of collisions occurring in a finite time, poses a difficulty for
numerical simulations\cite{McNamara1}. To properly resolve such finite
time singularities, we take collisions to be perfectly elastic ($r=1$)
when the relative velocity falls below a pre-specified threshold ${\bf
g}\cdot {\bf n}<\delta$. This scheme, which mimics actual granular
particles where $r\to 1$ for sufficiently small relative velocities,
can be applied as long as the cutoff velocity $\delta$ is much smaller
than the root mean square (rms) velocity, $v_{\rm rms}\equiv\langle 
v^2\rangle^{1/2}$ \cite{Ben-Naim,Bizon1}. In our simulations, $v_{\rm
rms}>10^{-3}$, and $\delta=10^{-5}$.

Given the initial conditions, the collision sequence is deterministic,
and the system evolves freely (velocities are not controlled).  In the
absence of energy input, the system ``cools'' infinitely.  The average
particle energy, $E(t)={1\over 2}\langle v^2\rangle$, as a function of
time is shown in Fig.~1.  First, we confirmed that for sufficiently
small times, when the spatial distribution of particles remains
roughly uniform, the energy decay follows Haff's cooling law
\cite{Haff}
\begin{eqnarray}
E(t) = E_0\left(1+t/t_*\right)^{-2}.
\label{haff}
\end{eqnarray}
Here, \hbox{$(t_0/t_*)^{-1} = 3\sqrt{\pi}(1-r^2)g(2R)/20$}, where $g(r)$ is
the pair correlation function \cite{Luding1,Luding2}.


However, this cooling law eventually breaks down due to the formation
of dense clusters
\cite{Goldhirsch2,Deltour,Brito,Luding2,Ben-Naim,Chen,Orza}.  The
deviation from this law occurs at a time scale that ultimately
diverges when collisions become elastic.  Interestingly, our
simulations show that a universal decay law
\begin{eqnarray}
E(t) \simeq A(r) t^{-1},
\label{energy}
\end{eqnarray}
holds beyond this time scale. This implies that the typical velocity
decays with time according to \hbox{$v\sim t^{-1/2}$}. The typical
length scale explored by a particle, \hbox{${\cal L}(t)\sim \int_0^t
dt' v(t')\sim t^{1/2}$}, remains small compared with the system size,
${\cal L}\ll L$, and thus, finite size effects were negligible
throughout the simulations. The inset to Fig.~1 shows that the energy
decay law is independent of the cutoff velocity as long as $\delta
\leq 10^{-4}$ (we also verified that numerical errors were
irrelevant).

\bigskip
{\psfig{file=fig1.eps,width=240pt}}
\noindent
{\small FIG.~1. The average particle energy $E(t)$ as a function of
time $t$ for three different restitution coefficients, $r=0.1$, $0.5$,
and $0.9$.  A line of slope $-1$ is shown for reference.  The inset
shows $E(t)$ as a function of $t$ for $r=0.1$ and different cutoffs
$\delta=10^{-4}$, $10^{-5}$.  Each line represents an average over
four independent realizations.}
\bigskip

Overall, this behavior is consistent with the $r$-independent energy
decay law, $E\sim t^{-2/3}$, found in one-dimension \cite{Ben-Naim},
although the exponent differs from the 1D value of 2/3.  This
asymptotic law (\ref{energy}) has not been observed in previous
numerical studies primarily due to insufficient temporal
range\cite{Goldhirsch2,Deltour,Brito,Luding2,Chen,Orza}.  For example,
power-law behavior was previously suggested by deformable sphere
molecular dynamics simulations in two- and
three-dimensions\cite{Chen}. However, the value of the scaling
exponent depended on parameters such as the system size and the
restitution coefficient.

The simulations show that as the volume fraction $\alpha$ was reduced,
the dependence of the prefactor $A(r)$ on the restitution coefficient
$r$ became weaker and weaker, suggesting completely universal behavior
in the dilute limit, $\alpha\to 0$. Further numerical studies with
smaller volume fractions are needed to fully resolve this issue.  If
indeed the prefactor $A$ is independent of the restitution coefficient
$r$, this implies a certain behavior of $t_c$, the time scale marking
the transition from the homogeneous to the clustering regime.
Matching the two asymptotic behaviors $E(t)\sim [(1-r)t]^{-2}$ for
$t\ll t_c$ with $E(t)\sim t^{-1}$ for $t\gg t_c$ shows that this
crossover time scale diverges according to $t_c\sim (1-r)^{-2}$ in the
quasi-elastic limit, $r\to 1$.



Next, we examined whether the entire velocity distribution, not merely
the typical velocity scale is independent of the restitution
coefficient. We studied $P(v,t)$, the probability distribution
function of the velocity magnitude $v\equiv |{\bf v}|$ at time 
$t$. Figure 2 clearly shows that a universal scaling function
underlies the velocity distribution when $t\gg t_c(r)$
\begin{equation}
P(v,t)\sim {1\over v^2_{\rm rms}}\Phi(z),\qquad z={v\over v_{\rm rms}}.
\end{equation}
This provides evidence that the asymptotic dynamics are universal.

The simulations suggest that the corresponding scaling function is
Gaussian, $\Phi(z)=\pi^{-1}\exp(-z^2)$, (see Fig.~2). We therefore
studied the kurtosis, particularly, the quantity $\kappa=\langle
v^4\rangle/\langle v^2\rangle^2-2$.  Initially, $\kappa = 0$, the
expected value for a Gaussian distribution in two-dimensions.  In the
intermediate regime, $\kappa $ is distinctly nonzero and strongly
depends on the value of the restitution coefficient. However, in the
asymptotic regime ($t\gg t_c$), $\kappa $ does not depend on
restitution coefficient and is very close to zero ($|\kappa|<0.1$), a
manifestation of the universal velocity statistics.  Given this
universality, we present data for a representative value of $r=0.1$ in
the rest of this paper.

\bigskip
{\psfig{file=fig2.eps,width=240pt}}
\noindent
{\small FIG.~2.  The scaling velocity distribution $\Phi(z)$ versus
the square of scaling variable $z=v/v_{\rm rms}$. The data corresponds
to three different restitution coefficients $r=0.1$, $0.5$, and $0.9$,
at three different times, $t=5 \times 10^2$ (excluding $r=0.9$), $ 2
\times 10^3$, $2 \times 10^4$, all in the clustering regime.  These
eight distributions follow a universal scaling function, a Gaussian
distribution.  The inset shows the exponential high-energy tail of
$P(v)$ at $t=10$ for $r=0.1$, where the velocity was rescaled by the
rms velocity.}
\bigskip

For completeness, we mention that in the intermediate homogeneous
regime, the distribution has an exponential high-energy tail, in
agreement with kinetic theory studies
\cite{Goldhirsch2,Noije,Esipov,Brey2,Brilliantov,Huthmann}. The
measured slope of $3.3$ is consistent with Direct Simulation Monte
Carlo (DSMC) \cite{Brey2} and MD simulations \cite{Huthmann} (see the
inset in Fig.~2).

\bigskip
{\psfig{file=fig3.eps,width=240pt}}
\noindent
{\small FIG.~3.  The adjusted kurtosis $\kappa$ as a function of time
for three restitution coefficients $r=0.1, 0.5$, and $0.9$.}
\bigskip

Both the Gaussian and the exponential extremal velocity statistics are
consistent with the following heuristic argument. The large-velocity
tails are dominated by the fastest particles ($v\approx 1$) initially
present in the system, still moving with their initial velocity at
time $t$. These particles managed to avoid any collisions with other
particles. The ``survival'' probability, $S(v,t)$, for such particles
decays exponentially with the volume they sweep till time $t$,
\hbox{$S(v=1,t) \propto \exp(-{\rm const.}\times t)$}.  Assuming a
stretched exponential large-velocity tail for the scaling function
\hbox{$\Phi(z)\sim \exp(-{\rm const.}\times |z|^{\gamma})$} as
$z\to\infty$ and the typical velocity decay $v \sim t^{-\beta}$ as
$t\to \infty$, and then matching the dominant behavior of the fastest
particles $P(1,t)=P(1,0)S(1,t)\propto \exp(-{\rm const.}\times t)$
with the expected scaling behavior $P(1,t)\propto
\Phi\left(t^{\beta}\right)\propto \exp\left(-{\rm const.}\times
t^{\beta\gamma}\right)$ yields the exponent relation $\beta\gamma=1$.
In the homogeneous regime, $\beta=1$ and thus $\gamma=1$, while in the
clustering regime, $\beta=1/2$ and thus, $\gamma=2$.

\bigskip
{\psfig{file=fig4.eps,width=240pt}}
\noindent
{\small FIG.~4. The ratio of the thermal energy $E_{th}$ to the
kinetic energy $E_k$ as a function of time for $r=0.1$. A line of
slope $-1/2$ is also shown as a reference.}
\bigskip


To quantify collectives motions and the corresponding velocity
fluctuations in the clustering regime, we also measured the local
kinetic energy, $E_k$, and the local thermal energy, $E_{th}$.  The
local kinetic energy is defined by $E_k=\frac{1}{2} \langle\rho u^2
\rangle$ with the density $\rho$ and the velocity ${\bf u}\equiv (u_x,
u_y)$ obtained by averaging the corresponding quantities over a small
region of space.  The local thermal energy $E_{th}=\frac {1}{2}\langle
(v_x-u_x)^2+(v_y-u_y)^2\rangle$ is obtained by subtracting the average
local velocity ${\bf u}$ from the particle velocity ${\bf v}$ . Space
was divided into $128\times 128$ small boxes. As shown in Fig.~4,
initially, the thermal energy is large compared with the kinetic
energy, indicating that particles essentially move independently and
there are no collective motions.  In contrast, in the asymptotic
regime, $t\gg t_c(r=0.1)\approx 10^2$ the thermal energy becomes much
smaller than the kinetic energy, indicating that coherent motion of
particles becomes dominant.  The appearance of the collective motions
was suggested by linear stability analysis of the hydrodynamic
equations\cite{Goldhirsch1,Deltour} and observed in MD simulations
\cite{Goldhirsch1,Goldhirsch2,McNamara2,Deltour,Brito,Luding1,Luding2,Ben-Naim,Chen}.



Interestingly, we find that the ratio of thermal to kinetic energy
decays algebraically, $E_{th}/E_{k}\sim t^{-0.5}$, in the long time
limit.  The thermal energy, $E_{th}=\frac {1}{2}\langle (\Delta
v)^2\rangle$ quantifies local velocity fluctuations, $\Delta v$, and
using Eq.~(\ref{energy}),
\begin{equation}
 E_{th} \sim (\Delta v)^2\sim t^{-3/2}. 
\end{equation} 
In other words, local velocity fluctuations, $\Delta v\sim t^{-3/4}$
are small compared with the typical velocity, $v\sim t^{-1/2}$, as
$\Delta v/v\sim t^{-1/4}$. Thus, at least two distinct velocity scales
are needed to characterize particle velocities in the clustering
regime. Moreover, the relatively small velocity fluctuations imply
strong velocity correlations and a well-defined average local
velocity.  Thus, a hydrodynamic description is plausible in the
clustering regime.

\bigskip
{\psfig{file=fig5.eps,width=240pt}}
\noindent
{\small FIG.~5. The ratio of the longitudinal energy $E_l$ to the
transverse energy $E_t$ as a function of time.}
\bigskip

To distinguish between acoustic and shear modes, we transformed the
averaged velocity field into Fourier space and decomposed the velocity
as ${\bf u}({\bf k}) = {\bf u}_l({\bf k}) + {\bf u}_t({\bf k}) $.
Here ${\bf u}_t({\bf k})$ is the transverse velocity satisfying ${\bf
k} \cdot {\bf u}_t = 0$, and ${\bf u}_l({\bf k})$ is the longitudinal
velocity with $ {\bf k} \times {\bf u}_l =0 $. The corresponding
representation in physical space is $ {\bf u}(x,y)={\bf u}_l(x,y)+{\bf
u}_t(x,y)$, with $\nabla \cdot {\bf u}_t=0$, and $\nabla \times {\bf
u}_l=0 $.  The averaged longitudinal energy, i.e., the energy of
acoustic modes, is defined as the spatial average $E_l=\frac{1}{2}
\langle \rho u_l^2\rangle$, and the transverse energy, i.e., the
energy of shear modes, is defined as the space average of
$E_t=\frac{1}{2} \langle\rho u_t^2\rangle$. Initially, the system is
isotropic and consequently, the longitudinal and the transverse
energies are equal (see Figure 5). Since collective motions are
negligible initially, statistical fluctuations are isotropic as
well. The ratio of longitudinal to transverse energies decreases
steadily and eventually, it saturates at a value of roughly
$E_l/E_t\to 0.4$ for $t>t_c$. Therefore, shear modes dominate over
acoustic modes in the clustering regime. This behavior is consistent
with the formation of dense, thin, elongated clusters, where particles
move coherently, parallel to the cluster orientation
\cite{Goldhirsch1}. Interestingly, we did not observe obvious
correlations between the cluster size and its velocity.

Collective motion of granular media can be described by mass, momentum
and energy balance equations \cite{Goldhirsch1,Deltour,Jekins,Bizon2}.
We have seen that the thermal energy, i.e., the temperature $T$, is
negligible compared with the kinetic energy, and therefore, we may
expand the system around zero temperature. In particular, we can
ignore the pressure term in the momentum equation due to the fact that
$p \sim T$ while keeping the viscosity term since $\nu \sim
T^{1/2}$\cite{Deltour}.  The resulting governing equations are

\begin{eqnarray}
&\partial_t \rho + \partial_{\alpha} (\rho u_{\alpha})  =  0, \\
\label{density}
&\partial_t (\rho u_{\alpha}) + \partial_{\beta} (\rho u_{\alpha} 
u_{\beta})
 =  \partial_{\beta}\pi_{\alpha \beta}, \\
\label{mom}
&\pi_{\alpha \beta} = \partial_{\alpha} (\rho \nu u_{\beta}) +
 \partial_{\beta} (\rho \nu u_{\alpha}) -
 \partial_{\gamma} (\rho  \nu u_{\gamma}) \delta _{\alpha \beta},
\label{eos}
\end{eqnarray}
taken in the limit of vanishing viscosity, $\nu \rightarrow 0$.  The
above equations are similar to the two-dimensional Burgers equation
supplemented by the continuity equation, used to model large-scale
formation of matter in the universe\cite{Burgers,Shandarin}. Here, the
formation of shocks corresponds to dense, thin, string-like clusters
where particles move parallel to the cluster orientation. At least
qualitatively, this is consistent with the above findings that shear
modes dominate over acoustic modes.

In conclusion, our simulations show that in the low density and large
system limits, the kinetic energy of freely evolving granular media
decays with time as $t^{-1}$ in the long time limit. We have also
observed that the particle velocity distribution is Gaussian in this
asymptotic time regime. The above behavior is independent of the
degree of dissipation. We have studied the origin of this universality
and found that the asymptotic dynamics are dominated by collective
motions of particles or alternatively, shear modes.  In the clustering
regime, strong spatial velocity correlations develop, as local
velocity fluctuations are much smaller compared with the typical
velocity. We have also argued that the governing continuum equations
are similar to the two-dimensional Burgers equation with infinitely
small viscosity.

{\bf Acknowledgments.}  The simulations were performed on the Johns
Hopkins University cluster computer supported by the US NSF
(CTS-0079674).  This research was supported by China's NSF (10128204),
and the US DOE (W-7405-ENG-36)

\begin{references}

\bibitem{Haff} P.~K.~Haff, J. Fluid Mech. {\bf 134}, 401 (1983).

\bibitem{Goldhirsch1} I.~Goldhirsch, and G.~Zanetti,
Phys. Rev. Lett. {\bf 70}, 1619 (1993).

\bibitem{Goldhirsch2} I.~Goldhirsch, M.~L.~Tan, and G.~Zanetti,
J. of Sci. Comp. {\bf 8}, 1 (1993).

\bibitem{McNamara1} S.~McNamara, and W.~R.~Young,
Phys. Rev. E {\bf 50} R28 (1994).

\bibitem{McNamara2} S.~McNamara, and W.~R.~Young,
Phys. Rev. E {\bf 53}, 5089 (1996).

\bibitem{Brey1} J.~J.~Brey, M.~J.~Ruiz-Montero, and D.~Cubero,
Phys. Rev. E {\bf 54}, 3664 (1996).

\bibitem{Deltour} P.~Deltour, and J.~L.~Barrat,
J. Phys. I {\bf 7}, 137 (1997).

\bibitem{Noije} T.~P.~C.~van~ Noije, M.~H.~Ernst, 
Gran. Matt. {\bf 1}, 57 (1998).

\bibitem{Brito} R.~Brito, and M.~H.~Ernst,
Erophys. Lett. {\bf 43}, 497 (1998).

\bibitem{Luding1} S.~Luding, M.~Huthmann, S.~McNamara, and
A.~Zippelius, Phys. Rev. E {\bf 58}, 3416 (1998). 

\bibitem{Luding2} S.~Luding and H.~J.~Herrmann,
Chaos {\bf 9}, 673 (1999).

\bibitem{Ben-Naim} E.~Ben-Naim, S.~Y.~Chen, G.~D.~Doolen, and S.~Redner,
    Phys.\  Rev.\  Lett.\ {\bf 83}, 4069 (1999).

\bibitem{Chen} S.~Y.~Chen, Y.~f.~Deng, X.~B.~Nie, and Y.~H.~ Tu,
Phys. Lett. A, {\bf 269}, 218 (2000).

\bibitem{Trizac} E.~Trizac, and A.~Barrat, Euro. Phys Jour. E {\bf 3} 
291 (2000).

\bibitem{Orza} J.~A.~G.~Orza, R.~Brito, and M.~H.~Ernst, con-mat/0002383.

\bibitem{Esipov} S.~E.~Esipov, T.~Poschel, J. Stat. Phys. {\bf 86}, 
1385 (1997).

\bibitem{Brey2} J.~J.~Brey, D.~Cubero, and M.~J.~Ruiz-Montero,
Phys. Rev. E {\bf 59}, 1256 (1999).

\bibitem{Brilliantov} N.~V.~Brilliantov and T.~Poschel,
Phys. Rev. E {\bf 61}, 2809 (2000).

\bibitem{Huthmann} M.~Huthmann, J.~A.~G.~Orza, and R.~Brito,
Gran. Matt. {\bf 2}, 189 (2000).

\bibitem{Alder} B.~J.~Alder and T.~E.~Wainright,
J. Chem. Phys. {\bf 27}, 1208 (1957).

\bibitem{Rapaport} D.~C.~Rapapaport,
J. of Comp. Phys. {\bf 30}, 1158 (1984). 

\bibitem{Marin} M.~Marin, D.~Risso and P.~Cordero,
J. of Comp. Phys. {\bf 109}, 306 (1993).

\bibitem{Bizon1} C.~Bizon, M.~D.~Shattuck, J.~B.~Swift, W.~D.~McComick, 
and Harry L. Swinney, Phys. Rev. Lett. {\bf 80}, 57 (1998).

\bibitem{Jekins} J.~T.~Jenkins and M.~W.~Richman, Phys. Fluids {\bf 28},
3585(1985).

\bibitem{Bizon2} C.~Bizon, M.~D.~Shattuck, J.~B.~Swift, and H.~L.~Swinney,
 Phys. Rev. E {\bf 60}, 4340 (1999).

\bibitem{Burgers}J.~Burgers, {\sl The Nonlinear Diffusion Equation},(Reidel, Dordrecht, 1974). 

\bibitem{Shandarin} S.~F.~Shandarin, and Ya.~B.~Zeldovich, Rev. of Mod. Phys.
{\bf 61}, 185 (1989). 

\end{references}
\end{multicols}
\end{document}



