\documentstyle[graphicx,multicol,aps,epsf]{revtex}
\begin{document}
\title{Dynamics of Vibrated Granular Monolayers}
\author{X.~Nie$^{1}$, E.~Ben-Naim$^{1}$, and S.~Y.~Chen$^{1,2}$}
\address{${}^{1}$Theoretical Division and Center for Nonlinear Studies,
Los Alamos National Laboratory, Los Alamos, NM 87545\\
${}^{2}$IBM Research Division, T. J. Watson Research Center,
P.O. Box 218, Yorktown Heights, NY 10598}
\maketitle

\begin{abstract}
  
  We study statistical properties of vibrated granular monolayers
  using molecular dynamics simulations. We show that at high
  excitation strengths, the system is in a gas state, particle motion
  is isotropic, and the velocity distributions are Gaussian. As the
  vibration strength is lowered, the system's dimensionality is
  effectively reduced from three to two.  Below a critical excitation 
  strength, a cluster phase occurs, and the velocity distribution
  becomes bimodal. In this phase, the system consists of clusters of
  immobile particles arranged in close-packed hexagonal arrays, and
  gas particles whose energy equals the first excited state of an
  isolated particle on a vibrated plate.

\end{abstract}
\smallskip
PACS: 81.05.Rm, 05.70.Fh, 02.70.Ns
\begin{multicols}{2}
  
  Granular media, i.e., ensembles of hard macroscopic particles
  exhibit rich, interesting, and only partially understood collective
  behavior \cite{jnb,lpk,dg}.  The dynamics of driven granular media
  is particularly important since the external energy source balances
  the energy loss due to the dissipative nature of the interactions.
  Understanding of the collective behavior of such systems is
  therefore important for establishing a more complete theoretical
  description of granular media. Here, we focus on monolayer
  geometries where collapse, clustering, and long range order have
  been observed recently
  \cite{olafsen1,olafsen2,losert1,losert2,dlk,gzb}. In particular,
  experimental studies reported that the velocity distribution
  function may exhibit both Gaussian and non-Gaussian behavior under
  different driving conditions.
    
  In this study, we carry out molecular dynamics simulations of
  vertically vibrated granular monolayers. To validate the simulation
  method, we verified the experimentally observed phase transition
  from a gas phase in high vibration strengths, to a cluster 
  phase in low vibration strengths \cite{olafsen1,olafsen2,losert2}.
  Additionally, we checked that several other details including the
  nature of the phase transition, the transition point, and the
  statistics of the horizontal velocities agree quantitatively with
  the experimental results. In particular, the horizontal energy
  vanishes linearly near the transition point, and the corresponding
  velocity distribution changes from a Gaussian to a non-Gaussian as
  the vibration strength is reduced.
  
  In contrast with the experimental studies, the simulations enables
  us to probe the vertical motion, an important characteristic of the
  dynamics. Our results show that as the system approaches the
  transition point, the vertical energy drops by several orders of
  magnitude. Furthermore, the cluster phase is characterized by a
  coexistence of clusters of immobile particles, and energetic gas
  particles, whose energy can be understood by considering an isolated
  particle on a vibrated plate. We also find that the deviation from
  the Gaussian behavior in the gas phase is directly related to the
  development of an anisotropy in the motion, i.e., significant
  differences between the horizontal and the vertical velocities.
  
  To study the dynamics of vibrated monolayers, we used the standard
  molecular dynamics simulation technique \cite{herrmann}. We
  considered an ensemble of $N$ identical weakly deformable spheres of
  mass $m$, radius $R$, and moment of inertia $I={2\over 5}mR^2$.  The
  simulation integrates the equations of motion for the linear and
  angular momentums, $m\ddot{\bf r}_i = \sum_{j\ne i}{\bf
    F}_{ij}+mg\hat z$, and $I\dot {\bf w}_i=\sum_{j\ne i}{\bf
    r}_{ij}\times {\bf F}_{ij}$, respectively. Here, ${\bf r}_i$ is
  the position of the $i$th particle, ${\bf w}_i$ is its angular
  velocity, and $g$ is the gravitational acceleration.  The contact
  force in the direction normal (tangential) to the vector ${\bf
    r}_{ij}={\bf r}_j-{\bf r}_i$ is denoted by ${\bf F}^n_{ij}$ (${\bf
    F}^t_{ij}$).  The force between two particles is nonzero only when
  they overlap, i.e., ${\bf F}_{ij}=0$ when $|{\bf r}_{ij}|>2R$.  When
  there is an overlap, the normal contact force ${\bf F}^n_{ij}={\bf
    F}^{\rm rest}_{ij}+{\bf F}^{\rm diss}_{ij}$ between the particles
  is a sum of the following forces: (a) A restoring force, ${\bf
    F}^{\rm rest}_{ij} = Y m_i(|{\bf r}_{ij}| - 2R) {\bf r}_{ij}/|{\bf
    r}_{ij}|$, with $Y$ the Young's modulus, and (b) An inelastic
  dissipative force, ${\bf F}^{\rm diss}_{ij} = - \gamma_n m_i{\bf
    v}_{ij}^{n}$.  The tangential force is the frictional force ${\bf
    F}^t_{ij}={\bf F}^{\rm shear}_{ij}= -\gamma_s m_i{\bf
    v}_{ij}^{t}$.  In the above, ${\bf v}_{ij}^{n}=({\bf
    v}_{ij}\cdot{\bf r}_{ij}){\bf r}_{ij}/|{\bf r}_{ij}|^2$, and ${\bf
    v}_{ij}^{t} = {\bf v}_{ij} - {\bf v}_{ij}^n $ are projections of
  the relative velocities (at the point of contact) in the normal and
  tangential directions, respectively. The coefficients $\gamma_n$ and
  $\gamma_s$ account for the dissipation due to the relative motion in
  the normal and tangential directions, respectively. Overall, the
  molecular dynamics method has the advantage that it is amenable for
  parallel implementation, and that it allows handling of collisions
  involving an arbitrary number of particles.  
  
  Initially, particles are randomly distributed on the vibrating
  plate, with a filling fraction $\rho$.  The velocities were drawn
  independently from an isotropic Gaussian distribution.  The plate
  undergoes harmonic oscillations in the vertical direction according
  to $z_p(t)=A(t) \cos(\omega t)$ with $A(t)$ the (slowly varying)
  vibration amplitude, and $\omega=2\pi\nu$ with $\nu$ the frequency.
  When a particle collides with the plate, it experiences the same
  force as if it were to collide with an infinite mass particle moving
  with the plate velocity. The simulation was carried out in a finite
  box with a height chosen to be large enough so that no collisions
  can occur with the box ceiling, and periodic boundary conditions
  were implemented horizontally.  Overall, the simulation parameters
  were chosen to be as compatible as possible with the experimental
  values \cite{olafsen1}: $R=0.595mm$, $g=9.8ms^{-2}$, $Y=10^7
  s^{-2}$, $\nu=70Hz$, $N=2000$, $\rho=0.463$, $\gamma_s=100s^{-1}$,
  and $\gamma_n=200s^{-1}$. Without loss of generality, the particle
  mass is set to unity, $m=1$.  The above parameters imply a
  restitution coefficient of $r=0.95$, and that contacts between
  particles last at most one thirteenth of an oscillation period. We
  verified that the results reported below did not depend sharply on
  the values of $g$, $\nu$, $\rho$, $\gamma_s$, and $\gamma_n$, as
  well as the nature of boundary conditions.
