The Forest Fire Model: The subtleties of criticality and scale invariance

Amongst the numerous models introduced with SOC, the Forest Fire Model (FFM) is particularly attractive for its close relationship to stochastic spreading, which is central to the study of systems as diverse as epidemics, rumours, or indeed, fires. However, since its introduction, the nature of the model's scale invariance has been controversial, and the lack of scaling observed in many studies diminished its theoretical attractiveness. In this study, we analyse the behaviour of the tree density, the average cluster size and the largest cluster and show that the model could be of high practical relevance for the activation dynamics seen in brain and rain studies. From this perspective, its peculiar scaling properties should be regarded as an asset rather than a limitation.


INTRODUCTION
Soon after the seminal paper introducing Self-Organised Criticality (SOC) [1], it was suggested that examples of SOC could include models describing the spread of activation in a manner reminiscent of forest fires or infectious diseases. The degree to which these models were examples of scale invariance and criticality instantly became a subject of intense debate, see, e.g. [2] and [3]. Despite the controversy and indications that the Drossel-Schwabl Forest Fire Model lacks scale invariance [2], the dynamics of the model is of a type that seems of direct relevance to many real systems such as the brain [4] and precipitation [5]. So if for no other reason, it is still worthwhile to develop a better understanding of the behaviour generated by this kind of stochastic spreading and relaxation dynamics and to develop ways to probe the dynamics which can be applied to data from real systems.
As the name implies, the approach of SOC is to focus on criticality in the sense of equilibrium statistical mechanics, and for this reason, one typically looks for scale invariance and dynamics that can tune sharply to a critical point. Consequently, broad crossover behaviour is not seen as properly belonging to the SOC paradigm. However, the relevance of a theoretical scientific framework is, in the end, determined by how useful it is for the description and analysis of real systems. Seen from this perspective, it is of importance to bear in mind that exact fine-tuning to a completely scale-invariant state is not always observed in systems that exhibit SOC-like behaviour. This is demonstrated, e.g. by the studies of the size distribution of rain showers [5], and the bursts of brain activity measured during fMRI scans [4]. Indeed, neither study finds sharp critical behaviour but, despite the systems being of totally different microscopic nature, both identify similar indications of critical behaviour in terms of approximate power laws and even features reminiscent of peaked, or perhaps diverging, fluctuations or susceptibilities. Moreover, both studies find that the dynamics pulls the systems into a crossover region of large fluctuations within which one most frequently finds configurations to reside. In other words, in both cases, the distribution of residence times, i.e. the amount of time spent at a certain value of the control parameter, is found to exhibit a broad peak centred about what appears to be a critical-like value of the control parameter. From our experience with equilibrium critical phenomena, one would expect a broad peak for small system sizes only, and that the width of the peak will shrink with increasing system size. Therefore, for systems as big as the atmosphere or with as many components as the brain, broad peaks would not be expected.
The experimental observations in [4,5] suggest that the dynamics couples the control parameter to the fluctuations in a way that makes the system move around in a critical region, rather than tuning to a critical point. This is similar to suggestions previously put forward, such as [6,7]. Intuitively, one may imagine something like the following in the case of precipitation: nucleation of drops can happen at a particular value of the vapour content. The vapour in the atmosphere over the ocean gradually builds up towards that value, and sometimes overshooting may even occur before nucleation is seeded. When precipitation events occur, a significant amount of vapour is removed from the atmosphere, and one observes oscillations between subcritical and supercritical regions. For some reason, the coupling between the driving (vapour formation) and the response (precipitation) produces such large fluctuations that a precise tuning to the critical value of the control parameter is excluded. A parallel situation may take place in the brain as neurones have to go through a refractory period before they are able to fire again after discharge.
A link between the behaviour observed in [4,5] and the Forest Fire Model was already suggested in [8]. Here we analyse in detail the residence time distribution (which corresponds to the tree density distribution in the FFM) and its relationship to the control parameters used in standard percolation [9] and [4], namely the mean cluster size and the size of the largest cluster. We find that there is indeed a critical-like region where the fluctuations peak up and where the system seems to spend most of the time, exactly like in the brain and rain observations. Furthermore, analysing the distribution of the largest cluster we find that it displays scale invariance for our reachable sizes of θ and that the average cluster size seemingly grows superlinearly with θ, in a way that reminds the superlinear growth of the instantaneous correlation length observed in [10].

