\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}