\vspace{-.1in} 
\begin{figure}
%\centerline{\epsfxsize=8cm \epsfbox{fig1.ps}} 
\centerline{\includegraphics[width=8cm]{fig1}}
{\small Fig.~1.  The gas versus the cluster phase.  Shown are
  top-projections of the instantaneous particle positions (left) and
  the cumulative positions over 100 consecutive oscillation cycles
  (right) for vibration intensities $\Gamma=1.0$ (top) and
  $\Gamma=0.6$ (bottom).}
\end{figure}
\vspace{-.15in} 
  
  We are primarily interested in statistical properties of the system
  in the steady state, and especially their dependence on the
  vibration strength, which can be quantified by the dimensionless
  acceleration $\Gamma = A \omega^2 /g $. The quantity $\Gamma$ can be
  tuned by varying either $\omega$ or $A$. We chose to fix the
  frequency and vary the vibration amplitude. This method should be
  valid as long as the time scale underlying the variation is larger
  than the systems' intrinsic relaxation time scales. Indeed, we
  verified that results obtained using slowly varying amplitudes
  and fixed amplitudes are in agreement. Each of our simulations was
  initially run at $\sim 10^3$ oscillation cycles at a constant
  amplitude $A_0=0.25mm$.  Then, the amplitude was slowly reduced in a
  linear fashion according to $A/A_0 = 1-t/\tau$, with the decay time
  $\tau\cong 10^3 s$ or alternatively $\sim 10^5$ cycles. Throughout
  this paper we report measurements of average quantities such as the
  temperature and the velocity distribution.  These were obtained by
  averaging over  $100$ consecutive oscillation cycles using $100$
  equispaced points per cycle.