MODEL DESCRIPTION
As reported in [3], the original forest fire model proposed in [11] was revised in [2] and became known in the literature as the Drossel-Schwable Forest Fire Model (DS-FFM), despite the resulting model coinciding with the one introduced a few years before in [12] as Self-organised percolation. In the following, we will analyse the model proposed in [2,12] and refer to it simply as forest fire model (FFM). The FFM consists of a dissipative dynamic that involves the occupation of empty sites (planting of trees) and the removing of clusters of trees (burning of a forest). In the following, we will restrict our simulations to a two-dimensional square grid with periodic boundary conditions and side length L. Our implementation follows [13][14][15][16], and is summarized by the following pseudo-code. At each time step, we select a site at random; if it is empty, we try to plant θ trees, if it is occupied, we burn the whole cluster connected to it, considering 4 neighbours for each site. In order to avoid finite-size effects, θ has to be tuned, taking into account the system size L 2 . In our simulations we keep fixed the ratio k = θ L 2 and use periodic boundary conditions as in [10]. The great majority of the simulations are done at k = 10 −3 , but we also present results for smaller values of k in Fig. 3 and in Fig. 10. Since the model is known for taking a long time to thermalize, we performed 5 * 10 6 burning steps for thermalisation and 10 6 to collect the statistics as in [16]. The only exception is for the heatmaps in Fig. 7 and Fig. 8, where 10 7 datapoints were used to produce the statistics.

The distribution of densities
To compute the distribution of the densities P (ρ), we assign one spatially averaged density ρ to each generated configuration and sample the distribution over the ensemble. In this way, we obtain an object that is equivalent to the distribution of residence times in [4,5]. In the early studies of the FFM, the average density of trees ρ was assumed to behave as The value of the power was estimated as b = 0.47 in [17,18] and b = 0.5 in [13]. However, it was noted in [19] that for θ > 10 4 the average density starts to deviate dramatically from the estimates made at lower values, ending up to be more than 100 standard deviations higher than the expected value for θ ∼ 10 5 , and seemingly growing as a pure power law for θ > 10 4 . Assuming that the power law behaviour holds asymptotically in θ, it was estimated in [13] that the tree density would reach the percolation value (ρ p ≃ 0.5927 . . .) for θ ∼ 10 40 .
To avoid finite-size effects, in [19] the values of θ and L were chosen as follows: for every value of θ, several simulations were carried out at different system sizes, and it was verified that the distribution of the rms radius, cluster sizes, and burning time were independent of L. The distribution of densities, on the other hand, depends strongly on the system size. For a fixed value of θ, increasing the system size would make the variance of the distribution decrease as L −2 in the absence of finite size effects [16]. Conversely, fixing L and increasing θ would lead to an increase in the standard deviation and ρ , as can be seen in Fig. 1. In the following, we will focus our attention on the behaviour of P (ρ) for increasing system sizes and fixed values of k.
In Fig. 2 and Fig. 3 we plot the distribution of densities rescaling the x-axis by θ ν , and tune ν to align the peaks of the distributions computed at different values of θ. Interestingly, in both cases the curves seem to collapse on the same shape for quite a wide range of θ, suggesting that the distribution has reached a stationary state and that increasing θ would imply only a shift of ρ . However, increasing θ even further clearly shows that the asymptotic behaviour suggested by the data collapse is only transient, as the curves start to deviate considerably from the shape observed at lower values of θ, as can be seen in Fig. 4. Rescaling ρ by θ 0.005 one can align the position of the peaks in Fig. 4 for θ > 10 4 , but it is clear that, in order to perform a data collapse, the y-axis must be rescaled as well, as the height of the peaks decreases with θ.
For a fixed θ, we found that P (ρ) is well fitted by a Gaussian distribution. In Fig. 5, we plot the mean and standard deviations of P (ρ) for k = 10 −3 and find that the average density seemingly increases as a power-law for θ > ∼ 10 4 , in agreement with the numerics presented in [19]. Although the average µ(θ) seems to enter an asymptotic regime after θ > ∼ 10 4 , it is less convincing whether the standard deviation σ(θ) has reached its asymptotic form.
If we assume a Gaussian behaviour for P (ρ) and that the mean and the standard deviation behave as µ = aθ νµ and σ = bθ νσ .
Using these estimates, we find that a data collapse seems to hold well for very large values of θ, as can be seen in Fig. 6. However, such extrapolations have to be taken with great care, firstly because of the small data available for the fit and the not so convincing behaviour of σ(θ), and secondly because the average density can never exceed 1, meaning that the apparent power-law growth has to stop at some point.
Since ν σ > ν µ , the fluctuations grow at a faster rate than the average, as one would expect for a system close to criticality. From the ratio σ(θ) µ(θ) ∼ 0.180 · θ 0.094 , one can obtain a crude estimate for when the fluctuations would become of the same order of magnitude as the average. This would happen at θ ∼ 10 17 and ρ ≃ 0.425. However, this argument can hardly hold, as it would imply that within two standard deviations, we would have values of the densities that exceed ρ = 1, which is of course, impossible. Therefore, we can conclude that for k = 10 −3 and within the Gaussian approximation, we can expect the standard deviation to be smaller than the average, despite growing at a faster rate.
As a reference, we computed a rough estimate of the value θ * that would correspond to an average density equal to the percolation threshold ρ(θ * ) = ρ p . This gives θ * ∼ 10 37 , which is consistent with the estimate made in [19] once the confidence bounds on b and ν µ are taken into account. However, there is no reason why the asymptotic density of the FFM should be the same as the percolation, as even for very large systems the clusters would still be correlated and therefore intrinsically different from a percolation process.
Like previous studies of the FFM, we are unable to settle the true asymptotic scaling behaviour of the model, which is still unreachable with today's computers. However, given that effective scaling behaviour in terms of highquality finite-size collapses can be obtained for very large system sizes and significant ranges of θ, we do think that a detailed phenomenological understanding of the model can be useful as a reference point for the discussion of real systems such as precipitation and brain dynamics. Given this, we present in the next section a careful numerical study of the effective onset region similar to the regions where the precipitation study [5] and the brain study [4] found a broad peak in the residence times.

