The role of intervention mechanisms on a self-organized system: dynamics of a sandpile with site reinforcement

The first and most straightforward manifestation of the possible effects of reinforcement mechanisms is the resulting statistics of the avalanches generated by the sandpile. Figure 2 shows the normalized probability density functions of the sandpile avalanche areas, P(A), upon the imposition of the different reinforcement parameter measures, f.

Zoom In Zoom Out Reset image size

The original sandpile f = 0 produces a stable power-law distribution across a minimum of three decades, $P(A) \sim A^{-\alpha}$, where the statistically obtained range for α from various trials is $\alpha = 1.39 \pm 0.03$. In particular, the best fit power-law curve for the specific run of f = 0 that is presented in figure 2 is the solid gray line following $P(A) \sim A^{-1.41}$, which slightly raised in the plot for clarity and as a guide to the eye. As the value of f is increased by orders of magnitude, deviations from the power-law scaling behavior start to emerge, shortening the tails of the distributions. To illustrate the extent of this deviation, the exponential distribution with the same mean value as the f = 1 data is plotted as the dashed gray curve, $P(A) \sim \text{exp}\left(-A/10.1\right)$. For the same rate and duration of external driving, the tails of the distribution significantly degrade. The power-law trend for f = 0 is characteristic of the long-range spatial correlations in the sandpile [24]. On the other hand, the exponential distribution is usually obtained for memoryless processes with a well-defined characteristic value and finite variance. The fact that the higher f grids produce statistics that lack the power-law trends and instead show decay trends closer to the exponential indicates that the local reinforcement mechanisms disrupt the spatial memory in the sandpile grid. Similar trends were obtained previously for processes that involve static spatial sinks [37], where the resulting distributions are parametrized by stretched exponential curves [38]. It should be noted, however, that the memory-based reinforcement is fundamentally different from static implementations; here, the reinforced sites, while more stable, are still active participants that may collapse again at some point. It is also worth noting that the decrease in occurrence of very large avalanches is compensated by the increase in the occurrence of the smallest ones as the fraction f is increased. This is most apparent from comparing the f = 0 and f = 1 plots in figure 2 for the regimes of the smallest avalanches A < 15.

The trends in P(A) for increasing f highlight the effect of the local stabilization mechanisms on the self-organized sandpile grid driven for finite iteration times. Introducing stabilizing intervention mechanisms for finite driving times lead to disruption in the original mechanisms, which, in turn, leads to a reduction of the probability of very large events and a corresponding increase in the occurrence of very small events. The inset of figure 2 further highlights this in the form of the decaying trend for the average avalanche size $\langle A \rangle$ vs. f.

Moreover, the observed decay trends in figure 2 is recovered for various grid sizes, as shown in figures 3(a)–(d). For finite iteration times, the memory-based reinforcement creates very high threshold sites that affect the internal cascade. In the presence of these reinforced sites, the avalanches in the grid fail to propagate to a wider area, impeding the long-range spatial correlations within the grid. The system therefore releases its accumulated stress in a multitude of small-area avalanches with smaller characteristic sizes.

Zoom In Zoom Out Reset image size

To show that the observed distributions are not due to transient effects, representative distributions for the f = 1 case (the largest reinforcement parameter used) that are run for longer (at least an order of magnitude greater) iteration times are presented in figure 3(e). The P(A) distributions for L = 128 gradually lose the truncation observed for the case of finite waiting times in figure 2. In fact, for smaller grids L = 32 run for $t_{\text{max}} = 10^9$, the P(A) shows a power-law behavior, albeit with steeper exponents: $P(A) \sim A^{-\alpha^{^{\prime}}}$ where $\alpha^{^{\prime}} = 1.76 \pm 0.06$. For comparison, the P(A) distributions for f = 0 (original sandpile) with L = 128 and $t_{\text{max}} = 10^8$, along with the solid line denoting the power law trend A^−1.41 (also shown in figure 2) is plotted together with the f = 1 distributions, with the dashed line A^−1.77 showing a good agreement, particularly with the L = 32 runs. The recovered steeper power-laws for higher f still indicate the increase is occurrence of the smaller avalanches and the corresponding decrease in the largest avalanches. However, the continued existence of power-law distributions appear to indicate the persistence of self-organization in the sandpile grid. By keeping the external driving and the internal cascade rules, and adding only a within-system adjustment, power-law avalanche area distributions still emerge without the need for a specific, finely-tuned set of parameters, consistent with the SOC nature of the sandpile.