\vspace{-.1in} 
\begin{figure}
%\centerline{\epsfxsize=7.5cm \epsfbox{fig2.eps}} 
\centerline{\includegraphics[width=7.5cm]{fig2}}
{\small {\bf Fig.~2} The vertical and horizontal temperatures versus the 
vibration strength.}
\end{figure}
\vspace{-.1in} 
  
  Qualitatively, we observe that above a critical vibration intensity,
  $\Gamma>\Gamma_c$, particles are in a gas phase, in which their
  motion is random, as seen in Fig.~1. When $\Gamma<\Gamma_c$ in
  addition to particles in the gas phase, hexagonal ordered clusters
  form, as shown in a snapshot of the system.  Furthermore, particles
  inside these hexagonal clusters are stationary, i.e., they are
  attached to the vibrating plate, while particles outside the
  clusters move appreciably.  This behavior is rather robust as it
  holds at least in the parameter range
  $100s^{-1}<\gamma_s<200s^{-1}$, $50s^{-1}<\gamma_n<200s^{-1}$, and
  $20Hz<\nu<100Hz$.


Experimental measurements of the horizontal temperature, defined by
$T_H ={1\over 2} \langle v_x^2+v_y^2\rangle$, indicate a linear dependence in
the vicinity of the critical point \cite{olafsen2} 
\begin{equation}
 T_H \propto (\Gamma-\Gamma_c).
\end{equation}
Our simulations confirm this linear behavior, as shown in Fig.~2 and
Fig.~3. We have also examined the temperatures corresponding to the
rotational motion. The horizontal rotational temperature is very close
to the horizontal temperature. In contrast, the vertical rotational
temperature is roughly $1/3$ the vertical temperature for
$\Gamma>\Gamma_c$, and it is negligible otherwise.  Furthermore, the
horizontal energy practically vanishes below the transition point, as
this quantity decreases by 3 orders of magnitude for
$\Gamma<\Gamma_c$.  This is reminiscent of a sharp phase transition
and it is therefore sensible to view $\Gamma_c$ as a critical point.
This linear behavior can be used to estimate the critical point, and a
linear least-square-fit yields $\Gamma_c=0.763$, in good agreement
with the experimental value $\Gamma_c=0.77$ \cite{olafsen1}. We
conclude that the near critical behavior observed numerically agrees
both qualitatively and quantitatively with the experiment. 


Simulations also allow measurements of the vertical velocities. We
find that the vertical energy $T_V=\langle{v_z}^2\rangle$ decreases
sharply near the transition point as well.  However, in contrast with
the horizontal energy, it does not vanish below the transition point.
Therefore, the velocities develop a strong anisotropy as the vertical
and the horizontal velocities behave quite differently. This reflects
the fact that the system is far from equilibrium.


At high accelerations the particle motion is nearly isotropic, i.e.,
the ratio of horizontal to vertical energy is of the order unity.
Indeed, this ratio approaches a value of roughly 0.5 (see Fig.~2). In
addition, the velocity distribution is an anisotropic Gaussian (see
Fig.~4), and the system is practically three-dimensional.  However, as
the acceleration is decreased, the two-dimensional geometry becomes
more and more pronounced. The vertical motion dominates over the
horizontal one, and the horizontal velocity distribution departs
strongly from a Gaussian distribution.  Near the phase transition
point the large velocity tail becomes nearly exponential (see Fig.~4).
Below the transition point, a significant fraction of the particles
have a nearly vanishing horizontal velocity, i.e., the distribution of
horizontal velocities is strongly enhanced near the origin.


