Numerical Study of the Correspondence Between the Dissipative and Fixed Energy Abelian Sandpile Models
aa r X i v : . [ c ond - m a t . s o f t ] A p r Numerical Study of the Correspondence Between the Dissipative and Fixed EnergyAbelian Sandpile Models
Su.S. Poghosyan , V.S. Poghosyan , V.B. Priezzhev and P. Ruelle Institute for Informatics and Automation Problems NAS of Armenia, 375044 Yerevan, Armenia Institut de Physique Th´eorique, Universit´e catholique de Louvain, B-1348 Louvain-La-Neuve, Belgium Laboratory of Theoretical Physics, Joint Institute for Nuclear Research, 141980 Dubna, Russia
We consider the Abelian sandpile model (ASM) on the large square lattice with a single dissipativesite (sink). Particles are added by one per unit time at random sites and the resulting density ofparticles is calculated as a function of time. We observe different scenarios of evolution dependingon the value of initial uniform density (height) h = 0 , , ,
3. During the first stage of the evolution,the density of particles increases linearly. Reaching a critical density ρ c ( h ), the system changesits behavior sharply and relaxes exponentially to the stationary state of the ASM with ρ s = 25 / ρ c (0) = ρ s and ρ c ( h > = ρ s . Our observations suggest that theequality ρ c = ρ s holds for more general initial conditions with non-positive heights. In parallel withthe ASM, we consider the conservative fixed-energy Abelian sandpile model (FES). The extensiveMonte-Carlo simulations for h = 0 , , , ρ c ( h )coincides with the threshold density ρ th ( h ) of FES. Therefore, ρ th ( h ) can be identified with ρ s ifthe FES starts its evolution with non-positive uniform height h ≤ Keywords : Self-organized criticality, driven dissipativesystems, fixed-energy sandpile, Abelian sandpile model.
I. INTRODUCTION
A long-standing discussion of the comparative criti-cal properties of the dissipative Abelian sandpile model(ASM) [1] and the conservative fixed-energy sandpile(FES) [2] has gained recently renewed impetus due tothe works by Fey et al [3]. The point of discussion in alaconic form can be reduced to the single question: giventhe same lattice with open and closed boundary condi-tions, whether the stationary density of the dissipativeASM ρ s coincides with the threshold density of the FES ρ th ? Using large-scale simulations on the square lattice,the authors of works [3, 4] gave a negative answer to thisquestion and supported their numerical findings by exactsolutions for some graphs of higher symmetry. A moredetailed answer contains a description of the average den-sity of grains ρ ( τ ) for the given density of added particles τ of the dissipative sandpile: ρ ( τ ) = tN , τ = tN . (1)where N is the linear size of the lattice. It was shownin [3] that ρ ( τ ) for the ASM exhibits a transition at thepoint ρ c which coincides with the threshold density ofthe FES: τ = ρ c = ρ th . For τ > ρ c , the system continuesevolution and its density tends to the stationary value ρ s of the ASM when τ → ∞ .In this article, we continue the analysis of the corre-spondence between critical points of the ASM and FES.Confirming the main result of [3] on the transition pointat ρ c = ρ th , we present arguments against the generalclaim ρ th = ρ s for the two-dimensional square lattice.The focus of the work will be the dependence of thetransition point on the initial conditions. Specifically, we study the nearly closed sandpile, namely the ASM onthe square lattice with periodic boundary conditions anda single dissipative site. One of the results suggested bythis study is that there exist initial conditions which leadduring the evolution to the transition point ρ c coincidingwith ρ s = 25 /
8, and the subsequent evolution does notchange this density. On the other hand, our simulationsof the ASM and FES show that the transition point ofthe nearly closed ASM coincides for large lattices withthe threshold density ρ th of the FES for all consideredinitial conditions. To fix our notations, we will specifythe two-dimensional ASM as in [7], i.e. all stable siteshave heights h i ≤ h i ≥ h i = 0) and negativeheights are not considered neither analytically nor nu-merically because they exist in transient states only andvanish when a system approaches criticality. Neverthe-less, these sites may affect the evolution of the ASM froman initial state to the transition point. The dynamics ofthe nearly closed ASM is very close to that of the FES asdissipation through the single site is strongly restricted.Non-ergodicity of the FES causing a dependence of thethreshold density upon the initial conditions has beendiscussed recently in [12]. A conclusion from this dis-cussion was: “the self-organized criticality exhibited bythe ASM has nothing to do with phase transitions in theFES” [12]. Being reasonable, this statement does not ex-clude however specific initial conditions whose evolutionleads the FES exactly to the stationary density of theASM ρ s = 25 / h = 0 , , , ρ as a function of the density ofadded particles τ behaves differently for different h andhas a characteristic kink at ρ c ( h ), also depending on h . For h = 1, the value ρ c (1) = 3 . . . . is close tothe threshold density of the FES obtained by Fey et al[3]. For h = 0, the density ρ increases linearly with τ ,reaches the value ρ s ≃ / ρ s = 25 / h = 0 , , ,
3. The special case h = 0is considered in more details.In parallel with the ASM, we present simulations ofthe FES for finite N × N square lattice with various N and demonstrate the critical behavior of the thresholddensity ρ th ( h ). Namely, we find numerically the finitesize correction 1 /N γ with critical index γ = γ ( h ) thatdepends on the initial conditions. II. THE MODEL AND CRITICAL DENSITY
We consider the standard Abelian sandpile model [1]on the N × N square lattice. A peculiarity of our con-sideration is the nearly closed boundary conditions: thelattice is a torus with a single selected site i . The height z i at any site beside i takes values 0 , , , , i serves as a sink where particles disappear.During a long evolution, the system gets eventually intoa set of recurrent configurations. The theory of recurrentstates of the ASM was developed by Dhar [7]. One ofimportant consequences of this theory is the existence ofa mapping of the set of recurrent configurations onto theset of spanning trees of the same lattice. The mappingallows one to find the height probabilities in the recurrentstate and therefore to find the stationary density ρ s . Ma-jumdar and Dhar [8] calculated the probability of height h = 1 in the thermodynamic limit P = 2 π − π . (2)The probabilities P , P and P = 1 − P − P − P werefound in [9]: P = 12 − π − π + 12 π + I , (3) P = 14 + 32 π + 1 π − π − I − I , (4) P = 14 − π + 4 π + I I . (5)Here I ν , ν = 1 , I ν = 1(2 π ) ZZZZ π i sin β det M ν dα dα dβ dβ D ( α , β ) D ( α , β ) D ( α + α , β + β ) , (6)where D ( α, β ) = 2 − cos( α ) − cos( β ) (7) and M , M are the matrices M = e iα e i ( β + β ) e i ( α − β ) e − iβ /π − e i ( α + α ) e − iα /π − e − i ( α + α ) e iα e iα , (8) M = e iβ e − i ( α + α ) − i ( β + β ) e iβ e − iα e − iα e iα e − i ( α + α ) e iα . (9)Later on, an exact relation between P and P wasproved in [5], thereby eliminating one of the two integrals.A conjecture on the value of the remaining one, based onits numerical evaluation to twelve decimal places, thenled to the following values for the P i [5], P = 14 − π − π + 12 π , (10) P = 38 + 1 π − π , (11)and for the stationary density, ρ s = P + 2 P + 3 P + 4 P = 258 . (12)This value of ρ s was first suggested by Grassberger [11]based on a high precision calculation of the integrals I , I . The present numerical accuracy reaches 10 − ,and forms the basis for the conjecture that ρ s = 25 / ρ s with the threshold density of theFES, ρ th [2]. The threshold density is defined as theaverage maximal number of randomly dropped particlesper site which allows the relaxation of the closed sand-pile to stop. The value ρ th marks the border of stabil-ity: for any ρ > ρ th , an avalanche process started byadding a particle never stops with probability one. Theconjecture ρ th = ρ s was supported by numerical calcula-tions [2]. However, recently, Fey et al [3] have obtained ρ th = ρ s + ∆, where ∆ = 0 . i . An idea is to retrace carefully the process ofadding particles from an initial state to the crirical stateof the ASM. At the initial stage of the process, the prob-ability of dissipation via the sink at i is negligibly smallfor the large lattices and the evolution of the system al-most coincides with the dynamics of the FES. When thedensity approaches a critical value, the large avalanchesappear and the probability of dissipation increases. If thetransition to the dissipative regime is sufficiently sharp,we can identify the transition point with ρ c . Accord-ing to the process described in [3], the system evolveswith further adding of particles and reaches eventuallythe stationary density of the ASM ρ s . The main dif-ference from previous studies is in the initial conditions.Beside the usual initial conditions 1 ≤ h ≤
3, we con-sider the empty lattice h = 0 and show that just thiscondition leads to the transition point ρ c = ρ s . III. SIMULATIONS FOR ASM WITH ONE SINK
Consider the Abelian sandpile model with one sink on a N × N square lattice L with periodic boundary conditionin both directions. The initial configuration is uniform, h i = h with h = 0 , , ,
3. We add particles by oneper unit time at a random site and the resulting den-sity of particles is measured as a function of time. Atthe first stage of the evolution, the system loses negligi-bly small fraction of particles to the sink, and the den-sity of particles increases with time linearly. Eventually,the system reaches some critical density ρ c ( h ) where itstarts losing a macroscopic amount of particles. At thispoint, the density profile changes sharply and the den-sity tends monotonically to the stationary density of theASM ρ s = 25 /
8. At the transition point, the density ρ ( τ ) has a kink which allows one to determine ρ c ( h )with reasonable accuracy. A. h = 3 During the initial stage of evolution, the number ofparticles in the system grows linearly with time τ . Raredissipative events occur only if the cluster of sites with h = 4 is situated in a close vicinity of i and the addedparticle hits there. When the density of particles reachesa threshold level, a kind of percolation appears whichcan be termed a ”weak percolation”. It implies that siteswith h = 4 do not percolate yet, but additional sitesincreasing their height during an avalanche, produce to-gether with them a percolation set. The point of the weakpercolation is rather low ( ρ c ≈ .
09 in our simulations).Therefore, the density of particles continues a slow in-crease despite the big avalanches spreading through thesystem. The weak percolation point is visible in Fig.2 asa kink following the linear part of the function ρ ( τ ). B. h = 2 Due to lowering of the initial density, the percolationpicture of the sites with h = 4 becomes more pronounced,getting close to the usual site percolation. The percola-tion dynamics versus the avalanche dynamics for differ-ent background densities in the FES has been discussedrecently by Park [12]. It was demonstrated that the crit-ical density of transition into the unstable state stronglydepends on the background density and may differ con-siderably both from accepted values of ρ s and ρ th . Our ΤΡ H Τ L FIG. 1: The evolution of the dissipative ASM with one sinkfor h = 3, number of samples – 1000, N = 500. ΤΡ H Τ L FIG. 2: The evolution of the dissipative ASM with one sinkfor h = 2, number of samples – 1000, N = 500. results for h = 2 confirm this conclusion. The den-sity of particles ρ ( τ ) grows strictly linearly in time up tothe transition point ρ c ≈ .
134 which can be associated,due to rare dissipation events, with the critical densityof the FES for the same initial condition. At the tran-sition point, the avalanche mechanism is activated andthe dissipation reduces density towards ρ s . It may seemstrange that lowering of the initial density from h = 3 to h = 2 leads to the increase of the transition point from ρ c ≈ .
09 to ρ c ≈ . h = 4 for h = 2 and the reduced density of parti-cles outside the percolation cluster. Of course, the valueof the transition point cannot be derived directly fromthe site-percolation critical probability P s = 0 . ... , be-cause the distribution of heights outside the percolationcluster remains unknown. C. h = 1 This case, basically, is very similar to the previous one.However, a tendency of increasing of the transition pointwith lowering of the initial density observed in the case h = 2 is reversed. Remarkably, the value ρ c ≈ . N = 500 is very close (af-ter reduction by one) to the threshold value of the FES ρ th = 2 . ± . N = 512. This supports the idea that the nearly closedsandpile behaves almost identically to the FES duringthe non-dissipative part of evolution. Indeed, the simu-lations in [3] were started with the initial conditions ζ = 0for the sandpile model with stable occupation numbers ζ = 0 , , , h = 1 , , , ζ = 0 for the FES in [3] cor-responds to h = 1 in our simulations. Analogously, our h = 0 corresponds to ζ = − D. h = 0 The initial condition h = 0 is special and turns outto be different from the other initial conditions. Thesimilarity between dynamics of the FES and the nearlyclosed ASM observed in the previous cases allows us touse an important statement proved for the FES: if a con- ΤΡ H Τ L FIG. 3: The evolution of the dissipative ASM with one sinkfor h = 1, number of samples – 1000, N = 500. Twohorizontal lines correspond to values ρ th (0) = 3 .
125 and ρ th (1) = 3 . . . . . ΤΡ H Τ L FIG. 4: The evolution of the dissipative ASM with one sinkfor h = 0, number of samples – 1000, N = 500. Twohorizontal lines correspond to values ρ th (0) = 3 .
125 and ρ th (1) = 3 . . . . . figuration stabilizes after an avalanche process, there issome site that has not toppled [13] (see a similar the-orem in [14, 15]). We may expect (and we do observethis in our simulations) that there is a set of sites withnon-zero density which never topple during the evolutionof the nearly closed ASM from the initial state up to thethreshold value where all sites topple ultimately. Then,the occupation number at these sites does not affect thedynamics, giving at the the same time a contribution tothe total density. If so, the threshold point ρ c ≈ . h = 1 can be lowered for h = 0 to the criticalpoint of the ASM. Our simulations confirm this expecta-tion. Adding particles to the empty lattice, we observe asabove the linear growth of density ρ ( τ ) which indicatesalmost non-dissipative character of evolution. The linearpart changes sharply at the threshold value ρ c ≈ . ρ ( τ ) = const .The coincidence of ρ c and ρ s implies that the percolationdynamics does not dominate anymore and the avalancheprocess approaches its standard critical behavior.Among the initial conditions that we have considered, h = 0 is the only one for which the threshold density ρ s appears to be close or equal to the stationary density25 /
8. We expect that this remains true for other initialconditions with non-positive heights.
IV. THE THRESHOLD DENSITY OF FES
Along with the nearly closed sandpile, we consideredthe FES for initial conditions h = 0 , , , ρ c = ρ th and to show the power-law conver-gence of the threshold density to its asymptotical value.To do that, we have considered N × N lattices for vari-ous N , and periodic boundary conditions in both direc-tions. On each of these lattices, we run 10 samples offixed energy sandpile model and measure average thresh-old density ρ th ( N, h ) and its standard deviation. Afterthe extrapolation of ρ th ( N, h ) we found numerically theasymptotic value of the threshold density and its finite-size corrections for large N . The accuracy δρ th ( N, h ) ofthe result ρ th ( N, h ) is estimated via its standard devia-tion Dρ th ( N, h ) by the following formula: δρ th ( N, h ) ∼ Dρ th ( N ) √ h = 0 we found ρ th ( N, ≃ . − . N − + 1 . N − . (14)For h = 1 we found ρ th ( N, ≃ . − . N − + 1 . N − . (15)For h = 2 we found ρ th ( N, ≃ . − . N − + 18 . N − . (16) (cid:144) N Ρ t h H N , L FIG. 5: The threshold densities of the fixed energy sandpilefor h = 0, number of samples – 10 and N = 10 , , . . . , ρ th (0) = 3 .
125 and ρ th (1) = 3 . . . . . (cid:144) N Ρ t h H N , L FIG. 6: The threshold densities of the fixed energy sandpilefor h = 1, number of samples – 10 and N = 10 , , . . . , ρ th (0) = 3 .
125 and ρ th (1) = 3 . . . . . (cid:144) N Ρ t h H N , L FIG. 7: The threshold densities of the fixed energy sandpilefor h = 2, number of samples – 10 and N = 10 , , . . . , For h = 3 we found ρ th ( N, ≃ . . N − + 1 . N − . (17) (cid:144) N Ρ t h H N , L FIG. 8: The threshold densities of the fixed energy sandpilefor h = 3, number of samples – 10 and N = 10 , , . . . , V. CONCLUSION
The results of our numerical simulations support thefollowing three statements: (1) the density of sand heldby the lattice with a single dissipative site (or with dis-sipation rate tending to zero), as a function of the quan-tity of sand deposed, shows a discontinuity at a value ρ c ,which can be identified with the threshold density of thecorresponding FES (with the same initial condition); (2)the critical value ρ c depends on the initial condition, andconsequently, there is presumably also a range of valuesfor the FES threshold density, also depending on the ini-tial condition and potentially on the way sand is dropped;(3) for the specific uniform initial configuration h = 0,and probably for others as well, the critical is close orequal to the stationary density, ρ c = ρ s = 25 /
8. Onecan expect that if the system starts its evolution fromsufficiently deep initial conditions, i.e. from initial con-dition with large negative heights, the stationary value25 / ρ c and ρ th exactly. ACKNOWLEDGMENTS [1] P. Bak, C. Tang and K. Wiesenfeld, Phys. Rev. Lett. ,381 (1987).[2] A. Vespignani, R. Dickman, M.A. Mu˜nnoz, and S. Zap-peri, Phys. Rev. E , 4564-4582 (2000).[3] A. Fey, L. Levine and D.B. Wilson, Phys. Rev. Lett. ,145703 (2010).[4] A. Fey, R. Meester and F. Redig, Ann. Prob. , 654-675(2009).[5] G. Piroux and P. Ruelle, Phys. Lett. B 607 , 188 (2005);M. Jeng, G. Piroux and P. Ruelle, J. Stat. Mech. P10015(2006).[6] E.V. Ivashkevich and V.B. Priezzhev, Physica
A 254 ,97-116 (1998); D. Dhar, Physica
A 369 , 29 (2006).[7] D. Dhar, Phys. Rev. Lett. , 1613 (1990).[8] S.N. Majumdar and D. Dhar, J. Phys. A: Math. Gen. , L357-L362 (1991).[9] V.B.Priezzhev, J. Stat. Phys , 955 (1994).[10] V.S. Poghosyan, S.Y. Grigorev, V.B. Priezzhev and P.Ruelle, Phys. Lett. B , 768 (2008);J. Stat. Mech. P07025 (2010).[11] P. Grassberger, private comunications.[12] Su-Chan Park, arXiv:1001.3359v1 [cond-mat.stat-mech]19 Jan 2010.[13] G. Tardos, SIAM J. Discrete Math. , 397398 (1988).[14] A. Fey, R. Meester and F. Redig, Ann. Prob. , 2, 654-675 (2009).[15] A. Fey, L. Levine and Y. Peres, J. Stat. Phys.138