Restricted sandpile revisited

a r X i v :c o n d -m a t /0601415v 1 [c o n d -m a t .s t a t -m e c h ] 18 J a n 2006Restricted sandpile revisited

Ronald Dickman ?Departamento de F′?sica,ICEx,Universidade Federal de Minas Gerais,Caixa Postal 702,30161-970Belo Horizonte,Minas Gerais,Brazil (Dated:February 6,2008)Abstract I report large-scale Monte Carlo studies of a one-dimensional height-restricted stochastic sand-pile using the quasistationary simulation method.Results for systems of up to 50000sites yield estimates for critical exponents that di?er signi?cantly from those obtained using smaller systems,but are consistent with recent predictions derived from a Langevin equation for stochastic sandpiles [Ramasco et al.,Phys.Rev.E 69,045105(R)(2004)].This suggests that apparent violations of universality in one-dimensional sandpiles are due to strong corrections to scaling and ?nite-size e?ects.?email:dickman@?sica.ufmg.br

I.INTRODUCTION

Sandpile models are the prime example of self-organized criticality(SOC)[1,2],a control mechanism that forces a system with an absorbing-state phase transition to its critical point[3,4],leading to scale invariance in the apparent absence of parameters[5].SOC in a slowly-driven sandpile corresponds to an absorbing-state phase transition in a model having the same local dynamics,but a?xed number of particles[3,6,7,8,9,10].The latter class of models is usually designated as?xed-energy sandpiles(FES)or conserved sandpiles.Continuous absorbing-state phase transitions characterized by a nonconserved order parameter(activity density)coupled to a conserved?eld that does not di?use in the absence of activity,are expected to de?ne a universality class[11].This class,referred to as C-DP(that is,a model-C version,in the sense of Halperin and Hohenberg[12],of directed percolation,or DP),appears to be distinct from that of directed percolation[13].

In recent years considerable progress has been made in characterizing the critical proper-ties of conserved stochastic sandpiles,although no complete,reliable theory is yet at hand. As is often the case in critical phenomena,theoretical understanding of scaling and univer-sality rests on the analysis of a continuum?eld theory or Langevin equation(a nonlinear stochastic partial di?erential equation)that reproduces the phase diagram and captures the fundamental symmetries and conservation laws of the system.Important steps in this direc-tion are the recent numerical studies of a Langevin equation[13,14]for C-DP.(The latter which appears to incorporate the essential aspects of stochastic sandpiles.)The critical ex-ponent values reported in Ref.[13]are in good agreement with simulations of conserved lattice gas(CLG)models[19,20],which exhibit the same symmetries and conservation laws as stochastic sandpiles.

The Langevin equation exponents are also consistent with the best available estimates for stochastic sandpiles in two dimensions[13],with the exception of the exponentθgoverning the initial decay of the order parameter.(The discrepancy regardingθlikely re?ects strong corrections to short-time scaling in sandpiles,due to long memory e?ects associated with initial density?uctuations[15].)Pending a better understanding of this question,it appears that stochastic sandpiles are consistent with C-DP in two dimensions.In the one-dimensional case,however,there is a signi?cant discrepancy between the Langevin equation results and those for sandpile models.

Speci?cally,analysis of the Langevin equation for C-DP yields,in one dimension,the order-parameter critical exponent valueβ=0.28(2),while previous studies[15,16,17,18] of stochastic sandpiles furnish values near0.40for this exponent.There are also smaller discrepancies for other critical exponents.If this discrepancy were to persist,one would be forced to conclude that the proposed Langevin equation misses some essential aspect of sandpiles(at least in the one-dimensional case),or that not all models with the same symmetries and conserved quantities belong to the same universality class.In an e?ort to clarify the situation,I apply the recently devised quasistationary simulation method [21,22,23]to the restricted-height sandpile introduced in Ref.[16].

The balance of this paper is organized as follows.In Sec.II I de?ne the model and summarize the simulation method.Numerical results are analyzed in Sec.III,and in Sec. IV I discuss the?ndings in the context of universality.

II.MODEL