The deviation from the Gaussian behavior can be quantified using the
kurtosis, defined via the second and fourth moments of the
distribution, $\kappa=\langle v^4\rangle/\langle v^2\rangle^2$.
Indeed, in the limit of high vibration intensities, $\Gamma\gg 1$,
this parameter approaches the Gaussian value $\kappa\to 3$.  On the
other hand, near the phase transition point, i.e., as
$\Gamma\to\Gamma_c$, this parameter approaches the exponential value
$\kappa\to 6$.  It proves useful to examine how the kurtosis depends
on $T_H/T_V$, the ratio between the horizontal and the vertical
energies. As shown in the inset to Fig.~4, the smaller the ratio (or
equivalently, the larger the anisotropy), the larger the deviation
from a Gaussian distribution. Hence, whether the velocity distribution
is Gaussian or not reflects the degree of anisotropy in the particle
motion.  Non-Gaussian distributions has been observed
experimentally\cite{olafsen1,olafsen2,losert2}, theoretically
\cite{godhirsch,puglisi,javier,bcdr}, and numerically
\cite{gzb,javier,isobe,taguchi} in one-, two-, and three-dimensional
geometries.  \vspace{-.1in}
\begin{figure}
%\centerline{\epsfxsize=7.5cm \epsfbox{fig3.eps}}
\centerline{\includegraphics[width=7.5cm]{fig3}}
{\small {\bf Fig.~3} Near critical behavior of the horizontal temperature.  
The critical acceleration $\Gamma_c=0.763$ was determined from the 
linear fit $T_H=B(\Gamma-\Gamma_c)$.} 
\end{figure}
\vspace{-.1in}

\begin{figure}
%\centerline{\epsfxsize=8cm \epsfbox{fig4.eps}} 
\centerline{\includegraphics[width=8cm]{fig4}}
{\small {\bf Fig.~4}
    The distribution of horizontal velocities in the gas state, and in
    the vicinity of the critical point. The velocities were normalized
    by the RMS velocity $v_{\rm rms}=\langle v_x^2\rangle^{1/2}$. A Gaussian
    distribution is plotted as a reference. The inset plots $\kappa$
    the kurtosis of the distribution of horizontal velocities versus
    $T_H/T_V$ the ratio of horizontal to vertical energies. The two
    extremal points correspond to the data sets plotted in this
    figure, i.e., $\Gamma=0.8$ and $\Gamma=6$.} 
\end{figure}
\vspace{-.1in}

The distribution of vertical velocities can be used to distinguish
between cluster and gas particles.  Indeed, the vertical velocity
distribution changes its character from a unimodal to a bimodal
distribution in the cluster phase (see Fig.~5) with the low
velocity peak corresponding to the cluster particles, and the high
velocity peak corresponding to the gas particles. Interestingly, the
location of the high velocity peak does not change as $\Gamma$
decreases. In fact, the nonvanishing energy level of Fig.~2 can be
understood by considering the first excited state energy of a single
particle bouncing on a vibrating plate, which can be calculated to be
\cite{losert1}
\begin{equation}
E_1 = {1\over 6}\left({\pi g \over\omega}\right)^2,  
\end{equation}
or in our case $E_1=8.16\ cm^2s^{-2}$.  We verified this result by
simulating the motion of a single particle on a vibrating plate.
Remarkably, the energy of the gas particles falls within less then
10\% of this value.  Therefore, below the phase transition point
particles residing in clusters are in the ground state, i.e., they are
moving with the plate.  Furthermore, the rest of the particles are in
the first excited state of an isolated ball on a vibrating surface.
Additionally, these gas particles collide mostly ($>95\%$) with the
plate, rather than with other particles.  This indicates that the gas
particles are essentially noninteracting.

The vertical velocity distribution can be used to study the fraction
of particles in each phase by simply integrating the area under the
respective energy peaks.  As shown in Fig.~6, $P_0$, the fraction of
particles in clusters is almost independent of the vibration intensity
below the transition point. As the transition point is approached, this
fraction rapidly decreases and ultimately vanishes for $\Gamma\gg
\Gamma_c$.  Although this quantity does not undergo a sharp
transition, its behavior is consistent with our previous estimate of
the transition point from the horizontal energy behavior.

\vspace{-.1in}
\begin{figure}
%\centerline{\epsfxsize=7.5cm \epsfbox{fig5.eps}} 
\centerline{\includegraphics[width=7.5cm]{fig5}}
{\small {\bf Fig.~5}
    The distribution of vertical velocities in the cluster phase for $\Gamma=0.6$.}
