A stochastic theory for temporal fluctuations in self-organized critical systems

A stochastic theory for the toppling activity in sandpile models is developed, based on a simple mean-field assumption about the toppling process. The theory describes the process as an anti-persistent Gaussian walk, where the diffusion coefficient is proportional to the activity. It is formulated as a generalization of the It\^{o} stochastic differential equation with an anti-persistent fractional Gaussian noise source. An essential element of the theory is re-scaling to obtain a proper thermodynamic limit, and it captures all temporal features of the toppling process obtained by numerical simulation of the Bak-Tang-Wiesenfeld sandpile in this limit.


Introduction
The existence of self-organized critical dynamics in complex systems has traditionally been demonstrated through numerical simulation of certain classes of cellular automata referred to as sandpile models [1]. Non-linear, spatio-temporal dynamics is always essential for the emergence of SOC behavior, but the details of this dynamics for a specific natural system is often poorly understood and/or not accessible to observation. In many cases the information available is in the form of time-series of spatially averaged data like stock-price indices, geomagnetic indices, or global temperature data. For scientists who deal with such data a natural question to ask is: are there specific signatures of SOC dynamics that can be detected from such data?
In this letter we shall report some results which provide a partial answer to such a question. Some important statistical features of the toppling activity are common to most weakly driven sandpile models described in the literature, and these are used to formulate a stochastic model for the toppling activity signal. A benchmark case against which our results are tested, is a numerical study of the Bak-Tang-Wiesenfeld (BTW) sandpile [2]. A crucial step in our work is a re-scaling of the dynamical variables which allows a natural passage to the thermodynamic (continuum) limit. We demonstrate that this leads to new results concerning SOC scaling laws. We find that the probability density function (pdf) for the toppling activity is a stretched exponential or close to the Bramwell-Holdsworth-Pinton distribution [3], depending on whether the sandpile is so slowly driven that avalanches are well separated, or it is driven so hard that several avalanches run simultaneously. The pdf for avalanche durations is unique in the thermodynamic limit, but is not a power law, unless we redefine the meaning of an avalanche to be the activity burst between successive times for which the activity rises above a positive threshold. Implementing such a threshold yields an exponent for the avalanche duration pdf of 1.63, in agreement with [4], but in contradiction to [7]. It also gives power-law quiet-time statistics as in [4] and thus refuting the claim in [8] that SOC implies power-law distributed avalanche durations, but Poisson-distributed quiet times.
The sandpile models considered in this short paper deal with a d ≥ 2-dimensional lattice of N d sites each of which are occupied by a certain integer number of quanta which we conveniently can think of as sand grains. The dynamics on the lattice is given by a toppling rule which implies that if the number of grains on a site exceeds a prescribed threshold, the grains on that site are distributed to its nearest neighbors. If the occupation number of some of these neighbors exceed the toppling threshold these sites will topple in the next time step, and the dynamics continues as an avalanche until all sites are stable. The details of this toppling rule can vary, but a useful theory for a broad class of natural phenomena should not be very sensitive to such detail.
In natural systems the SOC dynamics is usually driven by some weak random external forcing. In sandpile models this can be modeled by dropping of sand grains at randomly selected sites at widely separated times. In numerical algorithms this is often done by dropping sand grains only at those times when no avalanche is running. This ensures that the drive does not interfere with the avalanching process. Usually it will then only take a few time steps from one avalanche has stopped until a new starts, so for a large system the quiet times between avalanches will appear insignificant compared to their durations. A more physical drive would be to drop sand also during avalanches. If the dropping rate is slower than the typical duration of a system-size avalanche the drive would still not interfere with the avalanche dynamics, but the quiet times would depend on the statistical distribution of dropping times, which is typically a Poisson distribution. In many natural systems, however, avalanching occurs all the time, corresponding to a higher driving rate. In such cases, and also because there will always be noise in timeseries data, we cannot identify the start and termination of an avalanche from a zero condition of our observable. In practice we have to define avalanches as bursts in the time series identified by a threshold on the signal [4]. In a sandpile simulation such bursts are correlated and therefore the quiet times between the bursts are power-law distributed even if the dropping of sand grains is chosen to be a Poisson process. Hence if focus is on modeling features that can be detected in observational data we shall think of avalanches as activity bursts starting and terminating at a non-zero threshold value. Moreover, one of the main results of this work is that power-law shape of the pdf for avalanche duration is true only if one defines avalanches in this way.

