# Hysteresis in vibrated granular media

## I Introduction

The phenomenology of granular materials is very rich. One of the most characteristic dynamical behaviors is compaction. If a system of grains in a loosely packed configuration is submitted to vertical vibrations, or tappings, it approaches slowly a more compact state. The relaxation is well described by an inverse logarithmic law, and both the relaxation function and the density of the “stationary” state seem to depend only on a dimensionless parameter characterizing the vibration intensity .

More recently, it has also been shown that granular materials exhibit irreversible-reversible cycles if the vibration strength is increased and decreased alternately. When the system is tapped with increasing intensity starting from a loosely packed state, the density shows a non-monotonic behavior, with a maximum at a certain value of the intensity of vibration. If, afterwards, the granular medium is tapped with decreasing intensity, the density increases monotonically, and hysteresis effects show up. Interestingly, if the vibration intensity is again increased, the density follows a curve that is approximately equal to the evolution in the decreasing intensity process, being thus “reversible” . Of course, the rate of variation of the vibration intensity is the same for all the processes described above.

The static properties of granular materials have been studied by Edwards and co-workers . The starting point is the plausible idea that all the “microscopic” configurations of a powder having the same volume are equiprobable, provided that the powder has been prepared by extensive manipulation, i. e., by processes which do not act on individual grains. This has been called the “ergodic hypothesis” of powders. On this basis, it is possible to make an analogy between the variables of a molecular system and the parameters characterizing the state of a powder. The volume of a powder is analogous to the energy, while the entropy remains the same quantity, measuring the available number of configurations. The derivative of the volume with respect to the entropy is called the “compactivity” of the powder, and plays the role of the temperature. The lowest compactivity corresponds to the densest state (minimum volume), and the highest compactivity to the fluffiest stable configuration (maximum volume). The stationary value of the volume is an increasing function of the compactivity, like the energy increases with the temperature in a thermal system.

It seems interesting to try to establish a connection between the dynamical and the equilibrium properties of powders. It looks reasonable that, if the steady state is reached in a tapping process, there should be a relationship between the vibration intensity and the compactivity of a powder. The tapping process allows the system to explore the phase space of available configurations, and the same role is played by the temperature in a thermal system. In this way, it is tempting to explain the tendency of a granular system to move over an almost reversible curve, when the tapping strength varies in time, as the approach towards the stationary curve of the powder. If the above is true, the compactivity will be an increasing function of the vibration intensity, because the density over the “reversible” curve decreases with the tapping strength. In fact, this kind of behavior has been found in a simple model of a granular system described in terms of a Fokker-Planck equation .

Also, glass formers exhibit hysteresis effects when cooled and heated through the glass transition region. When a glass former is cooled, the system follows the equilibrium curve until a certain temperature $`T_g`$, at which it gets frozen due to the fast increase of the relaxation time. If this structural glass is reheated from its frozen state, the equilibrium curve is only approached for temperatures larger than $`T_g`$. This is due to the fact that the system starts from a configuration in which the structural rearrangements are very difficult . In this way, also hysteresis effects show up in glass formers, when they are submitted to thermal cycles.

The study of simple models have been very useful in order to understand, at least qualitatively, the behavior of structural glasses. In particular, hysteresis cycles are also shown by simple models when cooled and heated . The analytical approach to thermal cycles is a very difficult task, because it is necessary to solve the kinetic equations of the model with time-dependent coefficients. Nevertheless, some models whose dynamics is formulated by means of a simple master equation can be exactly solved . In those situations, the role played by the “normal solution” of the master equation, which is monotonically approached by all the other solutions , has been shown to play a fundamental role.

Here we will consider a simple model for granular media, which has been previously introduced . Its dynamical evolution is governed by a master equation, with the transition rates given as functions of the tapping strength. When the system is vibrated with a given intensity from a low density state, the density of the system increases monotonically until reaching an stationary value. The relaxation of the density is very slow, being very well fitted by an empirical inverse-logarithmic law. Interestingly, the stationary state of the system is described by Edward’s theory of powders. Then, it seems worth studying its behavior when the tapping strength changes in time for several reasons. Firstly, in order to verify if its behavior resembles that of real granular systems, showing the irreversible and reversible branches found in experiments . Secondly, to understand whether the hysteresis effects are related to the existence of a “normal solution” of the master equation with time-dependent tapping strength.

The paper is organized as follows. In Sec. II some general properties of models for granular media based on a master equation formulation of the dynamics are considered. The existence of a “normal solution” of the master equation is analyzed for the case of time-dependent transition rates. Also, the conditions to be verified by the law of variation of the tapping strength in order to guarantee that the system approaches the steady state curve are discussed. Section III is devoted to the analysis of the normal solution by means of Hilbert’s expansion of the master equation. A quite general expression for the first order deviation of the normal solution from the stationary curve is obtained. The specific model to be considered is presented is Sec. IV, as well as a brief discussion of its equilibrium state. Section V deals with the behavior of the model in processes with time-dependent vibration intensity. Firstly, processes in which the tapping strength decreases in time are studied. The existence of a phenomenon similar to the laboratory glass transition is analyzed. Secondly, we discuss processes with increasing vibration intensity. The role played by the normal solution turns out to be fundamental. Finally, the main conclusions of the paper are summarized in Sec. VI.

## II Some general dynamical properties

Let us consider a model system whose dynamics is described by means of the master equation 

$$\frac{dp_i(t)}{dt}=\underset{j}{}\left[W_{ij}(t)p_j(t)W_{ji}(t)p_i(t)\right].$$
(1)   
Here $`p_i(t)`$ is the probability of finding the system in state $`i`$ at time $`t`$, and the rates $`W_{ij}`$ for transition from state $`j`$ to state $`i`$ depend on time in a given way, independently of the state of the system. Let us define a function