As such, reinforcement paradigms, especially with the aim of possible applications to actual SOC systems, have inherent practical and fundamental limitations. At a practical level, not only do we need to introduce intervention mechanisms whose magnitudes are comparable to that of the relaxation event, we need to introduce it repreatedly, resulting in the existence of 'super' sites with very large thresholds of stability. Realistically speaking, there is a finite limit to our capability to reinforce any system, rendering this mechanism improbable, if not completely impossible. At the fundamental level, on the other hand, the system will drive itself into a critical state that produces power-law magnitude distributions despite the presence of memory-based reinforcements at the local site level. Adaptive local interventions within the system itself, in general, cannot circumvent the SOC mechanisms. This may explain the robustness of the observed power-law tails of the distributions of well-known complex systems.

The origins of the nonzero avalanches provide information regarding the space and time separations of the successive avalanches. In the following, we present the statistics of the successive-event (interevent) distances R and the successive-event (waiting) times T for the various f values used.

Because the sandpile grid is finite, the probability density functions of the separation distances of successive events, P(R), show a common unimodal distribution that has a characteristic value $\langle R \rangle = 67.3 \approx L/2$, i.e. half of the grid dimension. The random choice of trigger origins every time step in the application of equation (1) results in these common P(R) distributions across all f values, as shown in figure 4. In the figure, the normal distributions corresponding to the mean and standard deviation values of the R data are presented as the lines with the same color. The inset of figure 4 shows that the average separation distance $\langle R \rangle$ is consistent regardless of the stabilization fraction used.

Zoom In Zoom Out Reset image size

Interestingly, however, the waiting time statistics P(T) show different characteristic distributions across the range of reinforcement fractions f investigated. In figure 5(a), the P(T) distributions all show cutoff values at the tails in a double logarithmic scale. Using the f = 0 (original sandpile) as the reference, we observe the shorter tails for up to f = 0.001, and then a lengthening of the tails for f = 0.01 and above. To further highlight this behavior, the inset of figure 5(a) shows the average successive-event time, $\langle T \rangle$, showing a dip for f = 0.0001 and f = 0.001 before increasing again for higher f values.

Zoom In Zoom Out Reset image size

Previous works have shown that the waiting time distribution for the sandpile is best fitted with exponential decays [44]. As such, figures 5(b)–(e) shows comparisons between P(T) statistics for representative f values in (a) with the corresponding exponential distribution having the same characteristic (mean) value. While the original sandpile f = 0 closely follows the exponential trend, increasing f tends to make the tails heavier. It should be noted that the plots in figures 5(b)–(e) are presented in semilogarithmic scales, making the exponential curve behave like a straight line with slopes equal to $\langle T \rangle$ (whose values are shown in the inset of figure 5(a)). We observe that the P(T) statistics of $f = 0.000\,01$, f = 0.0001, and f = 0.001 (which have lower $\langle T \rangle$ than f = 0 as shown in the inset of figure 5(a)) have high statistical goodness-of-fit values with an exponential distribution, while those of f = 0.01, f = 0.1, and f = 1 (which have higher $\langle T \rangle$ than f = 0) have large deviations at the tails that degrade the exponential fitting.

The results observed in the P(T) statistics emphasize the dual nature of the stabilization rule equation (3) imposed on the sandpile grid. On one hand, the increase in the threshold helps the local site resist the local toppling and redistribution from equation (2). On the other hand, once the stabilized site gets toppled, the large amount of 'grains' that it contains may introduce a wider cascade into the neighborhood. In the observed regime of shorter (than the original sandpile) average waiting times, $f \in \left\{0.0001, 0.001\right\}$, the interplay between these two effects are deemed to be at a critical balance. It should be noted, however, that under these regimes, the grid is only more active (i.e. shorter average interevent times; see inset of figure 5(a)) but still less catastrophic (i.e. lower characteristic avalanche size than the original sandpile; see inset of figure 2). Even at the most active regimes, the probability of occurrence of the largest avalanches are suppressed, so the grid is still globally more stabilized.