The most frequently visited region
Even though the studies of precipitation [5] and brain activity [4] related their findings to critical behaviour, both observed a remarkable broad onset of the order parameter, which is certainly not the behaviour seen as one approaches the critical point of a second-order phase transition of an infinite system. The absence of a sharp onset of the order parameter is remarkable because, for system sizes like the earth atmosphere or the human brain, one would not expect any significant finite-size effect if the usual phenomenology of equilibrium phase transitions were any guidance.
While the activation dynamics characteristic of precipitation and brain phenomena are not at all similar to thermal equilibrium dynamics, both systems are at least at a schematic qualitative level similar to the dissipative dynamics of the FFM, with its cycles of loading and discharging. Here we want to investigate further the relationship between the distribution of residence times and the control parameter suggested in [8] and to determine to what extent the FFM exhibits an onset region similar to those observed in rain and brain. If that is the case, we may perhaps take this as indicative of a kind of "universality" different from the stringent universality definition we know from equilibrium systems and more of a pragmatic nature. Needless to say, could such a universality be established, it may be a great help in attempts to classify the behaviour of activation dynamics in complex systems. Furthermore, it could be taken as indicative of a level of emergent behaviour that is independent of the microscopic details, since rain and brain clearly do operate on totally different substrates.
We are of course just repeating the original hope of SOC research and suggesting that systems of entirely different nature may indeed exhibit similar emergent collective behaviour, even if the dynamics does not operate in a critical state, but rather is found to inhabit a broad region of approximately scale-free nature. Crucial for our suggestion is the observation that the true asymptotic behaviour of the FFM model appears to happen for such enormous systems sizes that they hardly are of relevance to real macroscopic systems. In contrast, the quasi-scale free behaviour observed in the FFM model for intermediate system sizes [16,19] may very well be helpful for the understanding of observations such as those presented in [4,5].
To investigate the presence of a critical-like region in the FFM, we focus on the onset of measures that characterise the clusters of trees and keep track of the frequencies at which the system visits different regions of the parameter space. The order parameter for precipitation [5] was taken to be the precipitation rate and for the brain [4], the size of the largest cluster of activated voxels in the fMRI scans -a choice that appears very natural in the light of ordinary percolation analysis. Fig. 7 and Fig. 8 we present the contour plot of the bivariate histograms of the average cluster size [9] S vs the density ρ, and of the largest cluster normalised to the number of trees S max vs ρ. The colour map represents the probability of observing a certain point in the parameter space, which is the same as the proportion of time spent by the system at that point. To create the heatmaps, we sampled 10 7 configurations and grouped points with similar probabilities for better visual representation. Therefore, the histograms are not perfectly normalised. The study of both the average cluster size S and the normalised largest cluster S max is inspired by the resemblance between the clusters of sites occupied by trees and the ordinary geometrical percolation transition. In percolation, either S and S max are used to construct an order parameter.