\end{figure}
\vspace{-.1in}




In summary, we have studied the dynamics of vibrated granular
monolayers using molecular dynamics simulations.  We find that the
transition between the gas and the cluster phases can be regarded as a
sharp phase transition, and that the horizontal energy decreases
linearly near the transition point. We have shown that at high
vibration strengths, the particle motion is isotropic, and the
velocity distributions are Gaussian. The deviations from a Gaussian
distribution were found to be closely related to the degree of
anisotropy in the motion. We have also shown that below the phase
transition point, the velocity distribution is bimodal. The cluster
particles move with the plate, while the gas particles are
noninteracting, as they collide primarily with the plate, and their
energy agrees with the first excited state of an isolated vibrated
particle.



Our results agree both qualitatively and quantitatively with the
experimental data. This shows that the underlying phenomena can be
explained solely by the simulated interactions, i.e., dissipative
contact forces.  Other mechanisms, possibly present in the experiment,
such as electrostatic forces, etc., are therefore not responsible for
the phase transition.  It will be interesting to use molecular
dynamics simulations to determine the full phase diagram of this
system by varying the density and the driving frequency.

\vspace{-.1in}
\begin{figure}
%\centerline{\epsfxsize=7.5cm \epsfbox{fig6.eps}} 
\centerline{\includegraphics[width=7.5cm]{fig6}}
{\small {\bf Fig.~6}
    The fraction of particles moving with the plate versus the
    vibration strength. The fraction $P_0$ was obtained from the
    horizontal velocity distribution by integrating the area enclosed
    under the low velocity peak (see Fig.~5).}
\end{figure}
\vspace{-.1in}

\begin{references}

\bibitem{jnb}H.~Jaeger and S.~Nagel, R.~P.~Behringer, 
{\sl Rev. Mod. Phys.} {\bf 68}, 1259 (1996).

\bibitem{lpk} L.~P.~Kadanoff, Rev. Mod. Phys. {\bf 71}, 435 (1999).

\bibitem{dg} P.~G.~deGennes, Rev. Mod. Phys. {\bf 71}, S374 (1999).

\bibitem{olafsen1} J.~S.~Olafsen and J.~S.~Urbach, 
Phys. Rev. Lett.  {\bf 81}, 4369 (1998).

\bibitem{olafsen2} J.~S.~Olafsen, and J.~S.~Urbach, 
Phys. Rev. E {\bf 60}, R2468 (1999). 


\bibitem{losert1} W.~Losert, D.~G.~W.~Cooper, and J.~P.~Gollub, 
Phys. Rev. E {\bf 59}, 5855 (1999).

\bibitem{losert2} W.~Losert, D.~G.~W.~Cooper, J.~Delour, A.~Kudrolli, and 
J.~P.~Gollub, Chaos {\bf 9}, 682 (1999). 

\bibitem{dlk} Y.~Du, H.~Li, and L.~P.~Kadanoff, Phys. Rev. Lett. {\bf 74}, 1268 (1995). 

\bibitem{gzb} E.~L.~Grossman, T.~Zhou, E.~Ben-Naim, Phys. Rev. E {\bf 
55}, 4200 (1997).


\bibitem{herrmann} H.~J.~Herrmann, 
{\em Computer simulation of granular media}, in {\em Disorder and Granular 
Media}, ed. D. Bideau and A. Hansen, 1993 Elsevier Science Publishers.

\bibitem{godhirsch} I.~Goldhirsch, and M.~L.~Tan, Phys. Fluids {\bf 8}, 
1752 (1996).

\bibitem{puglisi} A.~Puglisi, V.~Loreto, U.~M.~B.~Marconi, A.~Petri, 
and A.~Vulpiani, Phys. Rev. Lett. {\bf 81}, 3848 (1998).

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

\bibitem{bcdr}E.~Ben-Naim, S.~Y.~Chen, G.~D.~Doolen, and S.~Redner, 
      Phys.\  Rev.\  Lett.\ {\bf 83}, 4069 (1999). 
 
\bibitem{isobe} M. Isobe, and H. Nakanishi, J.\ Phys.\ Soc.\ Japan\ {\bf
    68}, 2882 (1999).

\bibitem{taguchi} Y.~H.~Taguchi, and H.~Takayasu, Europhys. Lett.
{\bf 30}, 499 (1995).


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


