\documentstyle[graphicx,multicol,aps,epsf]{revtex} \begin{document} \title{Slow relaxation in granular compaction} \author{E.~Ben-Naim$^a$, J.~B.~Knight$^b$, E.~R.~Nowak$^{c,d}$, H.~M.~Jaeger$^d$, and S.~R.~Nagel$^d$} \address{$^a$Theoretical Division and Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, NM 87545} \address{$^b$Department of Physics, Princeton University, Princeton, NJ 08540} \address{$^c$ Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801} \address{$^d$The James Franck Institute, The University of Chicago, Chicago, IL 60637} \maketitle \begin{abstract} Experimental studies show that the density of a vibrated granular material evolves from a low density initial state into a higher density final steady state. The relaxation towards the final density follows an inverse logarithmic law. As the system approaches its final state, a growing number of beads have to be rearranged to enable a local density increase. A free volume argument shows that this number grows as $N=\rho/(1-\rho)$. The time scale associated with such events increases exponentially $\sim e^N$, and as a result a logarithmically slow approach to the final state is found $\rho_{\infty}-\rho(t)\sim 1/\ln t$. Furthermore, a one dimensional toy model that captures this relaxation dynamics as well as the observed density fluctuations is discussed. \end{abstract} \vspace{.2in} \noindent{PACS numbers: 05.40.+j, 81.20.Ev, 82.65.My} \begin{multicols}{2} Systems consisting of many macroscopic particles such as sand and powders exhibit complex behavior despite their apparent simplicity\cite{Jaeger}. Vibrated sand may result in size segregation, rich pattern formation \cite{Evesque}, solitary waves \cite{Melo} or convection rolls\cite{Ehrichs}. Despite a growing interest, a comprehensive understanding of the basic principles underlying granular materials is lacking. Although individual grains are solid, it is inappropriate to classify their collective properties as entirely solid-like or liquid-like. Conventional thermodynamic theory is not applicable to sand as thermal fluctuations are negligible, i.e., $k_B T\equiv0$. A simple, yet fundamental property of granular materials is their densification under applied vibrations. Granular compaction is relevant to production, packing, and transportation of a wide array of products such as food, grains, chemicals, and drugs. Granular compaction can be viewed as a model system for non-thermal relaxation in a disordered medium. A granular assembly provides us with a practically uniform system where upon vibration, the well-defined bulk density evolves from a loosely packed mechanically stable initial state into a denser final state. The system explores available microscopic configurations, and slowly eliminates low-density metastable configurations. In a series of recent compaction experiments, monodisperse glass beads were confined to a long tube and were tapped vertically \cite{Knight,Nowak1,Nowak2}. The waiting time between successive taps was large enough to allow the beads to come to a rest before the next tap. Due to the vibration, the volume fraction increased from a loosely packed initial value of $\rho_0\cong 0.58$ to a a final density $\rho_{\infty}$ close to the random close packed limit $\rho_{\infty}\cong 0.64$. A large number of taps ($>10^5$) was necessary to to reach the final density. Moreover, the time dependence of the density was most consistent with an inverse logarithmic four parameter fit, $\rho(t)=\rho_f-\Delta\rho_{\infty}/[1+B\ln(1+t/\tau)]$. The parameters $\rho_f$, $\Delta\rho_{\infty}$, $B$, and $\tau$ depend only the tapping strength or the acceleration $\Gamma$. Several mechanisms were proposed to explain the temporal relaxation of the density of vibrated granular material \cite{Edwards,Mehta,Mehta1,Barker,Hong,Linz,Nicodemi}. However, these studies did not address the relevant excluded volume interaction between the particles. As the compaction progresses, individual particles move slowly, and when a void of the size of a particle is created, it is quickly filled by a new particle. Particles can not move into the space occupied by other particles. In other words, they interact with their neighbors via a hard core interaction. When the packing fraction is large, voids the size of a particle are rare and a large number of particles must be rearranged to accommodate an additional particle and a local density increase. Following this line of reasoning, we propose a simple heuristic picture based on free volume counting to explain the observed relaxation. We complement this picture with an analytically tractable one-dimensional model where the heuristic picture turns out to be exact. This model also proves very useful for studying density fluctuations around the steady state. \begin{figure}\vspace{-.1in} %\centerline{\epsfbox{ill2.ps}} %\centerline{\epsfxsize=5cm \epsfbox{fig1.ps}} \centerline{\includegraphics[width=5cm]{fig1}} \label{fig1} {\small {\bf Fig. 1} A rearrangement of a growing number of particles is necessary to enable a local density increase.} \end{figure} Consider an ensemble of identical rigid spherical particles of volume $V$ with average density $\rho$ (see Fig. 1). Denoting by $V_0$ the pore volume per particle, we have $\rho=V/(V+V_0)$, or alternatively, \begin{equation} V_0=V{1-\rho\over \rho}. \end{equation} Let us draw an imaginary box of size much larger then a particle diameter. Several particles have to move in a cooperative way to increase the number of the particles in the box by one. This number can be estimated by simply counting the amount of free volume. Assuming that the $N$ particles are rearranged in such a way that they contribute their entire free volume to create a void large enough to accommodate a particle, $NV_0=V$, or \begin{equation} \label{N} N={\rho\over 1-\rho}. \end{equation} Indeed, as the density vanishes no rearrangements are necessary, and when the density approaches its maximal value the number of rearrangements diverges. We further assume that the motion of particles is not correlated and therefore, the time associated with such a rearrangement should increase exponentially with $N$, $T\sim e^N$. Consequently, at large times the density increases according to the following rate equation, \begin{equation} \label{drhot} {d\rho\over dt}\propto (1-\rho){1\over T}=(1-\rho)e^{-\rho/(1-\rho)}. \end{equation} The rate at which the density increases is proportional to the void volume and inversely proportional to the rearrangement time. The latter exponential factor effectively reduces the density increase rate and it dominates as $\rho\to 1$. The solution of this equation is given asymptotically by \begin{equation} \label{rhot} \rho(t)\cong \rho_{\infty}-{1\over \ln t}, \end{equation} with $\rho_{\infty}=1$. This heuristic argument ignores the structure of the granular assembly and the fact that the density cannot exceed the close-packed value. In the physical case of three dimensions, the available free volume for rearrangements should vanish as the density approaches the close packing limiting value, thereby suggesting a more realistic form for $V_0$ such as $V(\rho_{\rm max}-\rho)/\rho$. Despite the simplifying assumptions, this heuristic picture is useful as it highlights the basic mechanism underlying granular compaction, i.e., the diverging number of rearrangements needed for the density to increase. In one-dimension, no geometrical complications occur and the exponential rate reduction factor is in fact exact. In the following, we present an analytically tractable model and discuss its relevance to density relaxation and steady state density fluctuations. Consider a stochastic adsorption-desorption process on a continuous one dimensional substrate. Identical particles of unit diameter adsorb uniformly from the bulk to a substrate with rate $k_+$ and desorb with rate $k_-$. In other words, $k_+$ adsorption attempts are made per unit time per unit length, and similarly, the probability that an adsorbed particle desorbs in an infinitesimal time interval $(t,t+dt)$ is $k_-dt$. While the desorption process is unrestricted, the adsorption process is subject to excluded volume constraints, {\it i.e.}, particles can not adsorb on top of previously adsorbed particles. The attempted adsorption event in Fig.~2 is thus rejected. This ``car-parking'' process was previously studied in the context of chemisorption\cite{Evans} and protein binding\cite{McGhee}. \begin{figure}\vspace{.2in} %\centerline{\epsfxsize=7.5cm \epsfbox{fig2.ps}} \centerline{\includegraphics[width=7.5cm]{fig2}} \vspace{.2in} \label{fig2} Fig 2 {\small The adsorption-desorption process.} \end{figure} Let us consider $P(x)$, the distribution of voids of size equal to $x$ between particles. This distribution satisfies two normalization conditions : $\rho=\int dx P(x)$ and $1=\int dx (x+1)P(x)$. The first condition states that there are as many voids as particles, while the second reflects conservation of the total length. Ignoring correlations between neighboring voids, this distribution evolves according to \begin{eqnarray} \label{pxt} &&{\partial P(x)\over \partial t}=2k_+\int_{x+1} dy P(y)-2k_-P(x)\\ &&+\theta(x-1)\!\left[{k_-\over\rho(t)}\int_0^{x-1}dy P(y)P(x-1-y)-k_+(x-1)P(x)\right].\nonumber \end{eqnarray} The first term represents gain due to adsorption and the second loss due to desorption. The last two terms apply only for voids larger than a particle. The above process satisfies detailed balance, and consequently, the system approaches its equilibrium state after waiting sufficiently long time. In this state, neighboring voids are uncorrelated and the above master equation holds. The equilibrium void distribution is exponential \cite{Krapivsky} \begin{equation} \label{px} P_{\rm \infty}(x)={\rho_{\infty}^2\over (1-\rho_{\infty})} \exp\left[-{\rho_{\infty}\over 1-\rho_{\infty}}x\right]. \end{equation} Using this equilibrium distribution, one can obtain the exact equilibrium density. Furthermore using a quasistatic (near-equilibrium) approximation time dependent properties can be studied as well. The density evolves according to the following rate equation \begin{equation} {d\rho\over dt}=-k_-\rho+k_+(1-\rho)e^{-{\rho/(1-\rho)}}. \label{rateq} \end{equation} Since desorption is unrestricted, the loss term is linear in the density. The second term is obtained by evaluating the total adsorption rate $\int_1^{\infty} dx(x-1)P(x)$ and substituting for $P(x)$ the equilibrium distribution of Eq.~(\ref{px}). Only voids larger than a particle contribute to adsorption, and there is a correction factor $(x-1)$ that vanishes at the minimal gap. This quasistatic approximation assumes that the system evolves slowly and thus is close to equilibrium. The total adsorption rate can be understood as follows: the factor $(1-\rho)$ accounts for the total void density, while the factor $s(\rho)=e^{-\rho/(1-\rho)}$ is the sticking probability, i.e., the probability that an adsorption event is successful. This factor is roughly unity in the dilute limit, and vanishes exponentially as $\rho\to 1$. The excluded volume interaction effectively reduces the adsorption rate, $k_+\to k_+(\rho)= k_+s(\rho)$. Interestingly, the sticking probability is identical to the one of Eq.~(3), and therefore the heuristic picture is shown to be exact in one dimension. We also note that the density increases with rate proportional to the density of voids the size of a particle, and it is quite possible that this holds in higher dimensions as well. Using Eq.~(\ref{rateq}), the {\it exact} equilibrium density is obtained from the following transcendental equation, $\alpha e^{\alpha}=k$, with $k\equiv k_+/k_-$ and $\alpha=\rho_{\infty}/(1-\rho_{\infty})$. The ratio of adsorption to desorption determines the equilibrium density in the model. Using Eq.~(\ref{rateq}), the equilibrium density is found. The leading behavior in the two limiting cases is \begin{equation} \rho_{\infty}(k)\cong\cases{k&$k\ll1$;\cr1-(\ln k)^{-1}&$k\gg1$.\cr} \label{rhoeq} \end{equation} While the behavior in the dilute limit is linear, the approach to the close packed state is very slow. The effect of the volume exclusion constraint is striking, a huge adsorption to desorption rate ratio, $k\cong10^9$, is necessary to achieve a $0.95$ steady-state occupancy. One can associate the parameter $k$ in the model with the vibration intensity in experiments, where a monotonic correspondence between the steady-state density and the acceleration $\Gamma$ \cite{Nowak1}. We now focus on the relaxation properties of the system. The granular compaction process corresponds to the high density limit, and we thus focus on the desorption-controlled case, $k\gg1$. Hence, we fix $k_+=1$ and consider the limit $k_-\to 0$ of Eq.~(\ref{rateq}). The early time behavior is dominated by adsorption and can be obtained by neglecting the desorption term. The system approaches complete coverage according to Eq.~(\ref{rhot}). This behavior was confirmed by numerical simulations\cite{Krapivsky,Jin}. It also occurs in a situation where diffusion plays the role of desorption\cite{Privman}. Use of Eq.~(\ref{rateq}) is justified {\it a posteriori} since the system evolves slowly and has enough time to equilibrate. The inverse logarithmic behavior is simply a reflection of the exponentially suppressed adsorption in the dense limit. Eq.~(\ref{rhot}) holds indefinitely only for the truly irreversible limit of the parking process, {\it i.e.}, for $k=\infty$. For large but finite rate ratios, the final density is given by Eq.~(\ref{rhoeq}). As the density approaches this steady state value, the loss term becomes significant and should be taken into account. The crossover time between the two different relaxation regimes, $t_0$, can be conveniently estimated by equating the time dependent density of Eq.~(\ref{rhot}) with the equilibrium density of Eq.(\ref{rhoeq}) $1-1/\ln t_0=1-1/\ln k$, and as a result $t_0\cong 1/k_-$. For $t\gg t_0$, the loss term is no longer negligible. By computing how a small perturbation from the steady state decays with time, an exponential relaxation towards the steady state is found $|\rho_{\infty}-\rho(t)|\propto e^{-t/\tau}$ for $t\gg t_0$. The relaxation time is related to $t_0$, however, an additional logarithmic correction occurs, $\tau=t_0(1-\rho_{\infty})^2\simeq t_0/(\ln k)^2$. In summary, the early time behavior of the system follows the irreversible limit of $k_-=0$. Once the system is sufficiently close to the steady-state, the density relaxes exponentially to its final value. The experimentally observed relaxation curves which corresponds to large compaction are indistinguishable over the observed time range\cite{Knight}. Also, steady state is achieved in a finite time $<10^6$ taps. Both features are consistent with our theory. The inverse logarithmic relaxation can be viewed as a sum of many exponential functions with growing decay times. In a finite system, however, a maximal relaxation scale should eventually dominate and the relaxation law of Eq.~(\ref{rhot}) should apply only up to this time scale. The parking process can be used to study density fluctuations in the steady state in finite systems. In a finite system, after the system relaxed to the steady state, the observed density can deviate from its average expected value. Since the equilibrium properties of the parking model are understood, the distribution of these deviations can be calculated. A useful quantity is $G(x_1,x_2\ldots, x_n)$, the probability of finding $n$ consecutive voids of sizes $x_1$, $x_2$, etc. In equilibrium voids are uncorrelated and this distribution is given by \begin{eqnarray} \label{gx} G_{\infty}(x_1,x_2\ldots, x_n)&=&\rho_{\infty}^{1-n}P_{\infty}(x_1) P_{\infty}(x_2)\cdots P_{\infty}(x_n)\nonumber \\ &=& \rho_{\infty} \left({\rho_{\infty}\over 1-\rho_{_\infty}}\right)^n \exp\left[-{\rho_{\infty} \over 1-\rho_{\infty}}V\right] \end{eqnarray} where $V=\sum x_i$ is the total void space of that configuration. The conditional probabilities $\rho_{\infty}^{-1}P_{\infty}(x_2)$ {\it etc.} has to be used to ensure proper normalization $\int dx_1\cdots\int dx_n G_{\infty}(x_1,x_2\ldots, x_n)=\rho_{\infty}$ (the normalization $\int dx P(x)=\rho_{\infty}$ is used here). The multiple void distribution depends exponentially on the total void space $V$, and a factor of $\rho{\infty}/(1-\rho_{\infty})$ is generated by each void. To calculate $P(\rho)$, the probability that the observed density in a system of size $L$ is $\rho$ we note that the total void length in such a case is $V=(1-\rho)L$ and the total number of cars is $n=\rho L$. Every configuration with $n$ voids $x_1$,$\ldots$, $x_n$, and $V=\sum x_i$ contributes to this probability and therefore \begin{equation} P(\rho)=\sum_{i=1}^n\int dx_i G_{\infty}(x_1,x_2\ldots, x_n) \end{equation} Evaluating the multiple integral using the steepest decent method we find a Gaussian density distribution \begin{equation} P(\rho)={1\over\sqrt{2\pi\sigma^2}} \exp\left[-{(\rho-\rho_{\infty})^2\over 2\sigma^2}\right], \end{equation} with the variance $\sigma^2=\rho_{\infty}(1-\rho_{\infty})^2/L$. Of course, in the limits of infinite systems, the variance vanishes. The fluctuations width, $\sigma$, and the relaxation time scale, $\tau$, are related via $L\sigma^2=k_-\rho_{\infty}\tau$, an analog of the fluctuation dissipation theorem. Surprisingly, this prediction is in agreement with the experimental observations. Most of the observed distributions of the density fluctuations are Gaussian. There are exceptions, however, especially near the bottom of the column, where a high density configuration is slightly preferred and the positive tail of the distribution is enhanced. This suggests that the system gets locked in a long lived high density metastable state. It is possible that averaging over longer time records will restore the Gaussian nature. The experimental variance also decreases with increasing density, as is the case for the theoretical variance that vanishes as $\rho_{\infty}\to 1$. Additionally, one can study the spectrum of the steady state density fluctuations. It is remarkable that despite the fact that the model is one dimensional, its density fluctuations power spectrum is very similar the experimental curves\cite{Nowak1}. In a realistic granular material, an individual particle can not penetrate its neighbors, and it is in contact with several other particles. Our model properly accounts for the hard core repulsion, but it ignores mechanical stability. We argue that in the long time limit mechanical stability can not play a significant role in determining the motion of individual grains during the compaction process. Instead, the motion is limited primarily by the presence of other beads. The tradeoff, however, is that such a simple theory can not make predictions about the final density. Nevertheless, it does elucidate the leading mechanism in granular compaction. It also motivates studying the form of the void distribution and specifically its large volume tail in higher dimensions . In conclusion, we have studied theoretically density relaxation towards steady state in granular compaction using heuristic arguments and a {\it microscopic} model in one dimension. Due to volume exclusion, exponentially growing time scales are associated with cooperative motion of grains. As a result the approach towards the steady state is an inverse logarithmic one. Since the argument leading to the logarithmic relaxation is quite general, the results should hold in a large class of physical situations such as horizontal shaking, aspherical particles, or even polydisperse distribution. Additionally, our toy model exhibits density fluctuations that agree remarkably well with the experiment. This provides an additional confirmation for the validity of this model. We thank L.~P.~Kadanoff, P.~L.~Krapivsky, A. Mehta, and T.~A.~Witten for useful discussions. \vspace{-.2in} \begin{thebibliography}{99} \bibitem{Jaeger}H.~Jaeger and S.~Nagel, R.~P.~Behringer, {\sl Rev. Mod. Phys.} {\bf 68}, 1259 (1996). \bibitem{Evesque}P.~Evesque and J.~Rajchenbach, {\sl Phys. Rev. Lett.} {\bf 62}, 44 (1989). \bibitem{Melo}F.~Melo, P.~B.~Umbanhowar, and H.~L.~Swinney, {\sl Phys. Rev. Lett.} {\bf 75}, 21 (1995). \bibitem{Ehrichs}E.~E.~Ehrichs, H.~M.~Jaeger, G. S. Karczmar, J. B. Knight, V. Y. Kuperman, and S.~R.~Nagel, {\sl Science} {\bf 267} 1632, (1995). \bibitem{Knight}J.~B.~Knight, C.~G.~Fandrich, C.~N.~Lau, H.~M.~Jaeger, and S.~R.~Nagel, {\sl Phys Rev. E} {\bf 51} 3957 (1995). \bibitem{Nowak1}E.~R.~Nowak, J.~B.~Knight, E.~Ben-Naim, H.~M.~Jaeger, and S.~R.~Nagel, {\sl Phys. Rev. E}, in press. \bibitem{Nowak2}E.~R.~Nowak, J.~B.~Knight, M.~Povinelli, H.~M.~Jaeger, and S.~R.~Nagel, {\sl Powder Technology}, in press. \bibitem{Edwards}S.~F.~Edwards in {\sl Granular Matter, An Interdisciplinary Approach}, edited by A. Mehta (Springer-Verlag, New York, 1994). \bibitem{Mehta}A.~Mehta in {\sl Granular Matter, An Interdisciplinary Approach}, edited by A. Mehta (Springer-Verlag, New York, 1994). \bibitem{Mehta1}A.~Mehta and S.~F.~Edwards, {\sl Physica A} {\bf 157}, 1093 (1989). \bibitem{Barker}G.~C.~Barker and Anita Mehta, {\sl Phys. Rev. E} {\bf 47}, 184 (1993). \bibitem{Hong}D.~C.~Hong, S.~Yue, J.~K.~Rudra, M.~Y.~Choi, and Y.~W.~Kim, {\sl Phys. Rev. E} {\bf 50}, 4123 (1994). \bibitem{Linz} S.~J.~Linz, {\sl Phys. Rev. E} {\bf 54}, 2925 (1996). \bibitem{Nicodemi}M.~Nicodemi, A.~Coniglio, and H.~Herrmann, {\sl Phys. Rev. E} {\bf 55}, 3962 (1997). \bibitem{Evans}J.~W.~Evans, {\sl Rev. Mod. Phys.} {\bf 65}, 1281 (1993). \bibitem{McGhee}J.~D.~McGhee and P.~H.~von Hippel, {\sl J. Mol. Bio} {\bf 86}, 469 (1974). \bibitem{Krapivsky}P.~L.~Krapivsky and E.~Ben-Naim, {\sl Jour. Chem. Phys.} {\bf 100}, 6778 (1994). \bibitem{Jin}X. Jin, G. Tarjus, and J. Talbot {\sl J. Phys. A} {\bf 27} L195 (1994). \bibitem{Privman}V.~Privman and M.~Barma, {\sl Jour. Chem. Phys.} {\bf 97}, 6714 (1992). \end{thebibliography} \end{multicols} %\centerline{FIGURE CAPTIONS} %{\bf Fig. 1} A rearrangement of a growing number of particles is necessary to %enable a local density increase. %{\bf Fig. 2} The adsorption-desorption process. %\centerline{\epsfxsize=5cm \epsfbox{ill1.ps}} %\centerline{\epsfxsize=5cm \epsfbox{ill2.ps}} \end{document}