I study the”independent”version of the model introduced in Ref.[16].The system,a continuous-time,restricted-height version of Manna’s stochastic sandpile[24],is de?ned on a ring of L sites.The con?guration is speci?ed by the number of particles,z i=0,1,or2,at each site i.Sites with z i=2are active,and have a toppling rate of unity.The continuous-time Markovian dynamics consists of a series of toppling events at individual sites.When site i topples,two particles attempt to move randomly(and independently)to either i?1 or i+1.(The two particles may both try to jump to the same neighbor.)Each particle transfer is accepted so long as it does not lead to a site having more than two particles.(If the target site is already doubly occupied the particle does not move.Thus an attempt to send two particles from site j to site k,with z k=1,results in z k=2and z j=1.)The next site to topple is chosen at random from a list of active sites,which is updated following each event.The time increment associated with each toppling is?t=1/N A,where N A is the number of active sites just prior to the event.

Any con?guration devoid of doubly occupied sites is absorbing.Although absorbing con?gurations exist for particle densities p=N/L≤1,the critical value p c(above which activity continues inde?nitely)appears to be strictly less than unity.In Ref.[16]the model was studied in the site and pair mean-?eld approximation(which yield a continuous

phase transition at p c=0.5and0.75,respectively,in one dimension),and via Monte Carlo simulation using system sizes of up to5000sites.The latter yield the estimates p c= 0.92965(3),β/ν⊥=0.247(2),z=ν||/ν⊥=1.45(3)andβ=0.412(4).A similar value,β=0.42(1),was obtained in Ref.[17]using a series of cluster approximations(of up to11 sites),combined with Suzuki’s coherent anomaly analysis[25].

The studies reported here employ the quasistationary(QS)simulation method,which, due to increased e?ciency in the critical region,permits a tenfold increase in the system size as compared to Ref.[16].The QS method,described in detail in[21],provides a just sampling of asymptotic(long-time)properties,conditioned on survival.In practice this is accomplished by maintaining(and gradually updating)a set of con?gurations visited during the evolution;when a transition to the absorbing state is imminent the system is instead placed in one of the saved con?gurations.Otherwise the evolution is exactly that of a “standard”simulation algorithm such as used in Ref.[16].

III.SIMULATION RESULTS

I performed two sets of studies using the QS method.The?rst is used to determine the QS order parameter(de?ned as the factionρof active sites),the moment ratio m= ρ2 /ρ2, and the mean lifetimeτof the quasistationary state,in the immediate vicinity of the critical point p c,for system sizes L=1000,2000,5000,10000,20000and50000.(The QS lifetime is taken as the mean number of time steps between successive attempts to visit the absorbing state.)A second set of simulations is used to study the supercritical regime(p>p c)for system sizes L=10000,20000and50000.(For p substantially larger than p c,the lifetime is much larger than the simulation time,so that the system never visits the absorbing state, and the QS method becomes identical to a standard simulation.)

Each realization of the process is run for109time steps;averages are taken in the QS regime,which necessitates discarding an initial transient that ranges from106time steps (for L=1000)to108time steps(for L=50000).The number of saved con?gurations ranges from1000(for L=1000)to400(for L=50000).The list updating probability p rep ranges from10?3(for L=1000)to5×10?6(for L=50000).During the initial relaxation period p rep is increased by a factor of ten to erase the memory of the initial con?guration.

I?rst discuss the studies focusing on the critical region.As in Ref.[16],I study,for each

system size,a series of particle number values N,chosen so that p=N/L lies immediately above or below p c.Since the particle density can only be varied in steps of1/L,estimates for properties at intermediate values of p are obtained via interpolation.The results of the QS simulations were found to agree,to within uncertainty,with the corresponding results of conventional simulations[16],for L=1000,2000and5000.The criterion for criticality is power-law dependence ofρandτon system size,i.e.,the familiar relationsρ～L?β/ν⊥and τ～L z,and constancy of the moment ratio m with L.The most sensitive indicator turns out to be the order parameterρ.Using the data for system sizes5000-50000,I rule out p values that yield a statistically signi?cant curvature of the graph of lnρversus ln L.This results in the estimate p c=0.929780(7).(For the remainder of the analysis p c is?xed at this value and is no longer available as an adjustable parameter.)The associated exponent isβ/ν⊥=0.213(6),where the uncertainty represents a contribution(±0.005)due to the uncertainty in p c and a small additional uncertainty in the linear?t to the data.Simulation results forρas a function of L,for various densities near p c,are shown in Fig.1;curvature of the plots for o?-critical values is evident in the inset.

