\documentstyle[aps,epsf,multicol,graphicx]{revtex}
\begin{document}
\title{Density fluctuations in vibrated granular materials}
\author{Edmund R.~Nowak, James B.~Knight, Eli Ben-Naim, † Heinrich
M.~Jaeger, and Sidney R.~Nagel}
\address{The James Franck Institute and the Department of Physics, The University of
Chicago, Chicago, Illinois 60637}
\maketitle
\begin{abstract} 
We report systematic measurements of the density of a vibrated
granular material as a function of time. Monodisperse spherical beads
were confined to a cylindrical container and shaken vertically. Under
vibrations, the density of the pile slowly reaches a final
steady-state value about which the density fluctuates. We have
investigated the frequency dependence and amplitude of these
fluctuations as a function of vibration intensity $\Gamma$. The
spectrum of density fluctuations around the steady state value
provides a probe of the internal relaxation dynamics of the system and
a link to recent thermodynamic theories for the settling of granular
material. In particular, we propose a method to evaluate the
compactivity of a powder, first put forth by Edwards and co-workers,
that is the analog to temperature for a quasistatic powder. We also
propose a stochastic model based on free volume considerations that
captures the essential mechanism underlying the slow relaxation. We
compare our experimental results with simulations of a one-dimensional
model for random adsorption and desorption.  
\\
\vspace{5pt} {PACS:} 81.05.Rm, 05.40.1j, 46.10.1z, 81.20.Ev
\end{abstract}
\begin{multicols}{2}
\section{Introduction}

One of the salient features of noncohesive granular materials is that
they can be packed over a range of densities and still retain their
resistance to shear. For example, a stable conglomeration of
monodisperse spheres can exist with a packing fraction r ranging from
$\rho\approx0.55$ (the random loose packed limit) to $\rho\approx0.64$
(the random close packed limit) and even to $\rho\approx0.74$ (the
crystalline state). Because thermal energies, $k_B T$, are
insignificant when compared to the energy it takes to rearrange a
single particle, each metastable configuration will persist
indefinitely until an external vibration comes along to knock it into
another state. Thus, no thermal averaging takes place to equilibrate
the system. The density of the material is determined both by its
initial preparation and by the manner in which it was handled or
processed, since such activities normally introduce some vibrations
into the material. The phase space for the granular medium is explored
not by fluctuations induced by ordinary temperature but by
fluctuations induced by external noise sources, such as vibrations. It
is the goal of this paper to provide an experimental foundation for
the use of such fluctuations as a probe of the dynamics as well as the
microstructure of granular media in the quasistatic, densely packed
limit.

Granular compaction involves the evolution from an initial low-density
packing state to one with higher final density and provides a model
system for nonthermal relaxation in a disordered medium. In a previous
study [1], we focused on the approach to a final steady-state density
as vibrations were applied to the system. In particular, we studied
the density of monodisperse spherical particles in a tall cylindrical
tube as a series of external excitations, consisting of discrete,
vertical shakes or ``taps,’’ were applied to the container. Such data
indicate that the compaction process is exceedingly slow: the density
approaches its final steady-state value approximately logarithmically
in the number of taps. A typical example of such behavior, in Fig. 1,
shows that in excess of $10^4$ taps may be required before the density
has relaxed to its steady-state value. However, if one vibrates for a
long enough time a steady-state density, depending on the intensity of
the taps, will be attained. Even after the density reaches the
steady-state value, one can discern fluctuations in the density about
that value: after each ``tap’’, the density will be slightly higher or
lower than it was before.  These fluctuations are reminiscent of
thermal fluctuations about an equilibrium state, yet such a connection
so far has not been investigated experimentally.

In statistical mechanics the study of fluctuations is of great
physical interest. The fluctuation-dissipation theorem relates the
dissipative response of a system to an external perturbation with the
microscopic dynamics of the system in a state of equilibrium.  Energy
fluctuations in thermal systems can be used to investigate the set of
distinct, microscopic states that are accessible to a system
maintained at a fixed temperature. Likewise, a study of density
fluctuations in granular media may provide a framework for
understanding the physical phenomenon of compaction, i.e., how a
vibrated powder, that is not in a steady state, finally approaches a
steady state.

In a granular system, density fluctuations from the steady state
represent the different volume configurations accessible to particles
subject to an external vibration. It is desirable to develop an
analogy between the role that vibrations play in nonthermal systems,
such as granular media, and the role of temperature in thermal
systems. Theoretically, this issue was addressed by Edwards and
co-workers [2-4] who introduced a statistical mechanics for
powders. The idea is based on the assertion that an analogy can be
drawn between the energy of a thermal system and the volume $V$
occupied by a powder. The entropy S of a powder is defined in the
usual sense, by the logarithm of the number of available
configurations. Edwards and co-workers then put forth the concept of
an effective temperature for a powder, called the compactivity $X$,
which is defined as $X\equiv\partial V/\partial S$. The significance
of this effective temperature is that it allows for the
characterization of a static granular system. This is distinct from
the case of rapid granular flows where a ``granular temperature’’
given by the mean-square value of the fluctuating component of the
particle velocities can be written down [5-7]. The compactivity is
then a measure of ``fluffiness’’ in the powder: when $X=0$, the powder
is in its most compact configuration, whereas for $X=\infty$ the
powder is the least dense.

Recently, another approach [8-10] that describes the static packing of
powders has adapted a statistical model that contains geometric
frustration as an essential ingredient. For granular materials,
frustration arises in the form of hard-core repulsive constraints and
the interlocking of grains of different shapes, which prevents local
rearrangements. Both the static and dynamic (in the presence of
vibration and gravity) properties of this model exhibit complex
behavior with features that are common to granular packing, such as
the logarithmic relaxation of density under tapping [1].

In this paper, we make contact with these ideas through a detailed
study of the process of granular compaction. In particular, we propose
a method for evaluating the compactivity of a vibrated powder through
a definition of a ``granular specific heat’’ and measurements of
density fluctuations observed in the reversible regime of steady-state
behavior. We also elaborate on a theoretical model [11,12], based on
the idea that the rate of increase in volume density is exponentially
reduced by the free volume, which captures many of the significant
features of our experiments. A model addressing the compaction of
binary mixtures consisting of grains with very different sizes was
recently proposed by de Gennes [13]. That model is similar to ours in
that it incorporates free volume constraints and also exhibits a
similar inverse logarithmic dependence for the density relaxation.

In the next section we will describe the experimental de-tails of the
system, review how to obtain reproducible and reversible densities,
and present our results for the density fluctuations.  In Sec.~III we
discuss several models in relation to our experimental results and
motivate the relevance of free volume constraints for granular
compaction. In Sec.~IV, we present the theoretical model and the
results of related simulations of compaction. Finally, in the last
section we discuss the central result of this paper, namely, how our
data can be related to thermodynamic approaches for understanding
granular media.

\section{EXPERIMENTAL RESULTS}

\subsection*{Experimental method}

