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