The stochastic model
We shall assume that the lattice has linear extent L = 1 with N d sites, so the thermodynamic limit N → ∞ can be thought of as a continuum limit. The sandpile evolves in discrete time steps labeled by k = 1, 2, 3 . . ., and the number of sites whose occupation number exceeds the toppling threshold at time k is called the toppling activity x N (k). The toppling increment is ∆x N (k) def = x N (k + 1) − x N (k). Let us define two active sites as dynamically connected if they have at least one common nearest neighbor, and define a connected cluster as a collection of active sites which are linked trough such connections. From numerical simulations of sandpiles we observe that such clusters never consist of more than a few elements and that the instantaneous number of clusters n N increases in proportion to x N . This implies that at each time k we can label the clusters by i = 1, . . . , cx N (k), where c < 1 is a constant depending on the specific toppling rule and the dimension d of the sandpile. We can then decompose the increment ∆x N (k) into a sum of local increment contributions ξ N,i (k) produced by each of the clusters, i.e. ∆x . We think of the local increment contributions as random variables which take values in a finite sample space. Indeed, if each cluster i only consists of a single overcritical site, then ξ i,N takes values in the set As a first step to a stochastic model we make a mean-field assumption [10,11], which impiles that ξ N,i (k) and ξ N,j (k) are statistically independent for i = j. Then the central limit theorem states that in the limit N → ∞, . This has been verified numerically in the two-dimensional BTW-model as shown in Fig. 1. The figure demonstrates the need to introduce a conditional probability: The conditional variance of the increments is proportional to x N and the conditional mean is not zero.
In fact, numerical simulations show that the the conditional mean increment, E[∆x N |x N ], is positive for small x N , reflecting the natural tendency for the activity to grow when it is small. On the other hand the mean increment decays exponentially to zero for moderate x N , and becomes negative when x N is comparable to the activity of a system-size avalanche, reflecting the limiting influence of the finite system size. These effects will be incorporated as a drift-term correction to the model, but for now we consider for simplicity of argument a Gaussian process with non-stationary increments and no drift term: where w(k) is a stationary Gaussian stochastic process with unit variance. From the numerical sandpile data (see Fig. 1) we observe that the normalized toppling process has the characteristics of a fractional Brownian walk with Hurst exponent H ≈ 0.37 on time scales shorter than the characteristic growth time for a system-size avalanche, consistent with a power spectrum which scales like f − Numerical simulations show that x N ∼ N D 1 ‡, where 0 < D 1 ≤ d can be interpreted as a fractal dimension of the set of active sites imbedded in the d-dimensional lattice space. This property is used to re-scale x N (k) such that it has a well-defined limit as N → ∞. We also have to re-scale the time variable by letting t = k∆t, where ∆t = N −D 2 . The value of D 2 will become apparent if we define the normalized activity variable X N (t) = N −D 1 x N (t/∆t), such that the corresponding increment becomes where ∆W H (t) = W H (t + ∆t) − W H (t). A well-defined thermodynamic limit N → ∞ requires D 2 = D 1 /2H, for which Eq. (3), by introduction of the limit function X(t) = lim N →∞ X N (t), reduces to the stochastic differential equation where we have heuristically added a drift term f (X) dt to account for the non-zero mean of the conditional increment. We take f (X) to be an exponentially decaying function based on the numerical results from the sandpile. In the 2-dimensional BTW model we find that D 1 ≈ 0.86 and hence D 2 = 1.16. This defines re-scaled coordinates X N = x N /N 0.86 and t N = k/N 1.16 .

