Emergence of quasi-units in the one dimensional Zhang model
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] N ov Emergence of quasi-units in the one dimensional Zhang model
Tridib Sadhu ∗ and Deepak Dhar † Department of Theoretical Physics,Tata Institute of Fundamental Research,Homi Bhaba Road, Mumbai 400005, India. (Dated: November 30, 2018)We study the Zhang model of sandpile on a one dimensional chain of length L , where a randomamount of energy is added at a randomly chosen site at each time step. We show that in spite ofthis randomness in the input energy, the probability distribution function of energy at a site in thesteady state is sharply peaked, and the width of the peak decreases as L − / for large L . We discusshow the energy added at one time is distributed among different sites by topplings with time. Werelate this distribution to the time-dependent probability distribution of the position of a markedgrain in the one dimensional Abelian model with discrete heights. We argue that in the large L limit, the variance of energy at site x has a scaling form L − g ( x/L ), where g ( ξ ) varies as log(1 /ξ )for small ξ , which agrees very well with the results from numerical simulations. PACS numbers: 64.60.av, 05.65+b
I. INTRODUCTION
After the pioneering work of Bak, Tang and Wiesen-feld (BTW) in 1987 [1], many different models for self-organized criticality have been studied in different con-texts; for review see [2-5]. Of these, models in the generalclass known as Abelian distributed processors have beenstudied a lot, as they share an Abelian property thatmakes their theoretical study simpler [2]. The originalsandpile model of Bak et al. [1], the Eulerian walkersmodel [8], and the Manna model [6] are all members ofthis class. Models which do not have the Abelian prop-erty have been studied mostly by numerical simulations.In this paper, we discuss the Zhang model [7], which doesnot have the Abelian property.In the Zhang model, the amount of energy added at arandomly chosen site at each time step is not fixed, butrandom. In spite of this, the model in one dimensionhas the remarkable property that the energy at a site inthe steady state has a very sharply peaked distributionin which the width of the peak is much less than thespread in the input amount per time step, and the widthdecreases with increasing system size L . This behaviorwas noticed by Zhang using numerical simulations in oneand two dimension [7], and he called it the ‘emergence ofquasi-units’ in the steady state of the model. He arguedthat for large systems, the behavior would be same asin the discrete model. Recently, A. Fey et al. [9] haveproved that in one dimension, the variance of energy doesgo to zero as the length of the chain L goes to infinity,but they did not study how fast it decreases with L .In this paper, we study this emergence of ‘quasi-units’in one dimensional Zhang sandpile by looking at how theadded energy is redistributed among different sites in the ∗ Electronic address: [email protected] † Electronic address: [email protected] avalanche process. We show that the distribution func-tion of the fraction of added energy at a site x ′ reaching asite x after t time steps following the addition is exactlyequal to the probability distribution that a marked grainin the one-dimensional height type BTW model added atsite x ′ reaches site x in time t . The latter problem hasbeen studied recently [10]. We use this to show that thevariance of energy asymptotically vanishes as 1 /L . Wealso discuss the spatial dependence of the variance alongthe system length. In the large L limit, the varianceat site x has a scaling form L − g ( x/L ). We determinean approximate form of the scaling function g ( ξ ), whichagrees very well with the results of our numerical simu-lations.There have been other studies of the Zhang model ear-lier. Blanchard et al. [11] have studied the steady stateof the model, and found that the distribution of energieseven for the two site problem is very complicated, andhas a multi-fractal character. In two dimensions, the dis-tribution of energy seems to sharpen for larger L , butthe rate of decrease of the width is very slow [12]. Mostother studies have dealt with the question as to whetherthe critical exponents of the avalanche distribution inthis model are the same as in the discrete Abelian model[13, 14] . A. Fey et al. ’s results imply that the asymptoticbehavior of the avalanche distribution in one dimension isindeed the same as in the discrete case, but the situationin higher dimension remains unclear [15, 16].The plan of the paper is as follows. In Section II, wedefine the model precisely. In Section III, we show thatthe calculation of the way the energy added at a site isdistributed among different sites by toppling is same asthe calculation of the time-dependent probability distri-bution of the position of a marked grain in the discreteAbelian sandpile model. This correspondence is used inSection IV to determine the qualitative dependence of thevariance of the energy variable at a site on its position x , and on the system size L . We propose a simple ex-trapolation form that incorporates this dependence. Wecheck our theoretical arguments with numerical simula-tions in Section V. Section VI contains a summary andconcluding remarks. A detailed calculation of the solu-tion of an equation, required in Section IV, is added asan Appendix. II. DEFINITION AND PRELIMINARIES
We consider our model on a linear chain of size L . Thesites are labelled by integers 1 to L and a real continuousenergy variable is assigned to each site. Let E ( x, t ) bethe energy variable at site x at the end of the time-step t .We define a threshold energy value E c , same for each site,such that sites with E ( x, t ) ≥ E c are called unstable, andthose with E ( x, t ) < E c are called stable. Starting froma configuration where all sites are stable, the dynamicsis defined as follows.(i) The system is driven by adding a random amount ofenergy at the beginning of every time-step at a randomlychosen site. Let the amount of energy added at time t be∆ t . We will assume that all ∆’s are independent, identi-cally distributed random variables, each picked randomlyfrom an uniform interval 1 − ǫ ≤ ∆ t ≤ ǫ . Let the siteof addition chosen at time t be denoted by a t .(ii) We make a list of all sites whose energy exceeds orbecomes equal to the critical value E c . All these sites arerelaxed in parallel by topplings. In a toppling, the energyof the site is equally distributed to its two neighbors andthe energy at that site is reset to zero. If there is topplingat a boundary site, half of the energy at that site beforetoppling is lost.(iii) We iterate Step (ii) until all topplings stop. Thiscompletes one time step.This is the slow driving limit, and we assume that allavalanche activity stops before the next addition event.In this limit, the model is characterized by two parame-ters ǫ and E c . In the limit ǫ = 0, and 1 < E c ≤
2, themodel reduces to the discrete case, where the behavioris well understood [17]. For non-zero but small ǫ , thebehavior does not depend on the precise value of E c . Infact, starting with a recurrent configuration of the pile,and adding energy at some chosen site, we get exactlythe same sequence of topplings for a range of values of E c [9]. To be precise, for any fixed initial configuration,and fixed driving sequence (of sites chosen for additionof energy), whether a site x topples at time t or not isindependent of E c , so long as we have 1+ ǫ < E c ≤ − ǫ .In the following, we assume for simplicity that E c = 3 / ≤ ǫ ≤ / E ( x, t ) = 0 and allother sites have energy in the range 1 − ǫ ≤ E ( x, t ) ≤ ǫ .The position of the empty site is equally distributedamong all the lattice points. There are also some re-current configurations in which all sites have energy E ( x, t ) ≥ − ǫ . In such cases, we shall say that the sitewith zero energy is the site L + 1. Then, in the steady state, there is exactly one site with energy equal to 0,and the L + 1 different positions of the site are equallylikely.If E c does not satisfy the inequality 1 + ǫ < E c ≤ − ǫ , this simple characterization of the steady stateis no longer valid. However, our treatment can be easilyextended to those cases. Since the qualitative behaviorof the model is the same in all cases, we restrict ourselvesto the simplest case here.It is easy to see that the toppling rules are in generalnot Abelian. For example, start with a two site model inconfiguration (1 . , .
0) and E c = 1 .
5. The final configu-ration would be (1 . , , . . , . et al. [9] have shownthat only in one dimension, for 1 + ǫ < E c , the Zhangmodel has a restricted Abelian character, namely, thatthe final state does not depend on the order of topplingswithin an avalanche. However, topplings in two differentavalanches do not commute. III. THE PROPAGATOR, AND ITS RELATIONTO THE DISCRETE ABELIAN MODEL
It is useful to look at the Zhang model as a perturba-tion about the ǫ = 0 limit. For sufficiently small ǫ , giventhe site of addition and initial configuration, the top-pling sequence is independent of ǫ . It is also independentof the amount of energy of addition ∆ t , and is same asthe model with ǫ = 0, which is the 1-dimensional Abeliansandpile model with integer heights (hereafter referred tosimply as ASM, without further qualifiers). We decom-pose the energy variables as E ( x, t ) = Nint[ E ( x, t )] + ǫη ( x, t ) , (1)where Nint refers to the nearest integer value. Then theinteger part of the energy evolves as in the ASM. Wewrite ∆ t = 1 + ǫu t , for all t. (2)Here u t is uniformly distributed in the interval [ − , +1].The linearity of energy transfer in toppling implies thatthe evolution of the variables η ( x, t ) is independent of ǫ . Thus, η ( x, t ) is a linear function of u t ; the precisefunction depends on the sequence of topplings that tookplace. These are determined by the sequence of additionsites { a t } up to the time t , and the initial configuration C . These together will be called the evolution historyof the system up to time t , and denoted by H t . Weassume that at the starting time t = 0, the variables η ( x, t = 0) are zero for all x , and the initial configurationis a recurrent configuration C of the ASM. Then, fromthe linearity of the toppling rules, we can write η ( x, t ) asa linear function of { u t ′ } for 1 ≤ t ′ ≤ t , and we can writefor a given history H t , η ( x, t |{ u t } , H t ) = t X t ′ =1 G ( x, t | a t ′ , t ′ , H t ) u t ′ . (3)This defines the matrix elements G ( x, t | a t ′ , t ′ , H t ). Thesecan be understood in terms of the probability distribu-tion of the position of a marked grain in the ASM asfollows. Consider the motion of a marked grain in theone dimensional height type BTW model. We start withconfiguration C and add grains at sites according to thesequence { a t } . All grains are identical except the oneadded at time t ′ , which is marked. In each toppling, themarked grain jumps to one of its two neighbors with equalprobability. Consider the probability that the markedgrain will be found at site x after a sequence of relax-ation processes at time t . We denote this probability asProb( x, t | a t ′ , t ′ , H t ). From the toppling rules in both themodels, it is easy to see that G ( x, t | a t ′ , t ′ , H t ) = Prob( x, t | a t ′ , t ′ , H t ) . (4)Averaging over different histories H t , we get the proba-bility that a marked grain added at x ′ = a t ′ at time t ′ is found at a position x at time t ≥ t ′ in the steadystate of the ASM. Denoting the latter probability byProb ASM ( x, t | x ′ , t ′ ), we get G ( x, t | x ′ = a t ′ , t ′ , H t ) = Prob ASM ( x, t | x ′ , t ′ ) , (5)where the over bar denotes averaging over different his-tories H t , consistent with the specified constraints. Here,the constraint is that H t must satisfy a t ′ = x ′ . At otherplaces, the constraints may be different, and will be spec-ified if not clear from the context.We shall denote the variance of a random variable ξ by V ar [ ξ ]. From the definition in Eq. (1), it is easy to showthat V ar [ E ( x, t )] = L/ ( L + 1) + ǫ V ar [ η ( x, t )] . (6)Different u t are independent random variables, also in-dependent of H t and have zero mean. Let V ar [ u t ] = σ .For the case when u t has a uniform distribution between − σ = 1 /
3. Then, from Eq. (3), weget
V ar [ η ( x, t )] = σ t X t ′ =1 G ( x, t | a t ′ , t ′ , H t ) . (7)As t → ∞ , the system tends to a steady state, and theaverage in the right hand side of Eq. (7) becomes a func-tion of t − t ′ . Also, for a given t ′ , all values of a t ′ areequally likely. We define F ( x, τ ) ≡ L lim t ′ →∞ X x ′ G ( x, t ′ + τ | x ′ , t ′ , H t ) . (8) Then, for large L , in the steady state ( t large), the vari-ance of energy at site x is 1 /L + ǫ Σ ( x ), whereΣ ( x ) = lim t →∞ V ar [ η ( x, t )] = σ ∞ X τ =0 F ( x, τ ) . (9)We define Σ to be the average of Σ ( x ) over x .Σ = 1 L X x Σ ( x ) . (10)Evaluation of G ( x, t | x ′ , t ′ , H t ) for a given history H t andaveraging over H t is quite tedious for t > G ,the problem has been studied in the context of residencetimes of grains in sand piles, and some exact results areknown in specific cases [10]. For G , the calculations aremuch more difficult. However, some simplifications occurin large L limit. We discuss these in the next section. IV. CALCULATION OF Σ ( x ) IN LARGE- L LIMIT
In order to find the quantity F ( x, τ ) in Eq. (8), wehave to average G ( x, t | x ′ , t ′ , H t ) over all possible histo-ries H t , which is quite difficult to evaluate exactly. How-ever, we can determine the leading behavior of F ( x, τ ) inthis limit.We use the fact that the path of a marked grain inthe ASM is a random walk [10]. Consider a particle thatstarts away from the boundaries, at x ′ = ξL , with L large, and 0 < ξ <
1. If it undergoes r ( H t ) topplingsbetween the time t ′ and t = t ′ + τ under some particularhistory H t , then its probability distribution is approxi-mately a Gaussian, centered at x ′ with width √ r . Then,we have G ( x, t | x ′ , t ′ , H t ) ≃ p πr ( H t ) exp (cid:18) − ( x − x ′ ) r ( H t ) (cid:19) . (11)Using this approximation for G , summing over x ′ , we get X x ′ G ( x, t | x ′ , t ′ , H t ) ≃ p πr ( H t ) . (12)Thus, we have to calculate the average of 1 / p r ( H t ) overdifferent histories. Here r ( H t ) was defined as the num-ber of topplings undergone by the marked grain. Differ-ent possible trajectories of a marked grain, for a givenhistory, do not have the same number of topplings. How-ever, if the typical displacement of the grain is muchsmaller than its distance from the end, differences be-tween these are small, and can be neglected. There aretypically O ( L ) topplings per grain per avalanche in themodel, and a grain moves a typical distance of O ( √ L )in one avalanche. Then, we can approximate r ( H t ) by N ( x ′ ), the number of topplings at x ′ .Let the number of topplings at x ′ at time steps τ =0 , , , . . . be denoted by N , N , N , . . . . Then, N ( x ′ ) = N + N + N + · · · . It can be shown that the num-ber of topplings in different avalanches in the one dimen-sional ASM are nearly uncorrelated (In fact the correla-tion function between N i and N j varies as (1 /L ) | i − j | .).By the central limit theorem for sum of weakly correlatedrandom variables, the mean value of N grows linearlywith τ , but the standard deviation increases only as √ τ .Then, for τ ≫
0, the distribution is sharply peaked aboutthe mean, and h / √ N i ≃ / p h N i .Clearly, for τ ≫ h N i = τ ¯ n ( x ′ ), where ¯ n ( x ′ ) is themean number of topplings per avalanche at x ′ in theASM, given by ¯ n ( x = ξL ) = Lξ (1 − ξ ) / . (13)The upper limit on τ for the validity of the above ar-gument comes from the requirement that the width ofthe Gaussian be much less than the distance from theboundary, (without any loss of generality, we can assumethat ξ < /
2, so that it is the left boundary ), else wecannot neglect events where the marked grain leaves thepile. This gives p τ ¯ n ( x ) ≪ ξL , or equivalently, τ ≪ ξL. Thus we get, F ( x, τ ) ≃ C L [ τ Lξ (1 − ξ )] − / , for 0 ≪ τ ≪ ξL, (14)where C is some constant.Also, we know that for τ ≫ L , the probability that thegrain stays in the pile decays exponentially as exp( − τ /L )[10]. Thus, G , and also G will decay exponentially with τ , for τ ≫ L . Thus, we have, for some constants C and a , F ( x, τ ) ≃ C L exp( − aτ /L ) , for τ ≫ L. (15)It only remains to determine the behavior of F ( x, τ ), for ξL ≪ τ ≪ L . In this case, in the ASM, there is asignificant probability that the marked grain leaves thepile from the end. This results in a faster decay of G , andhence of F with time. We argue below that the behaviorof the function F ( x, τ ) is given by F ( x, τ ) ∼ C Lτ , for ξL ≪ τ ≪ L, (16)where C is some constant. This can be seen as follows:Let us consider the special case when the particle starts ata site close to the boundary. Then ¯ n ( x ) is approximatelya linear function of x for small x . Its spatial variationcannot be neglected, and Eq. (12) is no longer valid. Wewill now argue that in this case G ( x, t ′ + τ | x ′ , t ′ ) ≃ x ′ τ − exp( − x/τ ) , for 0 ≪ τ ≪ L. (17)The time evolution of P rob
ASM ( x, t | x ′ , t ′ ) in Eq. (5) iswell described as a diffusion with diffusion coefficient pro-portional to ¯ n ( x ) which is the mean number of topplingsper avalanche at x in the ASM [10]. For understandingthe long-time survival probability in this problem, we can equivalently consider the problem in a continuous-timeversion: consider a random walk on a half line wheresites are labelled by positive integers, and the jump rateout of a site x is proportional to x . A particle starts atsite x = x at time t = 0. If P j ( t ) is the probability thatthe particle is at j at time t , then the equations for thetime-evolution of P j ( t ) are, for all j > ddt P j ( t ) = ( j +1) P j +1 ( t )+( j − P j − ( t ) − jP j ( t ) . (18)The long time solution starting with P j (0) = δ j,x is P j ( t ) ≃ x t − exp( − j/t ) (19)for t ≫ x and large j . The probability that the particlesurvives till time t decreases as 1 /t for large t . We havediscussed the calculation in the Appendix.Using Eq. (5), we see that G ( j, t ′ + τ | x , t ′ ) scales as x /τ . It seems reasonable to assume that G will scaleas G . Then, each term in the summation for F ( x, τ )in Eq. (8) scales as x /τ , and there are τ such terms,as the sum over x has an upper cutoff proportional to τ , and so F ( x, τ ) varies as 1 /τ for L ≫ τ ≫ x . Thisconcludes the argument.We can put these three limiting behaviors into a singlefunctional form that interpolates between these, as F ( x, τ ) ≃ L K exp( − aτ /L ) τ + B p τ Lξ (1 − ξ ) , (20)where K , a and B are some constants. In Section V,we will see that results from numerical simulation areconsistent with this phenomenological expression.Using this interpolation form in Eq. (9), and con-verting the sum over τ to an integration over a variable u = τ /L , we can writeΣ ( x = ξL ) ≃ σ L Z ∞ du K exp( − au ) u + B p uξ (1 − ξ ) . (21)This integral can be simplified by a change of variable au = z , givingΣ ( x = ξL ) ≃ Kσ L I (cid:16) B ′ p ξ (1 − ξ ) (cid:17) , (22)where K, B ′ are constants, and I ( y ) is a function definedby I ( y ) = 2 Z ∞ dz exp( − z ) z + y . (23)It is easy to verify that I ( y ) diverges as log(1 /y ) for small y . In particular, we note that the exponential term in theintegral expression for I ( y ) has a significant contributiononly for z near 1. We may approximate this by droppingthe exponential factor, and changing the upper limit of
0 1 2 3 4 -0.4 -0.2 0 0.2 0.4 P L ( E ) / √ L (E-1.0) √ L GaussianL=1000L=500L=200
FIG. 1: (Color online) Scaling collapse of the probability dis-tribution P L ( E ) of energy per site in the steady state fordifferent systems of size 200, 500 and 1000. The distributionis well described by a Gaussian of width 0 . the integral to 1. The resulting integral is easily done,givingΣ ( x = ξL ) ≃ K ′ σ L log B ′ p ξ (1 − ξ ) ! , (24)where K ′ is some constant. Averaging Σ ( x ) over x , weget a behavior Σ ( x ) ≃ /L . Of course, the answer isnot exact, and one could have constructed other inter-polation forms that have the same asymptotic behavior.We will see in the next Section that results from numer-ical simulations for Σ ( x ) can be fitted very well to thephenomenological expression in Eq. (24). V. NUMERICAL RESULTS
We have tested our non-rigorous theoretical argumentsagainst results obtained from numerical simulations. InFig. 1, we have plotted the probability distribution P L ( E ) of energy at a site, averaged over all sites. Weused L = 200, 500 and 1000, and averaged over 10 dif-ferent configurations in the steady state. We plot thescaled distribution function P L ( E ) / √ L versus the scaledenergy ( E − ¯ E ) √ L . A good collapse is seen, which veri-fies the fact that the width of the peak varies as L − / .The dependence of the variance of E ( x, t ) on x is plot-ted in Fig. 2 for systems of length 200, 300 and 400.The data was obtained by averaging over 10 avalanches.We plot ( L + λ )Σ ( x ) /σ versus x eff /L eff , where x eff differs from x by an amount δ to take into account thecorrections due to end effects. Then, for consistency, L is replaced by L eff = L + 2 δ . For the specific choice of λ = 5 ± δ = 1 . ± .
2, we get a good collapse of thecurves for different L . We also show a fit to the proposed ( L + λ ) Σ ( x ) / σ x eff /L eff L=400L=300L=200Theoretical function
FIG. 2: (Color online) Scaling collapse of Σ ( x ) /σ at site x for systems of different length L . ( L + λ ) Σ ( x ) / σ x eff /L eff L=400L=300L=200Theoretical function
FIG. 3: (Color online) The same plot in Fig. 2 resolved moreat the left boundary of the model and taking x axis in logscale. interpolation form in Eq. (24), with K ′ = 1 . ± . B ′ = 1 . ± .
2. We see that the fit is very good.In order to check the logarithmic dependence of Σ ( x )on x for small x , we re-plot the data in Fig. 3 usinglogarithmic scale for x . We get a good collapse of thedata for different L , supporting our proposed dependencein Eq. (24). VI. CONCLUDING REMARKS
To summarize, we have studied the emergence of quasi-units in the one-dimensional Zhang sandpile model. Thevariance of energy variables in the steady state is gov-erned by the balance between two competing processes.The randomness in the drive i.e., the energy of addition,tends to increase the variance in time. On the other hand,the topplings of energy variables tend to equalize the ex-cess energy by distributing it to the nearby sites. Thereare on an average O ( L ) topplings per avalanche. Hence,in one dimension there are, on an average, O ( L ) topplingsper site per avalanche. For large system size, the secondprocess dominates over the first and the variance becomeslow. We have shown that the variance vanishes as 1 /L with increasing system size and the probability distribu-tion of energy concentrates around a non-random valuewhich depends on the energy of addition. We have alsoproposed a functional form for the spatial dependence ofvariance of energy which incorporates the correct limit-ing behaviors, and matches very well with the numericaldata.An interesting question is whether one can extendthese arguments to the two-dimensional Zhang model.In this case, there are several peaks in the distribution ofenergies at a site, but there are some numerical evidencesfor the sharpening of the peaks as the system size is in-creased. However, as the number of topplings per sitevaries only as log L , the width is expected to decreasemuch more slowly with L , and the fluctuation effects canbe much stronger. This remains an open question forfurther study. APPENDIX
Here we discuss the solution of the Eq. (18) for thestarting values given by P j ( t = 0) = α j − . (A.1)We start with an ansatz P j ( t ) = b t exp( − a t j ), where both a t and b t are functions only of t . This form satisfies theEq. (18) for all j , t >
0, if a t and b t satisfy da t dt = 2 − e a t − e − a t , (A.2) db t dt = b t ( e − a t − e a t ) . (A.3)To solve the Eq. (A.2), we first make a change of variable z = e − a t . In terms of z , the equation becomes dz/dt =(1 − z ) , which can be easily integrated to give e − a t = t + A − t + A , (A.4)where A is an integration constant. To satisfy the initialcondition in Eq. (A.1), we choose A = (1 − α ) − . (A.5) Similarly, to solve the equation for b t , we use the form of e − a t given in Eq. (A4) and get db t dt = b t − t + A )( t + A )( t + A − . (A.6)This can be integrated to give b t = B ( t + A )( t + A − , (A.7)where B is an integration constant. Then the probabilitycan be written as P j ( t ) = B ( t + A − j − ( t + A ) j +1 . (A.8)To satisfy the initial condition at t = 0, we choose theintegration constant B = (1 − α ) − . Then, with thesevalues of A and B , we have the solution for all j , t > P j ( t ) = [(1 − α ) t + α ] j − [(1 − α ) t + 1] j +1 = φ j ( α, t ) , say . (A.9)Now, as φ j ( α, t ) satisfies the Eq. (18), ψ j,n ( α, t ) = 1( n − ∂ n − φ j ( t ) ∂α n − (A.10)will also satisfy the equation for any natural number n .In addition, ψ j,n ( α = 0 , t = 0) = δ j,n . (A.11)Hence, we see that the solution of the Eq. (18), startingwith P j ( t ) = δ j,n at t = 0 is P j ( t ) = ψ j,n ( α = 0 , t ) = 1( n − ∂ n − φ j ( α, t ) ∂α n − | α =0 , (A.12)for all j , t >
0, where φ j ( α, t ) is given in Eq. (A.9) and n is any natural number.It can be shown that for large t and j , the solutionasymptotically becomes P j ( t ) = nt − exp( − j/t ). [1] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. , 381 (1987); Phys. Rev. A , 364 (1988).[2] D. Dhar, Physica A , 29 (2006).[3] E. V. Ivashkevich and V. B. Priezzhev, Physica A , 97 (1998).[4] D. L. Turcotte, Rep. Prog. Phys. , 1377 (1999).[5] H. J. Jensen, Self-Organized Criticality - Emergent Com-plex Behavior in Physical and Biological Systems, Cam- bridge Lecture Notes in Physics 10 , (Cambridge Univer-sity Press, 1998 ).[6] S. S. Manna, J. Phys. A , L363 (1991).[7] Y. C. Zhang, Phys. Rev. Lett. , 470 (1989).[8] V. B. Priezzhev, A. Dhar, D. Dhar, S Krishnamurthy,Phys. Rev. Lett. , 5079 (1996).[9] A. Fey, R. Meestar, C. Quant, F. Redig,arXiv:math-ph/0701029.[10] D. Dhar and P. Pradhan, J. Stat. Mech.: Theor. Exp.(2004) P05002.[11] P. Blanchard, B. Cessac, T. Kr¨uger, J. Stat. Phys. , 307 (1997).[12] I. M. Janosi, Phys. Rev. A , 769 (1990).[13] S. Lubeck, Phys. Rev. E , 1590 (1997).[14] E. Milshtein, O. Biham, and S. Solomon, Phys. Rev. E , 303 (1998).[15] R. Pastor-Satorras and A. Vespignani, Eur. Phys. J. B , 197 (2000).[16] A. Giacometti and A. Diaz-Guilera, Phys. Rev. E , 247(1998).[17] P. Ruelle and S. Sen, J. Phys. A25