The data for the QS lifetimeτfurnish a similar but somewhat less precise estimate, p c=0.929777(17).Fitting the data for L=5000-50000,using the p c interval obtained from the analysis ofρ,I?nd z=1.50(4).The moment ratio m is also useful for setting limits on p c.As shown in Fig.2,this quantity appears to grow with system size for p

The present estimate for p c is signi?cantly greater than that found in Ref.[16],although the di?erence amounts to about0.01%.The results for the exponent z are consistent,but the present study yields a substantially(16%)lower estimate forβ/ν⊥than reported previously. The present result for m c is also substantially lower than the value1.1596(4)reported in Ref.[16].These di?erences highlight the strong?nite-size corrections a?ecting stochastic sandpiles.

I turn now to the results for the order parameter in the supercritical regime.Fig.3shows that the data for system sizes10000,20000and50000are well converged for?=p?p c≥10?3,that is,?nite-size e?ects are only present nearer the critical point.Evidently,the data are not consistent with a simple power law of the formρ～?β.Indeed this departure

from the familiar behavior of the order parameter was already noted(with data for smaller systems)in Ref.[16].In the latter work the power law was“restored”by introducing a size-dependent critical density p c(L)?p c,∞?Const/L1/ν⊥,leading to a series of estimates for the critical exponentβthat increase systematically with L,apparently converging to β=0.412(4).With the present data,which are converged over a broader range of?values, I?nd that shifting the critical value does not lead to an apparent power law.

One is therefore left to conclude that either the order parameter does not obey power-law scaling,or that there are unusually strong corrections to scaling.Including a correction to scaling term,one has

ρ～?β 1+A?β′ (1) so that there are now three adjustable parameters,β,β′and A.Even with a reasonably large number of data points(18for L=10000),this induces a huge range of variation in the exponentβ.Decent?ts can be obtained with values as low asβ=0.1and as large as 0.3.

To resolve this di?culty I return to the data in the immediate vicinity of p c.These data can be used to determine the correlation length exponentν⊥in the following manner. Finite-size scaling implies that for p?p c,the moment ratio obeys the relation

m(?,L)?F m(L1/ν⊥?).(2) where F m is a scaling function.This implies that

?m

?p

p c∝L1/ν⊥,(4)

and similarly for the derivative of lnτat the critical point.The derivatives are evaluated numerically as follows.For each value of L studied,data for?ve values of p clustered around p c are?t with a cubic polynomial;the derivative of the polynomial is then evaluated at p c. The resulting derivatives are plotted in Fig.4;clean power laws are observed,leading to

ν⊥=1.362(7),1.323(14)and1.372(21),using the data for lnρ,m and lnτ,respectively. Pooling these results yields the estimateν⊥=1.355(18).Then,using the values forβ/ν⊥and z reported above,I?ndβ=0.289(12)andν||=2.03(8).

Using this value forβ,the data for the order parameter in the supercritical regime can be ?t using the correction to scaling form,Eq.1,with parametersβ′=0.446and A=1.3505. For?=0.1,the correction term A?β′in Eq.1is0.48,showing that there are sizeable deviations from a pure power law.It is usual to verify scaling by seeking a data collapse, plottingρ?=Lβ/ν⊥ρversus??=L1/ν⊥?.For?>0.001the order parameter does not follow a pure power law and so the data cannot collapse.It is nevertheless of interest to construct such a scaling plot(Fig.5).Although the data do not collapse over most of the range,they do collapse in the interval?1≤?≤1.A linear?t to the data in this interval yields a slope of0.27(1).This is close to theβvalue obtained from the?nite-size scaling analysis,suggesting that simple scaling is restricted to a narrow interval very near the critical point.

IV.DISCUSSION

