Record-Breaking Avalanches in Driven Threshold Systems

Record-breaking avalanches generated by the dynamics of several driven nonlinear threshold models are studied. Such systems are characterized by intermittent behavior, where slow buildup of energy is punctuated by an abrupt release of energy through avalanche events which usually follow scale invariant statistics. From the simulations of these systems it is possible to extract sequences of record-breaking avalanches, where each subsequent record-breaking event is larger in magnitude than all previous events. In the present work, several cellular automata are analyzed among them the sandpile model, Manna model, Olami-Feder-Christensen (OFC) model, and the forest-fire model to investigate the record-breaking statistics of model avalanches which exhibit temporal and spatial correlations. Several statistical measures of record-breaking events are derived analytically and confirmed through numerical simulations. The statistics of record-breaking avalanches for the four models are compared to that of record-breaking events extracted from the sequences of independent identically distributed (\emph{i.i.d.}) random variables. It is found that the statistics of record-breaking avalanches for the above cellular automata exhibit behavior different from that observed for \emph{i.i.d.} random variables which in turn can be used to characterize complex spatio-temporal dynamics. The most pronounced deviations are observed in the case of the OFC model with a strong dependence on the conservation parameter of the model. This indicates that avalanches in the OFC model are not independent and exhibit spatio-temporal correlations.


