Analytic approach to stochastic cellular automata: exponential and inverse power distributions out of Random Domino Automaton
aa r X i v : . [ n li n . C G ] F e b Analytic approach to stochastic cellular automata:exponential and inverse power distributions out ofRandom Domino Automaton
Mariusz Bia lecki ∗ and Zbigniew Czechowski † Institute of Geophysics, Polish Academy of Sciencesul. Ks. Janusza 64, 01-452 Warszawa, Poland (Dated: September 12, 2018)We introduce the stochastic domino cellular automaton model exhibiting avalanches. Dependingof the choice of the parameters, the model covers wide range of properties: various types of expo-nential and long tail (up to inverse-power) distributions of avalanches are observed. The stationarystate of automaton is described by a set of nonlinear discrete equations derived in an exact wayfrom elementary combinatorial arguments. These equations allow to derive formulas explaining bothvarious exponential and inverse power distributions relating them to values of the parameters. Theexact relations between the state variable of the model (moments) are derived in two ways: fromdirect arguments and from the set of equations. Excellent agreement of the obtained analyticalresults with numerical simulations is observed.
PACS numbers: 45.70.Ht, 02.50.Ey, 05.65.+b, 91.30.PxKeywords: avalanches, discrete equations, exact solutions, solvable models, stochastic cellular automata, toymodels of earthquakes
I. INTRODUCTION
There is a significant interest in constructing simplestochastic models with avalanches reflecting propertiesof various natural phenomena, but only part of themhave the advantage of being solvable. One way to makemodels analytically tractable is to use the mean-field ap-proximation. However, requirements of exact results canforce further simplifications in a model. For example,in terms of self-organized criticality, in spite of the factthat the mean field theory of such critical phenomenawas already proposed [1], Dhar made further simplifica-tions for distinguished BTW model [2] (a unified viewof mean-field picture of stochastic self-organized criticalmodels is presented in [3]). We point out also that con-sidering structure of abelian algebras leads to analyticalresults in some sand-pile and related stochastic models[4, 5]. Solvable simple models can also be constructed inthe field of directed percolation; for an extremely simpleone, see [6].Another way towards solvability is to consider stochas-tic properties as non-essential and study deterministicmodels; see such an approach for sand-pile in [7]. In thecontext of deterministic cellular automata, an analyticapproach related to solvability is investigated as part ofthe theory of integrable systems and may employ sophis-ticated methods. [8–13]. We underline that constructionof integrable cellular automata is definitely not an easytask.Our aim here is to follow the ideas of analytic approachand apply them to construction of stochastic automata. ∗ Electronic address: [email protected] † Electronic address: [email protected]
In the article we propose and analyse in a direct elemen-tary way the Random Domino Automaton (RDA) - a newslowly driven systems exhibiting avalanche phenomena.We prefer here an elementary self-contained approach tothe description of the automaton and we make no useof applicable Markov processes terminology. We stresshere, all results, except of equation (28) (see [14] for de-tails), have elementary derivation in the text below. Onthe other hand, in spite of its simple formulation and be-ing analytically tractable, automaton covers wide rangeof behaviours depending on the choice of the parame-ters. The application of the RDA model for studying Itoequation is investigated in our parallel papers [15, 16].An inspiration for defining the rule for the domino au-tomaton comes from very simplified view of earthquakes.It corresponds to two tectonic plates moving with rel-ative constant velocity. The wedge may be irregularlyrough, and relative motion can be locked in some places,producing stress accumulation. Beyond some thresholdof stress, a relaxation took place. The size of relaxationdepends on the nearby accumulated stress. RDA is in-spired by earthquakes; however, a direct reproducing ofrealistic-like behaviour of such complicated phenomenais obviously beyond its scope in the present form. Never-theless, construction and analysis of models is one of pri-mary aims in geosciences [17, 18]. There are many verysimplified cellular automata models focusing on specificfeatures of the investigated behaviour in the field. Herewe point out a sequence of papers [19–21], where someinteresting cellular automata models were presented.Finally, as an unexpected property we mention an in-teresting link between stochastic cellular automata andinteger sequences (see [22] or The On-Line Encyclopediaof Integer Sequences). It is known, how to obtain Cata-lan numbers out of the bond directed percolation on asquare lattice [23]. Our cellular automaton in specificcase leads to Motzkin numbers (for details see [14]).Plan of the article is as follows. In Section II we in-troduce a definition of random domino automaton withrebound parameters. Section III contains full combina-torial derivation of equations for the distribution of clus-ters, which describes stationary state of the automaton.Several exact formulas - balance relations, average size ofclusters and avalanche, general formula for all moments- are displayed. We finish with introducing special formof rebound parameters, which leads to two distinguishedspecial cases studied in details in next Section IV. Thesecases correspond to exponential and inverse-power distri-butions as shown below. The obtained analytical resultsare compared with simulations. Section V gives a resumeof results presented in the article as well as directions forfuture work. Appendix A contains derivation of comple-mentary set of equations, describing empty clusters ofthe automaton.
II. DEFINITION OF RANDOM DOMINOAUTOMATON
In the random domino automaton model, the spaceconsists of discrete number N of cells on a line, and weassume periodic boundary conditions. Each cell may bein one of two states: empty state, when it is empty, or oc-cupied state, when contains a ball. The evolution of theautomaton is given by the following update procedureperformed in each discrete time step. A ball is addedto the system to the randomly chosen cell and we as-sume each cell to be equally possible. If the chosen cellis empty, there are two possibilities: it becomes occupiedwith probability ν or the ball is rebounded with probabil-ity (1 − ν ) leaving the state of the automaton unchanged.If the chosen place is already occupied, there are alsotwo possibilities: the ball is rebounded with probability(1 − µ ) or with probability µ the incoming ball triggers arelaxation. By relaxation we mean: balls from the chosencell and from all its adjacent occupied cells are removed.Thus, the relaxation produces an avalanche of size equalto the number of cells changing their state. Then theupdate procedure repeats in the next time step. All pos-sibilities are shown schematically on the diagram below. •ւ ց (1 − ρ ) ρ ւ ց . . . . . . . . . • . . . ւ ց ւ ց ν (1 − ν ) µ (1 − µ ) ւ ց ր տ ւ ց ր . . . • . . . . . . . . . . . . ↓ . . . . . . • . . .. . . • . . .An avalanche is represented by a symbol ↓ • . An exampleof relaxation of the size three is presented in the diagrambelow. ↓ • time = t · · · • • • • • • · · · time = t + 1 · · · • • ↓ ↓ ↓ • · · ·• • • The name ’domino automaton’ comes from the fol-lowing interpretation. Occupied cells are represented bystanding domino blocks and empty cells are representedby extra space between them. If the incoming ball hits anempty space, a domino block is added there with prob-ability ν . If the incoming ball strikes a domino block,with probability µ it falls down the chosen block and allits adjacent neighbours (on both sides) up to the gap(empty cell). Then thefallen dominoes are removed andthe procedure repeats.Presented defining rules of the domino automaton canbe easily modified leading to various extensions of thesystem. We mention some possibilities (like geometry,capacity of cells, many kinds of balls, different triggeringrules) in section V. However, the aim of this work is toprovide detailed description of the ”core” case describedabove. III. EQUATIONS OF RANDOM DOMINOAUTOMATONA. Notation
A state of the automaton is defined by altered se-quences of occupied and empty cells. A sequence of i consecutive occupied cells is called the cluster of length i (shortly i -cluster); a sequence of i consecutive empty cellsis called empty cluster of length i , where i = 1 , , . . . , N .The fixed size of the lattice N is assumed to be finite, butbig enough to make limit for the size of clusters negligi-ble. Denote the number of i -clusters by n i , the numberof empty i -clusters by n i , the total number of clusters by n , and the density (a number of occupied cells dividedby the number of all cells N ) by ρ . Then, it follows that n = X i ≥ n i = X i ≥ n i and ρ = 1 N X i ≥ n i i. (1)For higher densities ρ , the probability of relaxation isalso relatively higher, and the density is more likely todecrease. For smaller densities, triggering of an avalancheis relatively less probable and the density tends to growup. Hence, the variations of the density (see Fig.1) aresubjected to ’v-shape’ potential. The behaviour of theautomaton is described below under the assumption thatit is in a quasi-equilibrium state and the variables used inthe equations, like density and others, are average valuesand do not depend on time. The applicable description ofthe system as a Markov process is postponed to anotherpaper. i Ρ H i L i Ρ H i L FIG. 1: Examples of simulation results for time series of den-sity ρ ( i ) of the 1D domino automaton in case µ/ν = 1 withlattice size N = 500 and in case µ/ν = . i with lattice size N = 4000. The parameter i numbers avalanches. We consider the following settings for rebound param-eters. The coefficient µ = µ i does not depend on theposition of the chosen cell, but may depend on the size i of the chosen cluster. The parameter ν is fixed to be con-stant. There are two cases of special interest (discussedin detail in Section IV). One with µ = β = const refersto equal probability of provoking an avalanche for eachoccupied cell, while case µ = µ i = δi , where δ = const,refers to equal probability of triggering an avalanche foreach cluster, irrelevantly to its length. These two respec-tive cases are related to exponential and long-tail distri-butions of clusters. B. The balance equation for ρ . The stationarity condition requires the flow-in to beequal to flow-out, hence there must be a balance betweenlosses and gains in the number of occupied cells. In asingle time step, the number of occupied cells may stayunchanged (when ball is reflected), or may increase byone (an empty cell become occupied), or decrease by i (in a case of an avalanche of size i ). The expected valueof increase of density is equal to the probability that anempty space becomes occupied, namely ν (1 − ρ ). Theprobability of relaxation of a size i is µ i in i N . Since anypossible size i can trigger an avalanche, the stationaritygives ν (1 − ρ ) = 1 N ( X i ≥ µ i n i i ) or N = X i ≥ n i i ( µ i ν i + 1) . (2) C. The balance equation for n . In a single time step, the number of clusters may in-crease only when a new 1-cluster is created. The chanceis proportional to the number of interior cells in emptyclusters having the length three and bigger, as depictedin the diagram below. · · · | • | | ( i −
2) = interior z }| { | · · · · · · | | | {z } i | • | · · · Hence, the probability is ∼ X i ≥ ν ( i − i n i iN = ν (cid:18) (1 − ρ ) − nN + n N (cid:19) . (3)Losses in the total number of clusters may appear intwo ways: by joining of two clusters (separated by anempty 1-cluster) and by triggering an avalanche. In thefirst case, the probability is ∼ ν n N , in the second it is ∼ P i ≥ µ i n i iN . Hence, the balance equation for the numberof clusters is of the form(1 − ρ ) N − n = X i ≥ µ i ν n i i. (4) D. Balance equations for n i s. We underline the assumption used in the derivationof the equations below: clusters are distributed indepen-dently , by which we mean that the length of the ”next”cluster does not depend on the length of the ”previous”one. In other words, our investigations are done up tothe order of clusters (and order of empty clusters). Thepresented approach may be regarded as a generalizationof percolation, where subsequent cells are treated as in-dependent. Below, in subsection IV A, we compare bothapproaches numerically.To write down equations for the numbers of clustersof length i , i.e., for n i s, we consider all possibilities oflosses of such clusters as well as creation of them, andnext we claim that on the average the gains and losses ofrespective values compensate each other, as required bythe stationarity conditions. Losses.
There are two ways to destroy an i -cluster:by enlarging and by provoking the avalanche dependingon the cell where an incoming particle is thrown. · · · | | • | • | · · · | • | {z } i | | · · · (a) Enlarging. For any cluster there are two cells adja-cent to its ends, so the probability is ∼ ν n i N .
If the single empty cell is between two clusters of thelength i it is counted twice - it decrease the number of i -clusters by two. (b) Relaxation. In this case it is enough to knock outany of the occupied cell of the cluster, so the probabilityis ∼ µ i in i N .
Gains.
There are in general two possibilities to create i -cluster: enlarging ( i − (a) Enlarging. Case i = 1 was already considered and theprobability is given by formula (3). For i ≥ i − i is possible if the adjacentempty cluster is of a size bigger then one. Hence theprobability is ∼ ν n i − N P i ≥ n i n = 2 ν n i − N (cid:18) − n n (cid:19) . where the multiplier 2 counts left and right cases. (b) Joining two clusters. Two smaller clusters: one ofsize k ∈ { , , . . . , ( i − } and the other of size ( i − − k will be joined if the ball fills an empty cell between them. · · · | | k z }| { • | · · · | • | | ( i − − k ) z }| { • | • | · · · | • | {z } i | | · · · The probability is proportional to the number of empty1-clusters between k -cluster and ( i − − k )-cluster, hence ∼ ν n N i − X k =1 n k n · n i − − k n . The dot in the multiplication above underlines the inde-pendence assumption for the order of clusters. The lastformula introduces also an extra quadratic nonlinearityinto the system; so far, nonlinearity was present through n (and also n ).Finally the following set of equations for n i s is derived n = 1 µ ν + 2 (cid:0) (1 − ρ ) N − n + n (cid:1) , (5) n = 22 µ ν + 2 (cid:18) − n n (cid:19) n , (6) n i = 1 µ i ν i + 2 ×× n i − (cid:18) − n n (cid:19) + n i − X k =1 n k n i − − k n ! (7)for i ≥
3, where n = P i ≥ n i and ρ = N P i ≥ in i .Variable n can be cancelled from the above set byconsidering the balance for empty clusters and deriving, in analogous way, the respective set of equations for vari-ables n i (see Appendix A for details)). The result isgiven by equation (A1), namely n = 2 n (cid:16) νn P i ≥ µ i n i i (cid:17) . (8)The above, substituted into equations (5)-(7), form aclosed set for variables { n , n , . . . } . Then, equations(A1) and (A2) allow us to find n k for all k = 1 , , . . . .In conclusion, we obtained closed, nonlinear set ofequations for the distribution of the clusters in the cellu-lar automaton. E. Moments
The balance equations (2) and (4) (as well as equationsfor higher weighted moments of n i ) can be also obtainedfrom the above set (5)-(7).For n i , i = 1 , , . . . , we define a moment of order γ by m γ = 1 N X i ≥ n i i γ . (9)Then, the density ρ is equal to the first moment m , andthe normalized number of clusters nN is equal to the zeromoment m . We define also weighted moments b m γ asfollows b m γ = 1 N X i ≥ µ i ν n i i γ . (10)In this notation, equation (8) takes the following form n = 2 m N b m m . (11)Removing denominators from the set of equations (5)-(7),multiplying both sides by i z and performing summationwith respect to i , one obtains X i ≥ µ i ν i z +1 n i + 2 X i ≥ i z n i = (1 − ρ ) N − n + n + (12)+2(1 − n n ) X i ≥ ( i + 1) z n i + n n X i ≥ i − X k =1 i z n k n i − − k . Finally, after using formula (11), changing the order ofsummation (in indices i and k ), introducing variable j byformula i = k +1+ j , expanding binomials and performingsums, one obtains b m z +1 = 1 − m − m + 2 z − X k =0 (cid:18) zk (cid:19) m k + (13)+ 23 m + 2 b m l + p ≤ z X l,p =1 (cid:18) zl + p (cid:19)(cid:18) l + pl (cid:19) m l m p . The above equation for z = 0 and z = 1: b m = 1 − m − m , (14) b m = 1 − m , (15)are the balance of n equation (4) and the balance of ρ equation (2), respectively. Equations of moments for sub-sequent values of z allow to choose particular forms ofrebound parameters and isolate interesting cases as pre-sented below for formula (18).Definition of moments (9) leads to particularly neatform of expressions for the average size of a cluster < i > = P i ≥ n i i P i ≥ n i = m m , (16)and the average size of an avalanche < w > = P i ≥ µ i n i i P i ≥ µ i n i i = b m b m = 1 − m − m − m . (17)For rebound parameters in the form µ i = δi σ , ν = const., (18)where δ , σ and θ = δν are constants, equations (14) and(15) take the form m + θm − σ = 1 , (19)2 m + m + θm − σ = 1 . (20)Hence, cases σ = 0 and σ = 1, considered in detailsbelow, are very suitable for exact analysis. There is alsoan interesting behaviour when σ >
2, but we consider itin another paper.
IV. SPECIAL CASESA. Case µ = β = const : equal probability oftriggering relaxation for each occupied cell. To get the value of density from equation (2) in thejust considered case, extra relations are required. As thefirst step, consider the percolation approximation [24] n i = c (1 − ρ ) ρ i , (21)which means that cells are treated as independent. Fromequation (1) it follows that c = N . Then equation (2)gives ρ = 2 / (3 βν + q βν ) + 4(1 − βν )). Any value ofdensity ρ ∈ (0 ,
1) can be obtained for suitable βν ∈ (0 , ∞ ).The percolation approximation in the case βν = 1 givesthe value of ρ = , which differs by several percent fromthe value < ρ > ≃ . i n H i L FIG. 2: Values of n i in case µ = β = const from simulation(dots) compared with values computed from equations (solidline) and the approximation ke − γi (dashed line). Lattice size N = 500 and β/ν = 1.TABLE I: The average density < ρ > , the average size of acluster < i > , and the average size of an avalanche < w > from simulation results, equations of the model and percola-tion approximation for case µ/ν = 1.Simulation Equations Percolation a < ρ > < i > < w > a Based on the density balance. between the adjacent cells, since the relaxations takes outthe whole cluster, so treating cells as independent, likein formula (21), gives only very rough estimation of thedensity ρ .In the general case (not percolation approximation), µ = β = const, the balance equation for n gives thefollowing formula for an average size of the cluster < i > = N ρn = 2 ρ − ρ (1 + βν ) , (22)and the balance of ρ gives an average size of the avalanche < w > = νβ (cid:18) − ρρ (cid:19) . (23)A comparison of values of the average density < ρ > ,the average size of cluster < i > and the average sizeof the avalanche < w > in exemplary case µ/ν = 1 ob-tained from numerical solution of equations, from perco-lation approximation results (based on the density bal-ance), and from simulation is presented in Table I. Thelattice size is fixed as N = 500.Approximate solution for n i s in the case µ = β = const can be obtained by substituting in equation (7)the following formula n i = ke − γi for i = 3 , , . . . (24)where k and γ are some constants. For i >> e − γ = νβ kn n , or n i = k (cid:18) νβ kn n (cid:19) i . (25) i n H i L FIG. 3: Exact values of n i obtained from the model (line)in case µ i = δ/i compared with exemplary simulation data(dots). Lattice size N = 4000 and δ/ν = 1 / < ρ > and the average sizeof a cluster/avalanche < i > = < w > from simulation resultsand equations of the model for the case µ i = δi , for θ = 1 / < ρ > < i > = < w > The value of constant k can be found from equation (1).The above result indicates a close relation to the perco-lation dependence for n i (see equation (21)). Approxi-mation (25) works well even for n , as may be seen inFigure 2. B. Case µ = µ i = δ/i : equal probability of triggeringrelaxation for each cluster. For µ i = δi , where δ = const, the balance for ρ - equa-tion (2) - takes the form ρ = 1( θ + 1) (26)where θ = δν , which relates the density with the ratio ofcoefficients δ and ν only. Since θ ∈ (0 , ∞ ), any density ρ ∈ (0 ,
1) may be realized. The balance for n - equation(4) - is reduced to (1 − ρ ) N = (2+ θ ) n. Together equations(2) and (4) give n = N θ ( θ + 1)( θ + 2) . (27)The maximal number of clusters is (1+ √ − N ≈ . N ,as obtained for θ = √ < i > = < w > = 1 + 2 θ . Table II contains values of average density and averagecluster/avalanche size obtained from equations and sim-ulation in an exemplary case where θ = 0 .
25. The size oflattice is set to N = 4000 in order to avoid restrictions forsize of avalanches, which are relatively bigger comparingto the previous case.A simple form of n = 2 n/ (3 + 2 θ ) together with bal-ance equations (26) and (27) alows to reduce the set ofequations (5)-(7) to the solvable recurrence. The formu-las for n i s are rational functions of θ and N only. Theexplicit form of solutions and relation to the Motzkinnumbers is presented in [14].Figure 3 presents exact values of n i s compared in log-log scale with the simulation result in the exemplary casedescribed above. The slope of the distribution of n i ap-proximates inverse-power distribution in the middle partof a range of sizes i ; for bigger values of i , values of n i decrease quicker. This property agrees with convergenceof moments of n i we use above. For value of θ closer to 0the straight part of the plot extends for bigger values of i . In the limit case θ θN = const −→
0, the distribution is givenby the following power law [14] n i +1 ∼ i . (28) V. CONCLUSIONS.
We proposed and studied properties of simple stochas-tic cellular automaton with avalanches inspired by anextremely simplified model of earthquakes. The maingoal was to construct a simple model with ”transpar-ent” mathematical structure, which would enable ob-taining exact results yet covering a wide range of be-haviours. The rule of the automaton allows to derivefrom the first principles (using elementary combinatorics)the set of equations describing the average values of themodel parameters. These equations leads to equationsfor moments and to neat formulas for average clusterand avalanche sizes in terms of the zero and the firstmoments.After analysis of moment equations, we considered indetail two cases which differ by the form of parameter µ responsible for triggering of avalanches. We obtainedexponential type and inverse-power type distributions ofclusters within uniform framework. (The distribution ofavalanches w i is easily obtained from the formula w i = µ i n i i P i ≥ µ i n i i ∼ µ i n i i .) The quasi-equilibrium assumptionseems to be justified by simulation results.The exponential case µ = β = const is not fully analyt-ically solvable; nevertheless, the closed set of equationsleads to several exact relations. Moreover, an approxi-mate formula for n i was derived, and intimate relationto percolation was pointed out. Hence, the proposed ap-proach can be regarded as an extension of percolationapproximation results, and can provide basis for exacttreatment of other models. The automaton rule clearlystates that there is a ”coupling” between the adjacentcells, since the relaxations takes out the whole cluster, sotreating cluster – not cells – as independent is a substan-tial improvement in the presented case.The case µ i = δi is fully solvable, and all variablesare expressed as rational functions of ratio δν of reboundparameters only. The shape of the distribution of clus-ters n i approximate an inverse-power distribution in thevarious ranges of variable i , depending on chosen param-eters. In the limit case, distribution tends to ∼ n − .The critical density is equal to 1. The limit case of thepresented automaton reduces the set of equations to therecurrence, which leads to known integer sequence - theMotzkin numbers (see [14]). This result establishes anew, remarkable link between the combinatorial objectand the stochastic cellular automaton.The first application of RDA - related to Ito equa-tion - is already studied in our parallel papers [15, 16].The model serves as a fully controlled stochastic ”phe-nomenon” and an applicability of the reconstruction ofthe Ito equation from generated time series was tested.Due to its simplicity, it allows to derive exactly the suit-able Ito equation, and analytical results were comparedwith histogram method. The obtained results are part ofbroader studies of the privilege concept and its role forappearance of inverse-power distributions [25, 26].The model posses also some properties which may beused in seeking of applications to natural phenomena andchecking their adequacy. For example, in the contextof earthquakes, the discrete nature of the model makespossible an adjustment of positions and sizes of cells tothe geological structure of a fault. Bending and shrinkingof the 1-D grid have no influence on equations; the onlyimportant feature is an order of cells. In this context,relations describing average size of cluster < i > andavalanche < w > may be interpreted as relations givingtheir dependence on the density of energy on the wedge.The model establishes one-to-one correspondence be-tween distribution of avalanches (i.e. observable quan-tity) and respective parameter µ i ν responsible for trigger-ing avalanches, as it was illustrated in two cases consid-ered above. We emphasize that the relation works in bothdirections. Fixing rebound parameters leads to variousdistributions of avalanches, but also, for any distribution n i , one can easily find unique values of rebound parame-ters µ i ν by simple use of the set of equations (5)-(7). Thus,within the proposed model, one can infer about trigger-ing properties from ”observed” statistics of avalanches.For an earthquake interpretation, it would give an in-sight into properties of microscopic mechanism of releas-ing energy on a fault on the basis of already collecteddata.Last but not least, we propose few generalizations pos-sibly leading towards more realistic extensions of themodel. Each element of the presented automaton - in-cluding the incidence rule - can be subjected to variousmodifications. To be more specific, one can consider dif-ferent geometry of the array (for example, a tree shaped like Bethe lattice or any in bigger dimension), differentcapacities of cells and different kinds of balls distributedto the system. Also there are many kinds of dependenceof energy release threshold on other parameters (on spaceposition, on states of cells in a neighbourhood etc.). Weleave these topics for further investigations. Appendix A: Derivation of balance equations forempty clusters.
In analogy to the set of equations (5)-(7) for clusters,one can deduce the following set of equations for theempty clusters.
Losses.
There are two following possibilities. (a)
Occupation of any cell belonging to the empty k -cluster. The probability is just ∼ ν kn k N . (b)
Provoking an avalanche. An empty cluster of thelength k is lost when one of the two adjacent clusters isknocked out and forms an avalanche. · · · | | i z }| { • | • | • | k z }| { | | | | | | • | · · · The probability is ∼ n k n X i ≥ µ i in i N .
It comes from the probability of knocking out of i -clustertimes the probability that on the end of the cluster isempty k -cluster, times two ends, summed for all possiblevalues of i . Gains.
Again, there are two possibilities. (a)
Shortening an empty cluster. An empty cluster of alength k can be obtained from a bigger one of the length i ≥ k + 1 when an incoming ball hits its ( k + 1)th emptycell · · · | • | i z }| { | | | {z } k | × | | | | | | • | · · · or symmetrically, its ( i − k )th cell. Hence the relevantgain term is of the form ∼ ν X i ≥ k +1 n i N . (b)
Provoking an avalanche. If an incoming ball knockout the cluster of the length j ∈ { , , . . . , k − } sur-rounded by two empty clusters of the respective lengths l ∈ { , , . . . , k − − j } and ( k − j − l ), · · · | • | l z }| { | | j z }| { • | • | • | k − l − j z }| { | | | | {z } k | • | · · · then the empty cluster of the length k is created. There-fore, for k ≥ ∼ k − X j =1 k − − j X l =1 µ j jn j N n l n n k − j − l n . This is the probability of knocking out of the j -clustertimes probability the l -empty cluster is on the left timesprobability the empty ( k − l − j )-cluster is on the right,for all possible values of l . Balance equations.
The set of equations for the av-erage values in the quasi equilibrium is of the form n k k + 2 n X i ≥ µ i ν in i = 2 X i ≥ k +1 n i , k = 1 , , (A1) n k k + 2 n X i ≥ µ i ν in i = 2 X i ≥ k +1 n i ++ k − X j =1 k − − j X l =1 µ j ν jn j N n l n n k − j − l n , k ≥ . (A2)Equation (A1) may be written as n k = ( k + 2) + n P i ≥ µ i ν in i (cid:16) n − P k − i =1 n i (cid:17) for k = 1 , , and equation (A2), valid for k ≥ n k = ( k + 2) + n P i ≥ µ i ν in i n − P k − i =1 n i ) + P k − j =1 P k − − jl =1 µ j ν jn j N n l n n k − j − l n . Thus, on right hand sides there are terms n j with j < k only. Balance of n out of empty clusters. The set ofequation for empty clusters can be treated by methodswe use in Subsection III E in order to derive relations forrespective moments.Below we consider zero moment for n i , which leadsalso to equation (4). This is implied by the fact thatthe number of empty clusters is equal to the number ofclusters (see equation (1)). To perform sum of (A1) and(A2) for all values of k notice the following identities: X k ≥ n k k + 2 n X i ≥ µ i ν in i = (1 − ρ ) N + 2 X i ≥ µ i ν in i , X k ≥ X i ≥ k +1 n i = X k ≥ ( i − n i = (1 − ρ ) N − n, and X k ≥ k − X j =1 k − − j X l =1 µ j ν jn j N n l n n k − j − l n = 1 N X j ≥ µ j ν jn j . In the last identity we assumed that it is allowed tochange the order of summation with respect to indices j and k . Thus, a sum of (A1) and (A2) for all k leadsdirectly to (1 − ρ ) N − n = X i ≥ µ i ν n i i. It is the balance equation (4) already obtained.
Acknowledgement
This work was partially supported by the project IN-TAS 05-1000008-7889. [1] C. Tang and P. Bak, J. Stat. Phys.
51 No 5/6 , 797(1988).[2] D. Dhar and R. Ramaswamy, Phys. Rev. Lett. , 1659(1989).[3] A. Vespignani and S. Zapperi, Phys. Rev. E , 6345(1998).[4] F.C. Alcaraz and V. Rittenberg, Phys. Rev. E , 041126(2008).[5] D. Dhar, Physica A , 4 (1999).[6] B. Gaveau and L. Schulman, J. Phys. A: Math. Gen. ,L475 (1991).[7] S.C. Lee, N.Y. Liang, and W.J. Tzeng, Phys. Rev. Lett. , 1479 (1991).[8] D. Takahashi and J. Satsuma, J. Phys. Soc. Jpn.
59 No10 , 3514 (1990).[9] T. Tokihiro, D. Takahashi, J. Matsukidaira, and J. Sat-suma, Phys. Rev. Lett. , 3247 (1996).[10] T. Tokihiro, in Discrete Integrable Systems , edited by B. Grammaticos, T. Kosmann-Schwarzbach, andT. Tamizhmani (Springer, 2004), vol. 644 of
LectureNotes in Physics , pp. 383–424.[11] M. Bia lecki and A. Doliwa, Commun. Math. Phys. ,157 (2005).[12] M. Bialecki, RIMS Kokyuroku , 154 (2009).[13] A. Doliwa, M. Bia lecki, and P. Klimczewski, J. Phys. A:Math. Gen. , 4827 (2003).[14] M. Bia lecki (2011), arXiv:1102.0437 [math-ph] .[15] Z. Czechowski and M. Bia lecki (2010), arXiv:1012.5902[nlin.CG] .[16] Z. Czechowski and M. Bia lecki (2010), arXiv:1101.0098[nlin.CG] .[17] D. Sornette, Critical Phenomena in Natural Sciences,2nd edition (Springer, 2006).[18] P. Bhattacharyya and B. Chkrabarti,
Modelling Criticaland Catastrophic Phenomena in Geoscience , vol. 705 of
Lecture Notes in Physics (Springer, 2006). [19] A. Tejedor, S. Ambroj, J. B. Gomez, and A. F. Pacheco,J. Phys. A: Math. Theor , 375102 (16pp) (2008).[20] M. Vazquez-Prada, A. Gonzalez, J. B. Gomez, andA. F. Pacheco, Nonlinear Processes in Geophysics , 513(2002).[21] A. Gonzalez, M. Vazquez-Prada, J. B. Gomez, and A. F.Pacheco, Tectonophysics , 319 (2006).[22] N. Sloane and S. Plouffe, The Encyclopedia of IntegerSequences (Academic Press, 1995). [23] N. Inui and M. Katori, J. Phys. A: Math. Gen. , 4347(1996).[24] D. Stauffer and A. Aharony, Introduction to PercolationTheory (Taylor & Francis, 1992).[25] Z. Czechowski, Geophys. J. Int. , 754 (2003).[26] Z. Czechowski and A. Rozmarynowska, Physica A387