$$H(t)=\underset{i}{}p_i(t)\mathrm{ln}\frac{p_i(t)}{p_i^{}(t)},$$
(2)   
where $`p_i(t)`$ and $`p_i^{}(t)`$ are two solutions of Eq. (1) corresponding to different initial conditions. The above definition for $`H(t)`$ assumes that $`p_i^{}(t)`$ is positive for all the states $`i`$. This condition will be fulfilled after a transient period if the process defined by Eq. (1) is irreducible, even in the case that some initial probabilities vanish. The time variation of $`H(t)`$ is given by

$$\frac{dH(t)}{dt}=A(t)$$
(3)   
with $`A(t)`$ being a complicated functional of the two solutions $`p`$ and $`p^{}`$ , that can be written in the form

$$A(t)=\underset{ij}{}W_{ij}p_j^{}\left[\left(\frac{p_j}{p_j^{}}\frac{p_i}{p_i^{}}\right)\left(\mathrm{ln}\frac{p_i}{p_i^{}}+1\right)+\frac{p_i}{p_i^{}}\mathrm{ln}\frac{p_i}{p_i^{}}\frac{p_j}{p_j^{}}\mathrm{ln}\frac{p_j}{p_j^{}}\right].$$
(4)   
Its main property is that $`A(t)0`$. Besides, if the transition rates define an irreducible process at time $`t`$, the equality sign holds only when

$$p_i(t)=p_i^{}(t),$$
(5)   
for all the states $`i`$. As the function $`H(t)`$ is bounded below, $`H(t)0`$, it must tend to a limit and, therefore,

$$\underset{t\mathrm{}}{lim}p_i(t)=\underset{t\mathrm{}}{lim}p_i^{}(t).$$
(6)   
Thus all the solutions of the master equation converge toward the same behavior, if the long-time limit of the transition rates still define an irreducible process. This equation can be understood as showing the existence of a long-time regime, where the influence of the initial conditions has been lost, and the state of the system and its dynamics is fully determined by the law of variation of the transition rates. Therefore, there will be a special solution of the master equation such that all the other solutions approach it after an initial transient period. We will refer to this special solution as the “normal” solution of the master equation for the given time dependence of the transition rates. A more detailed and general discussion of the H-theorem leading to the existence of this “normal” solution can be found in Ref. .

Now suppose that the master equation models the dynamical behavior of a granular system submitted to vertical vibration. The transition rates $`W_{ij}`$ would be functions of the parameter $`\mathrm{\Gamma }`$ characterizing the strength of the vibration. If the granular pile is vibrated with sinusoidal pulses of amplitude $`A`$ and frequency $`\omega `$, the quantity $`\mathrm{\Gamma }=A\omega ^2/g`$ , where $`g`$ is the gravity, is usually defined. In the case of time dependent intensity $`\mathrm{\Gamma }`$, the equation determining the evolution of the granular pile will be of the type given by Eq. (1). Assuming that for arbitrary $`\mathrm{\Gamma }0`$ the tapping process allows the system to explore the whole configuration space, the stochastic process will be irreducible and the existence of a normal solution for a given program of variation of the intensity $`\mathrm{\Gamma }`$ follows, provided that $`\mathrm{\Gamma }`$ does not vanish in the long-time limit.

Let us assume that for every value of the intensity $`\mathrm{\Gamma }>0`$, the equation

$$\underset{j}{}W_{ij}(\mathrm{\Gamma })p_j^{(s)}(\mathrm{\Gamma })=\underset{j}{}W_{ji}(\mathrm{\Gamma })p_i^{(s)}(\mathrm{\Gamma })$$
(7)   
has a “canonical” solution, of the form

$$p_i^{(s)}(\mathrm{\Gamma })=\frac{1}{Z(X)}\mathrm{exp}\left[\frac{V_i}{\lambda X}\right],$$
(8)   
where $`V_i`$ is the volume of the system in state $`i`$,

$$Z(X)=\underset{i}{}\mathrm{exp}\left[\frac{V_i}{\lambda X}\right]$$
(9)   
is a partition function, $`\lambda `$ is a constant with the dimension of volume, and $`X`$ is a variable termed the compactivity of the granular system in the framework of Edward’s statistical mechanics theory of powders . Of course, in the context considered here the compactivity $`X`$ will be a function of the vibration intensity $`\mathrm{\Gamma }`$, $`X=f(\mathrm{\Gamma })`$, which has to be found for each particular model. Equation (8) defines the steady state reached by the system if the vibration intensity is constant in time. The macroscopic value of the volume in this state is

$$\overline{V}^{(s)}(X)=\underset{i}{}V_ip_i^{(s)}\left[\mathrm{\Gamma }(X)\right],$$
(10)   
and the configurational entropy $`S(X)`$ reads

$$S(X)=\lambda \underset{i}{}p_i^{(s)}(X)\mathrm{ln}p_i^{(s)}(X).$$
(11)   
Of course, the compactivity $`X`$ can be obtained from its usual definition,

$$X=\frac{d\overline{V}^{(s)}}{dS}.$$
(12)   
For the sake of simplicity, in the following we will take the unit of volume such that $`\lambda =1`$. The stationary volume $`\overline{V}^{(s)}(X)`$ is an increasing function of the compactivity $`X`$, since the powder “compressibility”

$$\kappa (X)\frac{d\overline{V}^{(s)}(X)}{dX}=\frac{1}{X^2}\underset{i}{}\left[V_i\overline{V}^{(s)}(X)\right]^2p_i^{(s)}(X),$$
(13)   
is proportional to the volume fluctuations over the ensemble of granular systems considered. Then, the compactivity $`X`$ should be an increasing function of the vibration intensity $`\mathrm{\Gamma }`$, because it is reasonable to expect that the powder became fluffier with higher vibration intensities. This is, for instance, the behavior found in the simple two-volume model considered in Ref. . In general, this property must be checked once the relationship between $`\mathrm{\Gamma }`$ and $`X`$, $`X=f(\mathrm{\Gamma })`$, has been derived for each specific model.