I. INTRODUCTION
A record-breaking event is an event that is larger (smaller) than all previous events in the observational sequence of some physical measurements. Sequences of record-breaking events can be extracted for various natural phenomena where records of the variability of different measurements exist. The best known example is weather, where daily observations of temperatures are documented and new record-breaking temperatures are often reported. Monitoring such record-breaking temperatures can be used to infer possible trends in the weather variability [1][2][3][4][5]. For example, the excess of maximum record-breaking temperatures over minimum record-breaking ones is taken to be direct evidence of global warming. Other natural phenomena for which record-breaking events play an important role include earthquakes [6][7][8][9][10], floods [11], forest-fires, volcanic eruptions, and solar events [12], to name only a few.
Record-breaking events or records form a subsequence of all recorded measurements and are only distinguished based on their sizes or magnitudes. Since a record-breaking event is defined as the largest event up to its time of occurrence, a subsequent event becomes the next record-breaking event only if it exceeds the size of the previous record and all events in between are smaller than both records [1,4,[13][14][15]. The mathematical framework to analyze record-breaking events is related to extreme value statistics and in most cases assumes events which are independent and identically distributed (i.i.d.) [13][14][15][16]. In this case several statistical properties of records are independent of their underlying distribution and yield to analytical treatment [17]. This standard treatment of records assumes that the underlying events occur equidistant in time. As shown in Ref. [18], one can generalize the analysis of record-breaking events to situations in which the occurrence of events grows stochastically in time. It was found that such a temporal modification preserves several fundamental properties of the original model, i.e. independence from the underlying distribution of the i.i.d. events.
In contrast the situation can be significantly different if one considers more complicated stochastic processes. In recent years some progress has been made in analyzing records and record-breaking events in non-i.i.d. random sequences including sequences with time-varying underlining distributions as well as sequences with strongly correlated events. For example, a theoretical and observational analysis was performed for daily record temperatures in Philadelphia to investigate the effects of temperature trends and correlations [3]. In that work the authors found that the small rate of increase in the mean temperatures does not significantly influence the statistics of record-breaking events over the time span of the observations. The authors also considered the role of temporal correlations between events by generating synthetic time series with correlations decaying as a power-law. It was reported that the frequency and magnitude of the record-breaking events was almost identical to those obtained using i.i.d. random variables and were not sensitive to the power-law exponent of the correlation function. Events drawn from independent random variables but with progressively broadening or sharpening distributions were considered in [5]. It was shown that for a power-law growth of the width of the distribution three regimes can be identified if one considers the growth of the mean number of records. These regimes correspond to the asymptotic behavior of the tail of the distribution according to Gumbel, Fréchet, and Weibull classes often studied in the context of extreme value statistics [19]. To further investigate the effects of correlations in time series, record-breaking events in sequences generated by random walks and Lévy flights were also considered [20]. It was found that the statistics of records are independent of the distribution of jumps. The average number of records and its standard deviation grow as a square root of the number of time steps which implies that fluctuations remain large and of the same order as the mean. For related stationary and non-stationary stochastic processes with long-range memory, the situation is more complicated but memory typically leads to significant deviations from the i.i.d. behavior [12,21].
Despite this progress for some non-i.i.d. processes, we are still far away from a general understanding of recordbreaking events and their properties in complex systems that are often observed in nature. To understand such systems and specifically the emergence of complex patterns and long-range spatio-temporal correlations in physical systems, various cellular automata have been studied in the past. Usually such models possess rather simple and intuitive rules of evolution which mimic real observed properties of dissipative complex systems. Of particular interest are driven nonlinear threshold systems exhibiting avalanche-like behavior with scale invariant statistics where effects of correlations and self-organization play a crucial role [22,23]. A plethora of models has been introduced to illustrate such behavior. Among them the sandpile model [24,25], Manna model [26], the Olami-Feder-Christensen model [27], and the forest-fire model [28,29] are the most prominent ones. They all exhibit scale-invariant behavior in the distribution of avalanches with respect to their size and duration when driven into a steady state. They also possess nontrivial spatio-temporal correlation patterns and serve as conceptual models to explain such diverse phenomena as earthquakes, forest-fires, landslides, etc.
In the present work, we investigate the influence of spatio-temporal correlations on records by studying the recordbreaking statistics of avalanches in the four cellular automata described above. For each of the models we study the record statistics both in magnitude and time domains by constructing various distributions and the corresponding averages. Our numerical simulations show that there are significant differences to what would be expected for i.i.d. random variables. The most pronounced deviations are observed in the case of the OFC model with a strong depen-dence on the conservation parameter of the model. In addition to numerical simulations and analysis of records of the four cellular automata, we also derive an exact analytical expression, Eq. (18), for the probability that the kth record is broken after m time steps in the case of i.i.d. random variables. This allows us a more rigorous comparison between the models and the i.i.d. case.
The paper is organized as follows. In the next section we introduce basic known facts concerning the statistics of record-breaking events extracted from sequences of i.i.d. random variables. We also present the analysis of recordbreaking events generated from the truncated power-law distribution which is relevant for our analysis of threshold cellular automata. In Section III we recall the definition of the four cellular automata which we have studied to infer the structure of their record-breaking avalanches. In Section IV we report on the results of computer simulations of those models. Section V concludes our analysis.

II. RECORD-BREAKING STATISTICS
Consider a sequence of measurements, {m(t i )}, of a particular outcome of a physical experiment or computer simulation of a model over a period of time, where t i are the times of measurements. In general, this represents a stochastic variable. Examples include temperature measurements, earthquake magnitudes, flood areas, forest-fire areas, etc. A record-breaking event x(t n ) up to time t n , is defined as the event with the largest value among all previous measurements, x(t n ) = max{m(t 1 ), m(t 2 ), ..., m(t n−1 )} [13]. Subsequent measurements can produce recordbreaking events if they exceed in value the current record-breaking event. Therefore, for any historical data or for ongoing observations the sequence of record-breaking events can be extracted and analyzed. In the analysis which follows, we use discrete time steps n = 1, 2, 3, ... which define the moments when measurements are done or time steps when random variables are generated. In this case we use a simplified notation for records: x(t n ) ≡ x(n). To enumerate the record-breaking events in a sequence we use a subscript k = 1, 2, 3, ... to indicate their record order. For example, x k indicates the value of the kth record-breaking event, therefore, for a given realization of observations one has the following sequence of record-breaking events: {x k } = x 1 , x 2 , x 3 , .... To specify in addition the time of the occurrence of the kth record we write: x k (n).
In the most simple case it is assumed that the records are extracted from the sequence of i.i.d. random variables drawn from a particular probability density function f (x). This density can be specified a priori from known distributions or it can be the outcome of the measurements of natural systems or the computer simulations of model systems. It can be bounded or unbounded depending on the system under consideration. For the bounded distributions the corresponding record-breaking events will be bounded as well. The probability that the records will not exceed x is given by where x min is the lower bound of the support over which the distribution is defined.
To quantify different aspects of record-breaking events several fundamental statistical measures can be estimated. The existing theory of record-breaking events is mostly developed for the i.i.d. random variables drawn from a particular distribution [14,15]. In the rest of the section we assume such kind of random variables.

A. The average number of record-breaking events
The average number of record-breaking events, N n , characterizes the occurrence of records in time. It can be estimated by averaging over ensemble realizations of sequences of random variables drawn from a particular distribution with the probability density function f (x). Consider a random variable N n which is the number of records that occurred up to a time step n. This number can be expressed as where I j is the indicator function given by , if a record occurs in the jth time step, 0, otherwise.
Assuming that the events are i.i.d. random variables the probability of the occurrence of the jth record-breaking event is simply P j = 1 j [14]. This represents the harmonic decrease in time of the probability of the occurrence of record-breaking events. Therefore, the mean value of I j can be computed using this probability: Using Eqs. (2) and (4) the average number of records, N n , can be evaluated where H n is a harmonic number, H n = n j=1 1 j , and γ ≈ 0.577215665... is the Euler-Mascheroni constant [30]. This result indicates that the average number of record-breaking events grows as N n ∼ ln n for large values of n and is independent of the underlying probability density function f (x). For record-breaking events generated from processes with memories or long range correlations the average number can deviate from Eq. (5).
The variance of the average number of record-breaking events N n is given by [5] Var For i.i.d. random variables the ratio of the variance, Eq. (6), to the mean, Eq. (5), tends to unity as n → ∞ and the distribution of the N n becomes Poisson distribution with the mean ln n. As a consequence the record times can be described by a log-Poisson process [5]. In this case the probability that k records have occurred up to a time step n is For large values of n the distribution for the number of records N n tends to the Gaussian distribution.

B. Distribution of magnitudes
The probability density function that the kth record-breaking event is equal to x has the form [3]: where F (x) is the distribution function for the random variable x given by Eq. (1). Eq. (8) represents a recursive formula to compute the distribution of magnitudes for the kth record-breaking event by knowing the distribution for the (k −1)st record. The p 1 is simply equal to the underlying distribution from which the random variables are drawn, . It is also possible to evaluate Eq. (8) for k = 2, 3, .... For example, by noticing that p 1 (x) = f (x) = dF dx , one has for k = 2 One can prove by induction that the general form of Eq. (8) for any arbitrarily k is Eq. (10) is valid for i.i.d. random variables drawn from any underlying distribution function F (x) [14].
The average magnitude, x k , of the kth record-breaking event can be estimated using Eq. (10) where S = [x min , x max ] is the support of the distribution function F (x). It is also possible to define the average magnitude x(n) of a record-breaking event at a time step n: where q(x, n) is the probability density function for the records at a time step n. So that the probability to find a record-breaking event to be between x and x + dx at time n is q(x, n)dx. This probability density function q(x, n) can be expressed in terms of F (x) as [4,14] q Noticing that f (x) = dF dx , Eq. (12) can be rewritten as follows It is important to note that the results, Eqs. (10)- (14), are valid only for i.i.d. random variables. In Section IV we consider simulations of the four cellular automata and compare the obtained results for record-breaking avalanches with the behavior given by Eqs. (10)- (14) in order to quantify the deviations from the i.i.d. case.

C. Distribution of interevent times
The interoccurrence or interevent times between two subsequent record-breaking events is an important measure which gives an estimate as for how long one has to wait for the next record to occur. It was shown that for i.i.d. random variables the distribution is independent of the underlying distribution of the record-breaking events and the non-normalized histogram follows a simple power law: G(m) = 1/m, where m is the time difference between kth and (k − 1)st record-breaking event, m = t k − t k−1 [4]. This distribution is constructed by incorporating all interevent times including those early in the sequence as well as those later in the sequence of records.
A more detailed structure of time intervals between record-breaking events is given by the distribution, w k (m), which is the probability that the kth record is broken after m time steps. The formal expression for this quantity in the case of i.i.d. random variables is [3,15] where F (x) is the underlying cumulative distribution function for the random variables from which records are drawn.
gives the probability that a new record, (k + 1)st, is set after m ≥ 1 time steps after the kth record of value x k . To obtain Eq. (15) one has to average this probability over all possible values of x. It is possible to evaluate Eq. (15) for the 1st record (k = 1) to be broken after m time steps, w 1 (m). This can be done explicitly by noticing that p 1 (x) = f (x) = dF dx . Substituting this into Eq. (15) and performing integration by parts one obtains Due to power-law nature of the distribution, Eq. (16), the average waiting time between the first record and the second record is infinite, m 1 = ∞ m=1 m w 1 (m) = ∞. Using Eq. (10) the Eq. (15) can be rewritten in terms of the cumulative distribution function F (x) as follows The integration can be carried out explicitly (Appendix A) with the result: For k = 2 and 3 we have where H m is a harmonic number and Ψ(m) is a digamma function [30]. These results show that the probability distribution w k (m) is independent of the underlying distribution F (x) as long as the random variables are i.i.d..
The derived equation for w k (m) can be used to compute the non-normalized histogram G(m) as follows where to evaluate the integral we exchanged the operation of summation and integration and used the known result for the expansion of the exponent into the infinite sum Using Eq. (22) the integral in Eq. (21) can be evaluated easily. The probability distribution, u k (n), for the time of occurrence of the kth record-breaking event at a time step n can also be analyzed. Using u k (n) one can compute the average time n k = ∞ n=1 n u k (n) of the occurrence of the kth record-breaking event. The first event is always a record-breaking event so n 1 = 1. The distribution for the second record-breaking event u 2 (n) coincides with the distribution w 1 (n) which defines the distribution of interevent times between the first record (k = 1) and the second record (k = 2).
The time of the occurrence of the kth record, n k , is a random variable and can be represented as a sum of two random variables, n k = n k−1 + m k , where n k−1 is the time of the (k − 1)st record and m k is the interevent time between records k and k − 1. Therefore, the density u k (n) can be obtained recursively from the discrete convolution of the two distributions u k−1 (n) and w k (m): The above defined distributions and averages can be used to describe statistical properties of record-breaking events generated by different models and processes and documented in observations of various natural systems. Based on their behavior it is possible to infer the structure of record-breaking events and underlying processes which generate them.

D. Power-law distributed random variables
Many natural phenomena with nonlinear threshold dynamics exhibit scale-invariant statistics in their distribution of output events (avalanches). These statistics can often be approximated by a power-law distribution. One prominent example is earthquakes where the frequency-magnitude statistics follows Gutenberg-Richter power-law scaling [31]. The probability density function for a power-law distribution can be written as where γ is an exponent. Usually observed physical processes exhibit power-law scaling only for several orders of magnitude and are truncated on both ends. This introduces finite-size effects into distributions and can lead to non-trivial scaling behavior. For this analysis we consider the power-law distribution with a bounded support [x min , x max ], where x min ≥ 0 and γ = 1. In this case the distribution is well defined and can be normalized. The probability density function f (x) and the distribution function F (x) are: For i.i.d. random variables the distribution function of the kth record-breaking event, p k (x), can be evaluated using Eqs. (8) and (26). Using Eq. (11) one can compute the mean value x k of the magnitude of the kth record-breaking event. The integration can be performed analytically for specific values of k but the expressions are rather cumbersome so one can evaluate them numerically. Similarly the average magnitude x(n) of the record-breaking events at a time step n can be evaluated numerically using Eqs. (14) and (26).

III. DRIVEN NONLINEAR THRESHOLD MODELS
Observations of many natural phenomena reveal that the behavior of these systems is governed by the threshold dynamics where a slow build up of energy is punctuated by the fast release of the accumulated energy through avalanches. The analysis of the frequency-magnitude statistics of these events shows that in many instances the distributions follow power-law statistics with anomalous exponents. These systems also exhibit non-trivial correlations in time and space. To replicate such behavior a plethora of cellular automata has been proposed. Among them the sandpile model [24,25], Manna model [26], the Olami-Feder-Christensen model [27], and the forest-fire model [28,29] play a prominent role. These models were also introduced as examples to illustrate the hypothesis of self-organized criticality (SOC) where systems evolve into a steady state which lack any characteristic length and time scales [22,24,25]. The initial assumption was that this state can be reached without fine tuning of any model parameters, i.e. the model dynamics organizes the system into a scale-free state. Subsequent analyses of the above and other models revealed that the tuning of parameters is present even in the case of "parameterless" models such as the sandpile model. In many cases the tuning parameter is the separation of slow and fast time scales present in the models though it is still debatable to which extent such a separation is necessary [32]. In the limit of large separation (slow driving and fast release of energy) some systems exhibit power-law tails in their avalanche size distributions.
Another feature which plays an important role in the dynamics of these models is a local conservation of energy. During avalanche events in some models the energy is conserved at each toppling site and only allowed to dissipate at the boundaries (the sandpile model, Manna model). Other models introduce a possibility of dissipation at interior sites (OFC model, fores-fire model). In some cases presence of local dissipation can alter the critical state.
A. The sandpile model The sandpile model was introduced by [24,25] as a paradigm of SOC to illustrate the emergence of scale invariant behavior in driven nonlinear threshold systems. The model is defined on a graph where each vertex has a height variable z ij ∈ {1, 2, ...} representing the amount of sand grains. The system is driven by slowly adding new sand grains at random into the system. Each vertex is also associated with a critical height z c . When a sand grain is redistributed to the vertex with z c amount of sand present the vertex becomes unstable and topples redistributing sand grains to its nearest-neighbors. This process can initiate an avalanche of toppling events. Starting from any initial distribution of heights after a transient period of time the model self-organizes into a steady state where the distribution of avalanches with respect to their sizes and durations follow a power-law behavior and long-range spatial correlations emerge [33]. This state has been identified as a SOC state in analogy to the second-order phase transitions encountered in statistical mechanics systems.

B. The Manna model
The Manna model is another example of a model exhibiting a SOC type behavior [26]. It was introduced as a model for particle collision and repulsion. The dynamics is pretty similar to the sandpile model except the height variable describes only two states: empty or occupied. If a new particle arrives at an empty site the site becomes occupied. Otherwise if the site is already occupied by a particle the new particle makes it unstable and two particles are transferred randomly to its nearest neighbors. This dynamics mimics the collision and repulsion of particles without conservation of momentum. After passing the transient regime the model enters into a steady state where distribution of avalanches display power-law behavior with well defined critical exponents [26].

C. The Olami-Feder-Christensen (OFC) model
The Olami-Feder-Christensen (OFC) model was introduced by [27] to analyze the behavior of massless blocks interconnected by springs and pulled over a surface with friction. In turn, the spring-block model [34] can be considered as a simple realization of a fault surface with threshold avalanche events representing the analog of earthquakes. The OFC model is locally nonconservative which results in different regimes of behavior depending on the degree of non-conservation specified by the parameter α.
In the model, each site (i, j) of a square lattice is assigned a continuous variable h ij ∈ [0, h c ] which defines a stress level. The model is driven uniformly with a constant rate until one of the sites reaches the threshold value h c = 1.0. This site becomes unstable and redistributes stress to its nearest neighbors and is reset to zero. The amount of stress transferred to each nearest neighbor site is αh ij , where 0 < α ≤ 0.25 defines a degree of non-conservation. This transfer may cause some of the neighbors to topple initiating an avalanche. The model is locally non-conservative for α < 0.25. The parameter α plays an important role in the dynamics of the model with smaller values of α leading to longer transient regimes [35] and also affecting its critical properties. However, the question of criticality of the OFC model still remains open [35].
The model produces rich behavior that can be explored in order to understand its dynamics. After passing the transient regime the model forms highly correlated spatial structures with power-law distributed patches of different heights. The existence of non-conservation in the model introduces several scales over which the distributions of avalanches are different at different scales. Large patches are mostly responsible for large avalanches while small avalanches are generated by an interplay between small patches [35].
Recent results concerning the spatio-temporal organization of avalanches in the OFC model show that the model generates, in addition to regular avalanches, near periodic avalanches associated with asperity structures formed as a highly concentrated stress state [36,37]. These asperity events form spatio-temporal clusters with pronounced characteristic features. The processes of synchronization and desynchronization in the OFC model with short-term correlations among events were also observed [38]. It was suggested that synchronization in the model drives it into a state with scale-invariant power-law distributed events while desynchronization results in the occurrence of aftershocks and foreshocks.
It was also suggested that the OFC model can exhibit the maximum in the entropy production rate for the range of the parameter α ∈ [0.15, 0.2]. This was associated with finite spatial memory of past avalanches associated with asperity events [39][40][41]. There is a tendency of similar large avalanches to reoccur as the model is driven forward. These avalanches rupture approximately the same area and occur quasi-periodically in time. The range of the parameter α where the maximum of the entropy production rate is achieved corresponds to the subcritical regime of the OFC model. A critical point was suggested to exist at α = 0.25 with a diverging correlation length and average avalanche sizes [41].

D. The forest-fire model
A model which resembles the occurrence and propagation of natural forest fires was proposed in [28,29]. The model is considered on a two-dimensional square lattice where each site can be occupied by a green tree or be empty. The dynamics of the model defines how trees are planted and how fires destroy the trees. Two parameters define the behavior of the model. A rate of planting of new trees, p, and a firing rate f p with which lightning can strike a site and start a fire. At each time step a random site is chosen. With the probability p a tree is planted and with the probability f lightning strikes a green tree. The fast time scale defines how model fires propagate in the system. A green tree adjacent to a nearest-neighbor burning tree will catch fire. This processes of spreading fires can initiate an avalanche which propagates until all connected green trees are burnt. It is assumed in the model that the process of the burning of fires is instantaneous. The control parameter in the model is the ratio θ = p/f [42].

IV. MODEL SIMULATIONS AND RESULTS
One of the main goals of this paper is to investigate the statistics of record-breaking avalanches in the above mentioned four cellular automata. Since these cellular automata are simple conceptual models mimicking real threshold systems with avalanche type dynamics, this will help us to answer questions related to the occurrence of recordbreaking events in models and real systems which display spatio-temporal correlations and non-trivial power-law scaling.
In order to perform numerical simulations of the above mentioned four cellular automata we used a square lattice of size L × L with varying L = 64, 128, 256 and performed in most cases 10 4 ensemble realizations of sequences of 10 5 − 10 6 time steps, where time steps were defined as the time of slow driving of the models not taking into account the time of propagation of avalanches. The statistics was collected after each model entered into a steady state to exclude transient effects. For each of the four models we constructed the various statistical measures of record-breaking events defined in Section II. The avalanche size for the sandpile, Manna, and forest-fire models was defined as the number of sites that toppled at least once during an avalanche. For these models the value of the avalanches were bounded by the lattice sizes. As a result, the record-breaking avalanches were bounded from above.
When simulating the OFC model for smaller values of the conservation parameter α the model generates a significant number of avalanches where only a single site relaxes. This directly affects the statistics of record-breaking events for early records since it leads to a non-negligible probability that ties among record values occur. Such ties do also occur for later records. To eliminate ties as well as the pathological limit states we introduced a low amplitude noise, η, into our simulations: during each toppling or relaxation event, a random number was drawn from a uniform distribution in the interval [−0.00005; 0.00005] and the toppled site was reduced by the value of its height plus/minus this noise η. This allowed us to eliminate ties by defining the size of an avalanche in the OFC model as the total amount of stress redistributed and dissipated during an avalanche for α < 0.25.

A. Results: The average number of record-breaking events
The average numbers of record-breaking events N n versus time step n obtained from the simulations of the four cellular automata are given in Fig. 1a. For the sandpile, Manna, and forest-fire models the behavior is close to the record-breaking events drawn from the i.i.d. random variables, Eq. To model the asymptotic behavior of the average number of record-breaking avalanches, N n , for large time steps we propose the following asymptotic equation where A, B and δ are constants. We fitted this equation to the average number of records obtained from the model simulations with the varying conservation parameter α and lattice sizes L = 64, 128, and 256. One such fit is given in Fig. 1a as a solid curve. The analysis indicates that the exponent δ increases with decreasing values of α. Then it reaches a plateau and starts to decrease moderately. To illustrate this the exponent δ is plotted versus the conservation parameter α in Fig. 3a. This figure also illustrates the dependence of the exponent δ on the lattice size L. We also plot the dependence of the parameter A on δ in Fig. 3b, that shows a similar trend. When the parameter α tends to the limiting value of 0.25, which corresponds to the full local conservation of energy in the OFC model, the average number of record-breaking events approaches the behavior corresponding to the records drawn from i.i.d. random variables and given by Eq. (5).
To confirm that the observed behavior for the OFC model was due to correlations among avalanche magnitudes we reshuffled the original sequences of avalanches which effectively destroyed any correlations among magnitudes. For these reshuffled sequences the average number of record-breaking avalanches followed exactly Eq. (5) for all considered values of α (not shown).

B. Results: Statistics of magnitudes of record-breaking avalanches
The avalanche frequency-magnitude statistics in the steady state for the most of the models follow power-law trend with finite-size effects. Both the sandpile and Manna models produce well defined power-law distributions for the sizes of avalanches with known power-law exponents [33,43]. The forest-fire model exhibits some variation in the distributions depending on the firing rate θ with presence of system wide fires for high values of θ and characteristic fires for small values of θ [42]. The behavior of the OFC model varies for different values of the conservation parameter α [35].
In this subsection we consider the distribution of record-breaking avalanche magnitudes. For the sandpile, Manna, and forest-fire models by magnitude we assume the area of toppled sites during an avalanche whereas for the OFC model it is the amount of stress transferred during toppling events in one avalanche as defined earlier. The distributions of the magnitudes of the kth record-breaking event, p k (x), are given for the OFC model with α = 0.18 in Fig. 4a. In the same figure we also plot as solid curves the theoretical distributions obtained from the record-breaking events drawn from the truncated power-law distribution, Eq. (25), and given by Eq. (10) with γ = 1.99, x min = 1.0, and x max = 256 2 . The model distributions deviate from the theoretical curves for several reasons. It is assumed that the magnitudes are not i.i.d. random variables. It is also evident from the distributions for avalanche magnitudes that they are not precisely power-law distributed. The average magnitudes x k of the kth record-breaking event for the four models are shown in Fig. 4b. These are compared for reference with the theoretical results applicable for i.i.d. random variables drawn from the truncated power-law distribution, Eq. (25), and given by Eq. (11). For the theoretical curves we used x max = 128 2 and the corresponding values of the exponent γ estimated from model simulations. These comparisons with the i.i.d. random variables drawn from the truncated power-law distribution are for illustrative purposes to show how model simulations deviate from them.
It is appropriate to note that when the record order k tends to infinity one expects that the distribution of magnitudes of the records has to converge to one of the three known classes (Gumbel, Fréchet, and Weibull) for the statistics of extreme events. Particularly, for the records drawn from the bounded distributions the distributions have to converge to a Weibull distribution. Another measure which characterizes the records is the average magnitude of record-breaking avalanches x(n) at a time step n. These have been computed for the four models and are plotted in Fig. 4c. It is observed that the average magnitude x(n) increases monotonically and approaches the upper bound set by the system size L 2 . For comparison we also plot as solid curves the average magnitudes obtained from Eq. (14) with the power-law distribution function Eq. (26). It is evident from the plots that the distributions and averages deviate from the theoretical curves. This can be attributed to the fact that the frequency magnitude distributions of the models on finite latices are not exactly given by Eq. (26) and the avalanches are not strictly i.i.d. random events. We also plot in Fig. 4d the average magnitude x(n) of record-breaking avalanches for the OFC model with the varying conservation parameter α and a fixed lattice size L = 128. The curves show a prominent trend where they experience a change in slope when going from small avalanches to large ones. This break in the behavior approximately corresponds to the transition observed in the average number of record-breaking events shown in Fig. 2a.
To confirm the presence of correlations between avalanche sizes in the OFC model we considered the size differences between subsequent avalanches, ∆a i = log 10 (a i+1 ) − log 10 (a i ), where i is the time step. A similar approach was used to study correlations between earthquake magnitudes [44][45][46]. In the absence of correlations the distribution of ∆a's should not deviate significantly from the distribution of the size differences constructed from the reshuffled catalogs, ∆a * i = log 10 (a i * ) − log 10 (a i ), where i * corresponds to the randomly picked index among N avalanches considered. In order to quantify the existence of correlations between avalanche sizes we considered the difference between cumulative distributions: ∆P (d) = P (∆a < d) − P (∆a * < d). This is illustrated in Fig. 5a for the four models. The dashed curve that fluctuates around 0 corresponds to the randomly reshuffled catalog where size correlations were destroyed. The error bars are given for 3σ confidence intervals. The OFC model displays strong deviation from the uncorrelated case indicating that avalanches in the OFC model are correlated. The deviations imply in particular that avalanches in the OFC model have the tendency to be followed by ones of similar size. On the other hand the avalanches in the forest-fire model show the opposite trend. This is expected for the forest-fire model as large fire avalanches are usually followed by small fires. This anti-correlation of avalanches for the forest-fire model does not change the statistics of records significantly [12]. The size differences between successive avalanches do not show strong correlations for the sandpile and Manna models. In Fig. 5b we also plot ∆P (d) for the OFC model for different values of the conservation parameter α. The strongest correlations are observed for the values of α below ≈ 0.16 which also corresponds to the higher values of δ observed in Fig. 3a.

C. Results: Interevent time distributions
Here we summarize the results obtained for the temporal structure of record-breaking avalanches of the four cellular automata. We computed the distribution that the kth record was broken after m time steps. This is shown in Fig. 6a for the OFC model with α = 0.17 for the first several orders of k = 1, 2, ..., 10. For comparison we plot the densities computed for the record-breaking events drawn from i.i.d. random variables given by Eqs. (19) which are independent of the underlying distribution F (x) of records. The comparison shows that the distributions for the OFC model deviate from the i.i.d. case for higher orders of k and large interevent time steps m. The other three models do not show such a departure and the corresponding distributions resemble the i.i.d. case.
The results of simulations for the interevent times between all subsequent record-breaking events for all four models are given in Fig. 6b. The normalized distributions, g(m), for the sandpile, Manna, and forest-fire models follow a simple power law: ∼ 1/m which is similar to the case of record-breaking events drawn from i.i.d. random variables [4]. These distributions were constructed by counting all interevent times for all record orders k. The OFC model shows departure from the i.i.d. case where one can observe a change in the behavior for small and large interevent times and strong dependence on the conservation parameter α (Fig. 6c). The probability density functions, u k (n), for the time of the occurrence of the kth record-breaking event versus a time step n are shown in Fig. 7a for the OFC model with α = 0.17 for the first several orders of k = 1, 2, ..., 10. In addition, the average time n k of the occurrence of the kth record-breaking avalanche is given in Fig. 7b. For comparison we also plot the average times obtained from the simulations of i.i.d. random variables for sequences with 10 6 time steps. This also illustrates the effect of the finiteness of the sequences on the average numbers. Fig. 7b demonstrates that the results for the four cellular automate are pretty close to the i.i.d. case with some deviations for the OFC model. On the other hand the finiteness of the sequence plays an important role as the average times of the occurrence are limited by this upper time limit that affects the average times for larger record orders k. In Fig. 7c we also plot the average time n k for the OFC model with the varying parameter α.

V. CONCLUDING REMARKS
Record-breaking events can be extracted from various natural and model systems. Their statistics can be used to analyze underlying trends and effects of correlations in the dynamics of those systems. Previous analyses of recordbreaking events drawn from i.i.d. random variables have revealed robust and universal behavior of temporal statistics of records independent of the underlying distribution functions from which the random variables were generated. The main goal of this paper is to find deviations from that universal behavior by studying record-breaking avalanches generated by driven nonlinear threshold systems. Particularly, we have considered four prominent cellular automata: the sandpile, Manna, OFC, and forest-fire models. In addition, we have also derived an exact formula, Eq. (18), for the probability that the kth record is broken after m time steps. The performed analysis reported in this paper revealed that among four models studied the OFC model displayed the most significant deviations in the statistics of record-breaking events compared to ones drawn from i.i.d. random variables. It was found for the OFC model that the average number of records for large avalanches can be described by Eq. (27) with the exponent δ being dependent on the conservation parameter α. This result signifies the departure from the case of i.i.d. random variables and is the direct effect of non-trivial correlations between avalanches in the OFC model. The other three cellular automata (the sandpile, Manna, and forest-fire models) do not show such a drastic departure indicating that their record-breaking avalanches are weakly correlated.
The distributions of record magnitudes for the four models also displayed deviations from the i.i.d. case. These distributions are controlled by the statistics of avalanche sizes which are not pure truncated power-laws with strong finite-size effects and dependence on model parameters. The correlations between avalanches also contribute to the observed behavior. This results in a more complex form of the distributions of record magnitudes and the corresponding averages.
The temporal structure of the occurrence of records also displays deviations compared to the i.i.d. case. Specifically, the distributions of the times of the occurrence of the kth record, u k (n), and the distribution for the kth record being broken after m time steps, w k (m), computed for the OFC model display departure from the corresponding distribution constructed from the i.i.d. random variables. On the other hand, the distribution of the interevent times computed for all record orders k do not show strong departure as well as the average time of the occurrence of the kth record, n k .
The finite-size effects both in time and magnitude domains play a prominent role in model simulations. They strongly affect statistics of records and have to be taken into account. The finiteness of the lattice introduces an upper bound for the magnitude of record-breaking avalanches. More thorough understanding of these effects requires development of more sophisticated simulation techniques in order to generate enough statistics of large records.
The OFC model can be considered as a simplification of a single fault system to simulate earthquake occurrences.
The avalanches in the model reproduce some important aspects of statistics of earthquakes such as a power-law frequency-magnitude distribution. Therefore, the analysis of correlations in the OFC model can shed light on the dynamics and triggering mechanisms of real earthquakes. The obtained results concerning the record-breaking avalanches in the OFC model can be used to constrain the effects of correlation and self-organization in seismogenic zones. Concerning real seismicity, the sequences of record-breaking earthquakes were extracted and analyzed for world-wide earthquakes [9,10]. The analysis of annual maximum earthquake magnitudes for the global earthquake catalog was also performed [47]. The results from those studies indicate that global earthquakes are independent and their magnitudes follow approximately the exponential distribution. The analysis of record-breaking earthquakes occurring on a single fault is problematic due to the limited records of existing catalogs. This makes it difficult to directly compare the correlations between avalanches we observe in the OFC model to real earthquake data.