A study of the one-dimensional restricted-height stochastic sandpile using quasistationary simulations permits study of systems an order of magnitude larger than previously studied, and yields critical properties di?erent than those obtained previously.In the case of the critical density,the small change(about0.01%)from the previous estimate may be attributed to?nite-size e?ects,which are known to a?ect sandpile models strongly.

Of greater concern are critical exponent values,since they de?ne the universality class of the model.Since there is every reason(based on symmetry considerations)to expect the re-stricted sandpile to belong to the same universality class as the unrestricted version(indeed, this seems well established in two dimensions[16]),I collect,in Table I,critical exponent val-ues from various studies of stochastic sandpiles,C-DP and the conserved threshold transfer process(CTTP),also expected to belong to the same class.

The overall conclusion from Table I is that studies using smaller lattices yield values in the range0.38-0.42for the exponentβ(Ref.[18]is however an exception),and that the large-scale simulation of Ref.[20],the numerical study of the C-DP?eld theory[13] and the present work yield a consistent set of results,withβ?0.29.(A similar value has

been found for a modi?ed conserved lattice gas model[27].)Although the system size(4000 sites)used in the?eld theory simulations is not large,one should note that each‘site’in such a simulation may represent a region comprising many lattice sites in the original model. Compared with the earlier sandpile simulations,the distinctive feature of the present work may not be system size,but the fact that here the exponentβis determined via?nite-size scaling at the critical point,rather than from the usual analysis of the order parameter in the supercritical regime.Indeed,it is easy to see from Fig.3that data for?=p?p c in the range10?3-10?1will yield larger estimates forβ.(The same observation applies to the CAM analysis[17],which essentially probes the shape of the functionρ(?)at some distance from the critical point?=0.)I observe a simple power-law behavior,and data collapse for various lattice sizes,only in a restricted range of the scaling variable??=L1/ν⊥?.

Also included in Table I are exponent values for one-dimensional directed percolation [28].The values obtained in Refs.[13]and[20],as well as in the present work,are not very di?erent from those of DP.A clear di?erence from DP scaling was however demonstrated in Ref.[14],where the initial decay exponent for one-dimensional C-DP is found to be θ=0.125(2),as opposed to0.1595(1)for DP.The rather substantial di?erences found here inβ/ν⊥,and in the moment ratio m(1.142(8)for the restricted sandpile compared with 1.1736(1)for DP[21]),lend further support to the conclusion that the C-DP/stochastic sandpile universality class is distinct from that of directed percolation,as is evidently the case in two dimensions.(This despite the result[29],that when suitably modi?ed to include ‘sticky grains’,sandpiles fall generically in the DP class.)

In summary,I have applied the quasistationary simulation method to a one-dimensional restricted-height stochastic sandpile,and obtained results consistent with recent studies of C-DP.This supports the assertion that the latter class includes stochastic sandpiles,as would be expected on the basis of symmetry and conservation laws.

Acknowledgements

I am grateful to Hugues Chat′e,M′a rio de Oliveira and Miguel Angel Mu?n oz for helpful discussions and comments on the manuscript.This work was supported by CNPq and Fapemig,Brazil.

[1]P.Bak,C.Tang and K.Wiesenfeld,Phys.Rev.Lett.59,381(1987);Phys.Rev.A38,364

(1988).

[2] D.Dhar,Physica A263(1999)4,and references therein.

[3]R.Dickman,M.A.Mu?n oz,A.Vespignani,and S.Zapperi,Braz.J.Phys.30,27(2000).

[4]M.A.Mu?n oz,R.Dickman,R.Pastor-Satorras,A.Vespignani,and S.Zapperi,in Modeling

Complex Systems,Proceedings of the6th Granada Seminar on Computational Physics,J.

Marro and P.L.Garrido,eds.,AIP Conference Proceedings v.574(2001).

[5]G.Grinstein,in Scale Invariance,Interfaces and Nonequilibrium Dynamics,NATO Advanced

Study Institute,Series B:Physics,vol.344,A.McKane et al.,Eds.(Plenum,New York,1995).

[6] C.Tang and P.Bak,Phys.Rev.Lett.60,2347(1988).