Of course, $`p_i^{(s)}`$ is not a solution of the master equation when the intensity $`\mathrm{\Gamma }`$ is time-dependent, and in general the system does not monotonically approach the steady distribution. Nevertheless, define \[compare with Eq. (2)\]

$$H^{(s)}(t)=\underset{i}{}p_i(t)\mathrm{ln}\frac{p_i(t)}{p_i^{(s)}(X)},$$
(14)   
where $`p_i(t)`$ is again one solution of Eq. (1), and $`p_i^{(s)}(X)`$ depend on time through the compactivity $`X=X(t)`$. If we define a statistical entropy as 

$$S^{}(t)=\lambda \underset{i}{}p_i(t)\mathrm{ln}p_i(t),$$
(15)   
and use the notation

$$\overline{V}(t)=\underset{i}{}V_ip_i(t),$$
(16)   
for the actual average volume at time $`t`$, after a very simple algebra it is found that

$$H^{(s)}(t)=\frac{1}{X}\left[\left(\overline{V}(t)XS^{}(t)\right)\left(\overline{V}^{(s)}(X)XS(X)\right)\right].$$
(17)   
Thus $`H^{(s)}(t)`$ is proportional to the deviation of the actual “effective” volume  at time $`t`$,

$$Y(t)=\overline{V}(t)XS^{}(t),$$
(18)   
from its stationary value.

The time variation of $`H^{(s)}`$ is easily obtained as 

$$\frac{dH^{(s)}}{dt}=A^{(s)}(t)\underset{i}{}\frac{p_i(t)}{p_i^{(s)}(X)}\frac{dp_i^{(s)}(X)}{dX}\frac{dX}{dt},$$
(19)   
where $`A^{(s)}(t)`$ is given by Eq. (4), but replacing $`p_i^{}(t)`$ by $`p_i^{(s)}(X)`$. Taking into account Eqs. (8-10), it is found

$$\frac{dH^{(s)}}{dt}=A^{(s)}(t)\frac{1}{X^2}\frac{dX}{dt}\left[\overline{V}(t)\overline{V}^{(s)}(X)\right].$$
(20)   
Equation (20) does not have a well-defined sign and, therefore, the stationary distribution is not monotonically approached, in general, when the vibration intensity is time-dependent. However, in those processes such that the compactivity (or, equivalently, the vibration intensity $`\mathrm{\Gamma }`$) increases monotonically in time, only the term

$$B(t)=\frac{1}{X^2}\frac{dX}{dt}\overline{V}^{(s)}(X)$$
(21)   
is positive. Using the analogy between temperature and compactivity, we will refer to those processes as “heating” processes. If in a given “heating” process it is verified that

$$\underset{t\mathrm{}}{lim}B(t)=0,$$
(22)   
we can conclude

$$\underset{t\mathrm{}}{lim}\frac{dH^{(s)}(t)}{dt}0.$$
(23)   
Since $`H^{(s)}(t)`$ is bounded below, the only possibility is in fact the equality, which implies

$$\underset{t\mathrm{}}{lim}A^{(s)}(t)=0$$
(24)   
and

$$\underset{t\mathrm{}}{lim}\frac{1}{X^2}\frac{dX}{dt}\overline{V}(t)=0.$$
(25)   
From Eq. (24) and the irreducibility of the stochastic process, it follows that in the long time limit,

$$\underset{t\mathrm{}}{lim}p_i(t)=\underset{t\mathrm{}}{lim}p_i^{(s)}(X),$$
(26)   
i. e., the system goes to the steady curve for long enough times, if the “heating” program verifies the condition given by Eq. (22).

The following picture emerges for the evolution of a granular system submitted to vibrations with increasing intensity: Starting from an arbitrary initial condition, in a first step the system tends to a behavior which is determined by the law of variation of the vibration intensity $`\mathrm{\Gamma }`$, and the initial condition has been forgotten. This implies the existence of a special solution of the master equation, called the “normal” solution, that is approached by all the other solution after an initial transient period. Besides, if the system is “heated” slowly, in the sense that Eq. (22) is verified, the normal solution tends afterwards to the stationary curve. This picture is similar to the one found in some models of structural glasses .

## III Hilbert’s method around the steady curve

In this Section we will use Hilbert’s method to derive a quite general form of the normal solution of Eq. (1) near the steady state curve. We will focus on the class of models for granular media considered in the previous section, but the results can be directly extended to any system with a “canonical” distribution describing the stationary state.

Let us consider Eq. (1), rewritten in the form 

$$\frac{d𝒑(t)}{dt}=\widehat{𝑾}(t)𝒑(t),$$
(27)   
where $`𝒑`$ is a vector (column matrix) whose elements are the probabilities $`p_i(t)`$ of the $`i`$-th state of the system at time $`t`$, and $`\widehat{𝑾}`$ is a square matrix with elements $`\widehat{W}_{ij}`$ given by

$$\widehat{W}_{ij}(t)=W_{ij}(t)\delta _{ij}\underset{k}{}W_{kj}(t).$$
(28)   
If the transition rates are time independent and the detailed balance condition is verified, solving Eq. (27) is equivalent to obtain the solution of the eigenvalue problem

$$\widehat{𝑾}𝝋(q)=\lambda (q)𝝋(q),$$
(29)   
with $`\lambda (q)>0`$ for all $`q`$. The eigenvectors $`𝝋(q)`$ are completed with the stationary distribution $`𝒑^{(s)}`$, which is an eigenvector of $`\widehat{𝑾}`$ corresponding to the null eigenvalue, i. e.,

$$\widehat{𝑾}𝒑^{(s)}=0.$$
(30)   
Besides, if the Markov process is irreducible there is only one stationary distribution with all its components positive . The matrix $`\widehat{𝑾}`$ is hermitian with the following definition for the scalar product of any two vectors $`𝒂`$ and $`𝒃`$,

