- Email: [email protected]
Thermodynamic description in a simple model for granular compaction
Physica A 275 (2000) 310–324
www.elsevier.com/locate/physa
Thermodynamic description in a simple model for granular compaction a Fsica
J. Javier Brey a; ∗ , A. Pradosa , B. Sanchez-Rey b
Teorica, Facultad de Fsica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, Spain b Fsica Aplicada, E.U. Politecnica, Universidad de Sevilla, Virgen de Africa 7, E-41011 Sevilla, Spain Received 18 July 1999
Abstract A simple model for the dynamics of a granular system under tapping is studied. The model can be considered as a particularization for short taps of a more general one-dimensional lattice model with facilitated dynamics. The steady state reached by the system is discussed and the results are shown to be consistent with the thermodynamic granular theory developed by Edwards and coworkers. In particular, the basic assumption of the theory, i.e., that the probability distribution c 2000 Elsevier Science B.V. depends only on the volume of the con guration, is veri ed. All rights reserved. PACS: 81.05.Rm; 05.50.+q; 81.20.Ev Keywords: Granular compaction; Thermodynamics; Tapping processes; Lattice models; Facilitated dynamics
1. Introduction A signi cant body of literature has emerged recently trying to apply the methods of kinetic theory and uid dynamics to the study of granular ows [1]. Also, some e orts have been devoted to formulate a statistical mechanics theory for steady granular media prepared under di erent conditions and, in particular, when submitted to external vibrations [2,3]. These studies have been motivated strongly by a series of experiments carried out to identify the mechanisms leading to compaction [4 – 6]. Granular compaction consists in the evolution from an initial low-density state to one with higher ∗
Corresponding author. Tel.: +3495-462-6558; fax: +3495-461-2097. E-mail address: [email protected] (J.J. Brey)
c 2000 Elsevier Science B.V. All rights reserved. 0378-4371/00/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 9 9 ) 0 0 3 7 5 - 1
J.J. Brey et al. / Physica A 275 (2000) 310–324
311
density as vibrations or, with more generality, some kind of series of external excitations are applied to the system. In actual experiments, the most usual perturbations are discrete, vertical shakes or taps applied to the container. The aspects considered in the studies of granular compaction include both the statics, i.e., the steady state reached by the system in the long time limit, and the dynamics, characterized by a time evolution of the density that is accurately described by an inverse logarithmic law. To investigate the latter aspect a series of models have been introduced [7–13], presenting most of the qualitative features of granular compaction observed in experiments. In particular, the density evolution is well tted by the inverse logarithmic law. Let us mention that most of the results from these models have been obtained numerically, while analytical results are scarce. In this paper we will deal with the description of the steady state reached by a granular system under tapping, by using a very simple one-dimensional lattice model, whose dynamics is a Markov process speci ed by means of a master equation. It was introduced in Ref. [13], where it was shown to exhibit the characteristic inverse logarithmic behavior mentioned above, and some properties of the density of the steady state were discussed. The basic idea here will be to identify an e ective Markov process for the description of the evolution of the system when submitted to tapping. In general, tapping processes are dicult to describe and analyze analytically since they are composed by an alternating series of two di erent dynamical evolutions: taps and freely relaxations. The nal state of every tap is the initial state for the next free evolution and vice versa. In the experiments and also in the computer simulations, it has been found useful to describe the evolution under tapping in terms of the number of taps. The values of the properties of the system are recorded at the end of the relaxation following every tap. The characteristic inverse logarithmic law mentioned above appears when the values of the density measured in this way are considered as a function of the number of taps. Thus, it is tempting to try to coarse grain the time description provided by a given model and consider each tap and its following relaxation as a single event. The derivation of the associated equivalent reduced dynamics is not a simple task, especially taking into account that the two processes that must be described as a single one are quite di erent in nature and that the transition between them is not smooth. Nevertheless, we have been able to do it for our model in the limit of very short or weak taps. The e ective tapping process is also formulated in terms of a master equation with well-de ned constant transition rates. Besides, the equation is expressed in terms of a continuous time variable, although it is only meaningful when a time equal to the duration of a given number of taps is considered. A direct consequence of having a master equation with time-independent transition rates to describe the relaxation under tapping is that it is possible to study the existence and properties of the steady con gurations. This allows, in particular, to check whether the model agrees with the thermodynamic description provided by the theory developed by Edwards and his collaborators [2,3], and whose main ingredient is the assumption that the volume of a granular media plays a role
312
J.J. Brey et al. / Physica A 275 (2000) 310–324
analogous to the energy of a statistical molecular system. We will show that our e ective model veri es this principle, which is a nontrivial result. The steady probability distribution only depends on the con guration of the lattice through the associated density. The paper is organized as follows. In Section 2 the model of Ref. [13] is brie y reviewed. This model will be referred to in the following as the original model. Next, it is applied to describe tapping processes in the limit of very short taps, and the relevant elementary events connecting the state of the system before the tap and after the combination of tapping plus free relaxation are identi ed. The resulting process always ends at the lowest order in the duration of the taps with all the sites (except one) occupied by a particle, due to the inability of the system to decrease its density in any elementary transformation. To describe steady states of di erent densities it is then necessary to take into account also some next order processes, namely those increasing the number of empty sites in the lattice. The form of the probability distribution in the steady state is discussed in Section 3. It can be written in the same way as the canonical distribution for a normal statistical system, the role of the energy being played by the volume and that of the temperature by a new parameter, the compactivity. This is the required result to connect with the thermodynamic theory of Ref. [2,3]. Finally, the last section contains some remarks and comments about the main results in the paper. 2. Description of the model The model for tapping processes we are going to formulate is based on a previous model with facilitated dynamics considered recently in Ref. [13]. Let us review the latter brie y. It is a one-dimensional lattice model in which each of the N sites can be either occupied by a particle or can be empty. A variable mi is assigned to each site i, taking the value mi = 1 if the site is empty, while mi = 0 if there is a particle on it. Then, a con guration of the system is fully speci ed by giving the values of the set of variables m ≡ {m1 ; m2 ; : : : ; mN }. The dynamics of the system is assumed to be a Markov process de ned by the following master equation for the transition probability p1=1 (m; t|m0 ; t 0 ) for nding the system with the con guration m at time t given it was in the con guration m0 at time t 0 ¡ t, @t p1=1 (m; t|m0 ; t 0 ) =
N X
[Wi (Ri m)p1=1 (Ri m; t|m0 ; t 0 ) − Wi (m)p1=1 (m; t|m0 ; t 0 )] :
i=1
(1) Here Ri m ≡ {m1 ; : : : ; Rmi ; : : : ; mN }, with Rmi = 1 − mi , and the expressions of the transition rates are Wi (m) = (mi−1 + mi+1 )[ + mi (1 − 2)] ; (2) 2
J.J. Brey et al. / Physica A 275 (2000) 310–324
313
where and are real parameters de ned in the intervals ¿ 0 and 0661, respectively. Thus the transition rate for the adsorption of a particle at site i is Wi+ (m) = (1 − )mi (mi−1 + mi+1 ) (3) 2 and the probability per unit of time for the desorption of a particle from the same site is (4) Wi− (m) = (1 − mi )(mi−1 + mi+1 ) : 2 For =1 there are no adsorption processes and for =0 particles cannot be desorbed. The rational for this system as a model for granular compaction has been discussed in Ref. [13] and will not be reproduced here. We are interested in the particularization of the above general dynamics to tapping processes. They are modeled in the following way. Starting from an arbitrary state or con guration, the system is allowed to relax with = 0, i.e., without desorption events, until it reaches a metastable con guration. Such con guration is characterized by all the holes or empty sites being isolated so that the transition rates given by Eq. (3) also vanish. Therefore, the system cannot evolve any longer and becomes trapped. This metastable con guration will be the initial state to describe the tapping process and will be referred to as state n = 0. Next, the value of is instantaneously increased to a nite value in order to simulate a pulse. Its duration t0 is taken very short as compared with the characteristic time of the transition attempts −1 , i.e., t0 .1. Of course, the tap has in general the e ect of taking out the system out of the metastable con guration. Afterwards, the system relaxes again with = 0 until a new metastable con guration n = 1 is reached. Clearly, evolution with = 0 corresponds to free relaxation in real tapping experiments. By repeating the same sequence of steps a series of metastable states n = 2; 3; 4; : : : is generated. The index n indicates the number of taps or short pulses applied to the system before that con guration has been obtained. Let us consider one of the metastable con gurations of the system after relaxing with = 0. We know that each hole is surrounded by two particles. When the value of is increased, the rst transition occurring in the system must be the desorption of a particle having at least one of its nearest-neighbor sites empty. Therefore, the possible rst elementary events are ··· 1 0 0 ··· → ··· 1 1 0 ··· ; ··· 0 0 1 ··· → ··· 0 1 1 ··· ; ··· 1 0 1 ··· → ··· 1 1 1 ··· :
(5)
We have employed the mi variables to specify the con gurations and have represented only the site whose state is changing and the two nearest-neighbor ones, since they determine the transition rate. For t0 small enough, only one transition can take place at the most in the system during a pulse. Afterwards if it is allowed to relax freely again to a new metastable
314
J.J. Brey et al. / Physica A 275 (2000) 310–324
con guration, we can compute the probability for each of the possible series of transitions connecting two metastable con gurations. For instance, let us consider the rst process shown in Eq. (5). The probability of the desorption event during the pulse is t0 =2, since from Eq. (4) its transition rate is equal to =2. Next, in the free relaxation a particle can be adsorbed with equal probability on any of the two empty sites. Then, this process leads to the same state or to the e ective di usion of the hole to the right. We will not consider those global transitions in which the nal state coincides with the initial one, because they are not necessary when writing a master equation. In a trivial notation, in which only the particles involved in the process are indicated explicitly, the probability for the e ective transition considered is 1
(1) (100 → 010) = t0 × = t0 : (6) 2 2 4 The same line of reasoning allows us to calculate the probabilities for the other e ective transitions between two metastable con gurations, namely
(1) (001 → 010) = t0 ; (7a) 4
(1) (101 → 001) =
t0 ; 8
(7b)
(1) (101 → 010) =
t0 ; 4
(7c)
t0 : (7d) 8 The key point in our approximation is introduced now. It will be assumed that in the limit of very short taps, the time evolution of the system under tapping can be accurately described by a continuous Markov process with transition rates given by
(1) (101 → 100) =
(1)
(m|m0 ) =
N X
(1) 0 i (m|m )
;
(8)
i=1 (1) 0 i (m|m ) =
(i) (1 − m0i ){2[m0i−1 (1 − m0i+1 ) + (1 − m0i−1 )m0i+1 + m0i−1 m0i+1 ] 8 ×(1 − mi−1 )mi (1 − mi+1 ) +m0i−1 m0i+1 [(1 − mi−1 )(1 − mi )mi+1 + mi−1 (1 − mi )(1 − mi+1 )]} ; (9)
where (i) re ects that the e ective transition only a ects three consecutive sites, (i) =
N Y j6=i−1; i; i+1
mj ; m0j :
(10)
J.J. Brey et al. / Physica A 275 (2000) 310–324
315
In summary, the master equation X @t p1=1 (m; t|m0 ; t 0 ) = [ (1) (m|m00 )p1=1 (m00 ; t|m0 ; t 0 ) m00
−
(1)
(m00 |m)p1=1 (m; t|m0 ; t 0 )]
(11)
provides the lowest-order approximation for the description of the compaction process of the original model under tapping provided that t0 .1 and the equation is applied for times t = nt0 , with n being a whole number. The free relaxation with = 0 is considered as instantaneous in this description and assumed to take place just after the transition events associated with the taps. Then, the above master equation does not contain any explicit reference to the unstable con gurations of the system during the tapping. They have been eliminated in the time coarse-graining implicit in Eq. (11). As a consequence, the time scale in the approximated model is just an auxiliary variable, providing the number of taps n by dividing it by t0 . A point to stress is that the approximation we are discussing is not the same as the continuous tapping limit considered by some authors. In the latter the system is not allowed to freely relax between consecutive taps, while here it always relaxes to a metastable con guration after every pulse, and thus the e ective dynamics connects the metastable states among themselves. In order to get some insight into the implications of Eq. (11), let us study the time evolution of the empty site concentration x1 (t), given by X x1 (t) = mi p(m; t) ; (12) m
where p(m; t) ≡
X
p1=1 (m; t|m0 ; 0)
(13)
m0
is the one-time probability distribution. For homogeneous conditions, the right-hand side of Eq. (12) does not depend on the site i considered. Then, a rather simple calculation leads to @x1 = − hmi−1 mi+1 it : (14) @t 2 Here we have used the notation X hf(m)it = f(m)p(m; t) (15) and taken into account that metastable con gurations cannot have three consecutive holes and, therefore, expressions like hmi−1 mi mi+1 it identically vanish in those states. Since hmi−1 mi+1 it ¿0, the time approach of the hole concentration towards its steady value is monotonic in this approximation. This is a direct consequence on the e ective dynamics we have introduced, in which particles can never be desorbed (see Eqs. (6) and (7)). In fact, it is easily seen that all the states with only one empty site are absorbent states [14] for the tapping, connected by di usion processes. It follows that in order to be able to describe tapping evolution leading to steady states
316
J.J. Brey et al. / Physica A 275 (2000) 310–324
Fig. 1. All the elementary trajectories leading to a decrease in the number of holes, with only one desorption event during the tap. We have indicated the processes taking place during the pulse ( ¿ 0) and in the free relaxation ( = 0). The trajectories leading to a nal state being the same as the initial one are not shown.
Fig. 2. All the elementary paths in which the number of holes increases. Two desorption events during the pulse ( ¿ 0) are necessary, while only one transition occurs during the free relaxation ( = 0). Note that the plotted trajectories are just the inverse of the transitions in Fig. 1.
as a consequence of the balance between adsorption and desorption processes, we have to consider higher-order terms in the duration t0 and strength of the individual taps. When considering the next order in t0 contributions, two quite qualitatively di erent kinds of terms appear. There are terms that represent just corrections to nonvanishing lowest-order transition rates. These will be not incorporated in the simpli ed model, since the lowest order is the dominant one in the limit we are considering. On the other hand, some second-order terms provide transition rates for elementary processes that are not present to lowest order. This happens, in particular, for all the transformations leading to the creation of a hole. As discussed above, the incorporation of these processes modi es substantially the physics of the tapping process and they must be retained. In Fig. 1 we have represented the processes leading to the decrease of the number of holes taken into account upon deriving Eq. (11). In each of the elementary events we have indicated whether they take place during the tap ( ¿ 0) or in the free relaxation ( = 0). The question now is whether the inverse processes indicated in Fig. 2 are possible and what are their transition rates. The rst point to realize is that the inverse of the elementary transformations with = 0 are only possible with ¿ 0, since they correspond to desorption processes. Therefore, we need at least two transitions to occur during the tap. Nevertheless, the total number of elementary events involved does not increase since the free relaxation reduces to only one transformation. Although, in principle, it could be thought that there are other second-order transformations, i.e.,
J.J. Brey et al. / Physica A 275 (2000) 310–324
317
with two motions during the tap, increasing the number of holes, it is easy to convince oneself that there are none when trying to construct them. In particular, the group of three particles involved in the transition must have holes on both sides. Then, we have to consider the following e ective elementary processes: 0 1 0 0 0→0 1 0 1 0; 0 0 1 0 0→0 1 0 1 0; 0 0 0 1 0→0 1 0 1 0:
(16)
Accordingly, with the previous discussion we have eliminated those processes in Fig. 2 conserving the number of empty sites. It is important to realize that the transformations kept are just the inverse of the rst-order transformations leading to the loss of a hole. The computation of the e ective transition rates for the processes in Eq. (16) requires to take into account that the two transformations with ¿ 0 must occur in the time interval t0 , the duration of a tap, with probability equal to (t0 )2 =8 in all of the cases shown in Fig. 2. Then the following expressions for the probabilities are obtained:
(2) (01000 → 01010) =
2 2 t02 ; 16
(17a)
(2) (00100 → 01010) =
2 2 t02 ; 8
(17b)
2 2 t02 : (17c) 16 The equality of probabilities (17a) and (17c) is a direct consequence of the isotropy of the dynamics of the original model. The inclusion of these elementary e ective events is a necessary condition to describe steady states with nonvanishing concentration of holes since, as we have already discussed, there are no processes decreasing the particle concentration to the lowest order. To each of the probabilities in Eqs. (17) a transition rate given by (2) = (2) =t0 can be associated, and the corresponding contributions must be incorporated into the expressions of the total transition rate (m|m0 ) in Eq. (8). Nevertheless, this leads to a rather complicated and not very illuminating expression that will not be written here. On the other hand, it is important to note that with the inclusion of the transition rates (2) all the metastable states of the system are connected through a possible chain of transitions. As a test of the validity of the e ective model, as an approximation of the original one in Fig. 3, we have plotted the time evolution, as a function of the number of taps n, of the particle density as predicted by both models. In the simulations periodic boundary conditions and a lattice with 104 sites have been used. In the e ective model n is obtained by dividing the time in units of t0 . The solutions of the equations can be characterized by the parameter
(2) (00010 → 01010) =
= t0 =2 :
(18)
For the curves in Fig. 3 it is = 0:015. It is seen that there is a quite good agreement for all the evolution. As expected, the agreement is not so good as the value of
318
J.J. Brey et al. / Physica A 275 (2000) 310–324
Fig. 3. Evolution of the density of particles with the number of taps n obtained from Monte Carlo simulation of the system, for = 0:015. The diamonds correspond to the e ective model introduced in this paper, while the crosses correspond to the original model of Ref. [7]. The dotted line is the theoretical prediction for the steady-state density of the e ective model, from Eq. (27).
Fig. 4. The same as in Fig. 3, but for = 0:5. Note that the deviations of the e ective model from the original one are larger, but the relative error is still quite small.
increases. An example is given in Fig. 4 obtained with = 0:5. Nevertheless, notice that the relative di erence is still very small. As discussed in Ref. [13], the plotted relaxation curves can be accurately tted by means of an inverse logarithmic law. The dotted lines in Figs. 3 and 4 are the steady values of the density as predicted by the e ective model. They will be discussed in the next section. 3. Steady state The Markov process de ned in the previous section is irreducible and, therefore, the master equation has a unique steady solution [14]. We are going to see that this solution
J.J. Brey et al. / Physica A 275 (2000) 310–324
319
veri es detailed balance and this property can be used to determine it explicitly. Thus, we a priori bet on a time-independent solution ps (m) of the e ective master equation for tapping verifying (m|m0 )ps (m0 ) = (m0 |m)ps (m) :
(19)
If such distribution exists, it will have the property that all the con gurations having the same number of empty sites or holes will have the same probability. This follows from the fact that they can be connected by di usive transitions and di usion is an isotropic process in our system, so that (m|m0 )= (m0 |m) if m and m0 have the same number of holes. Therefore, the value of ps (m) will depend only on the number of holes in the con guration m. In the following we will write m(k) to indicate that this con guration has k empty sites. For processes increasing (decreasing) the number of particles we have p s (m0(k+1) ) = p s (m(k) )
(m0(k+1) | m(k) ) : (m(k) | m0(k+1) )
(20)
A look at Eqs. (6); (7) and (17) shows that (m0(k+1) | m(k) ) t0 = (m(k) | m0(k+1) ) 2
(21)
for all the nonvanishing transition rates. Then, p s (m0(k+1) ) t0 = p s (m(k) ) 2
(22)
which does not depend on the particular con guration under consideration, in agreement with the conclusion reached above that p s (m(k) ) only depends on the value of k. Eq. (22) implies p s (m(k) ) =
k ; Z
(23)
where Z is a normalization constant and is de ned in Eq. (18). For a system with N sites we can formally write Z=
N=2 X k=1
k(N ) k ;
(24)
where k(N ) is the degeneracy of a lattice in a metastable state with k holes, i.e., the number of di erent con gurations having exactly k isolated holes. Note that the maximum number of holes is N=2 (assumed N is even) and the minimum one is unity. The steady hole concentration is N=2 N X X 1 1 X X s (k) @ 1 s s ln Z : (25) mj p (m) = k p (m ) = x1 = N m N @ N (k) j=1
k=1
m
320
J.J. Brey et al. / Physica A 275 (2000) 310–324
Fig. 5. Steady value of the density of particles s as a function of the parameter . The continuous line is Eq. (27), and the diamonds are the Monte Carlo values for the original model in Ref. [13].
The evaluation of the “partition function” Z, which requires previously the calculation of the number of states k(N ) , is sketched in the appendix. In the limit of a very large system (N → ∞), the result is 1 + (1 + 4 )1=2 : 2 Use of this expression into Eq. (25) yields ln Z = N ln
x1s =
(1 + 4 )1=2 − 1 : 2(1 + 4 )1=2
(26)
(27)
As discussed in Section 2, the e ective model we are considering is expected to be equivalent to the original model of Ref. [13] in the limit .1. In this limit it is x1 ' which agrees with the result from the original model [13]. The theoretical prediction given by Eq. (27) is compared with the numerical results obtained by solving numerically the original master equation, Eq. (1) in Fig. 5. The agreement can be considered as quite good for all the range of values of considered. Let us stress here that Eq. (27) does not present any singularity when considered as a function of . The steady distribution (5) can be rewritten in the usual form of a “canonical” probability distribution p s (m(k) ) = Z −1 e−k=X ;
Z=
N=2 X k=1
k(N ) e−k=X ;
(28)
where 1 (29) ln is the analogous of a temperature parameter. It was introduced by Edwards and coworkers in the context of a granular thermodynamic theory and called compactivity [2,3]. X =−
J.J. Brey et al. / Physica A 275 (2000) 310–324
321
The basic assumption of the theory is that the volume plays the same role in a vibrated granular material as the energy in a statistical molecular system. Then, all the con gurations of a given volume are equally probable. We have seen that this is veri ed in our model as a consequence of the dynamics. In the limit of small the steady concentration of particles s will be s ≡ 1 − x1s ' 1 − = 1 − e−1=X ;
(30)
1=X = ln(1 − s ) :
(31)
i.e.,
A relationship similar to this has been recently obtained from the analysis of experimental data [5, See Eq. (6)]. The most compact granular system, i.e., a system with x1s = 0 corresponds to = 0 and X = 0. When X increases also increases and, therefore, the equilibrium density of particles 1 − x s also decreases. This is in complete agreement with the picture on which the thermodynamic theory of powders is based [2,3]. In the limit X → ∞ it is x1s ' 0:28. In principle, negative compactivities must also be considered as we will see below. The limit X → 0− appears for → ∞, and then x1s → 12 , i.e., the least compact arrangement. Nevertheless, we must keep in mind that our model is expected to provide an accurate description of tapping processes in the original model for ¡ 1, and this clearly excludes negative values of X . The analogous of the free energy per site is =−
X ln Z ; N
(32)
so that Eq. (26) is equivalent to x1s =
d(=X ) : d(1=X )
It is also natural to de ne an “entropy” by X N p s (m)ln p s (m) = x1s + ln Z : S≡− X m
(33)
(34)
As expected it turns out to be an extensive property, i.e., proportional to N for large N , so that we can consider a speci c entropy per site = S=N , verifying the relationship = x1s − X :
(35)
Since we know the explicit expression for the partition function of our model, we can compute all the above quantities. In particular, a simple calculation gives =
(1 + 4e−1=X )1=2 − 1 1 + (1 + 4e−1=X )1=2 : + ln −1=X 1=2 2X (1 + 4e ) 2
(36)
Of course, this expression can be equivalently written in terms of the particle concentration s , that is the function plotted in Fig. 6. The entropy exhibits a maximum
322
J.J. Brey et al. / Physica A 275 (2000) 310–324
Fig. 6. Speci c entropy as a function of the density of particles s . The derivative of this function gives the inverse of the compactivity, except for a change of sign.
= 0:48 at a particle density s = 0:72, which corresponds to the in nite compactivity limit X → ∞. This is a result that could be foreseen, because 1 d = : dm0s X
(37)
In the context of the e ective tapping model, the whole entropy curve is physically meaningful and, the part of the curve to the left of the maximum, which corresponds to a negative compactivity X , cannot be discarded. This is, of course, a direct consequence of the existence of a limited number of con gurations. A similar behavior of the compactivity for a model of particle parking has been found by Monanson and Pouliquen [15].
4. Final remarks In this paper we have shown that a one-dimensional model for granular compaction, exhibiting the main characteristic features observed in real experiments, veri es the main assumption of the thermodynamic theory proposed by Edwards and his collaborators [2,3]. Because of the simplicity of the model, we have been able to compute explicitly the probability distribution for the con gurations in the steady state reached by the system in the limit of a great number of taps. This distribution depends only on the density of the system, being independent of other details of the con guration. Of course, a di erent point is that the possible metastable con gurations are limited by the stability conditions following from the rules de ning the dynamics of the system. In our simple model every empty site in the lattice must be surrounded by two particles. From the steady distribution, other steady properties of the system such as the “entropy” or the “compactivity” are also easily computed. This latter quantity, that is
J.J. Brey et al. / Physica A 275 (2000) 310–324
323
the analogous of the temperature for usual macroscopic thermodynamic systems, is negative in the small density region. As it is the case with the temperature in ordinary statistical systems, this is a consequence of the number of accessible states being bounded, and does not imply by itself any inconsistency or limitation of the theory. The results mentioned above refers to the e ective model, introduced as an approximation to another model referred to as the original one [13]. The former is expected to give an accurate description of the latter for tapping processes in the limit of very short taps. A question arising in a natural way is whether the thermodynamic theory is also valid for tapping processes in the original model with not so short taps. Although we have not been able to obtain an analytical proof, it seems on physical grounds that the answer is negative. A careful analysis of the reasoning leading to the form of the steady probability distribution, Eq. (5), indicates that an essential ingredient is that the elementary processes leading to the increase in the number of holes be the inverse of the elementary processes decreasing this number and that the relative probability of two given increasing processes be the same as that for their inverse. For long taps, the number of transitions taking place in the system during a tap can be quite large, and the associate probability rather involved. It is hard to believe that the structure of the probability distribution for processes increasing the density be conserved for processes creating holes. Perhaps a more detailed and de nite answer could be found by means of numerical measurements along the same lines as in Ref. [5]. A point we have not addressed here but deserves attention is the extension of the thermodynamic theory to describe the time relaxation of the density through a
uctuation-dissipation theorem [16] or some generalization for o -equilibrium situations [17]. The applicability of these ideas to our model will be discussed elsewhere.
Acknowledgements This research was partially supported by the Direccion General de Investigacion Cient ca y Tenica (Spain) through Grant No. PB98-1124.
Appendix A Calculation of the partition function Let us consider a lattice with N sites, from which k are empty. We need the number of possible con gurations k(N ) such that each empty site or hole is surrounded by two particles. A simple combinatorial calculation gives
k(N ) =
(N − k + 1)! k!(N − 2k + 1)!
(A.1)
324
J.J. Brey et al. / Physica A 275 (2000) 310–324
and substitution of this expression into Eq. (24) yields Z=
N=2 X
Zk ;
(A.2)
k=1
Zk =
(N − k + 1)! k : k!(N − 2k + 1)!
(A.3)
For large N we can evaluate the sum in Eq. (A.2) by means of the steepest descent method. The maximum of ln Zk appears for k = k˜ such that =
˜ − k) ˜ k(N ˜ 2 (N − 2 k)
(A.4)
and since does not depend on N =
(1 − ) (1 − 2)2
(A.5)
with k˜ = N. Solving this equation to get as a function of , and taking into account ˜ that k6N=2 the following expression is obtained: =
(1 + 4 )1=2 − 1 : 2(1 + 4 )1=2
(A.6)
Substitution of k˜ = N, with given by this expression, into ˜ ln(N − k) ˜ − k˜ ln k˜ − (N − 2 k) ˜ ln(N − 2 k) ˜ ; ln Z ' ln Z k˜ ' k˜ ln + (N − k) (A.7) valid for large N , gives Eq. (26). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
H.M. Jaeger, S.R. Nagel, R.P. Behringer, Rev. Mod. Phys. 68 (1996) 1259. S.F. Edwards, R.B.S. Oakeshott, Physica A 157 (1989) 1080. S.F. Edwards, C.C. Moun eld, Physica A 210 (1994) 279. J.B. Knight, C.G. Frandich, C.N. Lau, H.M. Jaeger, S.R. Nagel, Phys. Rev. E 51 (1995) 3957. E.R. Nowak, J.B. Knight, E. Ben-Naim, H.M. Jaeger, S.R. Nagel, Phys. Rev. E 57 (1998) 1971. E.R. Nowak, J.B. Knight, M. Povinelli, H.M. Jaeger, S.R. Nagel, Powder Technol. 94 (1997) 79. E. Caglioti, V. Loreto, H.J. Herrmann, M. Nicodemi, Phys. Rev. Lett. 79 (1997) 1575. M. Nicodemi, A. Coniglio, H.J. Herrmann, Phys. Rev. E 55 (1997) 3962. M. Nicodemi, A. Coniglio, H.J. Herrmann, Physica A 240 (1997) 405. M. Nicodemi, A. Coniglio, H.J. Herrmann, J. Phys. A 30 (1997) L397. M. Nicodemi, A. Coniglio, H.J. Herrmann, Phys. Rev. E 59 (1999) 6830. E. Ben-Naim, J.B. Knight, E.R. Nowak, H.M. Jaeger, S.R. Nagel, Physica D 123 (1998) 380. J.J. Brey, A. Prados, B. Sanchez-Rey, preprint cond-mat/9907144, Phys. Rev. E, to appear (Oct. 1999). N.G. van Kampen, Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam, 1993. R. Monasson, O. Pouliquen, Physica A 236 (1997) 395. S.F. Edwards, D.V. Grinev, Phys. Rev. E 58 (1998) 4758. M. Nicodemi, Phys. Rev. Lett. 82 (1999) 3734.
www.elsevier.com/locate/physa
Thermodynamic description in a simple model for granular compaction a Fsica
J. Javier Brey a; ∗ , A. Pradosa , B. Sanchez-Rey b
Teorica, Facultad de Fsica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, Spain b Fsica Aplicada, E.U. Politecnica, Universidad de Sevilla, Virgen de Africa 7, E-41011 Sevilla, Spain Received 18 July 1999
Abstract A simple model for the dynamics of a granular system under tapping is studied. The model can be considered as a particularization for short taps of a more general one-dimensional lattice model with facilitated dynamics. The steady state reached by the system is discussed and the results are shown to be consistent with the thermodynamic granular theory developed by Edwards and coworkers. In particular, the basic assumption of the theory, i.e., that the probability distribution c 2000 Elsevier Science B.V. depends only on the volume of the con guration, is veri ed. All rights reserved. PACS: 81.05.Rm; 05.50.+q; 81.20.Ev Keywords: Granular compaction; Thermodynamics; Tapping processes; Lattice models; Facilitated dynamics
1. Introduction A signi cant body of literature has emerged recently trying to apply the methods of kinetic theory and uid dynamics to the study of granular ows [1]. Also, some e orts have been devoted to formulate a statistical mechanics theory for steady granular media prepared under di erent conditions and, in particular, when submitted to external vibrations [2,3]. These studies have been motivated strongly by a series of experiments carried out to identify the mechanisms leading to compaction [4 – 6]. Granular compaction consists in the evolution from an initial low-density state to one with higher ∗
Corresponding author. Tel.: +3495-462-6558; fax: +3495-461-2097. E-mail address: [email protected] (J.J. Brey)
c 2000 Elsevier Science B.V. All rights reserved. 0378-4371/00/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 9 9 ) 0 0 3 7 5 - 1
J.J. Brey et al. / Physica A 275 (2000) 310–324
311
density as vibrations or, with more generality, some kind of series of external excitations are applied to the system. In actual experiments, the most usual perturbations are discrete, vertical shakes or taps applied to the container. The aspects considered in the studies of granular compaction include both the statics, i.e., the steady state reached by the system in the long time limit, and the dynamics, characterized by a time evolution of the density that is accurately described by an inverse logarithmic law. To investigate the latter aspect a series of models have been introduced [7–13], presenting most of the qualitative features of granular compaction observed in experiments. In particular, the density evolution is well tted by the inverse logarithmic law. Let us mention that most of the results from these models have been obtained numerically, while analytical results are scarce. In this paper we will deal with the description of the steady state reached by a granular system under tapping, by using a very simple one-dimensional lattice model, whose dynamics is a Markov process speci ed by means of a master equation. It was introduced in Ref. [13], where it was shown to exhibit the characteristic inverse logarithmic behavior mentioned above, and some properties of the density of the steady state were discussed. The basic idea here will be to identify an e ective Markov process for the description of the evolution of the system when submitted to tapping. In general, tapping processes are dicult to describe and analyze analytically since they are composed by an alternating series of two di erent dynamical evolutions: taps and freely relaxations. The nal state of every tap is the initial state for the next free evolution and vice versa. In the experiments and also in the computer simulations, it has been found useful to describe the evolution under tapping in terms of the number of taps. The values of the properties of the system are recorded at the end of the relaxation following every tap. The characteristic inverse logarithmic law mentioned above appears when the values of the density measured in this way are considered as a function of the number of taps. Thus, it is tempting to try to coarse grain the time description provided by a given model and consider each tap and its following relaxation as a single event. The derivation of the associated equivalent reduced dynamics is not a simple task, especially taking into account that the two processes that must be described as a single one are quite di erent in nature and that the transition between them is not smooth. Nevertheless, we have been able to do it for our model in the limit of very short or weak taps. The e ective tapping process is also formulated in terms of a master equation with well-de ned constant transition rates. Besides, the equation is expressed in terms of a continuous time variable, although it is only meaningful when a time equal to the duration of a given number of taps is considered. A direct consequence of having a master equation with time-independent transition rates to describe the relaxation under tapping is that it is possible to study the existence and properties of the steady con gurations. This allows, in particular, to check whether the model agrees with the thermodynamic description provided by the theory developed by Edwards and his collaborators [2,3], and whose main ingredient is the assumption that the volume of a granular media plays a role
312
J.J. Brey et al. / Physica A 275 (2000) 310–324
analogous to the energy of a statistical molecular system. We will show that our e ective model veri es this principle, which is a nontrivial result. The steady probability distribution only depends on the con guration of the lattice through the associated density. The paper is organized as follows. In Section 2 the model of Ref. [13] is brie y reviewed. This model will be referred to in the following as the original model. Next, it is applied to describe tapping processes in the limit of very short taps, and the relevant elementary events connecting the state of the system before the tap and after the combination of tapping plus free relaxation are identi ed. The resulting process always ends at the lowest order in the duration of the taps with all the sites (except one) occupied by a particle, due to the inability of the system to decrease its density in any elementary transformation. To describe steady states of di erent densities it is then necessary to take into account also some next order processes, namely those increasing the number of empty sites in the lattice. The form of the probability distribution in the steady state is discussed in Section 3. It can be written in the same way as the canonical distribution for a normal statistical system, the role of the energy being played by the volume and that of the temperature by a new parameter, the compactivity. This is the required result to connect with the thermodynamic theory of Ref. [2,3]. Finally, the last section contains some remarks and comments about the main results in the paper. 2. Description of the model The model for tapping processes we are going to formulate is based on a previous model with facilitated dynamics considered recently in Ref. [13]. Let us review the latter brie y. It is a one-dimensional lattice model in which each of the N sites can be either occupied by a particle or can be empty. A variable mi is assigned to each site i, taking the value mi = 1 if the site is empty, while mi = 0 if there is a particle on it. Then, a con guration of the system is fully speci ed by giving the values of the set of variables m ≡ {m1 ; m2 ; : : : ; mN }. The dynamics of the system is assumed to be a Markov process de ned by the following master equation for the transition probability p1=1 (m; t|m0 ; t 0 ) for nding the system with the con guration m at time t given it was in the con guration m0 at time t 0 ¡ t, @t p1=1 (m; t|m0 ; t 0 ) =
N X
[Wi (Ri m)p1=1 (Ri m; t|m0 ; t 0 ) − Wi (m)p1=1 (m; t|m0 ; t 0 )] :
i=1
(1) Here Ri m ≡ {m1 ; : : : ; Rmi ; : : : ; mN }, with Rmi = 1 − mi , and the expressions of the transition rates are Wi (m) = (mi−1 + mi+1 )[ + mi (1 − 2)] ; (2) 2
J.J. Brey et al. / Physica A 275 (2000) 310–324
313
where and are real parameters de ned in the intervals ¿ 0 and 0661, respectively. Thus the transition rate for the adsorption of a particle at site i is Wi+ (m) = (1 − )mi (mi−1 + mi+1 ) (3) 2 and the probability per unit of time for the desorption of a particle from the same site is (4) Wi− (m) = (1 − mi )(mi−1 + mi+1 ) : 2 For =1 there are no adsorption processes and for =0 particles cannot be desorbed. The rational for this system as a model for granular compaction has been discussed in Ref. [13] and will not be reproduced here. We are interested in the particularization of the above general dynamics to tapping processes. They are modeled in the following way. Starting from an arbitrary state or con guration, the system is allowed to relax with = 0, i.e., without desorption events, until it reaches a metastable con guration. Such con guration is characterized by all the holes or empty sites being isolated so that the transition rates given by Eq. (3) also vanish. Therefore, the system cannot evolve any longer and becomes trapped. This metastable con guration will be the initial state to describe the tapping process and will be referred to as state n = 0. Next, the value of is instantaneously increased to a nite value in order to simulate a pulse. Its duration t0 is taken very short as compared with the characteristic time of the transition attempts −1 , i.e., t0 .1. Of course, the tap has in general the e ect of taking out the system out of the metastable con guration. Afterwards, the system relaxes again with = 0 until a new metastable con guration n = 1 is reached. Clearly, evolution with = 0 corresponds to free relaxation in real tapping experiments. By repeating the same sequence of steps a series of metastable states n = 2; 3; 4; : : : is generated. The index n indicates the number of taps or short pulses applied to the system before that con guration has been obtained. Let us consider one of the metastable con gurations of the system after relaxing with = 0. We know that each hole is surrounded by two particles. When the value of is increased, the rst transition occurring in the system must be the desorption of a particle having at least one of its nearest-neighbor sites empty. Therefore, the possible rst elementary events are ··· 1 0 0 ··· → ··· 1 1 0 ··· ; ··· 0 0 1 ··· → ··· 0 1 1 ··· ; ··· 1 0 1 ··· → ··· 1 1 1 ··· :
(5)
We have employed the mi variables to specify the con gurations and have represented only the site whose state is changing and the two nearest-neighbor ones, since they determine the transition rate. For t0 small enough, only one transition can take place at the most in the system during a pulse. Afterwards if it is allowed to relax freely again to a new metastable
314
J.J. Brey et al. / Physica A 275 (2000) 310–324
con guration, we can compute the probability for each of the possible series of transitions connecting two metastable con gurations. For instance, let us consider the rst process shown in Eq. (5). The probability of the desorption event during the pulse is t0 =2, since from Eq. (4) its transition rate is equal to =2. Next, in the free relaxation a particle can be adsorbed with equal probability on any of the two empty sites. Then, this process leads to the same state or to the e ective di usion of the hole to the right. We will not consider those global transitions in which the nal state coincides with the initial one, because they are not necessary when writing a master equation. In a trivial notation, in which only the particles involved in the process are indicated explicitly, the probability for the e ective transition considered is 1
(1) (100 → 010) = t0 × = t0 : (6) 2 2 4 The same line of reasoning allows us to calculate the probabilities for the other e ective transitions between two metastable con gurations, namely
(1) (001 → 010) = t0 ; (7a) 4
(1) (101 → 001) =
t0 ; 8
(7b)
(1) (101 → 010) =
t0 ; 4
(7c)
t0 : (7d) 8 The key point in our approximation is introduced now. It will be assumed that in the limit of very short taps, the time evolution of the system under tapping can be accurately described by a continuous Markov process with transition rates given by
(1) (101 → 100) =
(1)
(m|m0 ) =
N X
(1) 0 i (m|m )
;
(8)
i=1 (1) 0 i (m|m ) =
(i) (1 − m0i ){2[m0i−1 (1 − m0i+1 ) + (1 − m0i−1 )m0i+1 + m0i−1 m0i+1 ] 8 ×(1 − mi−1 )mi (1 − mi+1 ) +m0i−1 m0i+1 [(1 − mi−1 )(1 − mi )mi+1 + mi−1 (1 − mi )(1 − mi+1 )]} ; (9)
where (i) re ects that the e ective transition only a ects three consecutive sites, (i) =
N Y j6=i−1; i; i+1
mj ; m0j :
(10)
J.J. Brey et al. / Physica A 275 (2000) 310–324
315
In summary, the master equation X @t p1=1 (m; t|m0 ; t 0 ) = [ (1) (m|m00 )p1=1 (m00 ; t|m0 ; t 0 ) m00
−
(1)
(m00 |m)p1=1 (m; t|m0 ; t 0 )]
(11)
provides the lowest-order approximation for the description of the compaction process of the original model under tapping provided that t0 .1 and the equation is applied for times t = nt0 , with n being a whole number. The free relaxation with = 0 is considered as instantaneous in this description and assumed to take place just after the transition events associated with the taps. Then, the above master equation does not contain any explicit reference to the unstable con gurations of the system during the tapping. They have been eliminated in the time coarse-graining implicit in Eq. (11). As a consequence, the time scale in the approximated model is just an auxiliary variable, providing the number of taps n by dividing it by t0 . A point to stress is that the approximation we are discussing is not the same as the continuous tapping limit considered by some authors. In the latter the system is not allowed to freely relax between consecutive taps, while here it always relaxes to a metastable con guration after every pulse, and thus the e ective dynamics connects the metastable states among themselves. In order to get some insight into the implications of Eq. (11), let us study the time evolution of the empty site concentration x1 (t), given by X x1 (t) = mi p(m; t) ; (12) m
where p(m; t) ≡
X
p1=1 (m; t|m0 ; 0)
(13)
m0
is the one-time probability distribution. For homogeneous conditions, the right-hand side of Eq. (12) does not depend on the site i considered. Then, a rather simple calculation leads to @x1 = − hmi−1 mi+1 it : (14) @t 2 Here we have used the notation X hf(m)it = f(m)p(m; t) (15) and taken into account that metastable con gurations cannot have three consecutive holes and, therefore, expressions like hmi−1 mi mi+1 it identically vanish in those states. Since hmi−1 mi+1 it ¿0, the time approach of the hole concentration towards its steady value is monotonic in this approximation. This is a direct consequence on the e ective dynamics we have introduced, in which particles can never be desorbed (see Eqs. (6) and (7)). In fact, it is easily seen that all the states with only one empty site are absorbent states [14] for the tapping, connected by di usion processes. It follows that in order to be able to describe tapping evolution leading to steady states
316
J.J. Brey et al. / Physica A 275 (2000) 310–324
Fig. 1. All the elementary trajectories leading to a decrease in the number of holes, with only one desorption event during the tap. We have indicated the processes taking place during the pulse ( ¿ 0) and in the free relaxation ( = 0). The trajectories leading to a nal state being the same as the initial one are not shown.
Fig. 2. All the elementary paths in which the number of holes increases. Two desorption events during the pulse ( ¿ 0) are necessary, while only one transition occurs during the free relaxation ( = 0). Note that the plotted trajectories are just the inverse of the transitions in Fig. 1.
as a consequence of the balance between adsorption and desorption processes, we have to consider higher-order terms in the duration t0 and strength of the individual taps. When considering the next order in t0 contributions, two quite qualitatively di erent kinds of terms appear. There are terms that represent just corrections to nonvanishing lowest-order transition rates. These will be not incorporated in the simpli ed model, since the lowest order is the dominant one in the limit we are considering. On the other hand, some second-order terms provide transition rates for elementary processes that are not present to lowest order. This happens, in particular, for all the transformations leading to the creation of a hole. As discussed above, the incorporation of these processes modi es substantially the physics of the tapping process and they must be retained. In Fig. 1 we have represented the processes leading to the decrease of the number of holes taken into account upon deriving Eq. (11). In each of the elementary events we have indicated whether they take place during the tap ( ¿ 0) or in the free relaxation ( = 0). The question now is whether the inverse processes indicated in Fig. 2 are possible and what are their transition rates. The rst point to realize is that the inverse of the elementary transformations with = 0 are only possible with ¿ 0, since they correspond to desorption processes. Therefore, we need at least two transitions to occur during the tap. Nevertheless, the total number of elementary events involved does not increase since the free relaxation reduces to only one transformation. Although, in principle, it could be thought that there are other second-order transformations, i.e.,
J.J. Brey et al. / Physica A 275 (2000) 310–324
317
with two motions during the tap, increasing the number of holes, it is easy to convince oneself that there are none when trying to construct them. In particular, the group of three particles involved in the transition must have holes on both sides. Then, we have to consider the following e ective elementary processes: 0 1 0 0 0→0 1 0 1 0; 0 0 1 0 0→0 1 0 1 0; 0 0 0 1 0→0 1 0 1 0:
(16)
Accordingly, with the previous discussion we have eliminated those processes in Fig. 2 conserving the number of empty sites. It is important to realize that the transformations kept are just the inverse of the rst-order transformations leading to the loss of a hole. The computation of the e ective transition rates for the processes in Eq. (16) requires to take into account that the two transformations with ¿ 0 must occur in the time interval t0 , the duration of a tap, with probability equal to (t0 )2 =8 in all of the cases shown in Fig. 2. Then the following expressions for the probabilities are obtained:
(2) (01000 → 01010) =
2 2 t02 ; 16
(17a)
(2) (00100 → 01010) =
2 2 t02 ; 8
(17b)
2 2 t02 : (17c) 16 The equality of probabilities (17a) and (17c) is a direct consequence of the isotropy of the dynamics of the original model. The inclusion of these elementary e ective events is a necessary condition to describe steady states with nonvanishing concentration of holes since, as we have already discussed, there are no processes decreasing the particle concentration to the lowest order. To each of the probabilities in Eqs. (17) a transition rate given by (2) = (2) =t0 can be associated, and the corresponding contributions must be incorporated into the expressions of the total transition rate (m|m0 ) in Eq. (8). Nevertheless, this leads to a rather complicated and not very illuminating expression that will not be written here. On the other hand, it is important to note that with the inclusion of the transition rates (2) all the metastable states of the system are connected through a possible chain of transitions. As a test of the validity of the e ective model, as an approximation of the original one in Fig. 3, we have plotted the time evolution, as a function of the number of taps n, of the particle density as predicted by both models. In the simulations periodic boundary conditions and a lattice with 104 sites have been used. In the e ective model n is obtained by dividing the time in units of t0 . The solutions of the equations can be characterized by the parameter
(2) (00010 → 01010) =
= t0 =2 :
(18)
For the curves in Fig. 3 it is = 0:015. It is seen that there is a quite good agreement for all the evolution. As expected, the agreement is not so good as the value of
318
J.J. Brey et al. / Physica A 275 (2000) 310–324
Fig. 3. Evolution of the density of particles with the number of taps n obtained from Monte Carlo simulation of the system, for = 0:015. The diamonds correspond to the e ective model introduced in this paper, while the crosses correspond to the original model of Ref. [7]. The dotted line is the theoretical prediction for the steady-state density of the e ective model, from Eq. (27).
Fig. 4. The same as in Fig. 3, but for = 0:5. Note that the deviations of the e ective model from the original one are larger, but the relative error is still quite small.
increases. An example is given in Fig. 4 obtained with = 0:5. Nevertheless, notice that the relative di erence is still very small. As discussed in Ref. [13], the plotted relaxation curves can be accurately tted by means of an inverse logarithmic law. The dotted lines in Figs. 3 and 4 are the steady values of the density as predicted by the e ective model. They will be discussed in the next section. 3. Steady state The Markov process de ned in the previous section is irreducible and, therefore, the master equation has a unique steady solution [14]. We are going to see that this solution
J.J. Brey et al. / Physica A 275 (2000) 310–324
319
veri es detailed balance and this property can be used to determine it explicitly. Thus, we a priori bet on a time-independent solution ps (m) of the e ective master equation for tapping verifying (m|m0 )ps (m0 ) = (m0 |m)ps (m) :
(19)
If such distribution exists, it will have the property that all the con gurations having the same number of empty sites or holes will have the same probability. This follows from the fact that they can be connected by di usive transitions and di usion is an isotropic process in our system, so that (m|m0 )= (m0 |m) if m and m0 have the same number of holes. Therefore, the value of ps (m) will depend only on the number of holes in the con guration m. In the following we will write m(k) to indicate that this con guration has k empty sites. For processes increasing (decreasing) the number of particles we have p s (m0(k+1) ) = p s (m(k) )
(m0(k+1) | m(k) ) : (m(k) | m0(k+1) )
(20)
A look at Eqs. (6); (7) and (17) shows that (m0(k+1) | m(k) ) t0 = (m(k) | m0(k+1) ) 2
(21)
for all the nonvanishing transition rates. Then, p s (m0(k+1) ) t0 = p s (m(k) ) 2
(22)
which does not depend on the particular con guration under consideration, in agreement with the conclusion reached above that p s (m(k) ) only depends on the value of k. Eq. (22) implies p s (m(k) ) =
k ; Z
(23)
where Z is a normalization constant and is de ned in Eq. (18). For a system with N sites we can formally write Z=
N=2 X k=1
k(N ) k ;
(24)
where k(N ) is the degeneracy of a lattice in a metastable state with k holes, i.e., the number of di erent con gurations having exactly k isolated holes. Note that the maximum number of holes is N=2 (assumed N is even) and the minimum one is unity. The steady hole concentration is N=2 N X X 1 1 X X s (k) @ 1 s s ln Z : (25) mj p (m) = k p (m ) = x1 = N m N @ N (k) j=1
k=1
m
320
J.J. Brey et al. / Physica A 275 (2000) 310–324
Fig. 5. Steady value of the density of particles s as a function of the parameter . The continuous line is Eq. (27), and the diamonds are the Monte Carlo values for the original model in Ref. [13].
The evaluation of the “partition function” Z, which requires previously the calculation of the number of states k(N ) , is sketched in the appendix. In the limit of a very large system (N → ∞), the result is 1 + (1 + 4 )1=2 : 2 Use of this expression into Eq. (25) yields ln Z = N ln
x1s =
(1 + 4 )1=2 − 1 : 2(1 + 4 )1=2
(26)
(27)
As discussed in Section 2, the e ective model we are considering is expected to be equivalent to the original model of Ref. [13] in the limit .1. In this limit it is x1 ' which agrees with the result from the original model [13]. The theoretical prediction given by Eq. (27) is compared with the numerical results obtained by solving numerically the original master equation, Eq. (1) in Fig. 5. The agreement can be considered as quite good for all the range of values of considered. Let us stress here that Eq. (27) does not present any singularity when considered as a function of . The steady distribution (5) can be rewritten in the usual form of a “canonical” probability distribution p s (m(k) ) = Z −1 e−k=X ;
Z=
N=2 X k=1
k(N ) e−k=X ;
(28)
where 1 (29) ln is the analogous of a temperature parameter. It was introduced by Edwards and coworkers in the context of a granular thermodynamic theory and called compactivity [2,3]. X =−
J.J. Brey et al. / Physica A 275 (2000) 310–324
321
The basic assumption of the theory is that the volume plays the same role in a vibrated granular material as the energy in a statistical molecular system. Then, all the con gurations of a given volume are equally probable. We have seen that this is veri ed in our model as a consequence of the dynamics. In the limit of small the steady concentration of particles s will be s ≡ 1 − x1s ' 1 − = 1 − e−1=X ;
(30)
1=X = ln(1 − s ) :
(31)
i.e.,
A relationship similar to this has been recently obtained from the analysis of experimental data [5, See Eq. (6)]. The most compact granular system, i.e., a system with x1s = 0 corresponds to = 0 and X = 0. When X increases also increases and, therefore, the equilibrium density of particles 1 − x s also decreases. This is in complete agreement with the picture on which the thermodynamic theory of powders is based [2,3]. In the limit X → ∞ it is x1s ' 0:28. In principle, negative compactivities must also be considered as we will see below. The limit X → 0− appears for → ∞, and then x1s → 12 , i.e., the least compact arrangement. Nevertheless, we must keep in mind that our model is expected to provide an accurate description of tapping processes in the original model for ¡ 1, and this clearly excludes negative values of X . The analogous of the free energy per site is =−
X ln Z ; N
(32)
so that Eq. (26) is equivalent to x1s =
d(=X ) : d(1=X )
It is also natural to de ne an “entropy” by X N p s (m)ln p s (m) = x1s + ln Z : S≡− X m
(33)
(34)
As expected it turns out to be an extensive property, i.e., proportional to N for large N , so that we can consider a speci c entropy per site = S=N , verifying the relationship = x1s − X :
(35)
Since we know the explicit expression for the partition function of our model, we can compute all the above quantities. In particular, a simple calculation gives =
(1 + 4e−1=X )1=2 − 1 1 + (1 + 4e−1=X )1=2 : + ln −1=X 1=2 2X (1 + 4e ) 2
(36)
Of course, this expression can be equivalently written in terms of the particle concentration s , that is the function plotted in Fig. 6. The entropy exhibits a maximum
322
J.J. Brey et al. / Physica A 275 (2000) 310–324
Fig. 6. Speci c entropy as a function of the density of particles s . The derivative of this function gives the inverse of the compactivity, except for a change of sign.
= 0:48 at a particle density s = 0:72, which corresponds to the in nite compactivity limit X → ∞. This is a result that could be foreseen, because 1 d = : dm0s X
(37)
In the context of the e ective tapping model, the whole entropy curve is physically meaningful and, the part of the curve to the left of the maximum, which corresponds to a negative compactivity X , cannot be discarded. This is, of course, a direct consequence of the existence of a limited number of con gurations. A similar behavior of the compactivity for a model of particle parking has been found by Monanson and Pouliquen [15].
4. Final remarks In this paper we have shown that a one-dimensional model for granular compaction, exhibiting the main characteristic features observed in real experiments, veri es the main assumption of the thermodynamic theory proposed by Edwards and his collaborators [2,3]. Because of the simplicity of the model, we have been able to compute explicitly the probability distribution for the con gurations in the steady state reached by the system in the limit of a great number of taps. This distribution depends only on the density of the system, being independent of other details of the con guration. Of course, a di erent point is that the possible metastable con gurations are limited by the stability conditions following from the rules de ning the dynamics of the system. In our simple model every empty site in the lattice must be surrounded by two particles. From the steady distribution, other steady properties of the system such as the “entropy” or the “compactivity” are also easily computed. This latter quantity, that is
J.J. Brey et al. / Physica A 275 (2000) 310–324
323
the analogous of the temperature for usual macroscopic thermodynamic systems, is negative in the small density region. As it is the case with the temperature in ordinary statistical systems, this is a consequence of the number of accessible states being bounded, and does not imply by itself any inconsistency or limitation of the theory. The results mentioned above refers to the e ective model, introduced as an approximation to another model referred to as the original one [13]. The former is expected to give an accurate description of the latter for tapping processes in the limit of very short taps. A question arising in a natural way is whether the thermodynamic theory is also valid for tapping processes in the original model with not so short taps. Although we have not been able to obtain an analytical proof, it seems on physical grounds that the answer is negative. A careful analysis of the reasoning leading to the form of the steady probability distribution, Eq. (5), indicates that an essential ingredient is that the elementary processes leading to the increase in the number of holes be the inverse of the elementary processes decreasing this number and that the relative probability of two given increasing processes be the same as that for their inverse. For long taps, the number of transitions taking place in the system during a tap can be quite large, and the associate probability rather involved. It is hard to believe that the structure of the probability distribution for processes increasing the density be conserved for processes creating holes. Perhaps a more detailed and de nite answer could be found by means of numerical measurements along the same lines as in Ref. [5]. A point we have not addressed here but deserves attention is the extension of the thermodynamic theory to describe the time relaxation of the density through a
uctuation-dissipation theorem [16] or some generalization for o -equilibrium situations [17]. The applicability of these ideas to our model will be discussed elsewhere.
Acknowledgements This research was partially supported by the Direccion General de Investigacion Cient ca y Tenica (Spain) through Grant No. PB98-1124.
Appendix A Calculation of the partition function Let us consider a lattice with N sites, from which k are empty. We need the number of possible con gurations k(N ) such that each empty site or hole is surrounded by two particles. A simple combinatorial calculation gives
k(N ) =
(N − k + 1)! k!(N − 2k + 1)!
(A.1)
324
J.J. Brey et al. / Physica A 275 (2000) 310–324
and substitution of this expression into Eq. (24) yields Z=
N=2 X
Zk ;
(A.2)
k=1
Zk =
(N − k + 1)! k : k!(N − 2k + 1)!
(A.3)
For large N we can evaluate the sum in Eq. (A.2) by means of the steepest descent method. The maximum of ln Zk appears for k = k˜ such that =
˜ − k) ˜ k(N ˜ 2 (N − 2 k)
(A.4)
and since does not depend on N =
(1 − ) (1 − 2)2
(A.5)
with k˜ = N. Solving this equation to get as a function of , and taking into account ˜ that k6N=2 the following expression is obtained: =
(1 + 4 )1=2 − 1 : 2(1 + 4 )1=2
(A.6)
Substitution of k˜ = N, with given by this expression, into ˜ ln(N − k) ˜ − k˜ ln k˜ − (N − 2 k) ˜ ln(N − 2 k) ˜ ; ln Z ' ln Z k˜ ' k˜ ln + (N − k) (A.7) valid for large N , gives Eq. (26). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
H.M. Jaeger, S.R. Nagel, R.P. Behringer, Rev. Mod. Phys. 68 (1996) 1259. S.F. Edwards, R.B.S. Oakeshott, Physica A 157 (1989) 1080. S.F. Edwards, C.C. Moun eld, Physica A 210 (1994) 279. J.B. Knight, C.G. Frandich, C.N. Lau, H.M. Jaeger, S.R. Nagel, Phys. Rev. E 51 (1995) 3957. E.R. Nowak, J.B. Knight, E. Ben-Naim, H.M. Jaeger, S.R. Nagel, Phys. Rev. E 57 (1998) 1971. E.R. Nowak, J.B. Knight, M. Povinelli, H.M. Jaeger, S.R. Nagel, Powder Technol. 94 (1997) 79. E. Caglioti, V. Loreto, H.J. Herrmann, M. Nicodemi, Phys. Rev. Lett. 79 (1997) 1575. M. Nicodemi, A. Coniglio, H.J. Herrmann, Phys. Rev. E 55 (1997) 3962. M. Nicodemi, A. Coniglio, H.J. Herrmann, Physica A 240 (1997) 405. M. Nicodemi, A. Coniglio, H.J. Herrmann, J. Phys. A 30 (1997) L397. M. Nicodemi, A. Coniglio, H.J. Herrmann, Phys. Rev. E 59 (1999) 6830. E. Ben-Naim, J.B. Knight, E.R. Nowak, H.M. Jaeger, S.R. Nagel, Physica D 123 (1998) 380. J.J. Brey, A. Prados, B. Sanchez-Rey, preprint cond-mat/9907144, Phys. Rev. E, to appear (Oct. 1999). N.G. van Kampen, Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam, 1993. R. Monasson, O. Pouliquen, Physica A 236 (1997) 395. S.F. Edwards, D.V. Grinev, Phys. Rev. E 58 (1998) 4758. M. Nicodemi, Phys. Rev. Lett. 82 (1999) 3734.