We perform extensive numerical simulations of different versions of the
sandpile model. We find that previous claims about universality classes are
unfounded, since the method previously employed to analyze the data suffered a
systematic bias. We identify the correct scaling behavior and conclude that
sandpiles with stochastic and deterministic toppling rules belong to the same
universality class.Comment: 4 pages, 4 ps figures; submitted to Phys. Rev. 

A. Ben-Hur

A. Chessa

A. Corral

A. Díaz-Guilera

A. Garciamartin

A. Vespignani

Alessandro Chessa

Alessandro Vespignani

D. Dhar

D. Spasojević

E. Milshtein

G. Durin

G. Gutenberg

H. Eugene Stanley

K. Christensen

L. A. N. Amaral

L. Pietronero

M. De Menech

P. Bak

P. Grassberger

R. Dickman

S. Ciliberto

S. Field

S. Lübeck

S. S. Manna

Stefano Zapperi

Y. C. Zhang

English

arXiv

We perform extensive numerical simulations of different versions of the sandpile model. We find that previous claims about universality classes are unfounded, since the method previously employed to analyze the data suffered from a systematic bias. We identify the correct scaling behavior and provide evidences suggesting that sandpiles with stochastic and deterministic toppling rules belong to the same universality class

Chessa, Alessandro

Stanley, H. Eugene

Vespignani, Alessandro

Zapperi, Stefano

IMT Institutional Repository

Universality in sandpiles

Crossref

S. Zapperi

AIR Universita degli studi di Milano