$$(𝒂,𝒃)=\underset{i}{}\frac{a_ib_i}{p_i^{(s)}}.$$
(31)   
If the Markov process is irreducible and the detailed balance condition is fulfilled for all times, Eqs. (29-31) remain valid for time dependent transition rates. Of course, the eigenvalues and eigenvectors will depend on time in general. The usual situation is that the transition rates depend on time through an externally controlled parameter like the temperature in a thermal system or the compactivity $`X`$ in a granular medium. Thus, we will write sometimes in the following $`\widehat{𝑾}(X)`$, $`\lambda (q,X)`$, $`𝝋(q,X)`$ and $`𝒑^{(s)}(X)`$.

Hilbert’s method consists in solving the master equation by means of the iterative process

$$\widehat{𝑾}(t)𝒑^{(0)}(t)=0,$$
(33)   
$$\widehat{𝑾}(t)𝒑^{(n)}(t)=\frac{d𝒑^{(n1)}(t)}{dt},n1.$$
(34)   
In this way we obtain a probability distribution

$$𝒑_H(t)=\underset{n=0}{\overset{\mathrm{}}{}}𝒑^{(n)}(t),$$
(35)   
which is a solution of Eq. (27). The first term $`𝒑^{(0)}(t)`$ gives the “stationary” distribution,

$$p_i^{(0)}(t)=p_i^{(s)}(X),$$
(36)   
where $`X`$ stands for the value of the compactivity at time $`t`$, $`X=X(t)`$. The Hilbert solution $`𝒑_H`$, constructed following the above rules, is a “normal” solution, because it only depends on the external law of variation of the transition rates, and it does not refer to any specific initial conditions. Nevertheless, the range of validity of $`𝒑_H(t)`$ is limited in general because of the divergence of Hilbert’s expansion . From a physical point of view, this divergence is connected to the fact that we are expanding $`𝒑(t)`$, which may describe a situation arbitrarily far from the steady state, around the stationary solution $`𝒑^{(s)}(X)`$.

In “heating” processes we have shown in Sec. II that a normal solution of the master equation exists and, under very general conditions, it tends to the stationary solution for very high compactivity. Thus in the high compactivity limit Hilbert’s method can be useful, since it is an expansion around the steady state. We will restrict ourselves to the first correction $`𝒑^{(1)}(t)`$, because the difference between the probability distributions $`𝒑_H(t)`$ and $`𝒑^{(s)}(X)`$ is expected to be small in the high compactivity limit. The normalization of $`𝒑^{(s)}(X)`$ implies that

$$(𝒑^{(s)}(X),\frac{d𝒑^{(s)}(X)}{dX})=\underset{i}{}\frac{dp_i^{(s)}(X)}{dX}=0,$$
(37)   
and, therefore,

$$𝒑^{(1)}(t)=\widehat{𝓣}(X)\frac{d𝒑^{(s)}(X)}{dX}\frac{dX}{dt},$$
(38)   
where $`\widehat{𝓣}(X)`$ is the inverse operator of $`\widehat{𝑾}(X)`$ in the space orthogonal to the equilibrium distribution $`𝒑^{\mathbf{(}𝒔\mathbf{)}}(X)`$. Using Eq. (8), we obtain

$$\frac{dp_i^{(s)}(X)}{dX}=\frac{1}{X^2}\left[V_i\overline{V}^{(s)}(X)\right]p_i^{(s)}.$$
(39)   
By introducing a function

$$\xi (q,X)=\underset{i}{}V_i\phi _i(q,X),$$
(40)   
where $`𝝋(q,X)`$ is the eigenvector defined in Eq. (29), we can rewrite Eq. (39) in a vectorial form as

$$\frac{d𝒑^{(s)}(X)}{dX}=\frac{1}{X^2}\underset{q}{}\xi (q,X)𝝋(q,X).$$
(41)   
The above expression follows from the completeness of the eigenvectors and the property

$$(𝝋(q,X),\frac{d𝒑^{(s)}(X)}{dX})=\frac{1}{X^2}\underset{i}{}\phi _i(q,X)\left[V_i\overline{V}^{(s)}(X)\right]=\frac{1}{X^2}\xi (q,X).$$
(42)   
The term proportional to $`\overline{V}^{(s)}(X)`$ vanishes as a consequence of the orthogonality of the eigenvectors $`𝒑^{(s)}(X)`$ and $`𝝋(q,X)`$ for all $`q`$. Then, to first order in deviations from the equilibrium curve we have

$$𝒑_H(t)=𝒑^{(s)}(X)\frac{1}{X^2}\frac{dX}{dt}\underset{q}{}\lambda ^1(q,X)\xi (q,X)𝝋(q,X).$$
(43)   
The notation used reflects the fact that the right hand side depends on time only through the compactivity $`X`$ of the system.

From Eq. (43) it is possible to evaluate the average values of the physical quantities of the system over Hilbert’s distribution in the first order approximation. For instance, the mean value of the volume is

$$\overline{V}_H(t)=\overline{V}^{(s)}(X)\frac{1}{X^2}\frac{dX}{dt}\underset{q}{}\lambda ^1(q,X)\xi ^2(q,X),$$
(44)   
This expression can be written in a more transparent way by using the following identity for the powder “compressibility” defined in Eq. (13),

$$\kappa (X)=\frac{1}{X^2}\underset{q}{}\xi ^2(q,X),$$
(45)   
and the expression for the average relaxation time of the volume found in linear relaxation theory (see the appendix)

$$\tau (X)=\frac{_q\lambda ^1(q,X)\xi ^2(q,X)}{_q\xi ^2(q,X)}.$$
(46)   
In this way, it is found

