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