In figures
It is clear from the heat-maps that the density around ρ * ≃ 0.4 stands out and that precisely like in the precipitation study [5], and even more so in the brain study [4], ρ * is indicative of an onset of the order parameter. Furthermore, the region over which the creation of large events happens is very broad, and the system spends a significant fraction of time in the critical region. We know from the Gaussianity of P (ρ), see Fig. 6, that this region stays broad for any reachable system size and values of θ. This very much suggests that the FFM model's quasi-tuning to a critical region is a stylistic feature of important relevance.

The distribution of the largest cluster
In the spirit of percolation, we now turn our attention to the non-normalized largest cluster, which we indicate with Λ to distinguish it from S max . For a fixed system size and θ, we find that the distribution of the largest clusters P (Λ) is very well fitted by a Fréchet distribution, which corresponds to the class of extreme value statistics with a power law tail. The good agreement of P (Λ) with a Fréchet distribution implies that the correlations between consecutive configurations generating the clusters are sufficiently weak to be ignored.
Although the distribution of cluster sizes does not obey simple scaling [16,19], it was found in [10] that the distribution of instantaneous correlation length is scale-invariant for k = 10 −3 and system sizes at least as big as θ = 9000. Given that P (λ) is fat-tailed, we now turn our attention to its scaling properties and see if simple scaling applies.
In order to measure Λ, one has to maintain and keep track of all clusters at all times, which makes the task more computationally intensive than just sampling the density ρ. For this reason, we analysed systems sizes that are smaller than those used for the analysis of P (ρ). In our simulations we used k = 10 −3 and values of θ up to 4000, but we also checked the scaling for k = 5 · 10 −4 .
Assuming that simple scaling holds, we can expect P (Λ) to follow for Λ 0 << Λ << Λ c , where Λ 0 is a constant lower cutoff and Λ c is an upper cutoff that diverges with θ. In Eq. 3, a is a constant metric factor, τ is the critical exponent, and G is a universal function that decays quickly for Λ >> Λ c . The form that is usually assumed for the upper cutoff is Λ c = bθ ν , where b is another constant metric factor and ν is the spatial dimensionality of the observable Λ [3].
If the data are consistent with the simple scaling ansatz, then it is possible to perform a data collapse by plotting P (Λ)Λ τ vs the rescaled variable Λ Λc for different values of θ. In Fig. 9, we performed a data collapse using ν = 0.055 and τ = 1.04. We also estimated the critical exponents analysing the first two moments of Eq. 3 and fitting the data with < Λ >= c 1 θ α1 and < Λ 2 >= c 2 θ α2 , obtaining c 1 = 10.8 ± 0.8 and α 1 = 1.04 ± 0.01 and c 2 = 111 ± 37 and α 2 = 2.11 ± 0.04 using 95% confidence bounds.
These estimates for ν and τ are consistent with the ones obtained via data collapse and that have been used in Fig. 9. Finally, we applied the whole procedure a second time using a smaller value of k, namely k = 5 · 10 −4 and θ up to 2000. From the data collapse in Fig. 10 we obtained an estimate of ν = 0.065 and τ = 1.04, while fitting the first two moments we obtained: Although we had to restrict our simulations to values of θ < 10 4 , we observe a very robust scaling for the largest cluster Λ, which is not observed in the distribution of cluster sizes P (S) over the same range of θ [16]. Interestingly, we found that Λ ∼ θ 1.04 , a super-linear growth which indicates that the correlations in the system increase rapidly with θ. This behaviour is consistent with the analysis of the distribution of instantaneous correlation lengths P (ξ) performed in [10], where it was found that the average instantaneous correlation length grows as ξ ∼ θ 0.56 over the same range of θ.