$$\overline{V}_H(t)=\overline{V}^{(s)}(X)\frac{dX}{dt}\kappa (X)\tau (X).$$
(47)   
It is important to note that the detailed balance condition is a key point in our derivation of Eq. (47), while the existence of the normal solution does not need detailed balance to be satisfied, but it follows if the Markov process is irreducible in the long time limit. On the other hand, Eq. (47) can also be applied to estimate the departure from the stationary curve in “cooling” processes, if the system is initially prepared in the stationary state. In other words, Eq. (47) should be applicable to any situation near the stationary curve. Therefore, it allows to calculate the first order in the deviation of the volume from its steady value over the normal curve when the system asymptotically tends to the stationary curve, as it is the case of “heating” processes. The same remark can be made about Eq. (43) for Hilbert’s distribution $`𝒑_H(t)`$. It contains all the information needed to evaluate one-time properties over the normal curve in situations close to the steady state.

## IV The model

In this Section we will briefly review a simple model for the vibrocompaction of a dry granular system which has been recently introduced . We consider a one-dimensional lattice with $`N`$ sites. Each site can be either occupied by a particle or empty, i. e., occupied by a hole, with the restriction that there cannot be two nearest neighbor holes (such a configuration would be unstable). A variable $`m_i`$ is assigned to each site $`i`$, taking the value $`m_i=1`$ if the site is empty, and $`m_i=0`$ if there is a particle on it. Then, a configuration of the system is fully specified by giving the values of the set of variables $`𝒎\{m_1,m_2,\mathrm{},m_N\}.`$

The dynamics of the system is defined as a Markov process, and formulated by means of a master equation for the probability $`p(𝒎,t)`$ of finding the system in the configuration $`𝒎`$ at time $`t`$ ,

$$\frac{d}{dt}p(𝒎,t)=\underset{𝒎^{}}{}\left[W(𝒎|𝒎^{})p(𝒎^{},t)W(𝒎^{}|𝒎)p(𝒎,t)\right],$$
(48)   
where $`W(𝒎|𝒎^{})`$ is the transition rate from state $`𝒎^{}`$ to state $`𝒎`$. The possible transitions in the system can be classified in three groups. Firstly, there are transitions conserving the number of particles, corresponding to purely diffusive processes. Their rates are given by

$$W(010|100)=W(010|001)=\frac{1}{2}\alpha ,$$
(49)   
with $`\alpha `$ being a constant. These transition rates must be understood as it is usually done, as the transition rates between states connected by the given rearrangement. Only the variables of the set of sites involved in the transition are indicated in the notation. Secondly, there are also transitions increasing the number of particles, with rates

$$W(010|101)=\frac{1}{2}\alpha ,$$
(51)   
$$W(001|101)=W(100|101)=\frac{1}{4}\alpha .$$
(52)   
Finally, the transition rates for those processes decreasing the number of particles are

$$W(01010|00100)=\frac{1}{2}\alpha ^2,$$
(54)   
$$W(01010|01000)=W(01010|00010)=\frac{1}{4}\alpha ^2.$$
(55)   
The transition rates in Eqs. (49)-(52) define an effective dynamics for tapping processes taking place in a previous more general model presented in Ref. . The parameter $`\alpha `$ characterizes the tapping process completely. Note that we have rescaled all the expressions of the transition rates in Ref.  in such a way that the time scale $`t`$ is a measure of the number of taps. For $`\alpha =0`$ no transition is possible in the system and, therefore, the parameter $`\alpha `$ measures the intensity of the vibration. The configuration with the highest density, i. e., no holes present, is completely isolated from the rest of the states, in the sense that no transition is possible from or towards it for any value of $`\alpha `$. Aside from this particular state, all the other possible configurations are connected through a chain of transitions with non-zero probability for $`\alpha 0`$, and the Markov process defining the dynamics is then irreducible.

When $`\alpha `$ does not depend on time, i. e., the intensity of vibration is constant, there is a unique stationary solution $`p^{(s)}(𝒎)`$ of the master equation (48). We will use the notation $`𝒎^{\mathbf{(}𝒌\mathbf{)}}`$ for a configuration of the system with $`k`$ holes. Clearly, it is verified that $`1k(N+1)/2`$. It is easily found that 

$$p^{(s)}(𝒎^{\mathbf{(}𝒌\mathbf{)}})=\frac{e^{k/X}}{Z},$$
(56)   
where $`Z`$ is the normalization constant, and $`X`$ is related to the vibration intensity $`\alpha `$ through the relation

$$\alpha =e^{1/X}.$$
(57)   
Comparison with Eq. (8) shows that in our model the number of holes $`k`$ plays the role of the volume of the system and $`X`$ is the compactivity of Edward’s theory of powders . More precisely, $`k`$ would be proportional to the excess volume from the densest state. Interestingly, the compactivity $`X`$ is an increasing function of $`\alpha `$, as it has been argued on general basis in Section II. Besides, the compactivity vanishes for the no vibrated case. The partition function can be analytically derived in the infinite system limit, with the result 

$$\mathrm{ln}\zeta \frac{1}{N}\mathrm{ln}Z=\mathrm{ln}2+\mathrm{ln}\left[1+\left(1+4\alpha \right)^{1/2}\right].$$
(58)   
From here, the stationary value of the density of holes is obtained in the standard way,

$$x_1^{(s)}=\frac{\overline{k}^{(s)}}{N}=\frac{d\mathrm{ln}\zeta }{d(1/X)}=\frac{1}{2}\left[1\left(1+4\alpha \right)^{1/2}\right].$$
(59)   
The stationary density of holes increases monotonically from the densest state $`x_1^{(s)}=0`$ to the fluffiest state $`x_1^{(s)}=1/2`$ as the vibration intensity increases monotonically from $`\alpha =0`$ to $`\alpha =\mathrm{}`$, since the analogous to the “compressibility”

$$\kappa =\frac{dx_1^{(s)}}{dX}=\frac{e^{1/X}}{X^2\left(1+4e^{1/X}\right)^{3/2}}$$
(60)   
is positive definite for any value of $`X`$. Nevertheless, the densest configuration cannot be actually reached, since for $`\alpha =0`$ no transition is possible and the system gets trapped in its initial state. As the number of holes $`k`$ is an upper bounded variable, the compactivity $`X`$ can take negative values, corresponding to configurations fluffier than those of positive compactivities. For instance, $`X0^{}`$ gives $`\alpha \mathrm{}`$ and $`x_1^{(s)}=1/2`$, the least dense state.

