Self-Organized Criticality on Self-Similar Lattice: Exponential Time Distribution between Extremes ? Dayana Mukhametshina1 , Alexsander Shapoval1,2[0000−0001−5340−1930] , and Mikhail Shnirman2 1 National Research University Higher School of Economics, 20 Myasnitskaya Ulitsa, 101000 Moscow, Russia hse@hse.ru https://www.hse.ru/en/ 2 Institute of Earthquake Prediction Theory and Mathematical Geophysics, Profsoyuznaya 84/32, 117997 Moscow, Russia Abstract. In 1987, Bak, Tang, and Wiesenfeld introduced a mechanism (hereafter, the BTW mechanism) that underlies self-organized critical systems. Extreme events generated by the BTW mechanism are be- lieved to exhibit an unpredictable occurrence. In spite of this general opinion, the largest events in the original BTW model are efficiently predictable by algorithms that exploit information that is hidden in ap- plications. Intending to relate the predictability of self-organized critical systems with the level of its asymmetry, we examine the inter-event dis- tribution of extreme avalanches generated by the BTW mechanism on symmetrical and asymmetrical self-similar lattices. Initially, we claim that the main part of the size-frequency relationship is power-law in- dependent of the asymmetry, but the asymmetry reduces the range of scale-free avalanches in the domain of small avalanches. Further, we turn to extremes and claim that they are located on the downward bend of the distribution of the avalanches over their sizes. Finally, we compare the probability distribution of waiting time between two successive extremes with the exponential distribution. The latter gives the reference point of the complete unpredictability naturally measured in terms of the sum of two rates related to type I and II statistical errors: the rate of the unpredicted avalanches and the alarm time rate. We posit that the devi- ations of the observed probability distribution from the exponential one do not affect the unpredictability of extremes drawn from the waiting time between them. Keywords: sandpile · power-law · waiting time distribution. 1 National Research University Higher School of Economics. * Copyright c 2019 for this paper by its authors. Use permitted under Creative Com- mons License Attribution 4.0 International (CC BY 4.0) 2 D. Mukhametshina , A. Shapoval, M. Shnirman 1 Introduction Extreme events constitute an essential phenomenon of complex systems. Their economic and social consequences are difficult to overestimate. Nevertheless, yet several decades ago, scholars primary focused on regular behavior of complex systems. Only recently the scenarios of extremes become better understood, but the prediction of extremes still remains a challenge for researchers, espe- cially if the system exhibits so called self-organized criticality associated with a power-law size-frequency distribution of “normal” events [12]. Bak et al. [2] introduced a simple mechanism (hereafter, the BTW mechanism) that gener- ates self-organized criticality. Two multi-scale processes: slow loading and quick stress release, which balance the stress on average, characterize the BTW mecha- nism. The stress release is modeled as an avalanche that redistributes the stress over the underlying system and carries the energy out at the boundary. The system attains a critical state observed through the power-law distribution of the avalanches over their size [9]. Typically, the existence of the power-law size- frequency distributions signals that the prediction of extremes are hardly pos- sible [1]. Under the BTW mechanism, a slight increase of loading affects the system alike its local redistribution. However, the BTW mechanism generates an out-of-equilibrium system whose observed variables oscillate around average values. Such oscillations, in general, could open a door for an effective prediction of extremes that topple the system from the super- to sub-critical state, if the oscillation exhibits a certain quasi-periodicity. Bak et al. [2] gave the first simple explanation of numerous critical phe- nomena reported by that time and hugely extended by nowadays [36]. The (truncated) power-laws describe the probability distributions of earthquakes [13], landslides [34], solar flares [7], sizes of large cities [10,15], city fires [33], and fi- nancial crashes [17]. Surprisingly, the power-law exponent observed in the BTW model is hardly tuned with the transformation of the BTW mechanism [3]. The random version of the model proposed by [16] results in a similar but distinct exponent [6,4]. However other numerous modifications leave the models inside the two classes of universality determined by these two exponents [29]. The BTW mechanism itself primary exhibits a concept of multiscale processes that exhibits self-organized criticality, but does not describe details of observed systems. The impossibility to adjust the exponent of the power-law and leave the (very) limited set of the universality classes reduces the range of the direct applicability of the BTW model and its adjacent generalizations. The exponent is tuned, if the BTW mechanism is realized on fractals [8,5], site-percolation lattice [20], and self-similar lattice [31]. Typically, a downward bend follows the scale-free range of the size-frequency relationship at the right part of the graph. In the BTW model, the transition is explained by the finite-size effect. The occurrence of the large events can exhibit time-clustering [28]. In this case, the largest events follow a foreshock activity or trigger aftershocks, or demonstrate the both phenomema [32]. Such time- clustering underlies the algorithms that predict the largest events. These events are characterized by specific waiting time probability distribution [22]. The BTW Self-Organized Criticality on Self-Similar Lattice 3 model is dissipative on the boundary, and large avalanches drive energy out of the system. Therefore, a large loading seems to be required for a consecutive extreme avalanche to occur. The existence of the quiescent episodes in the dynamics generates a certain quasi-periodicity of the largest avalanches. These arguments are confirmed for the BTW model [30] and the laboratory experiment imitating the model [25]. As far as the waiting time probability distribution deviates from the exponential one, the waiting time itself serves as a precursor of the main event in the models and their applications [26,11,27,35]. However, the extension of the set of the largest events to smaller ones destroys the efficiency of this precursor, tending the waiting time distribution to be exponential [23,11]. The avalanches that are located on the power-law part of the size-frequency relation graph are considered as unpredictable [24,1]. [19] constructed the framework based on the fraction of the unpredicted events and the alarm time rate in order to quantify the prediction algorithms. The fall of the prediction efficiency with the size of the forecasted events remains a general feature of the system [14] The goal of this paper is to extend the applicability of the BTW-like models to real-life processes by testing the existence of new universality classes gener- ated by asymmetry of the underlying system and inspecting the predictability of extremes within the universality classes. Even if the universality classes are insensitive to the asymmetry of the system, the asymmetry itself may affect the scenario of extremes. We examine the inter-event distribution of extreme avalanches generated by the BTW mechanism on the self-similar lattice. Ini- tially, we address the question of how the asymmetry of the lattice affects the size-frequency relationship of the avalanches, finding the scale-free part of the relationship and establishing that its exponent is insensitive to the asymmetry. Further, we turn to extremes and claim that they are located on the downward bend of the distribution of the avalanches over their sizes. Finally, we compare the probability distribution of waiting time between two successive extremes with the exponential distribution. The latter gives the reference point of the complete unpredictability naturally measured in terms of type I and II errors. We posit that the deviations of the observed probability distribution from the exponential one do not affect the unpredictability of extremes drawn from the waiting time between them. The rest of the paper is organized in the following way. We define the model in Section 2, discussing the size-frequency relationship in details. The waiting time probability distribution is studied in Section 3. The last section concludes. 2 Model 2.1 Self-Similar Lattice: In order to define a model we introduce the following notation. Consider d × d lattice. Let a pattern be an arbitrary partition of the lattice cells into two (marked and unmarked) sets. We note that it is enough to specify the marked set to determine the partition. The case d = 3 is considered. We deal with the two patterns that consist of four cells illustrated by Figure 1. The first pattern is 4 D. Mukhametshina , A. Shapoval, M. Shnirman symmetrical. It consists of four corner cells, colored in grey in Figure 1a. The second pattern consists of three cells located along a side, whereas the forth cell adjoins the corner one, Figure 1b. This pattern can be treated as an example of the largest asymmetry. As the pattern repeats the move of the chess knight, we further refer to it as a knight Fig. 1. Symmetrical (a) and chess knight (b) patterns We define a self-similar lattice that consists of cells of different sizes. The smallest cells determine the unit of measurements. The (linear) length of the lattice is L = 3n , where n is the depth of self-similarity realized in computations. In theoretical constructions, n tends to +∞. The recursive algorithm, constructing the self-similar lattice, at each step takes all cells that does not belong to the pattern, divide them into 9 equal squares, and specify the pattern among these new, smaller squares. The patterns are pre-defined; they are identical at each step. More precisely, let a pattern be fixed. We denote C0,L a lattice with L cells on the side. At every step r = 0, 1, . . . , n − 2 of the recursion, the algorithm repeats the following sub-steps: 1. Split each cell c ∈ Cr,L into 9 equal squares with 3n−r−1 cells on the side. Four out of nine cells form the pattern Pr+1 (c). 2. Ĉr+1,L denotes the set of all cells appeared at sub-step 1 that form the patterns Pr+1 (c), c ∈ Cr,L . 3. Cr+1,L denotes the set of all cells appeared at sub-step 1 that does not belong to the patterns Pr+1 (c), c ∈ Cr,L . Note that step 1 is not applied to the cells from Ĉr+1,L . The final step r = n − 1 differs from the previous ones: 1. Split Cn−1,L into 9 equal squares with 1 cell on the side. 2. Cells corresponding to the pattern will be denoted as Ĉn,L . Self-Organized Criticality on Self-Similar Lattice 5 The set ∪nr=1 Ĉr,L is further referred to as self-similar lattice. Figure 2 illustrates the self-similar lattices obtained with n = 3 and d = 3 Fig. 2. Self similar lattice with (a) symmetrical and (b) knight patterns; n = 3, d = 3. The number of neighbors is written inside each cell. 2.2 Dynamics The cells, i. e., c ∈ Ĉr,L , where r = 1, . . . , n, are assumed to be enumerated in some way c1 , c2 , . . .. We use the following notation. – hi is a number associated with the cell ci . Following tradition [2], the quantity hi is called a grain of sand, whereas the model is referred to as sandpile. – N (i) is the set of the cells that are adjacent to the cell ci by a common edge; these cells are called neighbors. |N (i)| denotes the number of the neighbors illustrated in Figure 2. – If the cell ci is located inside the lattice C0,L (i. e., neither edge belong to the lattice boundary) then we put Hi = |N (i)|. Otherwise, Hi is the sum of |N (i)| and the number of the ci ’s edges located on the lattice boundary. The cell ci is called stable if hi < Hi . At the beginning all cells are stable. For the sake of simplicity we put hi = 0 for all i. Our computer simulations give evidence that what follows is insensitive to the initial conditions. The model dynamics consists of two stages: accumulation and avalanche. Accumulation. A grain falls on a randomly chosen cell: 1. A cell c of the self-similar lattice is chosen at random with the probability being proportional to the cell’s area. If c ∈ Ĉr,L , r = 1, . . . , n, then the probability is p = 32(n−r) /32n = 3−2r . 2. Let i be the number of the just chosen cell c, i. e., c = ci . Then a grain of sand “falls” on ci : hi → hi + 1. 6 D. Mukhametshina , A. Shapoval, M. Shnirman 3. We perform a stability check. If hi < Hi , we repeat accumulation process. Otherwise, the avalanche starts. Avalanche. The unstable cell ci topples transferring grains to the neighbors equally: hi → hi − Hi (1) hj → hj + 1, ∀j ∈ N (i) (2) The transfer (1) and (2) conserves the total amount of grains in the lattice when an inner cell topples. If at least one neighbor becomes unstable, the transfers continue. In other words, while there are unstable cells ci , the rules (1) and (2) are applied. Each avalanche is finite because the transfers at the boundary are dissipative. We call the number of sand grains that are displaced during the avalanche the size of the avalanche. 2.3 Power-Law Size-Frequency Relationship First, we investigate the probability distribution of the avalanches over their sizes. Let N (s) be the number of avalanches that have the size s. Then we define 1 Q(s) = (N (s) + N (s − 1) + N (s − 2) + N (s − 3)) 4 We substitute Q(s) for N (s) to stabilize the graph in the domain of small sizes. Figure 3 displays the graphs of Q(s) for the symmetrical (a) and knight (b) patterns. The power-law part of the graphs turns to a downward bend at the right. Fig. 3. Symmetrical (a) and knight (b) patterns, n = 5, 6, 7, d = 3 Self-Organized Criticality on Self-Similar Lattice 7 The power-law part and the structure of the downward bend can be studied in more details when the exponential binning is applied. Let f (s) be the fraction of the avalanches in the catalogue that have the size s; X F (s) = f (σ), (3) σ∈[s/∆s,s∆s) where ∆s is a parameter. We stress that if f (s) follows a power-law function, then F (s) also does, but the exponent of the power-law is increased by 1. Indeed, the summation in equation (3) serves as the integral that transforms 1/sβ into 1/sβ−1 . For instance, if β = 1, the integral over the exponential bin is Z s∆s s∆s s−1 = ln s s/∆s = 2 ln ∆s = 2s0 ln ∆s. s/∆s Fig. 4. Symmetrical (a) and knight (b) patterns, n = 5, 6, 7, d = 3, ∆s = 1.6 Figures 4 a and b exhibit the size-frequency relation for the both patterns. The both graphs have a power-law part followed by a bend down that might be characterized by a steeper line in the log-log scale. The asymmetry affects the left part of the graph that corresponds to small avalanches. On the asymmetrical lattice, a significant share of the avalanches reaches the boundary and stops almost immediately. 3 Waiting Time Distribution In this section, we investigate the waiting time probability distribution and relate it to the prediction of rare events. Rare events are defined parametrically. We will 8 D. Mukhametshina , A. Shapoval, M. Shnirman call an avalanche the rare event if its size is bigger than some S ∗ . The inter-event probability distribution rather accurately follows the exponential function while S ∗ belongs to the power-law range of sizes (we do not support this statement by a graph). If S ∗ is located on the right part of the downward bend, this probability distribution deviates from the exponential function; see Figure 5, where the complement cumulative distribution function, 1 − cdf is displayed. Under symmetrical pattern, the deviation is in the domain of a large waiting time: a large time gap between successive rare events occurs more frequency then the exponential random probability distribution predicts. In the case of the knight pattern, the observed waiting time distribution is convex in the linear-log scale, but the convexity is small. Fig. 5. Symmetrical pattern, S ∗ = 32000 (a) and knight pattern, S ∗ = 20000 (b) , n=5 We introduce a simple algorithm predicting the next rare event based on the record of the previous one. If a rare event occurs at the time moment t0 , then an alarm is raised for the interval (t0 , t0 + T ], where the parameter T > 0 will be adjusted later. The alarm means that the algorithm predicts the occurrence of the consecutive rare event during time window (t0 , t0 + T ]. The unit of time is associated with the grain falling. It falls one grain per the time unit. Each alarm continues either up to the occurrence of a new rare event (and then a new alarm is raised) or T units of time. The rare events that occurs when the alarm is raised are predicted. The other rare events are unpredicted. Let η be the fraction of the unpredicted rare events. The first rare event is unpredicted by the construction of the algorithm. Therefore, it is ignored when η is computed. Further, let τ be the alarm rate (the sum of the intervals with the raised alarm divided by length of the whole Self-Organized Criticality on Self-Similar Lattice 9 time interval) and ε = η + τ . The quantities η and τ are related to type I and II statistical errors [18]. Their sum ε describes the efficiency of the prediction algorithm. The equation ε ≈ 1 characterizes unpredictability. If the underlying probability distribution of the waiting time is exponential then our algorithm results in ε = 1 (this is the result of an elementary mathematics skipped here). The both errors as functions of the parameter T are displayed in Figure 6 whereas their sum is shown in 7. In our computer simulations, S ∗ is chosen as 32000 (symmetrical pattern) and 20000 (knight pattern), that leads to the amount of rare events equalled to 193 and 139 respectively. Fig. 6. Symmetrical pattern, S ∗ = 32000 (a) and knight pattern, S ∗ = 20000 (b) , n=5 Fig. 7. Symmetrical pattern, S ∗ = 32000 (a) and knight pattern, S ∗ = 20000 (b) , n=5 10 D. Mukhametshina , A. Shapoval, M. Shnirman The quantity ε oscillates around 1. This suggests that the waiting time proba- bility distribution deviates from the power one insignificantly. 4 Conclusion We constructed the BTW mechanism on the self-similar lattice. The existence of the asymmetry in the geometry of the self-similar lattice does not affect the self- organized criticality of the underlying system. The size-frequency relationship follows the power-law and the power-law exponents are identical. Their value was earlier predicted by [31]. The power-law part of the size-frequency relation graphs is transformed to a downward bend located at the right. There are traces of another power-law in this bend whose exponent is (very) roughly estimated as 3. The best approximation to the downward bend deserves a separate study. We found that the waiting time distribution between rare events deviates from the exponential one, if the low border of these events is sufficiently large. However we did not detect a power-law inter-event distribution as [22] identified for the classical BTW sandpile. In the sandpile on the self-similar lattice studied here, the prediction of the rare avalanches drawn from their waiting time distri- bution fails completely. We construct the algorithm that expects the occurrence of a new extreme within T time units after the previous extreme. This algorithm is ineffective for any value of the parameter T . This result does not give evidence against the prediction of extremes in the sandpile model on the self-similar lattice in general. We believe that the pre- dictability of large avalanches found by [28] is translated onto the case of the BTW-mechanism on the self-similar lattice. At least, it is worth trying to con- struct a prediction algorithm that switches on the alarm when the surplus of sand grains on the lattice is observed. This algorithm could be effective since large avalanches are expected upon the system becomes supercritical. The BTW-like models considered in the paper can be used when study- ing earthquake formation processes, as the exponent in the power law is tuned through an appropriate choice of the self-similarity determinant. Narkunskaya & Shnirman [21] scrutinized long-term dynamics of the distribution of the earth- quakes over their magnitude and constructed a precursor of strong earthquakes drawn from an upward bend of this distribution. The existence of a similar precursor in the model would further justify its applicability to the earthquake formation process. Then the prediction algorithms adjusted on our artificial sys- tem can be further tested when predicting strong earthquakes. References 1. Bak, P., Paczuski, M.: Complexity, contingency, and criticality. Proceedings of the National Academy of Sciences 92(15), 6689–6696 (1995) 2. Bak, T.J., Tang, C., Wiesenfeld, K.: Self-organized criticality: an explanation of 1/f noise. Phys. Rev. Lett. 59, 381–383 (1987). https://doi.org/10.1103/PhysRevLett.59.381 Self-Organized Criticality on Self-Similar Lattice 11 3. Ben-Hur, A., Biham, O.: Universality in sandpile models. Physical Review E 53(2), R1317 (1996) 4. Biham, O., Milshtein, E., Malcai, O.: Evidence for universality within the classes of deterministic and stochastic sandpile models. Physical Review E 63(6), 061309 (2001) 5. Chen, J.P., Kudler-Flam, J.: Laplacian growth & sandpiles on the sierpinski gas- ket: limit shape universality and exact solutions. arXiv preprint arXiv:1807.08748 (2018) 6. Chessa, A., Stanley, E.H., Vespignani, A., Zapperi, S.: Universality in Sandpiles. Phys. Rev. E 59, R12–R15 (1999) 7. Crosby, N.B., Aschwanden, M.J., Dennis, B.R.: Frequency distributions and cor- relations of solar x-ray flare parameters. Solar Physics 143, 275–299 (1993) 8. Daerden, F., Priezzhev, V.B., Vanderzande, C.: Waves in the sandpile model on fractal lattices. Physica A: Statistical Mechanics and its Applications 292(1-4), 43–54 (2001) 9. Dhar, D.: Self-organized critical state of sandpile automaton models. Physical Re- view Letters 64(14), 1613 (1990) 10. Gabaix, X.: Zipf’s Law for Cities: An Explanation. Quarterly Journal of Economics 114, 739–767 (1999) 11. Geist, E.L., Parsons, T.: Distribution of tsunami interevent times. Geophysical Research Letters 35(2) (2008) 12. Ghil, M., Yiou, P., Hallegatte, S., Malamud, B., Naveau, P., Soloviev, A., Friederichs, P., Keilis-Borok, V., Kondrashov, D., Kossobokov, V., et al.: Extreme events: dynamics, statistics and prediction. Nonlinear Processes in Geophysics 18(3), 295–350 (2011) 13. Gutenberg, B., Richter, R.F.: Frequency of earthquakes in California. Bulletin of the Seismological Society of America 34, 185–188 (1944) 14. Hallerberg, S., Kantz, H.: Influence of the event magnitude on the predictability of an extreme event. Physical Review E 77(1), 011108 (2008) 15. Levy, M.: Gibrat’s Law for (All) Cities: Comment. American Economic Review 99, 1672–1675 (2009). https://doi.org/10.1257/aer.99.4.1672 16. Manna, S.: Two-state model of self-organized criticality. Journal of Physics A: Mathematical and General 24(7), L363 (1991) 17. Mantegna, R.N., Stanley, H.E.: An Introduction to Econophysics: Correlations and Complexity in Finance. Cambridge University Press, Cambrdige, England (2000) 18. Molchan, G.M.: Structure of optimal strategies in earthquake prediction. Tectono- physics 193(4), 267–276 (1991) 19. Molchan, G.M.: Earthquake prediction as a decision-making problem. Pure and Applied Geophysics 149(1), 233–247 (1997) 20. Najafi, M.: Bak–tang–wiesenfeld model on the square site-percolation lattice. Jour- nal of Physics A: Mathematical and Theoretical 49(33), 335003 (2016) 21. Narkunskaya, G., Shnirman, M.: Hierarchical model of defect development and seismicity. Physics of the Earth and Planetary Interiors 61, 29–35 (1990) 22. Paczuski, M., Boettcher, S., Baiesi, M.: Interoccurrence times in the bak-tang- wiesenfeld sandpile model: A comparison with the observed statistics of solar flares. Physical review letters 95(18), 181102 (2005) 23. Paczuski, M., Maslov, S., Bak, P.: Avalanche dynamics in evolution, growth, and depinning models. Physical Review E 53(1), 414 (1996) 24. Pepke, S., Carlson, J.: Predictability of self-organizing systems. Physical Review E 50(1), 236 (1994) 12 D. Mukhametshina , A. Shapoval, M. Shnirman 25. Rosendahl, J., Vekić, M., Rutledge, J.: Predictability of large avalanches on a sandpile. Physical review letters 73(4), 537 (1994) 26. Saichev, A., Sornette, D.: Theory of earthquake recurrence times. Journal of Geo- physical Research: Solid Earth 112(B4) (2007) 27. Shapoval, A.: Prediction problem for target events based on the inter-event waiting time. Physica A: Statistical Mechanics and its Applications 389(22), 5145–5154 (2010) 28. Shapoval, A.B., Shnirman, M.G.: Strong events in the sand-pile model. Interna- tional Journal of Modern Physics C 15(02), 279–288 (2004) 29. Shapoval, A.B., Shnirman, M.G.: Crossover phenomenon and universality: From random walk to deterministic sand-piles through random sand-piles. International Journal of Modern Physics C 16(12), 1893–1907 (2005) 30. Shapoval, A.B., Shnirman, M.G.: Scenarios of large events in the sandpile model. Selected Papers From Volumes 33 and 34 of Vychislitel’naya Seysmologiya 8, 179– 183 (2008) 31. Shapoval, A.B., Shnirman, M.G.: The btw mechanism on a self-similar image of a square: A path to unexpected exponents. Physica A: Statistical Mechanics and its Applications 391(1-2), 15–20 (2012) 32. Shnirman, M., Shapoval, A.: Variable predictability in deterministic dissipative sandpile. Nonlinear Processes in Geophysics 17(1), 85–91 (2010) 33. Song, W.G., Zhang, H.P., Chen, R., Fan, W.C.: Power law distribution of city fires. Fire Safety J. 38, 453–465 (2002). https://doi.org/10.1016/j.physa.2011.08.020 34. Stark, C.P., Hovius, N.: The characterization of landslide size distributions. Geo- phys. Res. Lett. 28, 1091–1094 (2001) 35. Talbi, A., Bellalem, F., Mobarki, M.: Turkey and adjacent area seismicity forecasts from earthquake inter-event time mean ratio statistics. Journal of Seismology pp. 1–13 (2019) 36. Watkins, N.W., Pruessner, G., Chapman, S.C., Crosby, N.B., Jensen, H.J.: 25 years of self-organized criticality: concepts and controversies. Space Science Re- views 198(1-4), 3–44 (2016)