Driving sandpiles to criticality and beyond
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J un Driving Sandpiles to Criticality and Beyond
Anne Fey, Lionel Levine, and David B. Wilson Delft Institute of Applied Mathematics, Delft University of Technology, The Netherlands Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Microsoft Research, Redmond, WA 98052, USA (Dated: December 16, 2009; revised March 4, 2010)A popular theory of self-organized criticality relates driven dissipative systems to systems withconservation. This theory predicts that the stationary density of the abelian sandpile model equalsthe threshold density of the fixed-energy sandpile. We refute this prediction for a wide variety ofunderlying graphs, including the square grid. Driven dissipative sandpiles continue to evolve evenafter reaching criticality. This result casts doubt on the validity of using fixed-energy sandpiles toexplore the critical behavior of the abelian sandpile model at stationarity.
PACS numbers: 64.60.av, 45.70.Cc
In a widely cited series of papers [1–5], Dickman,Mu˜noz, Vespignani, and Zapperi (DMVZ) developed atheory of self-organized criticality as a relationship be-tween driven dissipative systems and systems with con-servation. This theory predicts a specific relationshipbetween the abelian sandpile model of Bak, Tang, andWiesenfeld [6], a driven system in which particles addedat random dissipate across the boundary, and the cor-responding “fixed-energy sandpile,” a closed system inwhich the total number of particles is conserved.After defining these two models and explaining theconjectured relationship between them in the DMVZparadigm of self-organized criticality, we present datafrom large-scale simulations which strongly indicate thatthis conjecture is false on the two-dimensional square lat-tice. We then examine the conjecture on some simplerfamilies of graphs in which we can provably refute it.Early experiments [7] already identified a discrepancy,at least in dimensions 4 and higher, but later work fo-cused on dimension 2 and missed this discrepancy (it isvery small). Some recent papers (e.g., [8]) restrict theirstudy to stochastic sandpiles because deterministic sand-piles belong to a different universality class, but thereremains a widespread belief in the DMVZ paradigm forboth deterministic and stochastic sandpiles [9, 10].Despite our contrary findings, we believe that the cen-tral idea of the DMVZ paradigm is a good one: the dy-namics of a driven dissipative system should in some wayreflect the dynamics of the corresponding conservativesystem. Our results point to a somewhat different re-lationship than that posited in the DMVZ series of pa-pers: the driven dissipative model exhibits a second-orderphase transition at the threshold density of the conser-vative model.Bak, Tang, and Wiesenfeld [6] introduced the abeliansandpile as a model of self-organized criticality; for math-ematical background, see [11]. The model begins with acollection of particles on the vertices of a finite graph.A vertex having at least as many particles as its degree topples by sending one particle along each incident edge. A subset of the vertices are distinguished as sinks: theyabsorb particles but never topple. A single time step con-sists of adding one particle at a random site, and thenperforming topplings until each non-sink vertex has fewerparticles than its degree. The order of topplings does notaffect the outcome [12]. The set of topplings caused byaddition of a particle is called an avalanche.Avalanches can be decomposed into a sequence of“waves” so that each site topples at most once duringeach wave. Over time, sandpiles evolve toward a station-ary state in which the waves exhibit power-law statistics[13] (though the full avalanches seem to exhibit multifrac-tal behavior [14, 15]). Power-law behavior is a hallmarkof criticality, and since the stationary state is reachedapparently without tuning of a parameter, the model issaid to be self-organized critical .To explain how the sandpile model self-organizes toreach the critical state, Dickman et al. [1, 3] introducedan argument which soon became widely accepted: see, forexample, [16, Ch. 15.4.5] and [17–19]. Despite the appar-ent lack of a free parameter, they argued, the dynamicsimplicitly involve the tuning of a parameter to a valuewhere a phase transition takes place. The phase tran-sition is between an active state, where topplings takeplace, and a quiescent “absorbing” state. The param-eter is the density , the average number of particles persite. When the system is quiescent, addition of new par-ticles increases the density. When the system is active,particles are lost to the sinks via toppling, decreasingthe density. The dynamical rule “add a particle whenall activity has died out” ensures that these two densitychanging mechanisms balance one another out, drivingthe system to the threshold of instability.To explore this idea, DMVZ introduced the fixed-energy sandpile model (FES), which involves an explicitfree parameter ζ , the density of particles. On a graphwith N vertices, the system starts with ζN particles atvertices chosen independently and uniformly at random.Unlike the driven dissipative sandpile described above,there are no sinks and no addition of particles. Subse-quently the system evolves through toppling of unstablesites. Usually the parallel toppling order is chosen: ateach time step, all unstable sites topple simultaneously.Toppling may persist forever, or it may stop after somefinite time. In the latter case, we say that the system stabilizes ; in the terminology of DMVZ, it reaches an“absorbing state.”A common choice of underlying graph for FES is the n × n square grid with periodic boundary conditions. Itis believed, and supported by simulations [20], that thereis a threshold density ζ c , such that for ζ < ζ c , the systemstabilizes with probability tending to 1 as n → ∞ ; andfor ζ > ζ c , with probability tending to 1 the system doesnot stabilize. THE DENSITY CONJECTURE
For the driven dissipative sandpile on the n × n squaregrid with sinks at the boundary, as n → ∞ the stationarymeasure has an infinite-volume limit [21], which is a mea-sure on sandpiles on the infinite grid Z . One gets thesame limiting measure whether the grid has periodic oropen boundary conditions, and whether there is one sinkvertex or the whole boundary serves as a sink [21] (seealso [22] for the corresponding result on random span-ning trees). The statistical properties of this limitingmeasure have been much studied [23–25]. Grassbergerconjectured that the expected number of particles at afixed site is 17 /
8, and it is now known to be 17 / ± − [25]. We call this value the stationary density ζ s of Z .DMVZ believed that the combination of driving anddissipation in the classical abelian sandpile model shouldpush it toward the threshold density ζ c of the fixed-energy sandpile. This leads to a specific testable pre-diction, which we call the Density Conjecture. Density Conjecture [4].
On the square grid, ζ c =17 /
8. More generally, ζ c = ζ s .Vespignani et al. [4] write of FES on the square grid,“the system turns out to be critical only for a partic-ular value of the energy density equal to that of thestationary, slowly driven sandpile.” They add that thethreshold density ζ c of the fixed energy sandpile is “theonly possible stationary value for the energy density” ofthe driven dissipative model. In simulations they find ζ c = 2 . ζ c also found the value very close to 17 / ζ c ( Z ) is 2 . / n × n torus. After each ad- n trials estimate of ζ c ( Z n )64 2 . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± .
64 128 256 512 1024 2048 4096 8192 16384
PSfrag replacements nζ c ( Z n ) 2 . . TABLE I: Fixed-energy sandpile simulations on n × n tori Z n .The third column gives our empirical estimate of the thresh-old density ζ c ( Z n ) for Z n . The standard deviation in eachof our estimates of ζ c ( Z n ) is 4 × − . To six decimals,the values of ζ c ( Z ) , . . . , ζ c ( Z ) are all the same. Thedata from n = 64 to n = 16384 are well approximated by ζ c ( Z n ) = 2 . ± × − − (0 . ± . n − . , asshown in the graph. (The error bars are too small to be visi-ble, so the data are shown as points.) The rapid convergenceis due in part to periodic boundary conditions. We concludethat the asymptotic threshold density ζ c ( Z ) is 2 . ζ s ( Z ) is2 . dition, we performed topplings until either all sites werestable, or every site toppled at least once. For determin-istic sandpiles on a connected graph, if every site topplesat least once, the system will never stabilize [26–28]. Werecorded m/n as an empirical estimate of the thresholddensity ζ c ( Z n ), where m was the maximum number ofparticles for which the system stabilized. We averagedthese empirical estimates over many independent trials.We used a random number generator based on the Ad-vanced Encryption Standard (AES-256), which has beenfound to exhibit excellent statistical properties [29, 30].Our simulations were conducted on a High PerformanceComputing (HPC) cluster of computers. PHASE TRANSITION AT ζ c We consider the density conjecture on several otherfamilies of graphs, including some for which we can de-termine the exact values ζ c and ζ s analytically.Dhar [12] defined recurrent sandpile configurations andshowed that they form an abelian group. A consequenceof his result is that the stationary measure for the drivendissipative sandpile on a finite graph G with sinks is theuniform measure on recurrent configurations. The sta-tionary density ζ s ( G ) is the expected total number ofparticles in a uniform random recurrent configuration,divided by the number of non-sink vertices in G .The threshold density ζ c and stationary density ζ s fordifferent graphs is summarized in Table II. The onlygraph on which the two densities are known to be equal graph ζ s ζ c Z Z / = 2 .
125 2 . . . . bracelet / = 2 . . . . . flower graph / = 1 . . . . . . . . ladder graph − √ = 1 . . . . . . . . complete graph / × n + O ( √ n ) × n − O ( √ n log n )3-regular tree / / . . . . TABLE II: Stationary and threshold densities for differentgraphs. Exact values are in bold. is Z [17, 18, 27]. On all other graphs we examined, withthe possible exception of the 3-regular Cayley tree, itappears that ζ c = ζ s .Each row of Table II represents an infinite family ofgraphs G n indexed by an integer n ≥
1. For example,for Z we take G n to be the n × n square grid, and for theregular trees we take G n to be a finite tree of depth n . Assinks in G n we take the set of boundary sites G n \ G n − (note that on trees this corresponds to wired boundaryconditions). The value of ζ s reported is lim n →∞ ζ s ( G n ).The exact values of ζ s for regular trees (Bethe lattices)were calculated by Dhar and Majumdar [31]. The cor-responding values of ζ c we report come from simulations[32]. We derive or simulate the values of ζ s and ζ c for thebracelet, flower, ladder, and complete graphs in [32].As an example, consider the bracelet graph B n , whichis a cycle of n vertices, with each edge doubled (see Fig-ure 1). A site topples by sending out 4 particles: 2 toeach of its two neighbors. One site serves as the sink. Tocompare the densities ζ c and ζ s , we consider the drivendissipative sandpile before it reaches stationarity, by run-ning it for time λ . More precisely, we place λn particlesuniformly at random, stabilize the resulting sandpile, andlet ρ n ( λ ) denote the expected density of the resulting sta-ble configuration. In the long version of this paper [32]we prove Theorem 1 ([32]) . For the bracelet graph B n , in thelimit as n → ∞ ,1. The threshold density ζ c is the unique positive rootof ζ = − e − ζ (numerically, ζ c = 2 . ).2. The stationary density ζ s is / .3. The final density ρ n ( λ ) , as a function of initial den-sity λ , converges pointwise to a limit ρ ( λ ) , where ρ ( λ ) = min (cid:18) λ, − e − λ (cid:19) = ( λ, λ ≤ ζ c − e − λ , λ > ζ c . FIG. 1: The graphs on which we compare ζ s and ζ c : thegrid (upper left), bracelet graph (upper right), flower graph(2 nd row left), complete graph (2 nd row right), Cayley trees(Bethe lattices) of degree d = 3 , , rd row), and laddergraph (bottom). Part 3 of this theorem shows that despite the inequality ζ s = ζ c , a connection remains between the driven dissi-pative dynamics used to define ζ s and the conservativedynamics used to define ζ c : since the derivative ρ ′ ( λ )is discontinuous at λ = ζ c , the driven sandpile under-goes a second-order phase transition at density ζ c . For λ < ζ c , the driven sandpile loses very few particles tothe sink, and the final density equals the initial density λ ; for λ > ζ c it loses a macroscopic proportion of par-ticles to the sink, so the final density is strictly smallerthan λ . As Figure 2 shows, the sandpile continues toevolve as λ increases beyond ζ c ; in particular, its densitykeeps increasing.We are also able to prove that a similar phase transi-tion occurs on the flower graph , shown in Figure 1. In-terestingly, the final density ρ ( λ ) for the flower graph isa decreasing function of λ > ζ c (Figure 2 bottom).Our proofs make use of local toppling invariants onthese graphs. On the bracelet graph, since particles al-ways travel in pairs, the parity of the number of particleson a single vertex is conserved. On the flower graph, thedifference modulo 3 of the number of particles on the twovertices in a single “petal” is conserved.One might guess that the failure of the density con-jecture is due only to the existence of local toppling in- PSfrag replacements ρρ λλ ζ c ζ c ζ c ζ s = 5 / ζ s = 5 / − e − λ − e − λ PSfrag replacements ρρ λλ ζ c ζ c ζ c ζ s = 5 / − e − λ ζ s = 5 / ζ s = 5 / e − λ e − λ FIG. 2: Density ρ ( λ ) of the final stable configuration as afunction of initial density λ on the bracelet graph (top row)and flower graph (bottom row) as the graph size tends toinfinity. A phase transition occurs at λ = ζ c . At first glance(left panels) it appears that the driven sandpile reaches itsstationary density ζ s at this point, but closer inspection (rightpanels) reveals that the final density ρ ( λ ) continues to changeas λ increases beyond ζ c . These graphs are exact. variants, or else to boundary effects from the sinks. Theladder graph (Figure 1) has no local toppling invariants;moreover, it is essentially one-dimensional, so the bulk ofthe graph is well insulated from the sinks at the bound-ary. Nevertheless, we find [32] a small but undeniabledifference between ζ s and ζ c on the ladder graph. CONCLUSIONS
The conclusion of [5] that “FES are shown to exhibitan absorbing state transition with critical properties co-inciding with those of the corresponding sandpile model”should be re-evaluated.In response to this article, several researchers have sug-gested to us that perhaps the density conjecture holds forstochastic sandpiles even if not for deterministic ones.This hypothesis deserves some scrutiny.For the driven dissipative sandpile, there is a transitionpoint at the threshold density of the FES, beyond whicha macroscopic amount of sand begins to dissipate. Thecontinued evolution of the sandpile beyond ζ c shows thatdriven sandpiles have (at least) a one-parameter familyof distinct critical states. While the stationary state hasrightly been the object of intense study, we suggest thatthese additional critical states deserve further attention. [1] R. Dickman, A. Vespignani, and S. Zapperi, Phys. Rev. E , 5095 (1998), cond-mat/9712115.[2] A. Vespignani, R. Dickman, M. A. Mu˜noz, andS. Zapperi, Phys. Rev. Lett. , 5676 (1998),cond-mat/9806249.[3] R. Dickman, M. Mu˜noz, A. Vespignani, and S. Zapperi,Brazilian J. Phys. , 27 (2000), cond-mat/9910454.[4] A. Vespignani, R. Dickman, M. A. Mu˜noz, and S. Zap-peri, Phys. Rev. E , 4564 (2000).[5] M. A. Mu˜noz, R. Dickman, R. Pastor-Satorras,A. Vespignani, and S. Zapperi, in Proc. 6th GranadaSem. on Comput. Physics (2001) cond-mat/0011447.[6] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. ,381 (1987).[7] P. Grassberger and S. S. Manna, J. Phys. France , 1077(1990).[8] J. A. Bonachela and M. A. Mu˜noz, Phys. Rev. E ,041102 (2008), arXiv:0806.4079.[9] R. Vidigal and R. Dickman, J. Stat. Phys. (2005),cond-mat/0403280v1.[10] S. D. da Cunha, R. R. Vidigal, L. R. da Silva, andR. Dickman, (2009), arXiv:0906.1392.[11] F. Redig, Les Houches lecture notes (2005).[12] D. Dhar, Phys. Rev. Lett. , 1613 (1990).[13] D. V. Ktitarev, S. L¨ubeck, P. Grassberger, andV. B. Priezzhev, Phys. Rev. E , 81 (2000),cond-mat/9907157.[14] M. De Menech, A. L. Stella, and C. Tebaldi, Phys. Rev. E , R2677 (1998), cond-mat/9805045.[15] R. Karmakar, S. S. Manna, and A. L. Stella,Phys. Rev. Lett. , 088002 (2005), cond-mat/0312127.[16] D. Sornette, Critical phenomena in natural sciences , 2nded. (Springer-Verlag, 2004) pp. xxii+528.[17] R. Meester and C. Quant, Markov Process. RelatedFields , 355 (2005).[18] A. Fey-den Boer and F. Redig, Markov Process. RelatedFields , 425 (2005), math-ph/0510060.[19] L. T. Rolla and V. Sidoravicius (2009), arXiv:0908.1152.[20] F. Bagnoli, F. Cecconi, A. Flammini, and A. Vespignani,Europhysics Letters , 512 (2003), cond-mat/0207674.[21] S. R. Athreya and A. A. J´arai, Comm. Math. Phys. ,197 (2004).[22] R. Pemantle, Ann. Probab. , 1559 (1991).[23] S. N. Majumdar and D. Dhar, J. Phys. A , L357(1991).[24] V. B. Priezzhev, J. Stat. Phys. , 955 (1994).[25] M. Jeng, G. Piroux, and P. Ruelle, J. Stat. Mech.,P10015 (2006), cond-mat/0609284.[26] G. Tardos, SIAM J. Discrete Math. , 397 (1998) .[27] A. Fey-den Boer, R. Meester, and F. Redig, Ann. Probab. , 654 (2009), arXiv:0710.0939.[28] A. Fey, L. Levine, and Y. Peres, J. Stat. Phys. (2010), arXiv:0901.3805.[29] P. Hellekalek and S. Wegenkittl, ACM Trans. Model.Comput. Simul. , 322 (2003).[30] P. L’Ecuyer and R. Simard, ACM Trans. Math. Software , Art. 22, 40 (2007).[31] D. Dhar and S. N. Majumdar, J. Phys. A23