Finally, let us qualitatively study the dynamics of the system for time-independent transition rates in the low vibration limit, $`\alpha 1`$ or $`X0^+`$. In that limit, the stationary value of the density of holes $`x_1^{(s)}`$ is very small, and a typical configuration of the system consist of a few holes, separated by long arrays of particles. Since $`x_1^{(s)}\alpha `$, the mean distance between holes is of order $`\alpha ^1`$. Moreover, at least in linear relaxation theory and at low compactivities, the evolution of the system will be mainly associated to diffusive processes. The characteristic time $`\widehat{\tau }`$ of the relaxation process would be the square of the mean distance between holes divided by an effective diffusion coefficient, which can be estimated from Eq. (49) as of the order $`\alpha `$. Therefore,

$$\widehat{\tau }\alpha ^3.$$
(61)   
This result will be important in the following, since an estimate of the relaxation time in the linear relaxation regime is necessary when studying the time-dependent rates case, as it follows from Eq. (47).

## V Processes with time-dependent vibration intensity

Next we will study the dynamical behavior of the model described in the previous section, when submitted to processes in which the vibration intensity $`\alpha `$ changes in time in a given way. Due to the relationship between $`\alpha `$ and the compactivity $`X`$, Eq. (57), this is equivalent to consider that the compactivity $`X`$ varies in time following a given law. As already mentioned, we will refer to a process as a “heating” (“cooling”) one when the vibration intensity is monotonically increased (decreased). Such kind of processes have been already considered in the literature, both in real granular systems  and in simple models .

The general results of Sec. II are applicable to this particular model. The only limitation is due to the loss of irreducibility of the dynamics for $`\alpha =0`$. Thus the existence of a special solution of the master equation, such that all the others approach it in the long time limit, applies to any program of variation of $`\alpha `$ except for “cooling” processes up to $`X=0`$. As a consequence, for “heating” processes there will be a special “normal” curve, to which all the other solutions tend at a first stage. Later on, the system will approach the “stationary” distribution $`p^{(s)}[\alpha (t)]`$ in the long time limit, provided that the condition in Eq. (22) is verified, i. e.,

$$\underset{t\mathrm{}}{lim}\frac{1}{X^2(t)}\frac{dX(t)}{dt}x_1^{(s)}\left[\alpha (t)\right]=0.$$
(62)   
Taking into account that $`x_1^{(s)}`$ is upper bounded by $`1/2`$, and the relationship between $`\alpha `$ and $`X`$, Eq. (57), the above condition can also be written as

$$\underset{t\mathrm{}}{lim}\frac{d\mathrm{ln}\alpha (t)}{dt}=0.$$
(63)   
Equation (63) expresses a restriction for the “heating” programs driving the system to the stationary curve in the long time limit, but it does not affect at all the existence of the normal solution, which only depends on the ergodicity of the process.

Application of Eq. (47) to the present model yields

$$x_1(t)=x_1^{(s)}[X(t)]\frac{dX(t)}{dt}\kappa [X(t)]\tau [X(t)]+\mathrm{},$$
(64)   
where $`\tau (X)`$ is the mean relaxation time of the density in the linear relaxation approximation. Upon writing the above expression we have taken into account that in our model the number of holes plays the role of the volume. For high compactivities the second term on the right hand side of Eq. (64) is negligible against the first one, since the compressibility $`\kappa 0`$ when $`X\mathrm{}`$. Then, for high vibration intensities the system remains over the stationary curve. As the compactivity decreases, the system departs from that curve, as a consequence of the increase of the mean relaxation time $`\tau `$, which is expected to be proportional to the characteristic time $`\widehat{\tau }`$ estimated in Eq. (61).

There is some freedom when choosing the law of variation of the vibration intensity $`\alpha `$. Our choice will be motivated by simplicity, but also by the analogies with a glass-like behavior previously found in granular systems . In supercooled liquids the temperature is usually varied at a constant rate , so we have considered processes in which the compactivity changes linearly in time,

$$\frac{dX}{dt}=\pm r,$$
(65)   
with $`r>0`$, which is equivalent to

$$\frac{d\alpha }{dt}=\pm r\alpha \left(\mathrm{ln}\alpha \right)^2,$$
(66)   
the plus sign corresponding to “heating” processes and the minus sign to “cooling” programs.

The rest of this section is organized as follows. First, we study “cooling” processes. The system is initially put in the stationary state corresponding to a given value $`\alpha _0`$ of the vibration intensity. Then, the compactivity $`X=1/\mathrm{ln}\alpha `$ is decreased following Eq. (65). The existence of a phenomenon analogous to the laboratory glass transition of supercooled liquids arises. Afterwards, “heating” processes are considered, paying special attention to the appearance of hysteresis effects, and relating them to the trend of the system to approach the normal curve.

### A “Cooling” processes

We consider the continuous decreasing of the compactivity of the system from a given initial value $`X_0`$ down to $`X=0`$. The latter corresponds to $`\alpha =0`$, i. e., no vibration. The system is initially placed in the stationary state corresponding to the value $`\alpha _0=\mathrm{exp}(1/X_0)`$ of the vibration intensity. Then, the compactivity is decreased following the law

$$\frac{dX}{dt}=r_c.$$
(67)   
Our starting point will be Eq. (64), particularized for the “cooling” process we are considering, i. e.,