[7]M.Paczuski,S.Maslov,and P.Bak,Phys.Rev.E53,414(1996).

[8] A.Vespignani and S.Zapperi,Phys.Rev.Lett.78,4793(1997);Phys.Rev.E57,6345(1998).

[9]R.Dickman,A.Vespignani and S.Zapperi,Phys.Rev.E57,5095(1998).

[10] A.Vespignani,R.Dickman,M.A.Mu?n oz,and Stefano Zapperi,Phys.Rev.Lett.81,5676

(1998).

[11]M.Rossi,R.Pastor-Satorras,and A.Vespignani,Phys.Rev.Lett.85(2000)1803.

[12]P.C.Hohenberg and B.J.Halperin,Rev.Mod.Phys.49,435(1977).

[13]J.J.Ramasco,M.A.Mu?n oz,and C.A.da Silva Santos,Phys.Rev.E69,045105(R)(2004).

[14]I.Dornic,H.Chat′e,and M.A.Mu?n oz,Phys.Rev.Lett.94,100601(2005).

[15]R.Dickman,M.Alava,M.A.Mu?n oz,J.Peltola,A.Vespignani,and S.Zapperi,Phys Rev.

E64,056104(2001).

[16]R.Dickman,T.Tom′e,and M.J.de Oliveira,Phys.Rev.E66,016111(2002).

[17]R.Dickman,Phys.Rev.E66,036122(2002).

[18]S.L¨u beck,Phys Rev.E66,046114(2002).

[19]R.Pastor-Satorras and A.Vespignani,Phys.Rev.E62,R5875(2000).

[20]J.Kockelkoren and H.Chat′e,e-print cond-mat/0306039.

[21]M.M.de Oliveira and R.Dickman,Phys.Rev.E71,016129(2005).

[22]R.Dickman and M.M.de Oliveira,Physica A357,134(2005).

[23]M.M.de Oliveira and R.Dickman,eprint:cond-mat/0601163.

[24]S.S.Manna,J.Stat.Phys.59,509(1990);S.S.Manna,J.Phys.A24,L363(1991).

[25]M.Suzuki and M.Katori,J.Phys.Soc.Jpn.55,1(1986);M.Suzuki,M Katori,and X.Hu,

J.Phys.Soc.Jpn.56,3092(1987).

[26]S.L¨u beck and P.C.Heger,Phys.Rev.E68,056102(2003).

[27] C.E.Fiore and M.J.de Oliveira,Braz.J.Phys.,in press.

[28]I.Jensen,J.Phys.A32,5233(1999).

[29]P.K.Mohanty and D.Dhar,Phys Rev.Lett.89,104303(2002).

TABLE I:Summary of exponent values for one-dimensional models in the C-DP universality class. L max denotes the largest system size studied.Abbreviations:CAM:coherent anomaly method;FT:?eld theory.

Modelβz

Manna[15]0.42(2) 1.66(7)

Manna[26] 1.39(11)

CTTP[18]0.38(2) 1.66(7)

Rest.Manna[16]0.416(4) 1.50(9)

Rest.Manna CAM[17]0.41(1)

C-DP[20]0.29(2) 1.55(3)

C-DP FT[13]0.28(2) 1.47(4)

Rest.Manna(present work)0.289(12) 1.50(4)

0.2521

FIGURE CAPTIONS

FIG.1.Stationary order parameter versus system size for particle densities(bottom to top) p=0.92977,0.92978and0.92979.Inset:ln L0.213ρversus ln L for the same set of particle densities.

FIG.2.Moment ratio m versus system size for particle densities(top to bottom)p=0.92976, 0.92978and0.92980.

FIG.3.Stationary order parameter versus?=p?p c for system sizes(top to bottom) L=104,2×104and5×104.

FIG.4.Derivatives of(lower to upper)lnτ,lnρand m with respect to particle density, evaluated at p c,versus system size.The slope of the straight line is0.734.

FIG.5.Scaled densityρ?versus scaled distance from critical point??,as de?ned in text. System sizes:104(open squares);2×104(?lled squares);5×104(diamonds).