DISCUSSION
SOC was very much inspired by the successes of the renormalisation group studies of equilibrium critical phenomena of the 1970-ties and its phenomenal understanding of the origin of universality classes. Initially, the discussions of SOC focused on accurately establishing the scaling behaviour and related scaling exponents of the various models thought to represent the SOC phenomenology. The original sandpile model [1] was disappointingly far from simple scaling, and also the FFM model turned out to behave very differently from the familiar scaling of equilibrium models such as the Ising model or geometrical percolation. Though there is at least one class of SOC models that exhibits clear scaling, namely the class represented by the Manna model [20], the beauty of strict scaling and universality classes defined by scaling exponents seems not to really capture the less than ideal critical behaviour frequently observed in real systems, such as our two examples from rain and brain.
There is no doubt that the studies of the emergent dynamics of real complex systems from biology, geophysics, astrophysics, economics and more [21][22][23] keep identifying behaviour which is qualitatively in the spirit of the hopes and dreams behind SOC, namely lack of one characteristic scale in time and space, large fluctuations and no need for specific external tuning. So if the beauty of strict scaling and exact power laws does not carry over from equilibrium critical phenomena, the question is how we establish a systematic classification of the emergent phenomena observed in completely different systems. The study of the FFM model and its comparison with the behaviour of real systems suggests it is possible to establish a useful phenomenological understanding and classification reaching beyond the usual strict classification of universality classes defined in terms of shared scaling exponents.
The study presented here confirms earlier investigations of the behaviour of the density ρ, and the change in the behaviour of ρ for very large system sizes. Although we observed that the Gaussian behaviour of P (ρ) holds at least until θ = 10 5 , we also showed that it is possible to obtain very good but deceitful data collapses for P (ρ) at different ranges of θ. If, on the one hand, this should serve as a warning for the analysis of the asymptotic scaling behaviour of the FFM and other out of equilibrium systems, on the other hand, it shows that such effective scaling could be of guidance in the understanding and analysis of real systems that show similar dynamics.
Analysing the behaviour of the average cluster size and the largest cluster, we found that the FFM exhibits very similar behaviour to the one observed in experimental studies of rain and brain [4,5]. In particular, both studies and the FFM display a critical region over which both the residence times distribution and the order parameter peak up, meaning that the system spends most of the time in a highly susceptible state.
From the study of the largest cluster Λ it emerges that, although the FFM displays broken scaling in the distribution of cluster sizes, the distribution P (Λ) appears to be scale-invariant at least for θ < 10 4 . Furthermore, P (Λ) displays a super-linear growth of the first moment similar to the one observed for the average instantaneous correlation length in [10]. It is clear from the anomalous behaviour of the density ρ that such results should be taken with care, as the model seems to enter a new regime when θ > 10 4 . However, the robust scaling observed both in the largest cluster and in the instantaneous correlation length suggests that at least for system sizes below θ = 10 4 there is a fast and scale-free growth of the correlations with the activity θ.
Although the FFM does not display the reassuring scaling observed in equilibrium models, its phenomenology appears to summarise elegantly and robustly the emergent dynamics of spreading and recharging seen in such disparate phenomena as rain and brain. Moreover, examples of broken scaling and non-exact powerlaws are ubiquitous in nature and in the scientific literature [24] and, for this reason, we believe that the characteristic behaviour of the FFM should be seen as an asset rather than a limitation of the model.

CONFLICT OF INTEREST STATEMENT
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.