RAPID COMMUNICATIONS
PHYSICAL REVIEW E JANUARY 1999VOLUME 59, NUMBER 1Universality in sandpiles
Alessandro Chessa,1 H. Eugene Stanley,2 Alessandro Vespignani,3 and Stefano Zapperi4
1Dipartimento di Fisica and Unita´ INFM, Universita´ di Cagliari, Via Ospedale 72, 09124 Cagliari, Italy
2Center for Polymer Studies and Department of Physics, Boston University, Boston, Massachusetts 02215
3The Abdus Salam International Centre for Theoretical Physics (ICTP), P.O. Box 586, 34100 Trieste, Italy
4PMMH-ESPCI, 10 Rue Vauquelin, 75231 Paris Cedex 05, France
~Received 24 August 1998!
We perform extensive numerical simulations of different versions of the sandpile model. We find that
previous claims about universality classes are unfounded, since the method previously employed to analyze the
data suffered from a systematic bias. We identify the correct scaling behavior and provide evidences suggest-
ing that sandpiles with stochastic and deterministic toppling rules belong to the same universality class.
@S1063-651X~99!50701-0#
PACS number~s!: 05.65.1b, 05.40.2aSandpile automata @1# are among the simplest models to
describe avalanche propagation, a phenomenon of upsurging
experimental interest in a wide range of fields @2#. In the
stationary state, after suitable tuning of the driving fields @3#,
these models display critical behavior in the avalanche sta-
tistics. As for ordinary critical phenomena, it is possible to
define a set of scaling exponents to characterize the large
scale behavior of the system @3#.
The precise identification of universality classes in sand-
pile models @1# is an unresolved issue. From a theoretical
standpoint, it would be unusual that small modifications in
the dynamical rules of the model could lead to different uni-
versality classes. Real-space renormalization group calcula-
tions @5# suggest that different sandpile models, such as the
Bak, Tang, and Wiesenfeld ~BTW! @1#, and the Manna @4#
models, all belong to the same universality class. This result
is also confirmed by a recently proposed field theory ap-
proach @6# that shows that all sandpile models @7# are de-
scribed by the same effective field theory at the coarse
grained level. Universality is also found between BTW ~dis-
crete! and Zhang @8# ~continuous! models in the dynamical
renormalization group calculations of Ref. @9#.
The results obtained by numerical simulations are unclear.
Early large scale numerical simulations of the Manna @4# and
BTW models @10# show that the avalanche distributions are
described by the same exponents for the power law decay
and the scaling of the cutoffs. These results were questioned
by Ben-Hur and Biham @11# who analyzed the scaling of
conditional expectation values @12# of various quantities.
They found significant differences in the exponents for the
two models and therefore proposed a classification of univer-
sality in sandpile models, in which models with stochastic
update rules, such as the Manna model, fall into a universal-
ity class different from that of Abelian models, such as the
BTW @13#. The method was later applied to the Zhang
model, which was declared ‘‘nonuniversal’’ @14#. These re-
sults pose a puzzling problem since they contradict all of the
existing theories and do not agree with the scaling predicted
analyzing avalanche distributions @4,10#.
Here we present large scale numerical simulations of the
BTW and Manna sandpile models with the goal of settling
the issue of universality. First we show that the method ofPRE 591063-651X/99/59~1!/12~4!/$15.00conditional expectation values, introduced in Ref. @12# and
used in Ref. @11#, is systematically biased by nonuniversal
corrections and does not provide indications about universal-
ity classes. By removing the bias, we provide evidence that
the BTW and Manna models are universal. This conclusion
appears to be consistent with data collapse and moment
analysis of the distributions @15#.
Sandpile models are defined on a d-dimensional hypercu-
bic lattice. On each site i of the lattice we define an integer
variable zi , which we call ‘‘energy.’’ At each time step an
energy grain is added on a randomly chosen site (zi!zi
11). When one of the sites reaches or exceeds a threshold
zc a ‘‘toppling’’ occurs: zi5zi2zc and z j5z j11, where j
represents the nearest-neighbor sites of site i . In the BTW
model zc52d and each nearest-neighbor site receives a grain
after the toppling of the site i . In the Manna model zc52
and, therefore, only two randomly chosen neighboring sites
receive a grain. A toppling can induce nearest-neighbor sites
to topple on their turn and so on, until all of the lattice sites
are below the critical threshold. This process is called an
avalanche. A slow driving is usually imposed so that grains
are added only when all of the sites are below the threshold.
The model is conservative and energy is dissipated only at
the boundary sites @1#. Here we perform numerical simula-
tions of two-dimensional Manna and BTW models with open
boundary conditions and conservative dynamics. The lattice
size ranges from L5128 to L52048 in both models. In each
case, statistical distributions are obtained averaging over 107
nonzero avalanches.
Avalanches in sandpile models are usually characterized
by three variables: the number of topplings s , the area a
affected by the avalanche, and the avalanche duration T . The
probability distribution of each of these variables is usually
described as power law with a cutoff
P~x !5x2txG~x/xc!, ~1!
where x5s ,a ,T . When the system size L goes to infinity,
the cutoff xc diverges as xc;Lbx. Under the finite size scal-
ing ~FSS! assumption of Eq. ~1!, the set of exponents
$tx ,bx% defines the universality class of the model.R12 ©1999 The American Physical Society
RAPID COMMUNICATIONS
PRE 59 R13UNIVERSALITY IN SANDPILESIn two dimensions, an accurate numerical determination
of the power law exponents in Eq. ~1! proved to be a difficult
task @4,10,16,17# due to the large deviations at the lower and
upper cutoffs. For this reason, Christensen et al. @12#, in or-
der to distinguish among universality classes, proposed a
more refined numerical analysis based on the evaluation of
the expectation value E(xuy) of the variable x restricted to
all the avalanches with variable Y5y , where $X ,Y %
5$s ,a ,T% @12#. It is assumed that E(xuy);ygxy and the ex-
ponents gxy are used to distinguish among universality
classes @11#. These exponents satisfy the scaling relations
gxy5gyx
21 and gxz5gxygyz .
As stated in Ref. @12#, if the conditional probability dis-
tribution p(xuy) is sufficiently peaked, then gxy is well de-
fined, and to each value of the variable x we can unambigu-
ously associate a value of the variable y ~i.e., x;ygxy!. In
particular, the cutoff of the distributions should be related by
the same exponents ~i.e., xc;yc
gxy!, which implies gxy
5bx /by . For instance, we have gsa5bs/2 since, in two
dimensions, avalanches are compact for both the BTW @10#
and Manna models @11#, so that ba52. The data collapse
analysis shows that the BTW and Manna models both share
the same exponent, bs.2.7 @4,10,16#, which implies gsa
.1.35. On the contrary, Refs. @11,14# found gsa.1.06 for
the BTW model and gsa51.24 for the Manna model, which
would yield two different universality classes for the two
models. Less marked differences were also observed for the
other exponents gxy @11,14#.
In order to resolve this paradox, we return to the hypoth-
esis underlying the use of conditional expectation values:
p(xuy) must be symmetrical and strongly peaked around the
average value. We checked numerically that this assumption
is not fulfilled; in the BTW model the distribution p(sua) is
maximum for s5a and decreases for s.a , with a charac-
teristic value s* scaling as abs/2 ~see Fig. 1!. The distribution
is not symmetric ~see also Ref. @17#!, consistent with the
constraint s>a ~the avalanche area cannot be greater than
the number of topplings!. Similar considerations apply, as
well, to other quantities ~i.e., a>T , s>T!, in which condi-
tional probability distributions show asymmetry, although
less marked.
To understand the effect of nonsymmetric distributions on
conditional expectation values, consider a distribution of the
form
p~xuy !5u~x2y ! f @~x2y !/x*#/x*, ~2!
where the characteristic value scales as x*(y);ygxy, and
u(x) is the step function. The factor 1/x* ensures normaliza-
tion for any y ,
E dxp~xuy !51, ~3!
so that the conditional expectation value is given by
E~xuy ![E
y
`
dx
x
x*
f @~x2y !/x*# . ~4!
Performing the substitution z5x2y , we obtainE~xuy !5y1E
0
`
dz
z
x*
p~z/x*!5y1Cygxy, ~5!
where C is a nonuniversal constant.
In the BTW model p(sua) has the form of Eq. ~2!, as
shown by performing data collapse analysis @see Fig. 1 ~in-
set!#. Thus, we can easily subtract the linear bias from the
expectation value in order to obtain the correct scaling be-
havior to be compared with that of the Manna model ~Fig. 2!,
for which conditional distributions appear to be symmetric.
Data from avalanche areas up to a.106 provide striking
evidence that both models share the same asymptotic behav-
ior with an exponent gsa51.3560.05, in agreement with
other published results @4,10,16,17#. The scaling of the other
expectation values is also biased, as is apparent from the
bending in the curves reported in Refs. @11,14#. The correc-
tion of the bias is not so straightforward as in the case we
FIG. 1. Probability distribution of having an avalanche size s
given its area a for the BTW model. The inset shows that all data
collapse onto the universal scaling function p(sua)5a2gsa f @(s
2a)/agsa# , with gsa.1.35.
FIG. 2. Conditional expectation value E(sua) for the BTW and
Manna models ~after bias subtraction!. The slope is given by gsa
51.3560.05 for both curves.
RAPID COMMUNICATIONS
R14 PRE 59CHESSA, STANLEY, VESPIGNANI, AND ZAPPERIhave discussed, but can be obtained from the analysis of
p(xuy). This discussion clearly shows that conditional ex-
pectation values are not a reliable method to determine the
universality class of a model, unless a systematic analysis of
the bias is performed.
To provide further evidence about a single universality
class, we perform the moment analysis on the distributions
P(x ,L), in close analogy with the recent work of De Menech
et al. @15# on the two-dimensional BTW model. Here we
apply the moments analysis on both the BTW and Manna
models, taking advantage of the large sizes reached in our
numerical simulations. We define the q moment of x on a
lattice of size L as ^xq&L5*xqP(x)dx . If the FSS hypothesis
@Eq. ~1!# is valid, at least in the asymptotic limit (x!`), we
can transform z5x/Lbx and obtain
^xq&L5Lbx~q112t!E zq1tG~z !dz;Lbx~q112t!, ~6!
or, in general, ^xq&L;Lsx(q). The exponents sx(q) can be
obtained as the slope of the log-log plot of ^xq&L versus L .
Using Eq. ~6!, we obtain ^xq11&L /^xq&L;Lbx or sx(q11)
2sx(q)5bx , so that the slope of sx(q) as a function of q is
the cutoff exponent bx5]sx(q)/]q . This is, in general, not
true for small q because the integral in Eq. ~6! is dominated
by the lower cutoff. In particular, corrections to scaling of
the type ^xq&L;Lsx(q)F(L) are important for q<tx21. For
instance, when q.tx21, logarithmic corrections give rise to
effective exponents up to very large lattice sizes. Finally,
normalization imposes sx(0)50.
In Fig. 3 we show the results obtained from the moment
analysis of the distribution P(s) for the Manna and the BTW
models. In this case we can use the exact result ^s&;L2,
which implies ss(1)52, as a test for the convergence of our
simulations to the asymptotic scaling regime. This relation is
fulfilled and the ss(q) of the two models are indistinguish-
able for q>1, indicating universal scaling behavior. We ob-
serve small deviations for small q that are due to the non-
universal lower cutoff. By measuring the slope of ss(q), we
obtain bs.2.7. This value is larger than the value reported
in Ref. @15# ~i.e., D.2.5!, where small lattice sizes have
been used. We have repeated the same analysis for the
P(T ,L) and the P(a ,L), and the measured cutoff exponents
FIG. 3. Plot of ss(q) for the BTW and Manna models. The
linear part has slope 2.74. Note the nonuniversal corrections to the
linear behavior expected for q.t21.0.3.b t and ba are reported in Table I. Also in this case the
exponents for the two models share the same values within
error bars, confirming the presence of a single universality
class.
As a final consistency test, we use the data collapse
method in order the check the FSS hypothesis, which states
that rescaling qx[x/Lbx and Pqx[P(x ,L)L
bxtx, the data for
different L must collapse onto universal curves. If FSS is
verified, we can compute the exponent tx from the scaling
relation (22tx)bx5sx(1), which should be satisfied for
enough large sizes. Using the values of bx reported in Table
I and the values obtained for sx(1), we find the exponents tx
to be inserted into the data collapse. For instance, using the
exact result ss(1)52 and the estimated bs52.7, we obtain
ts51.27. The data collapse with these values is satisfactory
for both models ~see Fig. 4!.
In the same way, we obtain very good data collapse for
the Manna model P(a) and P(t) distributions, yielding t t
51.5 and ta51.35. On the other hand, we find that the
BTW data collapses for time and area distributions are not
compatible with the FSS hypothesis. The linear behavior of
the moments analysis, however, ensures that for large sizes
the FSS form must be approached. This result can be ex-
plained if we assume that the scaling in the BTW model
displays subdominant corrections of the form P(x)
5(C1x2t11C2x2t21fl)G(x/xc), where Ci are nonuniver-
sal constants. These corrections are compatible with the lin-
ear behavior at large q , but the decay of the P(x) is not a
simple power law for small x and thus FSS is not obeyed. It
TABLE I. Values of the critical exponents describing the scal-
ing of the cutoff of the distributions for different models in d52.
The results are obtained from the moments analysis ~see text!. Note
that the exponents bs , b t , and ba are usually reported in the lit-
erature as D , z , and d f , respectively.
Model bs b t ba ts
Manna 2.7460.02 1.5060.02 2.0260.02 1.2760.01
BTW 2.7360.02 1.5260.02 2.0160.02 1.2760.01
FIG. 4. Data collapse analysis of the avalanche size distribution
for the Manna and BTW ~inset! models. The values used for the
critical exponents are ts51.27 and bs52.7. Lattice sizes used are
reported in the figure.
RAPID COMMUNICATIONS
PRE 59 R15UNIVERSALITY IN SANDPILESis worth it to remark that the time and area distributions span
over much less order of magnitude than the size distribution,
which could explain why subdominant corrections are more
relevant in the first two cases. Subdominant corrections are
due to higher order operators in the dynamics and do not
determine the universality class, since the asymptotic scaling
behavior is ruled by the leading power.
In summary, we have presented numerical evidence point-
ing toward a single universality class for the Manna and the
BTW models. In particular, we show that previous analyses
@11,14# are not reliable because of systematic biases intro-
duced by the method employed. Further work is needed inorder to quantify the extent of subdominant corrections to
scaling in the BTW model.
We thank A. Barrat, D. Dhar, R. Dickman, S. Krishna-
murthy, E. Marinari, M. A. Mun˜oz, and D. Stauffer for use-
ful discussions and suggestions. The main part of the nu-
merical simulations have been run on the KALIX parallel
computer @18# ~a Beowulf project at Cagliari Physics Depart-
ment!. We thank G. Mula for leading the effort toward orga-
nizing this computer facility. A.V. and S.Z. acknowledge
partial support from the European Network under Contract
No. ERBFMRXCT980183. The Center for Polymer Studies
is supported by the NSF. A.C. acknowledges the hospitality
of the CPS, where this work was initiated.@1# P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 381
~1987!; Phys. Rev. A 38, 364 ~1988!.
@2# Avalanches were recorded in frictional sliding, S. Ciliberto and
C. Laroche, J. Phys. I 4, 223 ~1994!; type II superconductors,
S. Field, J. Witt, F. Nori, and X. Ling, Phys. Rev. Lett. 74,
1206 ~1995!; magnetic systems, G. Durin, G. Bertotti, and A.
Magni, Fractals 3, 351 ~1995!; D. Spasojevic´, S. Bukvic´, S.
Milosevic´, and H. E. Stanley, Phys. Rev. E 54, 2531 ~1996!;
fractures, A. Garciamartin, A. Guarino, L. Bellon, and S. Cil-
iberto, Phys. Rev. Lett. 79, 3202 ~1997!; earthquakes, G.
Gutenberg and C. F. Richter, Ann. Geofis. 9, 1 ~1956!.
@3# A. Vespignani and S. Zapperi, Phys. Rev. Lett. 78, 4793
~1997!; Phys. Rev. E 57, 6345 ~1998!.
@4# S. S. Manna, J. Phys. A 24, L363 ~1991!.
@5# L. Pietronero, A. Vespignani, and S. Zapperi, Phys. Rev. Lett.
72, 1690 ~1994!; A. Vespignani, S. Zapperi, and L. Pietronero,
Phys. Rev. E 51, 1711 ~1995!.
@6# R. Dickman, A. Vespignani and S. Zapperi, Phys. Rev. E 57,
5095 ~1998!; A. Vespignani, R. Dickman, M. A. Mun˜oz, and
S. Zapperi, e-print cond-mat/9806249.
@7# We note that directed models @D. Dhar and R. Ramaswami,
Phys. Rev. Lett. 63, 1659 ~1989!# belong to a different univer-
sality class. The classification is different in one dimension@L. A. N. Amaral and K. B. Lauritsen, Phys. Rev. E 56, 231
~1997!#.
@8# Y. C. Zhang, Phys. Rev. Lett. 63, 470 ~1989!.
@9# A. Dı´az-Guilera, Europhys. Lett. 26, 177 ~1994!; A. Corral and
A. Dı´az-Guilera, Phys. Rev. E 55, 2434 ~1997!.
@10# P. Grassberger and S. S. Manna, J. Phys. ~France! 51, 1077
~1990!; S. S. Manna, J. Stat. Phys. 59, 509 ~1990!; Physica A
179, 249 ~1991!.
@11# A. Ben-Hur and O. Biham, Phys. Rev. E 53, R1317 ~1996!.
@12# K. Christensen, H. C. Fogedby, and H. J. Jensen, J. Stat. Phys.
63, 653 ~1991!.
@13# This classification is not correct, since D. Dhar ~e-print cond-
mat/9808047! has recently shown that the Manna model is
Abelian.
@14# E. Milshtein, O. Biham, and S. Solomon, Phys. Rev. E 58, 303
~1998!; e-print cond-mat/9805206.
@15# M. De Menech, A. L. Stella, and C. Tebaldi, Phys. Rev. E 58,
R2677 ~1998!.
@16# A. Chessa, E. Marinari, A. Vespignani, and S. Zapperi, Phys.
Rev. E 57, R6241 ~1998!.
@17# S. Lu¨beck and K. D. Usadel, Phys. Rev. E 55, 4095 ~1997!; 56,
5138 ~1997!; S. Lu¨beck, ibid. 56, 1590 ~1997!.
@18# For information see the web site http://kalix.dsf.unica.it/


https://air.unimi.it/retrieve/handle/2434/657103/1259910/PhysRevE.59.R12.pdf

Universality in sandpiles

Abstract

Similar works

Full text

Available Versions

IMT Institutional Repository

Crossref

AIR Universita degli studi di Milano