The details of the experimental apparatus and measurement technique
were published elsewhere [1]. Monodisperse, spherical soda-lime glass
beads (of $2$ mm diameter) were confined to a $1.88$ cm diameter Pyrex
tube measuring 1 m in height. The tube was subjected to discrete
vertical shakes or taps’’! each consisting of one complete cycle of a
$30$ Hz sine wave. The vibration intensity was parametrized by
$\Gamma$, which is the ratio of the peak acceleration $A$ that occurs
during a single tap to the gravitational acceleration $g=9.8\ {\rm
m}/{\rm s}^2$: $\Gamma=A/g$.  The beads were baked prior to loading in
the tube and precautions were taken to minimize complications
resulting from electrostatic charging, convection, and external
humidity fluctuations. The column of beads was prepared in a low
density initial state by flowing high pressure, dry nitrogen gas from
the bottom to the top of the tube. The top layer of the beads was free
to move, i.e., there was no load or dead-weight surcharge applied to
the column of beads. The density, or equivalently the packing fraction
$\rho$, which is the percentage of volume occupied by the beads, was
determined either by a measurement of the total height of the beads
within the tube or using capacitors that were mounted on the outside
wall of the tube. For the latter, the capacitance was found to vary
linearly with packing fraction. Each capacitor averaged the density
over sections containing approximately $6000$ beads. Measurements of
$\rho$ were taken as a function of time, i.e., number of taps $t$ and
as a function of the intensity of the vibrations,
$\Gamma$. Corrections for instrumental drift were made by using
simultaneously acquired data from a second, stationary tube
(identically prepared with the same type of beads and connected to the
same vacuum system). Our instrumentation allowed shaking intensities
up to $\Gamma\approx 7$ and provided a resolution $\Delta\rho=0.0006$
in measured packing fraction changes.