How are the R and T statistics affected by the space and time scales used? To investigate this, we consider the longer runs and different grid sizes for the f = 1 case, which are also used to generate figure 3(e). In figure 6(a), we observe, as expected, that the P(R) statistics are only related to the spatial scales of the system. A simple rescaling by $\langle R \rangle$ that preserves the unit area under the probability density functions resulted in the data collapse for all the L and $t_{\text{max}}$ values obtained. On the other hand, the waiting time statistics is naturally affected by the area, i.e. the number of grid sites. In producing an avalanche event, especially the largest ones, bigger grids need to be triggered for longer periods compared to smaller ones. Incidentally, rescaling the T axis by the number of sites L² while preserving the unit areas under the probability density function curve also resulted in the data collapse, as shown in figure 6(b).

Zoom In Zoom Out Reset image size

To further explain the avalanching behavior of the grid, the distributions of the individual states of the grid sites, $P\left(z_{xy}\right)$, and the corresponding distribution of their thresholds, $P\left(Z_{xy}\right)$, are plotted side by side in figure 7 for representative f values. To accomplish this, we allowed for the system to evolve for up to $t = t_{\text{max}} = 10^8$ iterations. At that point, all the states z_xy and thresholds Z_xy of the L² individual grid sites are taken and analyzed for their statistics. The plots in figure 7 are therefore just snapshots of the grid, and may thus be different if other time instances are considered. Here, however, by fixing the time to $t = 10^8$ for all f values, we know that the system has already experienced repeated avalanching and stabilization mechanism; additionally, using a common finite driving time for analyses will isolate the effect of f on the grid.

Zoom In Zoom Out Reset image size

The inset of figure 7(a) shows that the original sandpile f = 0 evolves to a stable shape that peaks around z ≈ 0 (corresponding to the newly-collapsed sites) and $z = Z_0/4$, $z = Z_0/2$, and $z \approx Z_0$, i.e. the quartile values due to the redistribution rule to the four-nearest neighbors given by equation (2). Without stabilization, this shape of the distribution will persist, with the peaks just shifting near the multiples of $Z_0 / 4$. The scaling of the inset is boxed and shaded in the main plot for reference. Everything that is beyond this shaded region is due to the effect of stabilization, allowing for sites to have states greater than Z₀. Interestingly, in the main plot of figure 7(a), the distributions $P\left(z_{xy}\right)$ show tails whose extent differ for the various f values used. The $f = 0.000\,01$ has comparable statistics to the original sandpile; however, as f is increased, the tails become heavier, until the broadest tails at f = 0.01 is obtained. Interestingly, increasing f beyond 0.01 resulted in the reversal of the trend, i.e. a shortening of the tails, reminiscent of the trends observed in the waiting times in figure 5(a).

The mean and the standard deviation of z_xy for increasing f values are also shown in the lower plots of figure 7(a). While there is a monotonic increase for the mean, the standard deviation shows a prominent peak for f = 0.01, another manifestation of the heavy-tailed statistics for this reinforcement value. This also gives us an idea regarding the avalanching behavior of the grid. For the case of f = 0.001, the activity of the grid resulted in the dominance of small-neighborhood collapses, which resulted in the sharp peaks at z ≈ 0. As a result, the plot for f = 0.01 shows a preponderance of almost empty, recently collapsed sites and a very long tail corresponding to several very critical sites, which will eventually redistribute to its local neighborhood.

On the other hand, the distributions of the local site thresholds $P\left(Z_{xy}\right)$ in figure 7(b) show an expected lengthening of the tails, this time governed by the reinforcement rule equation (3). From an expected Dirac delta distribution centered at Z₀ = 1 for f = 0, the $P(Z_{xy})$ begin to show a lengthening of the tails as time progresses, as shown in figure 7(b). The mean and the standard deviation of the thresholds shown in the lower plots of figure 7(b) show generally increasing trends as well. Regardless of the size and frequency of the avalanches in the grid, the memory-based reinforcement equation (3) will tend to increase the thresholds of the sites of origin.