$$x_1(t)=x_1^{(s)}[X(t)]+r_c\kappa [X(t)]\tau [X(t)]+O(r_c^2).$$
(68)   
As we have already discussed, from the above equation follows that the system remains in equilibrium at “high” compactivities, since the second term in the right hand side vanishes for $`X\mathrm{}`$. Nevertheless, as the compactivity becomes smaller this term grows, due to the increase of the relaxation time $`\tau `$. This means that there will exist a range of values of the compactivity in which the second term is comparable to the first one. A rough estimate of the value of the compactivity $`X_g`$ at which the system would depart from the stationary curve can be obtained by equalling both terms. If we consider that the system is slowly “cooled”, $`r_c1`$, it is also $`\alpha _g=e^{1/X_g}1`$, and we can approximate both terms for their leading behaviors in that limit, namely

$$x_1^{(s)}(X)\alpha ,$$
(70)   
$$\kappa (X)=\frac{dx_1^{(s)}}{dX}\alpha \left(\mathrm{ln}\alpha \right)^2$$
(71)   
$$\tau (X)\tau _0\alpha ^3$$
(72)   
where $`\tau _0`$ is a constant of the order of unity. Therefore, we get

$$r_c\tau _0\left(\mathrm{ln}\alpha _g\right)^2=\alpha _g^3.$$
(73)   
For $`r_c1`$ it is

$$\alpha _gr_c^{1/3}\left|\mathrm{ln}r_c\right|^{2/3}.$$
(74)   
We have omitted the factors containing $`\tau _0`$, as well as any other factor of the order of unity.

For $`\alpha <\alpha _g`$ the system is effectively “frozen”, due to the divergent tendency of the relaxation time no more transitions are possible, and the density of holes will be approximately constant in this region. A measure of the effective number of transitions left to the system before reaching $`\alpha =0`$ from a given time $`t`$ is given by the scale 

$$s(t)=_t^{t_0}𝑑t^{}\frac{1}{\tau [X(t^{})]},$$
(75)   
where $`t_0`$ is the time instant for which the “cooling” program finishes, i.e. $`\alpha (t_0)=0`$. In that way, the system would get frozen for $`t>t_f`$ such that $`s(t_f)=1`$, being easily obtained that $`\alpha (t_f)\alpha _g`$. Using Eq. (74) it is possible to estimate the leading order value of the compactivity at which the system gets frozen,

$$X_g=\frac{1}{\mathrm{ln}\alpha _g}\frac{3}{|\mathrm{ln}r_c|}.$$
(76)   
This kind of behavior has been also found numerically in a granular system model . The inverse logarithmic dependence of the cooling rate is typical for the laboratory glass transition temperature of supercooled liquids , and it has also been analytically derived in some simple models of structural glasses .

Since the density of holes remains nearly constant for $`\alpha <\alpha _g`$, and the two first terms of Hilbert’s expansion (68) are of the same order for $`\alpha \alpha _g`$, it is reasonable to expect that

$$x_{1,\text{res}}=\underset{tt_0}{lim}x_1(t)x_1^{(s)}(\alpha _g)=\alpha _gr_c^{1/3}\left|\mathrm{ln}r_c\right|^{2/3},$$
(77)   
where $`x_{1,\text{res}}`$ is the residual value of the density of holes, extending again to this model of granular system the terminology of structural glasses.

In order to check the above results, Fig. 1 shows the residual value of the density of holes as a function of the “cooling” rate, measured by the parameter

$$\delta =r_c\left(\mathrm{ln}r_c\right)^2,$$
(78)   
in a log-log scale. The numerical result agrees with the theoretical prediction, Eq. (77), since the curve is well fitted by a straight line with a slope approximately equal to $`1/3`$. In Fig. 2 the evolution of the density of particles, $`\rho =1x_1`$, in a “cooling” process with rate $`r_c=10^5`$ is plotted. For comparison, the equilibrium curve, given by Eq. (59), is also shown. The estimate of the freezing compactivity from Eq. (76) is $`X_g0.26`$, which is seen to be in good agreement with the region in which the Monte Carlo density is approximately constant. Similar behaviors have been observed for other small values of the cooling rate.

### B “Heating” processes and hysteresis effects

Let us analyze “heating” processes, i. e., processes in which the vibration intensity is monotonically increased. In this kind of processes, the dynamics of the system is irreducible. Thus there is a “normal” solution of the master equation, such that any other solution tends to it. Besides, if the heating program verifies Eq. (63), the normal solution approaches the stationary curve for large enough times. We are going to discuss how these results can be applied to our model, in order to understand its behavior when “heated” from $`\alpha =0`$.

The compactivity will be increased according to the law

$$\frac{dX}{dt}=r_h,$$
(79)   
where $`r_h`$ is the rate for this process. According to Eq. (64), the system will approach the stationary curve in the high compactivity (long time) limit. In the vicinity of the stationary curve, the evolution of the system is given by

$$x_1(t)=x_1^{(s)}[X(t)]r_h\kappa [X(t)]\tau [X(t)]+O(r_h^2).$$
(80)   
This equation explains why, in “heating processes”, the system tends to the stationary curve following a curve different from the “cooling” one, even when $`r_c=r_h`$. It follows directly from the comparison of Eqs. (80) and Eq. (68), by noting that the deviation from the stationary behavior is of opposite signs in “heating” and “cooling” processes. Therefore, hysteresis effects show up.

However, perhaps the main result for “heating” processes is the existence of the normal solution. Fig. 3 shows the evolution of the density of particles in a heating process with $`r_h=10^5`$. The shape of the curve depends on the initial condition for $`X=0`$. Two different initial preparations of the system have been considered. The system was previously cooled down to $`X=0`$, following two different linear programs with $`r_c=10^5`$ and $`r_c=10^3`$, respectively. For the sake of clarity, these “cooling” curves are not shown. From Fig. 3 it is seen that both “heating” curves tend to a common behavior and, afterwards, they approach the stationary curve for high compactivities. Also plotted is the normal curve, which was obtained by starting from the loosest packing state, $`x_1=0.5`$ .