The desired outcome of a shake cycle is to provide clearly defined
periods of uniform dilation of the bead assembly. During these periods
of dilation the beads have some freedom to rearrange their positions
relative to their neighbors and thereby replace one stable
close-packed configuration by another. Previously [1,14], we have
shown that the overall behavior of the compaction process is
qualitatively similar at different depths into the container (see also
Fig.~(1). Spurious effects from continuous vibrations, such as period
doubling or surface waves [12], were avoided by spacing the taps
sufficiently far apart in time to allow the system to come to complete
rest between taps. Also, by using a tall container with smooth,
low-friction interior walls shear-induced dilation and granular
convection were suppressed [15]. Although friction between beads and
with the tube walls can affect the mechanical stability of a bead
configuration, we argue below that the motion of beads is limited
primarily by geometric constraints imposed by the presence of other
beads, particularly at the high densities investigated here.



The ratio of the container diameter to the bead diameter can also
influence the compaction process. For small values of this ratio,
ordering (crystallization) induced by the container walls [16] will
increase the measured packing fraction over its bulk value, leading to
densities that can exceed the random close-packed limit. This may be
responsible for the high maximum packing fractions seen in
Fig.~(2). Previous studies [1,14] indicate that the qualitative
behavior of the compaction process is similar for varying bead
sizes. The container walls can also place constraints on the density
fluctuations. Since it is our aim to investigate these density
fluctuations, the choice of bead size was a compromise between
maximizing the container-to-bead diameter ratio and not having the
amplitude of the density fluctuations be obscured by statistical
averaging over a large number of particles.
\begin{figure}
\centerline{\epsfxsize=9cm\epsfbox{fig1.eps}}
\caption{The time evolution of the volume density $\rho$ at three different
depths near the top, middle, and bottom of the pile of beads. The
curves represent a single run (no ensemble averaging!  at a vibration
intensity $\Gamma=6.8$. The pile settles slowly from its initial low
density configuration toward a higher steady-state density at long
times, $t>10^4$ taps. The dashed lines are fits to Eq.~(1) with
typical values of parameters: $0.637<\rho_\infty <0.647$,
$0.036<\Delta\rho_\infty<0.044$, $0.20<B<0.40$, $10<\tau<18$.}
\end{figure}

\subsection*{Reaching the steady state}

At a high acceleration $\Gamma$ the steady-state density, $\rho_{ss}$ can be
approached by simply applying a very large number of taps (often
greater than $10^4-10^5$). An example is shown in Fig.~(1) for
$\Gamma=6.8$. The three panels correspond to the capacitively measured
density near the top, middle, and bottom sections of the pile of
beads. (The tap number $t$ is offset by $1$ tap so that the initial
density can be included on the logarithmic axis.) Note that these
curves represent a single run, and separate runs starting from the
same initial density differ in the details of the density fluctuations
but show a similar overall behavior. The behavior of $\rho(t)$,
obtained by averaging many runs of this kind, appears to be
homogeneous throughout the pile at these high accelerations. As
discussed in Ref.~[1], the time evolution of this ensemble averaged
density is well fitted by the expression
\begin{equation}
\rho(t)= \rho_\infty-{\Delta\rho_\infty\over 1+B\ln (1+t/\tau)}
\end{equation}
where the parameters $\rho_\infty$ , $\Delta\rho_\infty$, $B$, and
$\tau$ depend only on the acceleration $\Gamma$. Equation (1) was
found to fit the ensemble averaged density over the whole range
$0<\Gamma<7$ better than other functional forms that were tried i.e.,
exponential, stretched exponential, or algebraic forms, see
Ref.~[1]. The dashed lines in Fig.~(1) show a fit to Eq.~(1).  Here,
the value of the final density, $\rho$ , is approximately equal to the
observed steady-state density $\rho_{ss}$ .


For small values of $\Gamma$, however, $\rho_\infty$ corresponds to a
metastable state and not the steady-state density. In particular, for
values of the applied acceleration $\Gamma<3$, it is difficult, if not
experimentally impossible, to reach the steady-state by merely
applying a sufficiently large number of taps of identical
intensity. In this case, the steady state can be reached by
``annealing’’ [14] the system. The annealing is controlled by the ramp
rate, $\Delta \Gamma/\Delta t$, at which the vibration intensity is
varied over time. Experimentally, we slowly raise the value of
$\Gamma$ from $0$ to a value beyond $\Gamma^*$ in increments of
$\Delta \Gamma \approx 0.5$. At each intermediate value of $\Gamma$ we
apply $\Delta t=10^5$ taps.  $\Gamma^*$ defines an ``irreversibility
point’’ in the sense that, once it has been exceeded, subsequent
increases as well as decreases in $\Gamma$ at a sufficiently slow rate
$\Delta\Gamma/\Delta t$ lead to reversible, steady-state behavior.  We
found that $\Gamma^*\approx 3$ for $1$, $2$, and $3$ mm beads [14]. A
typical run is shown in Fig.~2. Here we have used $2$ mm beads, and
started with an initial density of $\rho\approx 0.59$.  The highest
densities are achieved by annealing the system, i.e., decreasing
$\Gamma$ slowly from $\Gamma^*$ back down to $\Gamma=50$. If $\Gamma$
is rapidly reduced to $0$ (large $\Delta\Gamma/\Delta t$) then the
system falls out of ``equilibrium.’’ This leads to lower final
densities and a curve for $\rho(\Gamma)$ that is not reversible. A
crucial result emerging from data such as in Fig.~2 is that along the
reversible branch, the density is monotonically related to the
acceleration.  We note that in 3D simulations of granular compaction
by Mehta and Barker [17] a similar monotonic decrease in steady-state
volume fraction as a function of shaking intensity was found. Thus,
only once the steady-state has been reached is there a single-valued
correspondence between the average density and the applied 
acceleration.
\begin{figure}
\centerline{\epsfxsize=9.0cm\epsfbox{fig2.eps}}
\vspace{-1in}
\caption{The dependence of $\rho$ on the vibration history. The
beads were prepared in a low density initial configuration and
then the acceleration amplitude $\Gamma$ was slowly first
increased (solid symbols) and then decreased (open symbols). At
each value of $\Gamma$ the system was tapped $10^5$ times after
which the density was recorded and $\Gamma$ was subsequently
incremented by $\Delta\Gamma\cong 0.5$. The upper branch that has
the higher density is reversible to changes in $\Gamma$, see 
square symbols. The $\Gamma^*$ denotes the irreversibility point.}
\end{figure}

\subsection*{Density fluctuations about the steady state}


After the granular material has been vibrated for a sufficiently long
time, it reaches a steady-state density $\rho_{ss}$.  Although there
is a well-defined average density, Fig.~1 already hints that there are
large fluctuations about this value. The magnitude of the fluctuations
depends on the vibration intensity and depth within the
container. Figure 3 shows in more detail an example of these
fluctuations as a function of time,
$\delta\rho(t)=\rho(t)-\rho_{ss}$. In Fig. 3(a) we plot
$\delta\rho(t)$ for a fixed value of acceleration, $\Gamma=5.9$, but
measured at different depths in the container. Note that the rate at
which the density varies in time decreases with depth into the
pile. That is, the top of the pile has more high frequency noise than
the bottom. The curve marked ``reference’’ is the reference capacitor
to which no vibrations are applied. This last curve is essential to
compensate for drifts that could occur in the electronics over the
very long period of our measurements. Each record shown here is 4096
taps long and up to 132 successive such records were assembled to
produce one very long time sequence. Figure 3(b) shows the
fluctuations in the density measured at the bottom capacitor as a
function of acceleration $\Gamma$. As $\Gamma$ is increased both the
magnitude of the fluctuations and the amount of high-frequency noise
increase.

From data as in Fig.~3 we can obtain the shape of the distribution
function for the fluctuation amplitudes. We plot in Fig.~4 the
logarithm of the relative probability of occurrence $D(\delta\rho)$
versus $\Psi^2=(\rho-\rho_{ss})^2{\rm sgn}(\rho-\rho_{ss})$ so that
a Gaussian random process will have a triangular shape. In that figure
we plot $D(\delta\rho)$ for the entire range of accelerations,
$4<\Gamma<7$ for which fluctuations could be reliably measured with
our equipment. All data records were corrected for instrumental drifts
using the reference capacitor. As can be seen in Fig. 4, the majority
of data shows Gaussian character. For a small fraction of runs e.g.,
($\Gamma=5.9$), however, we find significant deviations from Gaussian
behavior, particularly near the middle and bottom of the pile. When
such deviations are present they tend to preferentially occur for
positive values of $\Psi^2$ , i.e., higher densities. The deviations
could be due to a metastable state, away from the mean, in which the
system gets trapped.  Fluctuations about this metastable state may
even be distributed in a Gaussian fashion. The reason why such
metastable states favor the lower portion of the column and why they
are prominent at certain values of $\Gamma$ is unclear.
\begin{figure}
\centerline{\epsfxsize=8.5cm\epsfbox{fig3.eps}}
%\vspace{-.0in}
\caption{Fluctuations in the volume density $\delta\rho(t)=\rho(t)-
\rho_{ss}$ after the system has had sufficient time to relax to a
steady-state density $\rho_{ss}$. In (a) the fluctuations at three
different depths are shown for $\Gamma=5.9$. The reference capacitor
is used to correct for any instrumental drift. The dependence of the
fluctuations on $\Gamma$ is shown in (b) for the beads near the bottom
of the pile. Fluctuations over a broad range of time scales are
evident.}
\end{figure}


\end{multicols}
\begin{figure}[t]
\centerline{\epsfxsize=13.0cm\epsfbox{fig4.eps}}
\vspace{-.5in}
\caption{The distribution functions $D(\delta\rho)$ for the occurrence
of fluctuation amplitudes in the steady state are shown (solid
circles) for the three depths at various $G$. Plotted as a function of
$\Psi^2=\delta\rho^2 {\rm sgn}(\delta\rho)$ a Gaussian distribution
has a triangular shape. For selected panels, the time dependence of
the distributions is shown by plotting the distribution functions for
only the first (open squares) or second (open triangles) half of the
time record. The majority of data appears stationary even when
significant non Gaussian deviations are observed, e.g., at
$\Gamma=5.9$ near the middle and bottom of the pile.}
\end{figure}
\begin{multicols}{2}

We can qualitatively check whether the distribution functions
correspond to a stationary random process or whether they conceal a
slow drift away from an originally well-defined mean density.
Strictly speaking, a stationary Gaussian process is one for which
correlation functions of order higher than second are zero, see
Ref.~[18]).  This is done by dividing each time record into two equal
length halves and then determining the distribution functions for each
half separately, as shown for selected values of $\Gamma$ and depths
by the open symbols in Fig.~4. We find that in practically all cases
the fluctuations do appear to be stationary and, moreover, that in the
very few nonstationary cases observed, the Gaussian character is
recovered at {\it later} times (i.e., in the second half of the record).

By assembling 132 successive time traces of the type shown in Fig.~3,
we can obtain continuous time records containing 540 672 data
points. From such records we calculate both the density
autocorrelation function and the power spectrum for the density
fluctuations, $S_\rho(\omega)$, where the frequency v is measured in
units of inverse taps. In Fig.~5 we plot $S_\rho(\omega)$ versus
$\omega$ for the three depths at various values of acceleration,
$\Gamma=4.3$, $5.1$, $5.9$, and $6.8$. We note several distinctive
features to these power spectra. In particular, three characteristic
regimes emerge: (i) a white noise regime,
$S_\rho(\omega)\propto\omega^0$ below a low-frequency corner
$\omega_L$ , (ii) an intermediate-frequency regime with nontrivial
power-law behavior, and (iii) a simple roll-off $S_\rho(\omega)\propto
\omega^{-2}$ above a high-frequency corner, $\omega_H$. This
classification appears to apply to all traces shown in Fig.~5. It is
most pronounced for the spectrum in the lower right hand panel. For
spectra where $\omega_L$ and $\omega_H$ are sufficiently separated in
frequency, the data show that the spectral dependence between
$\omega_L$ and $\omega_H$ cannot be approximated by just a simple
superposition of two separate Lorentzians each having a frequency
dependence $S\propto \tau/(1+\omega^2\tau^2)$ but different
characteristic times $\tau$. A comprehensive analysis of this region
reveals that the most consistent description for all the data is
obtained with a Lorentzian tail, $S_\rho(\omega)\sim \omega^{-2}$ just
above $\omega_L$ , followed by a region with $S_\rho(\omega)\propto
\omega^{-\delta}$ (with $\delta\approx 0.9\pm 0.2$) stretching up to
$\omega_H$ , the high-frequency corner.

\end{multicols}
\begin{figure}[t]
\vspace{-.3in}
\centerline{\epsfxsize=13cm\epsfbox{fig5.eps}}
\vspace{-.3in}
\caption{The power spectral density $S_\rho$ (taps) of the
fluctuations as a function of frequency $\omega$ in the steady state
is shown for the three depths at various $\Gamma$. For most spectra,
two characteristic corner frequencies, $\omega_L$ and $\omega_H$ , are
discernible which shift to higher frequencies for increasing $\Gamma$
and decreasing depth. The characteristic regimes of behavior are
denoted by the dashed lines in the lower right hand panel, which are
guides to the eye.}
\end{figure}
\begin{multicols}{2}


One result from the data in Fig.~5 is the dependence of both corner
frequencies on the acceleration $\Gamma$. To determine these
frequencies we used a combination of two methods, which we illustrate
here for the simple case of a Lorentzian spectrum.  First, for any
Lorentzian, the product $\omega S_\rho$ has a maximum precisely at
$\omega=1/\tau$ so that $\omega_L$ and $\omega_H$ can be associated
with the frequencies at which $\omega S_\rho$ exhibits peaks. Second,
even though we were using extremely long time records they are still
of finite length.  Figure 5 clearly indicates cases where $\omega_L$
is difficult to obtain because of the large statistical variance
($\approx 25\%$) in $S_\rho$ throughout the lowest decade in
frequencies. In these instances we employed the additional information
contained in the one-sided sine transform of the density-density
autocorrelation function. For example, for a single Lorentzian for
which the autocorrelation function is simply $\approx e^{-t/\tau}$ ,
the ratio of sine to cosine transform of the autocorrelation function
is given by $\omega t$, which depends only on $t$.  A plot of this
ratio versus v then allows one to obtain $\omega_L=1/\tau$ even if
this frequency falls {\it outside} the experimentally accessible frequency
window. A detailed discussion of the more general case, where the
signal consists of a superposition of independent fluctuators with a
distribution of relaxation times $\tau$ will be presented elsewhere
[19].

\begin{figure}
\centerline{\epsfxsize=8.5cm\epsfbox{fig6.eps}}
\caption{The characteristic frequencies, $\omega_L$ (open symbols) and
$\omega_H$ (solid symbols), in the power spectra plotted as a function
of $1/\Gamma$. The general trend is for both $\omega_L$ and $\omega_H$
to increase with increasing $\Gamma$ and decreasing depth into the
pile. The dashed line in (b) is a guide to the eye, indicating that
the trend is consistent with an activated-like behavior $\omega\approx
\exp(-\Gamma_0/\Gamma)$, with $\Gamma_0$. For comparison, (a) shows
the dependence of $\omega_H$ on $\Gamma$.}
\end{figure}


Figure 6 plots the resulting corner frequencies as a function of
applied acceleration. The trend is for both $\omega_L$ and $\omega_H$
to increase as a function of increasing $\Gamma$ and with decreasing
depth into the pile, see Fig.~6(a). We note that over the relatively
small available range of $\Gamma$, the variation of $\omega_H$ is
consistent with behavior reminiscent of thermal activation: $\omega_H=
\omega_0 \exp(-\Gamma_0/\Gamma)$. In this context, $\Gamma_0$ would
represent an energy barrier and $\omega_0$ would be an attempt
frequency. We find that a value of $\Gamma_0\approx 15$ is consistent
with all the data, and that the greatest variation is in the parameter
$\omega_0$, which varies from $2\times 10^{-3}$ to $7\times 10^{-2}$
for $\omega_L$ and $1$ to $15$ for $\omega_H$ .



\section{DISCUSSION}


Several mechanisms [17,20-22] have been proposed to explain the
kinetics of compaction. Although the proposed mechanisms are
compelling, their quantitative predictions fail to describe the time
dependence observed experimentally [1]. In light of our experimental
results, we pay special attention to models based on free volume
considerations as it appears that they not only capture the
experimentally observed slow relaxation towards the steady state, but
may also provide a valid framework for understanding the fluctuation
spectrum. Such models [8-10,12,13] include strong nearest-neighbor
repulsive interactions between particles that effectively block the
occupation of adjacent sites.  On very general grounds it is
reasonable to assume that for the case of granular compaction, the
rate of increase in volume density is exponentially reduced by free
volume [23-25]. One alternative approach to the compaction problem is
due to Linz [26] who proposes a phenomenological decay law for
successive inverse packing fractions. Moreover, recent models based on
the dynamics of crystalline clusters in the material have been
proposed by Gavrilov [27] and by Head and Rogers [28]. These
approaches lead to a time evolution that is essentially equivalent to
Eq.~(1). 

A simple heuristic argument [12,24] for the compaction
process illustrates how the effects of free volume can lead to the
observed inverse logarithmic behavior. If each grain has a volume
$V_g$ and we start with a number $n$ of grains per unit volume, then
the volume fraction is given by $\rho=nV_g$. In general there exists a
maximum possible volume fraction, $\rho_{\rm max}$ , corresponding to
the configuration of grains that occupies the least amount of volume.
Then, at some volume fraction $\rho$, the average free volume
available to each grain for rearrangements is $V_f=
V_g(1/\rho-1/\rho_{\rm max})$. During compaction, individual hardcore
grains move, and when a void large enough to contain a grain is
created, it is quickly filled by a new particle. When the volume
fraction is large, voids the size of a particle are rare and a large
number of voids must rearrange to accommodate an additional
particle. We can estimate the rate of compaction by assuming that $N$
grains must rearrange in such a way that they contribute their entire
free volume to create a grain-sized void, $NV_f=V_g$. We find that
this number increases as $N=\rho\rho_{\rm max}/(\rho_{\rm max}-\rho)$
during the compaction process. For independent, random grain motion
during a tap, the probability for $N$ grains to rearrange and open up
a grain-size void is then $e^{-N}$.  Consequently,
\begin{equation}
{d\rho/dt}\propto (1-\rho)e^{-\rho\rho_{\rm max}/ (\rho_{\rm max}-\rho)}.
\end{equation}
The rate at which the density increases is proportional to the void
volume and the probability for such a rearrangement. The latter
exponential factor reduces the rate and dominates for large $\rho$. In
the limit $\rho\to\rho_{\rm max}$ we have $N\approx \rho_{\rm
max}^2/(\rho_{\rm max}-\rho)$ and the solution of this equation is
given asymptotically by $\rho(t)=\rho_{\rm max} -A/[B+\ln(t)]$, where
$A$ and $B$ are constants [10,13,24,25]. This result closely
approximates our experimentally based fitting form, Eq.~(1), for the
ensemble-averaged $\rho(t)$ as it approaches the steady state.

The solution to Eq.~(2) always approaches the maximum density
$\rho_{\rm max}$ and does not allow for a lower steady-state
density. The reason that this model leads to jamming is the absence of
any void-creating mechanism that would be represented by a competing
term on the right-hand side of Eq. (2). (The consequences of including
a mechanism for the generation of voids are discussed in the next
section for the ``parking lot model’’ and in Appendix A.) The
competition between void annihilation and creation during tapping
naturally leads to density fluctuations.  We can also examine our data
for the dependence of the corner frequencies on the acceleration
$\Gamma$ [Fig. 6(b)], where we found that $\omega_H=\omega_0
\exp(-\Gamma_0/\Gamma)$. We use the fact that $\rho$ is a monotonic
function of $\Gamma$ in the reversible steady-state regime (Fig.~2)
and write to first order $\rho(\Gamma)\approx \rho_{\rm max}-m\Gamma$,
where $m$ is a positive constant, locally approximating the slope
$\Delta\rho/\Delta\Gamma$ (see Ref. [14] for data on other bead
sizes). Substituting in for $\Gamma$, the expression for $\omega_H$
can then be rewritten as $\omega_H=\omega_0\exp[-m\Gamma_0/(\rho_{\rm
max}-\rho)]$. This has the same form as the right-hand side of Eq.
(2) (in the limit that $\rho\to \rho_{\rm max}$ and indicates that the
kinetics depend sensitively on the available free volume, $1/\rho_{\rm
max}-1/\rho$. As $\Gamma$ is reduced and the density approaches the
maximum density, the kinetics slow down rapidly. The manner in which
the kinetics slow down is reminiscent of the Vogel-Fulcher form used
to describe another class of disordered metastable materials, namely,
glasses [29]. Similarities to glasses have recently been found in
another approach to the compaction process [8-10].


\section{THEORETICAL MODEL AND SIMULATIONS OF COMPACTION}

\subsection*{The parking lot model}


In an attempt to explicitly work out some of the consequences of the
free volume approach to granular compaction, we next discuss a
simplified model. The model was previously studied in the context of
chemisorption [23-25,30,31] and protein binding [32]. Despite its
simple nature, it gives remarkably good qualitative agreement with the
experimental data, both for the approach to the steady state and for
the spectrum of fluctuations in the steady state.  This model has the
advantage that it readily lends itself to computer simulations; we
restrict ourselves to the one-dimensional (1D) case, but extensions to
higher dimensions are straightforward. Moreover, much is known about
its low-density limit, for which mean-field equations exist that are
amenable to analytic treatment (see the Appendix). In 1D, the model
can be compared to parallel curbside parking where there are no marked
parking spaces. For the person wishing to park a vehicle, the familiar
situation is that there exist large, but not quite large enough,
spaces between previously parked cars. The analogous question to the
one we have been asking is ``How many other cars have to be moved just
a bit for the additional one to fit in?’’ 

\begin{figure}
\centerline{\epsfxsize=7.5cm\epsfbox{fig7.eps}}
\vspace{-.8in}
\caption{The adsorption-desorption process. Adsorption is successful
only in spaces large enough to accommodate a particle. Desorption of a
particle, on the other hand, is unrestricted.}
\end{figure}


The model is defined as follows: identical particles of unit length
adsorb uniformly from the bulk onto 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
between $t$ and $t+dt$ is $k_-dt$.  While the desorption process is
unrestricted, the adsorption process is subject to free volume
constraints, i.e., particles cannot adsorb on top of previously
adsorbed particles; see Fig.~7.  This stochastic process is
well-defined in arbitrary dimension and clearly satisfies detailed
balance so that the system eventually reaches a steady-state
density. In one dimension, $\rho_{\rm max}=1$.



Mapping the model on to the experiment, we associate an adsorption
event with the annihilation or filling of a void within the pile of
beads, whereas a desorption event is associated with the creation of a
void. The ratio of adsorption to desorption rates, $k=k_+ / k_-$,
determines the final steady-state density in the model (see also the
Appendix). Thus one can associate $k$ in the model with the magnitude of
the acceleration $\Gamma$ in our experiment.


\subsection*{Simulation of compaction based on the parking lot model}

In this section we compare the experimental results with Monte Carlo
simulations of the 1D parking process. The details of the simulation
are described elsewhere [24]. Here we report our results for a system
size of $00$ similar results were found for a system size of $25$.


The simulations were started from a zero density initial state and
allowed to evolve to various steady-state densities by varying $k_-$
at a fixed value of $k_+=1$. In Fig.~8 we show the time evolution of
the density as it approaches a steady-state density $\rho_{ss}=0.84$.
The steady-state densities obtained after equilibration coincide with
those predicted by Eq.~(A2a) in the Appendix. We find that the
simulations reproduce the slow logarithmic relaxation towards the
steady state in agreement with Eq. (1).

\begin{figure}
\centerline{\epsfxsize=7.5cm\epsfbox{fig8.eps}}
\vspace{-.5in}
\caption{The time evolution of the density in the simulation for
$k=10^3$ ($\rho_{ss}=0.84$). Time is in units of Monte Carlo steps
(MCS). The solid line represents a fit to Eq.~(1).}
\end{figure}

We now turn to the density fluctuations. In this case, the simulations
ran long enough to ensure that a steady-state density was attained
before density fluctuations were recorded. For low densities
$\rho_{ss}<0.37$ we find that in the simulation the power spectra of
the fluctuations are best described by a Lorentzian having a single
characteristic time scale, as expected from the mean-field analysis
[see Eq.~(A4) in the Appendix]. However, at higher densities (higher
$k$) local fluctuations dominate the dynamics and the power spectra
show the emergence of two distinct corner frequencies, which become
progressively more separated. This is shown in Fig.~9 for a wide range
of ratios $k$, where $k_+$ was fixed at a value of 1. Most notable is
the low-frequency corner, which shifts rapidly to lower frequencies
for small increments in density. By comparison, the high-frequency
corner decreases much more slowly.

Our simulations find power spectra, $S_\rho(\omega)$, strikingly
similar in shape to those obtained experimentally when $\rho_{ss}>
0.50$. Again we see three distinct regimes (Fig.~9). Below a corner
frequency, $\omega_L$ , there is white noise [$S_\rho(\omega)\propto
\omega_0$].  Above a high frequency corner, $\omega_H$ ,
$S_\rho(\omega)\propto \omega^{-2}$ . The simulations offer the
advantage of allowing the separation between $\omega_L$ and $\omega_H$
to be tuned by increasing the value of $k$ or, equivalently,
increasing the density. This allows the systematic investigation of
the spectral dependence in the intermediate regime between the two
corner frequencies. As in the experimental data, we find that there is
a Lorentzian tail, $S_\rho(\omega)\propto \omega^{-2}$ just above
$\omega_L$ . At higher frequencies, stretching up to the 
high-frequency corner $\omega_H$ , we find a power-law regime
$S_\rho(\omega)\propto \omega^{-\delta}$.  The exponent $\delta$
appears to depend slightly on the separation between the two corner
frequencies. For the largest separations that span nearly 5 decades in
frequency we find $\delta\approx 0.5$. This value is smaller than that
found in the experimental data ($\delta\approx0.9$), but again is
inconsistent with a simple superposition of two Lorentzians having
characteristic time scales $\omega_L^{-1}$ and $\omega_H^{-1}$.

\begin{figure}
\centerline{\epsfxsize=7.5cm\epsfbox{fig9.eps}}
\vspace{-.5in}
\caption{Power spectra, $S_\rho$ (MCS), of the density fluctuations in
the simulation of the one-dimensional parking process. The evolution
of the spectral dependence is shown for values of the ratio $k=33$,
$10^2$, $10^3$, $10^4$, corresponding to final steady-state densities
$\rho_{ss}=0.72$, $0.77$, $0.84$, and $0.88$, respectively (see
text). The strongest dependence on k is for $\omega_L$ , which
decreases rapidly as the density increases. Such spectra are similar
to those found in the experiment, see Fig.~5.}
\end{figure}



\section{ANALOGY WITH THERMAL FLUCTUATIONS}

In ordinary statistical mechanics, the fluctuation-dissipation theorem
allows the determination of the response of a system to a small
perturbation from its thermal fluctuations about equilibrium. In this
section, we will explore the possibility that we can derive similar
information about the granular system from its fluctuations about its
steady-state density. In the granular thermodynamics theory developed
by Edwards and coworkers, an analogy is made between granular and
thermal systems. The basic assumption is that the {\sl volume} $V$ of
a powder is analogous to the {\sl energy} of a statistical system (we
note that $V$ here refers to the total volume and not just to the free
volume). Instead of a Hamiltonian, there is a function that specifies
the volume of the system in terms of the positions of the individual
grains. The ``entropy’’ is thus the logarithm of the number of
configurations: $S=\lambda \ln \int d$ (all configurations) where
$\lambda$ is the analog of Boltzmann’s constant. Using this they
defined a quantity analogous to a temperature in a thermal system,
which they call the ``compactivity’’ $X$: $X= \partial V/ \partial
S$. In contrast to the notion of ``granular temperature,’’ which
depends on the random motion of the particles, the compactivity
characterizes the static system after it has reached a steady-state
density via some preparation algorithm. Such an algorithm would be one
as we have described above, where we have vibrated the granular system
until it has reached the reversible steady-state density. If this
theory is valid, then we should be able to define an equilibrium such
that two systems in equilibrium with a third system are also in
equilibrium with each other. That is, no net volume will be
transferred between the two systems when they are placed in contact
with each other if they have the same value of $X$.

In a thermal system we can write the specific heat in two ways as
follows:
\begin{equation}
C_V=dE_0/dT|_V=\langle(E-E_0)^2/k_BT^2,
\end{equation}
where $E_0$ is the equilibrium average of the energy $E$ of the
system, $k_B$ is Boltzmann’s constant, $T$ is the temperature, and
$\langle\cdots\rangle$ represents the time average. In Edwards’ theory
for a powder the analogous quantity to the specific heat of a thermal
system given in Eq.~(3) becomes
\begin{equation}
C=dV_{ss}/dX=\langle(V-V_{ss})^2\rangle/ \lambda X^2,
\end{equation}
where $V_{ss}$ is the steady-state volume. Since we have measured the
density fluctuations in the steady state (Figs.~3-5), we are in a
position to explicitly calculate the variance,
$\langle(V-V_{ss})^2\rangle$, of volume fluctuations for a given
steady-state volume $V_{ss}$ defined here as $V_{ss}=1/\rho_{ss}$. We
can then write
\begin{equation}
\int_{V_1}^{V_2}\!\!\! dV_{ss} / \langle(V-V_{ss})^2\rangle
=\int_{X_1}^{X_2}\!\!\!dX/\lambda X^2 =1/\lambda X_1-1/\lambda X_2.
\end{equation}
Equation (5) allows us to measure the difference in compactivities for
any two volumes as long as we know the fluctuations of the volumes
(i.e., densities) as a function of the average volume.  This is
equivalent to obtaining the difference in temperatures for a thermal
system between any two energies. Clearly, as $V_{ss}$ increases $X$ is
expected to increase as well. Equation (5) allows the determination of
an absolute value for the compactivity only once a suitable point of
reference can be found.  In Fig.~10 we show the experimentally
obtained values of $\langle(V- V_{ss})^2\rangle$ for several
steady-state volumes along the reversible branch of Fig.~2 over the
range $4<\Gamma<7$ [33].  The solid curve through the data for the top
of the pile represents a linear fit to the function:
\begin{equation}
\langle(V-V_{ss})^2\rangle=a+bV_{ss},
\end{equation}
where $a=-7.2\times 10^{-4}$ and $b=4.9\times 10^{-4}$. This implies
that the magnitude of the fluctuations goes to zero at $\rho=0.68$,
that is, near the close-packed density. Using this form for the
dependence of the fluctuations on $V_{ss}$ in Eq.~(5) we find that
\begin{equation}
1/\lambda X\propto \ln(a+bV_{ss}).
\end{equation}
This functional dependence is valid only over the limited range of
experimental data, and may not be an adequate description of the
general behavior. Below, we discuss a similar analysis for the
simulation data for which a broader range of volumes can be
explored. Using Eq. (6), we can evaluate the difference in inverse
compactivities between any two steady-state volumes. We find that
$1/\lambda X_1-1/\lambda X_2=0.04$ where the subscripts 1 and 2 refer
to the smallest and largest volumes for which we have data.  This
result explicitly demonstrates how the compactivity increases for
larger volumes (smaller densities).

It is also interesting to consider the size of the fluctuating volumes
that give rise to the observed variance. This can be estimated by
assuming that the fractional fluctuations scale as
$\langle\delta\rho^2\rangle / \rho_{ss}^2=\langle \delta V^2\rangle /
V_{ss}=\kappa^2/N$, where $\delta V=V-V_{ss}$ . This is the usual
$N^{-1/2}$ classical self-averaging property of $N$ in dependently
fluctuating variables. The parameter $\kappa$ accounts for the fact
that there exists a maximum range of density changes for each grain
that does not compromise the mechanical stability of the granular
assembly. Figure 10 indicates that $\langle\delta
V^2\rangle/V_{ss}^2>1/40,000$, and as an upper bound we let
$\kappa=\rho_{\rm xtal}-\rho_{\rm RLP}=0.74-0.55\approx 0.2$. We then
find that $N>1600$. Since we know that each capacitor averages over a
volume corresponding to $6000$ beads, this suggests that there are
roughly $1600$ independently fluctuating clusters each consisting of
$\sim 4$ beads (lower bound).

\begin{figure}
\centerline{\epsfxsize=8.0cm\epsfbox{fig10.eps}}
\vspace{-.5in}
\caption{The average variance of the experimental volume fluctuations
(open symbols) as a function of the steady-state volume. The trend is
for the variance to increase with increasing volume and {\it depth}
into the pile. The solid symbols represent the variance as determined
from the distribution of fluctuation amplitudes in Fig.~4 (see
text). The dashed lines are linear fits to the solid symbols.}
\end{figure}


An important feature that can be seen in Fig. 10 is that the variance
becomes systematically larger the deeper into the pile one goes [see
also $D(\delta\rho)$ for $\Gamma=5.1$ in Fig.~4].  For the middle and
bottom sections of the column, the variance appears nonmonotonic with
a peak near $V_{ss}=1.592$, see open symbols in Fig.~10. Examination
of the corresponding distributions of fluctuation amplitudes
($\Gamma=5.9$ in Fig.~4) indicates that the increased variance is due
to a non-Gaussian tail in the distributions (see the description of
Fig.~4, above). For comparison, we have also determined the variance
from the slopes of the distribution functions in Fig.~4. We used the
slopes corresponding to the low-density side of the distributions
because these were most consistent with a Gaussian form over all
accelerations and depths. In this way, the effect of the non-Gaussian
tails in some of the distributions can be avoided.  These results are
shown as solid symbols in Fig.~10 and the dashed lines correspond to
linear fits through the data. Here too, it is evident that the
variance is larger for larger depths. A larger variance implies a
correspondingly larger phase space. At first this seems
counterintuitive be-cause the time records [Fig. 3(a)] and power
spectra (Fig.~5) indicate that density fluctuations are slower at the
bottom of the pile. Although the kinetics near the bottom of the pile
may be slower, there is a greater number of configurations with
different volumes that are accessible to those beads.

\begin{figure}
\centerline{\epsfxsize=8.0cm\epsfbox{fig11.eps}}
\vspace{-.5in}
\caption{The average variance of the volume fluctuations in the
simulation as a function of the steady-state volume, $V_{ss}$.  The
fractional variance as a function of $V_{ss}-1$ is plotted on
logarithmic axes in the inset. The solid line is a power-law fit given
by $0.0124(V_{ss}-1)^{1.37}$.}
\end{figure}

A depth dependence to the variance also suggests the presence of a
gradient in the compactivity. In Fig.~10, we used the average
steady-state volumes $V_{ss}$ , obtained from optical measurements of
the total column height. One possibility is that this average volume
density does not accurately represent the density in the different
sections of the pile. If so, the larger compactivity near the bottom
of the pile then implies that the bottom beads are actually in a less
compact state than those at the top. However, from the trend in
$\langle \delta V^2\rangle$ versus $V_{ss}$ in Fig.~10 the difference
in packing fraction be-tween the top and bottom of the pile that would
be necessary to have the variances be equal would be
$\Delta\rho\approx 0.035$. Since this difference is huge on the scale
of $\rho(\Gamma)$ for the reversible branch in Fig.~2 we do not
believe this to be a plausible explanation.  Rather, it appears that
there is another variable, such as pressure, in addition to the
volume, that controls the depth dependence of the
fluctuations. Indeed, supporting evidence to this effect can also be
seen in Fig.~6, which shows that the high-frequency corner $\omega_H$
decreases with increasing depth into the pile. Nevertheless, we expect
that the system is entirely jammed $\langle\delta V^2\rangle\to 0$ at
the same density (i.e., $\rho_{\rm max}$) for all depths in the pile.

\begin{figure}
\centerline{\epsfxsize=8.5cm\epsfbox{fig12.eps}}
\vspace{-.9in}
\caption{The left hand side of Eq.~(5) is numerically evaluated and
plotted as a function of steady-state volume for the simulation
data. As plotted, the difference in inverse compactivities between the
highest density configuration ($V_{ss}-1$) and a low density
configuration (higher $V_{ss}$) can be read off directly.}
\end{figure}

With the simulation described above, a broader range of densities can
be explored than that which is experimentally accessible.  Figure 11
shows the dependence of the variance in volume fluctuations as a
function of steady-state volume for the 1D parking lot model. The
rapid decrease in variance near $V_{ss}-1$ suggests that there may be
a diverging length or time scale as the system approaches its most
compact state. Indeed, plotting the normalized variance as a function
of the free volume ($V_{ss}-1$) does reveal power-law-like behavior
$\langle\delta V^2\rangle\propto (V_{ss}-1)^\beta$ with $\beta\approx
1.4$. This is shown in the inset of Fig.~11. Proceeding with the
compactivity analysis, the data in Fig.~11 was numerically integrated
to yield the left-hand side of Eq.~(5) to within a constant. An
absolute value cannot be established with just our data. In Fig.~12 we
plot the difference in inverse compactivity as a function of volume.
This difference is with respect to the state $V_{ss}=1.1$. Figure 12
indicates a nontrivial functional dependence to the increase in
compactivity with system volume.

For comparison, the 3D experimental results correspond to relatively
high densities in the 1D simulation because in 3D the available void
volume is with respect to the random close packed limit ($\approx
0.64$) while the corresponding limit in 1D is 1. Taking this into
account, the 30\% increase in $\langle\delta V^2\rangle$ in the
experimental results shown in Fig.~10 compares well with the
simulation data in Fig.~11 over a similarly restricted range in
$V_{ss}$.



\section{CONCLUSIONS}

In this paper, we have examined the volume fluctuations about a
steady-state density for a granular system. For these measurements to
provide a useful analogy with a thermal system, it is essential that
the fluctuations be measured in steady-state conditions. For this
reason, we have explicitly taken data on the reversible density line
as shown in Fig.~2. From these measurements, we have been able to
determine experimentally the compactivity, which is the quantity
analogous to the temperature in the theory of Edwards et al.

Theories based on free volume seem particularly well suited for
describing the data. As the system approaches its final state, a
growing number of particles have to be rear-ranged in order for the
density to be increased locally. The rate of increase in the density
is exponentially reduced by this number leading to a logarithmically
slow approach to the steady-state density as observed
experimentally. Monte Carlo simulations of a one-dimensional
adsorption-desorption process based on these ideas show fluctuations
about the steady state density that are strikingly similar to those
observed experimentally. These results attest to the importance of
volume exclusion for granular relaxation and steady-state dynamics
under vibration.



Despite this model’s simplicity and obvious shortcomings, it appears
to capture an essential mechanism underlying the remarkably slow
relaxation and the nature of the density fluctuations. This mechanism
is associated with the reduction of free volume available for particle
motion as the density increases. Although our simple model cannot
predict the experimental values of the fitting parameters in Eq.~(1),
the inverse logarithmic density relaxation towards the steady state is
the same one observed for granular compaction (see Fig.~1).

In the simulation model, our treatment was restricted to
one-dimensional processes, but we expect that the results hold in
higher dimensions as well. There are other important distinctions
between the model and real granular media. One difference is that in
the structure of a granular assembly the particles form contact
networks. The creation of a void (a desorption event in the model)
therefore requires the rearrangement of several particles and is thus
restricted just as is the annihilation of voids. Another difference is
the mechanical stability of a granular assembly.  This property will
place limits on the free volume available for large void creation. For
instance, for spherical particles in three dimensions and in the
presence of gravity the available free volume is deter-mined by the
restricted range of accessible volume fractions, namely, between the
random close packed $\rho\approx0.64$ and random loose packed
$\rho\cong0.55$ configurations.


It is interesting to speculate whether the reduction in free volume
leads to a crossover from a simple independent particle picture for
compaction to a more complex process at higher densities, presumably
involving correlations over increasingly longer length and time
scales. In this regard, we have demonstrated that density fluctuations
are an important probe of the underlying microscopic dynamics. Indeed,
the study of fluctuations may elucidate the physics of independent and
cooperative-particle motions, which lead to the macroscopic response
of a powder subject to vertical vibrations. For instance, it is
interesting to note that from both the experimental and simulation
data there appear to be two characteristic time scales, related to the
corner frequencies $\omega_L$ and $\omega_H$ in the power spectra,
that characterize the steady-state dynamics. This behavior may be
related to the results found in 3D simulations of vibrated powders by
Mehta and Barker [17,20,34] and in simulations of a frustrated lattice
gas [8-10]. Those results suggested the existence of two exponential
relaxation mechanisms: the faster of the two involves the motions of
independent particles while the slower involves collective particle
motions, which were found to be diffusive. However, we emphasize again
that both our experimental and simulation data are not consistent with
a superposition of two independent exponentially decaying processes.

In this paper we have presented results for monodisperse spherical
particles subject to vertical shaking. Realistic powders are far more
complicated with properties that depend on factors like cohesive
forces, polydispersity in size, and ir-regularity or anisotropy in
shape.  Nevertheless, our results can provide a valuable benchmark for
evaluating the predictions of theoretical models and simulations. The
applicability of concepts such as compactivity or ``granular
temperature’’ in the description of quasistatic granular media
requires further exploration. In particular, it would be interesting
to ex-amine the properties of granular systems comprised of particles
with shape anisotropies and subject to isotropic shaking. Furthermore,
our experimental data suggest that the steady-state properties of a
granular assembly cannot be fully described by a single state
variable, i.e., the volume. Rather, another variable is required to
account for the depth dependence of the volume fluctuations.

{\it Note added in proof.} The width of the density fluctuations in
the parking lot model can be calculated in the mean-field approach
described in the Appendix. For details see E.~Ben-Naim et al., Physica
D (to be published). Such calculation predicts a power law as seen in
the inset to Fig.~11 for the simulation data, but with an exponent
$\beta=2$.


\section{ACKNOWLEDGMENTS}
We are grateful to M.~L.~Povinelli and S.~Tseng for assistance on
certain aspects of this work. It is a pleasure to acknowledge
stimulating discussions with S.~Coppersmith and T.~Witten. This work
was supported by the NSF through MRSEC Grant DMR-9400379 and through
Grant No. CTS-9710991. We acknowledge additional support from the
David and Lucile Packard Foundation, and from the Research
Corporation.


\appendix
\section*{}
For the one-dimensional parking lot model an analytical mean-field
description exists. On the continuum, an approximate rate equation for
the density evolution was constructed from the exact steady-state void
distribution.  This equation yields an approach to the steady state
that is essentially identical to that found in the experiment [i.e.,
Eq.~(1)]. However, it is less successful in capturing the fluctuation
behavior, particularly for the high densities relevant here. We
summarize the salient analytic results for the model. Details can be
found in Refs. [23-25]. A modified Langmuir equation can be written
for the rate of change in density [24]:
\begin{equation}
{d\rho/dt}=k_+ (1-\rho)e^{-\rho/(1-\rho)} -k_-\rho.
\end{equation}
The gain term is proportional to the fraction of unoccupied space,
which is modified by an ``excluded volume’’ constraint. It was
previously shown that in steady state the prob-ability $s(\rho)$
that an adsorption event is successful is given by
$s(\rho)=e^{-\rho/(1-\rho)}$ [24]. This so-called ``sticking
coefficient’’ vanishes exponentially as $\rho\to1$. This
effectively reduces the sticking rate, $k_+\to
k_+(\rho)=k_+s(\rho)$. The desorption process, on the other hand,
is unrestricted and so the loss term is proportional to the
density itself. The steady-state density $\rho_{ss}$ , which is
obtained by imposing $d\rho/dt=0$, can be determined as a function
of the adsorption to desorption rate ratio, $k=k_+/k_-$, from the
following transcendental equation:
\begin{displaymath}
\alpha e^{\alpha}=k,\qquad{\rm with}\qquad {\alpha=\rho_{ss}/(1-\rho_{ss})}.
\hfill(A2a)
\end{displaymath}
The following leading behavior in the two limiting cases is found
\begin{displaymath}
\rho_{ss}(k)\cong\cases{k, &$k\ll 1$;\cr 1-(\ln k)^{-1},&$k\gg1$.} \hfill(A2b)
\end{displaymath}

The effect of the volume exclusion constraint is striking, a huge
adsorption to desorption rate ratio, $k>10^9$, is necessary to achieve
a $0.95$ steady-state density. We now focus on the relaxation
properties of the system. The granular compaction process corresponds
to the high density limit, and we thus consider the
desorption-controlled case, $k\gg1$. Hence, let us fix $k_+=1$ and
consider the limit $k_-\to 0$ of Eq. (A1). For $t\gg 1/k_+$ , it can
be shown that the system approaches complete coverage,
$\rho_\infty=1$, according to [24,25]
\begin{equation}
\rho(t)\cong \rho_\infty-1/(\ln k_+ t).
\end{equation}
This is confirmed by numerical simulations in one dimension (see
Ref. [24] and Sec. IV). We conclude that the excluded volume
constraint gives rise to a slow relaxation.

Equation (A3) holds indefinitely only for the truly irreversible
limit of the parking process, i.e., for $k=\infty$. For large but finite rate
ratios, the final density is given by Eq.~(A2b). By computing how a
small perturbation from the steady state decays with time, an
exponential relaxation towards the steady state is found
$|\rho_{ss}-\rho(t)|\propto e^{-1/T}$ for $t\gg 1/k_-$. The relaxation time is
\begin{equation}
T=(1-\rho_{ss})^2/k_- .
\end{equation}
The above results can be simply understood: 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.

\begin{thebibliography}{99}

\bibitem{} J.~B.~Knight et al., Phys. Rev. E {\bf 51}, 3957 (1995).
\bibitem{} S.~F.~Edwards and R.~B.~S.~Oakeshott, Physica A {\bf 157}, 1080 (1989).
\bibitem{} A.~Mehta and S.~F.~Edwards, Physica A {\bf 157}, 1091 (1989).
\bibitem{} S.~F.~Edwards and C.~C.~Mounfield, Physica A {\bf 210}, 290 (1994).
\bibitem{} S.~Ogawa, in Proceedings of the U.S.-Japan Symposium on
Continuum Mechanics and Statistical Approaches in the Me-chanics of
Granular Materials, edited by S. ~C.~Cowin and M.~Satake (Gakujutsu
Bunken Fukyukai, Tokyo, 1979).
\bibitem{} C.~S.~Campbell, Annu. Rev. Fluid Mech. {\bf 22},57 (1990).
\bibitem{} I.~Ippolito et al., Phys. Rev. E {\bf 52}, 2072 (1995).
\bibitem{} A.~Coniglio and H.~J.~Herrmann, Physica A {\bf 225}, 1 (1996).
\bibitem{} M.~Nicodemi, A.~Coniglio, and H.~J.~Herrmann, J. Phys. A
{\bf 30}, L379 (1997).
\bibitem{} M.~Nicodemi, A.~Coniglio, and H.~J.~Herrmann, Phys. Rev. E {\bf 55}, 3962
(1997).
\bibitem{} H.~M.~Jaeger, S.~R.~Nagel, and R.~P.~Behringer, Phys.
Today {\bf 49}, 32 (1996).
\bibitem{} H.~M.~Jaeger, S.~R.~Nagel, and R.~P.~Behringer,
Rev. Mod. Phys. {\bf 68}, 1259 (1996).
\bibitem{} P.~G.~de Gennes (unpublished);
\bibitem{} T.~Boutreux and P.~G.~de Gennes, Physica A {\bf 224},59 (1997).
\bibitem{} E.~R.~Nowak et al., Powder Technol.  {\bf 94}, 79 (1997).
\bibitem{} J.~B.~Knight, H.~M.~Jaeger, and S.~R.~Nagel, Phys. Rev. Lett.
{\bf 70}, 3728 (1993).
\bibitem{} K.~Shinohara, in Handbook of Powder Science and Technology, edited by
M.~E.~Fayed and L.~Otten (van Nostrand Rein-hold Co., New York, 1984).
\bibitem{} A.~Mehta and G.~C.~Barker, Phys. Rev. Lett. {\bf 67}, 394 (1991).
\bibitem{} M.~Nelkin and A.-M.~S.~Tremblay, J. Stat. Phys. {\bf 25}, 253 (1981).
\bibitem{} E.~R.~Nowak, S.~R.~Nagel, and H.~M.~Jaeger (unpublished).
\bibitem{} G.~C.~~Barker and A.~Mehta, Phys. Rev. E {\bf 47}, 184 (1993).
\bibitem{} H.~S.~Caram and D.~C.~Hong, Phys. Rev. Lett. {\bf 67}, 828 (1991).
\bibitem{} D.~C.~Hong et al., Phys. Rev. E {\bf 50}, 4123 (1994).
\bibitem{} G.~Tarjus and P.~Viot, Phys. Rev. Lett. {\bf 68}, 2354 (1992).
\bibitem{} P.~L.~Krapivsky and E.~Ben-Naim, J. Chem. Phys. {\bf 100}, 6778 (1994).
\bibitem{} X.~Jin, G.~Tarjus, and J.~Talbot, J. Phys. A 27, L195
(1994).
\bibitem{} S.~J.~Linz, Phys.  Rev. E {\bf 54}, 2925 (1996).
\bibitem{} K.~Gavrilov (unpublished).
\bibitem{} D.~A.~Head and G.~J.~Rodgers (unpublished).
\bibitem{} M.~D.~Ediger, C.~A.~Angell, and S.~R.~Nagel,
J. Phys. Chem. {\bf 100}, 13200 (1996).
\bibitem{} V.~Privman and M.~Barma, J. Chem. Phys. {\bf 97}, 6714 (1992).
\bibitem{} J.~W.~ Evans, Rev. Mod. Phys. {\bf 65}, 1281 (1993).
\bibitem{} J.~D.~McGhee and P.~H.~von Hippel, J. Mol. Biol. {\bf 86}, 469 (1974).
\bibitem{} The variance is also equal to the integrated power spectrum
over all frequencies.  Because the power spectral dependence is weaker
than $\omega^{-2}$ at lowest frequencies and tends to fall off as $\omega^{-1}$ at the
highest frequencies we can have confidence that the variance is
finite.
\bibitem{} G.~C.~Barker and A.~Mehta, Phys. Rev. A {\bf 45}, 3435 (1992).
\end{thebibliography}


\end{multicols}
\end{document}