Analysis of avalanches
A time series X(t) ≥ 0, representing a succession of avalanches with zero quiet times, can be constructed numerically from the discrete-time version of Eq. (4) by integrating the equation using realizations of the fractional Gaussian noise process ∆W H (t). At those times when X(t) drops below zero we consider the avalanche as terminated, and a new, independent realization of ∆W H (t) is generated and used to produce the next avalanche. From long, stationary time-series generated from the stochastic model and from the sandpile model this way, we can construct pdfs P(X) which turn out to give almost identical results for the two models (see Fig. 2). The shape of this pdf is universal in the thermodynamic limit: a stretched exponential P(X) ∼ exp (−aX µ ) with µ ≈ 0.5.
A different pdf appears if the time-series are constructed by launching the avalanches at random times (Poisson-distributed) with characteristic time between launches shorter than the growth time of a system size avalanche. In this case several avalanches may run simultaneously, and P(X) from both models are close to the Bramwell-Holdsworth-Pinton distribution, which was claimed to be valid for the toppling-activity in the BTWmodel in [3]. Consider a solution of Eq. (4) with initial condition X(0) = Y > 0, and let P (X, t) be the evolution of the density distribution in X-space of an ensemble of realizations of the stochastic process X(t) all launched at activity X = Y at time t = 0. Every realization X(t) will sooner or later terminate at a finite time t = τ for which X(τ − 1) > 0 and X(τ ) ≤ 0, and then we remove it from the ensemble. P (X, t) contains information about all commonly considered avalanche characteristics. For example, it is easily found from from Eq. (4) that, on time scales shorter than the growth time of a system-size avalanche, X(t) is a self-similar process with non-stationary increments and self-similarity exponent h = 2H [6]. Hence the variance of X(t) with respect to P (X, t) will scale as ∼ t 2h . That this relation holds for the 2-dimensional BTW model can easily be verified through numerical simulation (Fig. 3(a)).
We can also compute the survival probability ρ(τ ) = ∞ 0 P (X, τ ) dX, which is the probability that a realization of an avalanche has not terminated at the time τ . This function is related to the pdf for avalanche durations by p dur (τ ) = −ρ ′ (τ ), so that p dur (τ ) is a power law if and only if ρ(τ ) is a power law. Fig. 3(b) shows the function ρ(τ ) for numerical simulations of the BTW sandpile in the re-scaled coordinates X N and t N , demonstrating that the pdf for avalanche durations does not represent a power law. The power-law form ρ(τ ) ∼ τ 0.5 proposed in [7] can only be obtained as a tangent to the log-log plot of ρ(τ ) at a given duration time τ , and the slope of this tangent depends crucially on the duration time τ for which this tangent is drawn.
The situation changes if we let all avalanches terminate when X drops below a small threshold X c > 0 as proposed in [4]. In this case avalanche durations are the return times to the line X = X c , and by changing coordinates to Y = X − X c we see that this corresponds to the return times to the time axis of the process given by the stochastic differential equation For small avalanches where X(t)−X c ≪ X c we can approximate this expression with dY (t) = σ √ X c dW H (t), i.e. can approximate Y (t) by a fractional Brownian motion with Hurst exponent H. Using the result of Ding and Yang [9] on the return times of a fractional Brownian motion we get p dur (τ ) ∼ τ 2−H = τ −1.63 .
Numerical simulations of the BTW model verify this result: The survival function ρ(τ ) becomes a power law on time scales shorter than a system-size avalanche (see Fig. 4(a)), and the slope of the graph in a log-log plot is approximately −0.63, which corresponds to a scaling of the pdf for duration times on the form p dur (τ ) ∼ τ −1.63 . The result is also reproduced by simulations of Eq. (4) with an exponentially decaying drift term. Fig. 4(b) shows the log-log plot of the pdf for duration times in the stochastic differential equation and a line with slope −1.63, demonstrating that the avalanche statistics in the BTW sandpiles is captured by the stochastic differential equation.  From the scaling ρ(τ ) ∼ τ −α we can deduce an exponent for the pdf of avalanche size as well. On the time scales where the toppling activity can be approximated by a fractional Brownian motion W H (t), the signal disperses with time as X ∼ t H , the size of an avalanche of duration τ scales like S(τ ) ∼ τ 0 t H dt ∼ τ H+1 . Assuming that the pdf for avalanche size is on the form p size (S) ∼ S −ν , the relation p size (S) dS = p dur (τ ) dτ yields τ −ν(H+1)+H ∼ τ −α−1 , so With H = 0.37 we obtain ν = 1.46. We also remark that if we omit the drift term and let H = 1/2 and X c = 0 we obtain the so-called mean-field theory of sandpiles. In this case the stochastic differential equation has a corresponding Fokker-Planck equation ∂P ∂t = σ 2 2 ∂ 2 ∂X 2 (XP ) . If we solve this equation on the interval [0, ∞) with absorbing boundary conditions in X = 0 we can obtain an analytical expression for P (X, t), and from some straightforward algebra we find for large τ that p dur (τ ) ∼ τ −2 [6]. Since X c = 0 we can not approximate the toppling activity by a Brownian motion on any scale and thus X(t) diverges like ∼ t h , where h = 2H. By replacing H with h = 2H in Eq. (5) we get p size (S) ∼ S −3/2 , in agreement with previous mean field approaches [10,11].

Concluding remarks
We point out that the validity of Eq. 4 is not restricted to the BTW model. For instance, the equation has been verified for the Zhang model [6,12], though with a different Hurst exponent H. Time series of global quantities derived from numerical simulation of different sandpile and turbulent fluid systems can be shown to be adequately described by Eq. 4, where H and the specific form of the drift term depend on the system at hand [6].