Figure 4 depicts a particular cycle of “cooling” and “heating” with the same rate, namely $`r_c=r_h=10^5`$, as well as the normal curve of the “heating” process. A behavior similar to the one found in real granular systems , and also in the “Tetris” model , is observed. When starting the heating process from the loosest packing state the normal curve is obtained, which tends to the stationary behavior in the high vibration intensity limit. Afterwards, cooling and reheating with the same rate leads to the other two curves of the figure. These are approximately “reversible” for very small rates, since the deviation from the stationary curve is smaller the smaller the rate. Nevertheless, they cannot be used to obtain the stationary values of the density for low compactivities, due to the glass-like kinetic transition. On the other hand, at high compactivities the deviations of the system from the stationary curve for the “cooling” and the “heating“ processes are symmetric, as predicted by Eqs. (68) and (80).

## VI Final remarks

In the first part of the paper we have considered a wide class of models of granular systems, namely those models whose dynamics under tapping is described by means of a master equation. It has been assumed that the tapping process is such that the system is able to explore the whole space of metastable states, when it is tapped with any vibration intensity $`\mathrm{\Gamma }0`$. Then, if $`\mathrm{\Gamma }`$ changes in time, and the Markov process remains irreducible in the long time limit, all the solutions of the master equation tend to approach a special solution, called the “normal” solution for the given vibration program. Besides, if the stationary distribution for constant $`\mathrm{\Gamma }`$ is consistent with Edward’s statistical mechanics theory of powders, the normal solution tends to the stationary curve in the high vibration intensity regime, provided that the “heating” program is not too fast, in the sense given by Eq. (22). It has also been argued that the compactivity $`X`$ of Edward’s theory should be an increasing function of the vibration intensity $`\mathrm{\Gamma }`$. A quite general prediction for the behavior of the normal solution in situations near the steady state has been found, by using Hilbert’s method to solve the master equation. This result can be applied to any system having a “canonical” steady state distribution.

Very recently, a simple model for a granular system submitted to vertical vibration with an stationary state described by Edward’s theory of powders has been introduced . The model belongs to the class of systems described in the previous paragraph, and here we have investigated its behavior in processes in which the vibration intensity depends on time. Such kind of processes have been carried out in real granular systems , and also investigated in some models . For the sake of concreteness we have taken the compactivity as a linear function of time.

The behavior of the model under “cooling” processes up to zero vibration intensity exhibits a phenomenon similar to the laboratory glass transition. The “cooling” evolution departs from the stationary curve, and freezes, in a narrow region around the value of the compactivity $`X_g`$ such that the effective number of transitions left to the system until reaching $`X=0`$ becomes smaller than unity. This allows us to estimate the dependence of $`X_g`$ and the residual values of the density upon the “cooling” rate. The results are analogous to those obtained previously in simple models of structural glasses .

In “heating” processes, a crucial role is played by the “normal” solution of the master equation, which is completely determined by the program of variation of the compactivity, and attracts any other solution, independently of the initial condition. The hysteresis effects found when the system is “cooled” and “reheated” are related to the trend of the system to approach the normal curve. The normal curve corresponds to “heating” the system from the loosest state. In a first stage, the system tends to approach monotonically the “normal” curve and, for longer times, corresponding to high enough compactivities, the stationary curve is reached if the “heating” program verifies Eq. (63).

The work reported here suggests a very close relationship between structural glasses and the behavior of a granular system under vibrations, supporting previous results in the same direction . It also raises some questions, both from a theoretical and experimental point of view. We think that it is worth checking the existence of the normal solution in real granular systems, and whether it can be constructed starting from the loosest packing state. Also, it seems interesting to know if the glass analogy also extends to the linear relaxation regime, giving a stretched exponential behavior of the response functions, or if the inverse logarithmic relaxation law  remains valid in such a region. The latter possibility would lead to the need of investigating the reason for such a behavior.

###### Acknowledgements.

This research was partially supported by the Dirección General de Investigación Científica y Técnica (Spain) through Grant No. PB98–1124.    

## A Average relaxation time of the volume

Linear relaxation describes the evolution of the system at constant value of the compactivity $`X`$, when starting from the steady state corresponding to a compactivity $`X+\mathrm{\Delta }X`$, with $`\mathrm{\Delta }X`$ very small, in the sense that we can approximate

$$𝒑^{(s)}(X+\mathrm{\Delta }X)=𝒑^{(s)}(X)+\frac{d𝒑^{(s)}(X)}{dX}\mathrm{\Delta }X.$$
(A1)   
Taking into account Eq. (41), the evolution of the probability distribution from this initial condition, $`𝒑(t=0)=𝒑^{(s)}(X+\mathrm{\Delta }X)`$, is

$$𝒑(t)=𝒑^{(s)}(X)+\frac{\mathrm{\Delta }X}{X^2}\underset{q}{}\xi (q,X)𝝋(q,X)e^{t\lambda (q,X)},$$
(A2)   
because $`𝝋(q,X)`$ is the eigenvector of the evolution operator corresponding to the eigenvalue $`\lambda (q,X)`$. From the above expression it is straightforward to calculate the time evolution of the average volume,

$$\overline{V}(t)\overline{V}^{(s)}(X)=\frac{\mathrm{\Delta }X}{X^2}\underset{q}{}\xi ^2(q,X)e^{t\lambda (q,X)},$$
(A3)   
where we have made use of the definition of the function $`\xi (q,X)`$, Eq. (40).

The volume linear relaxation function corresponding to a value of the compactivity $`X`$ is defined in the usual way,

$$\varphi _V(t;X)\frac{\overline{V}(t)\overline{V}^{(s)}(X)}{\overline{V}(t=0)\overline{V}^{(s)}(X)},$$
(A4)   
yielding

$$\varphi _V(t;X)=\frac{_q\xi ^2(q,X)e^{t\lambda (q,X)}}{_q\xi ^2(q,X)}.$$
(A5)   
On the other hand, the linear relaxation time for the volume $`\tau (X)`$ is given by definition by the area under the curve $`\varphi _V(t;X)`$ as a function of $`t`$. This leads to Eq. (46).