Critical behaviour of the Drossel-Schwabl forest fire model

We present high statistics Monte Carlo results for the Drossel-Schwabl forest fire model in two dimensions. They extend to much larger lattices (up to 65 536 × 65 536) than previous simulations and are much closer to the critical point (up to θ≡p/f = 256 000). They are incompatible with all previous conjectures for the (extrapolated) critical behaviour, although in general they agree well with previous simulations wherever they can be directly compared. Instead, they suggest that scaling laws observed in previous simulations are spurious, and that the density ρ of trees in the critical state was grossly underestimated. While previous simulations gave ρ≈0.408, we conjecture that ρ is actually equal to the critical threshold pc = 0.592... for site percolation in d = 2. This is still, however, far from the densities reachable with present day computers, and we estimate that we would need many orders of magnitude higher CPU times and storage capacities to reach the true critical behaviour - which might or might not be that of ordinary percolation.


Introduction
Empirical analyses suggest that power laws with anomalous exponents (i.e. exponents not obtainable from naive dimensional or mean-field arguments) are ubiquitous in nature, ranging from 1/f noise and earthquake distributions to fractal coast lines, species extinction rates, weather records and the statistics of DNA [1]- [4]. On the other hand, it is well known that such scaling laws-most clearly seen by linear relationships in log-log plots-can be spurious. Log-log plots have a notorious tendency to suggest linear curves, even if there are no real power laws.
It would thus be extremely useful if these empirical observations could be backed by theoretical models where power laws can be either proven exactly or at least verified beyond doubt by high statistics simulations. Unfortunately, equilibrium systems in general show anomalous 17.2 scaling laws only at critical points which are codimension one phenomena: one has to fine tune some parameter (e.g. temperature) to reach them, otherwise no power laws are obtained. Thus, they cannot be used to justify why such power laws should be seen in nature.
A possible solution to this puzzle was indicated in [2,5] where it was suggested that many non-equilibrium systems could be driven by their dynamics into a critical state. The main ingredients of this self-organized criticality (SOC) is slow driving towards some instability and a mechanism to relax the tensions built up by the drive locally and partially. Since the tensions are not relaxed completely (in many models they are just redistributed), the state becomes marginally stable and apt to relaxation events on wider and wider length scales, which then lead to anomalous scaling laws. The paradigmatic model is the 'sand pile' of [5] which does not describe real sand piles but which was proven exactly to show anomalous scaling laws at least for some of its observables [6]- [8].
Another model which was proposed to show SOC is the forest fire model introduced independently by Henley [9] and Drossel and Schwabl [10]. In this lattice model with discrete time each site can be either occupied by a tree or empty. New trees are grown at a small fixed rate p on empty sites, and at a rate f p sites are hit by lightning strokes which then burn the entire cluster of trees connected to this site. Burning happens infinitely fast, so the only relevant parameter is the ratio θ = p/f which also sets the scale for the average fire size (the number of trees burnt after one lightning stroke). Criticality is observed in the limit θ → ∞.
The Drossel-Schwabl model (called the DS model in the following) is different from other SOC models in two ways.
• It involves not only the separation between two time scales (the slow build-up of stress and the fast relaxation), but involves three time scales: the fast burning of connected clusters of trees, the slow re-growth, and the even slower rate of lightning strokes. • The growth of trees does not lead to a state which is inherently unstable (as the addition of sand grains to the top of a pile does), but only to a state susceptible to being burnt. Without lightning, the tree density in any patch of forest can go far beyond criticality. When the lightning strikes finally, the surviving trees have a density far above critical.
Indeed, it was observed in [11] that the stationary steady state of the model is not 'critical' in most regions, in the sense that its tree density is not marginal for the spreading of fire. Rather, it is composed of large patches of roughly uniform density, most of which are either far below or far above the critical density for spreading. Nevertheless, power laws were observed for several observables, partly because these patches occur at all sizes, so that the fires also had a broad spectrum of sizes.
While normal scaling with mean-field exponents had been seen in [10], all subsequent simulations [9,11,12,13,14] showed clear signs of anomalous scaling: • The fire size distribution scaled, for small s (s is the number of trees burnt in one fire) and in the limit θ → ∞, as P (s) ∼ s 1−τ with τ ≈ 2.15. • For finite θ, P (s) is essentially cut off at s max ∼ θ λ with λ ≈ 1.1.
• The average rms radius R(θ) = R 2 1/2 of all fires scaled as θ ν with ν ≈ 0.58. Here, R 2 is the Euclidean distance of a burning tree from the site of lightning, and the average is taken over all trees in all fires at fixed θ. • The rms radius R(s, θ) of fires of fixed size s scaled as s 1/D where the fractal dimension D of fires is D ≈ 1.96 [13] to 2.0 [11].

17.3
Finally, the average density of trees was found to be ρ = 0.408 − const/ √ θ [11,14]. There were seen already at that time, however, large corrections to this scaling law and deviations from conventional ansatzes. Thus, • The determination of τ was subject to large systematic uncertainties [11]; • The scaling R ∼ s 1/D was observed in [11] only for fires of typical size (s ∼ θ). For large s, R was significantly larger; • For finite θ, the fire distribution P (s) did not follow the normal scaling ansatz P = s 1−τ φ(s/s max ) [11]; • There are (at least) two differently divergent length scales [14]: the correlation length evaluated from all pairs of sites scales differently with θ than R(θ) 1 ; • Finite size behaviour is abnormal [15].
In two recent publications, these problems were taken up again. In [16] it was claimed that they are due to non-leading corrections to scaling. In [17] a more radical solution was proposed with two distinct classes of fires which both contribute to the scaling limit. In the latter paper a connection to ordinary percolation was also proposed, and the conclusion was backed by a 'coarse grained' model which supposedly could mimic the DS model at extremely large θ and on extremely large lattices.
It is the purpose of this paper to present very large simulations which show quite unambiguously that none of these really describe the true critical behaviour of the DS model. While our simulations agree perfectly with previous ones for the lattice sizes and θ values used there, we shall see that the supposed scaling laws are just transients. Even the present simulations do not reach the true asymptotic regime, but some suggestive features of the true asymptotics do emerge.
We describe our simulations and their results in the next section. In the last section we draw our conclusions.

The simulations
Our simulations are straightforward, with a few subtleties. They follow [11] in taking p → 0 by making no attempts to grow new trees while a cluster burns. As in [11] we made exactly θ growth attempts on randomly chosen sites between two successive lightnings, instead of letting this number fluctuate. In the large L limit this should not make any difference. But in contrast to [11], where a depth-first algorithm had been used, we used a breadth-first algorithm to burn the cluster 2 .
In order to simulate very large lattices with minimal storage, we use multi-spin coding, i.e. we use only one bit to store the status of a site. In this way we could simulate lattices of size 1 The length ξ c defined in [14] is closely related to R(θ). In the notation of [14] (with just a minimal change) it is given by is the occupation of site x), while R 2 (θ) = y∈Z 2 y 2 K(y)/ y∈Z 2 K(y). 2 In a depth-first approach, trees to be burnt are put onto a stack, i.e. the neighbours of the last ignited trees are handled first. In a breadth-first approach, items are put into a first-in first-out queue, so that distant neighbours of the site of lightning are burnt only after all closer neighbours have already been handled. The burning times quoted in [11] were obtained from the stack heights in a recursive (depth-first) implementation. They do not agree with the more realistic burning times obtained with breadth-first algorithms and should thus be considered as obsolete.

17.4
L × L, L = 65 536 on computers with 1 GB main memory. Notice that we do not need to store for every tree whether it is burning or not, since the burning trees are stored in a separate list which is updated at each time step. Boundary conditions were helical, i.e. sites are indexed by one scalar index i, with neighbours having indices i ± 1 and i ± L, and with i + L 2 ≡ i. The largest previous simulations [16] had used L = 19 000.
We were careful to discard sufficiently long transients (between 1.6 × 10 6 and 1.2 × 10 7 lightnings; this is up to one order of magnitude longer than those in [14,16]) before taking statistics. This is needed since excessively large fires occur during these transients, and thus the tail of the distribution P (s) is heavily influenced by them. We believe that previous analyses were affected by this problem, which is easily overlooked since bulk properties (such as ρ or P (s) for typical s) show much faster convergence. The total number of fires in each run used for averaging was between 10 9 (for θ ≤ 250) and 9.3 × 10 6 (for θ = 256 000; see table 1). Previous authors [16] went only to θ = 32 768 with 2.5 × 10 7 fires. Compared to that, our statistics is larger by roughly one order of magnitude. All simulations were done on fast Alpha workstations. The CPU times per run varied between 15 h and 4 weeks. As a random number generator we used Ziff's four-tap generator with period 2 9689 − 1 [18].
We tested for finite size effects by making several runs for the same θ on lattices of different sizes. Previously, finite size effects had been studied systematically in [15], but only for much smaller lattices (L ≤ 2000; the authors of [16] called their analysis a finite size analysis, but they actually made a conventional scaling analysis without checking the finite size behaviour). For θ = 64 000, for example, we made runs with L = 2 14 , 2 15 and 2 16 . We verified that distributions like P (s), R(s, θ) or P (t) (t is the burning time of a fire) were independent of L within the statistical accuracy, i.e. any systematic L dependence was masked by statistical fluctuations.
Systematic dependencies were seen only for averaged quantities like s , ρ or R(θ). Indeed, ρ can be measured either immediately before or immediately after a lightning. Since each lightning burns on average s trees (here s is averaged over all lightnings, whether they hit a tree or an empty site), we have ρ before − ρ after = s /L 2 . We found that ρ before decreased with L, while ρ after increased with it (see table 1). More precisely, we have Notice that s is given for L → ∞ exactly by s = (1 − ρ)θ [10]. For finite L, one has to replace just ρ by (ρ before + ρ after )/2, to leading order in 1/L. This was verified in the simulations, but it just tests stationarity and the absence of gross mistakes. The rms radius of fires depended on L and θ in a more complicated way, see table 1. The data were less clear in that case, but they could be fitted by Results for the fire size distributions P (s) are shown in figure 1. We see the approximate power decay s 1−τ for s < s max with s max roughly proportional to θ. But we also see the strong deviations from this power law first observed in [11]. Due to these deviations, P (s) decreases with θ in the scaling region (making the effective value of τ increase with θ), but has a growing bump at s ≈ s max . In [17] it was conjectured that the asymptotic value of τ , estimated from the scaling region 1 s θ in the limit θ → ∞, is τ = 2.45. This was based on heuristics and on simulations for small θ on very small lattices (L = 1300). In order to test it, we plotted P (s) in figure 1(b) after multiplying it by a suitable power of s. It is seen that s 1.19 P (s) becomes flat in a wider and wider region of s as θ increases. For very small fires (s < 300) P (s) decreases faster than s −1. 19 , but there is no indication that the slope increases with θ for Table 1. Statistics and main results: N is the number of lightnings used for averaging, N trans that of lightnings discarded during the transients. The density is measured after fires became extinct and before new trees are grown. Notice that the density has much larger errors in some runs than in others of comparable statistics. This results from the fact that there are important longranged negative autocorrelations in the density time series, and I had not written out the information needed to take them into account in these runs.  θ > 8000. Based on this evidence we would thus conclude that τ = 2.19 ± 0.01, thereby ruling out the value of [17]-provided, of course, that there is no change of behaviour as θ becomes even larger.
In [11] it was conjectured that the bump near s ≈ s max is due to the cut-off. If the integrated distribution P int (s) = ∞ s =s P (s ) were just a power multiplied by a sharp cut-off, its derivative P (s) would have a bump where the cut-off sets in. This bump would consist of those events which would have been in the tail which is cut off. If this were right and P int (s) would indeed show normal scaling, we should expect the height of the bump in figure 1(b) to be independent of θ. But this is obviously not the case. Instead, it increases with θ, suggesting a different scaling law s 1−τ for the envelope of the curves in figure 1(a), with τ = 2.111 ± 0.006. A similar conclusion is reached by looking directly at P int (s). Indeed, log-log plots of P int (s) versus s are much straighter in the scaling region ( figure 2(a)). But again a blow-up after multiplication by a suitable power of s shows that this is misleading ( figure 2(b)). Even P int (s) is not convex for θ > 1000, and it develops an increasingly sharp shoulder near s = s max as θ → ∞. We can try to obtain an alternative estimate of τ by fitting straight lines such that they touch both maxima in figure 2(b). Results of this are shown in figure 3. In contrast to the previous estimate of τ they would indicate that the effective τ decreases with θ and is clearly less than 2.19 for θ → ∞ (also the value 2.159 of [14] seems hardly compatible with the extrapolation of figure 3 to θ → ∞). But a convergence to 2.11 seems quite possible, as we would expect if the bumps in figure 1(b) have a width independent of θ (which is not excluded by the data).
This obviously means that these scaling violations are not corrections which disappear in the limit θ → ∞, as was claimed in [16]. It rather indicates that the conventional scaling picture, is basically wrong, as was already conjectured in [11,17]. Previous analyses indicated that s max actually increased faster than θ, roughly as with λ ≈ 1.08. We verified this qualitatively, but verified also the finding of [11] that s max and thus also λ are not well defined since the sharpness of the cut-off increases with θ. This is already obvious from figure 2(b), but it persists for larger values of s not shown in this plot. It can also be seen by making copies of figure 2(a) on transparencies and overlaying them. Better scaling than for P (s) was seen in [11] for R(θ). Our present data are fully compatible with those of [11,13,14] 3 , but involve much higher statistics and cover a wider range of θ. Thus it might not be too surprising that we now do not see perfect scaling any more. But the observed deviations from a pure power law (see figure 4) are much larger than expected from subasymptotic corrections 4 . They clearly show that the previously seen power law was spurious and does not describe the asymptotic behaviour. A power law fit through the last two points would give R(θ) ∼ θ 0.563 , but this is obviously not yet the asymptotic behaviour. Again, as for P (s) and P int (s), we cannot yet say what the correct asymptotic behaviour will be. In any case, the claim of [14] that there are two different diverging correlation lengths becomes obsolete.
An indication of the origin of all these puzzles comes from looking at R(s, θ). This is defined analogously to R(θ), except that the averaging is now done over all fires of fixed size s. We expect R(s, θ) ∼ s 1/D if fires are fractal with dimension D. Previous analyses had given hit either very small regions with supercritical tree density which burn off completely, or regions of subcritical density on which the fires form subcritical percolation clusters. Both mechanisms give compact clusters with fuzzy boundaries. In any case, region I becomes less and less important as θ increases. • Region II. This is roughly equal to the scaling region in figure 1(b). Here R(s) is nearly proportional to √ s, i.e. D is very close to 2-but not quite. Also, we see a clear decrease of the minimal slope in figure 5 with θ, but it seems not to be sufficient to give D = 2 in the limit θ → ∞ (see figure 6). According to figure 6 the critical exponent ν defined by R(s, θ) ∼ s ν in region II seems to converge to 0.505 for θ → ∞, corresponding to D = 1.98. This is halfway between the best previous estimates D = 1.96 [13] and D = 2.0 [11]. It indicates that clusters in region II are fractal, but more compact than critical percolation clusters (D = 1.89). • Region III: In the region of very large fires, corresponding to the increasing bumps in figure 1(b), the apparent fractal dimension decreases again. Unfortunately, due to the strong curvatures of the curves in figure 5 in this region, we cannot quote any definite dimension value. But by plotting R(s, θ)/s ν with suitable values of ν on log-log scales, we see that (i) the maximal values of ν increase slowly with θ (except for the very largest θ, where statistics is poor); and (ii) for our largest θ and s, we have ν ≈ 0.64. Thus the largest fires are definitely more fractal than critical percolation clusters! In terms of the scenario with roughly homogeneous patches with constant tree density [11,13,17], fires in region II correspond to single patches of typical size which are hit by lightning just when their density has reached about the critical percolation threshold. If lightning would always strike exactly at criticality, and if all trees in the patch would have burnt during the last fire so that all trees now are placed randomly, this would give D = 1.89. Fires in region II are more compact, mainly because it will take some time until a lightning strikes by chance the region after it has reached the critical density.  (1), but for the largest L the extrapolation shifted them by less than the error bars (which are indicated inside the squares). The broken line is a fit of the small θ values (θ ≤ 8000) to an ansatz The larger a patch is, the bigger will be its perimeter and thus also the chance that a fire 'spills over' and also burns a neighbouring patch, and from this also a next patch, etc. Since this will happen only along parts of the perimeter, the resulting fires will be rather fuzzy, with effective dimension <2. We propose that such fires dominate in region III. Although they are rather few in number, they are very important since they burn large parts of the entire lattice, and they lead to rearrangements of the global pattern of patches. Notice that fires in region III burn only small parts of the entire region they cover, leaving behind more unburnt trees than fires of type II. Since they dominate more and more as θ increases, one might suspect that this leads to an increase of ρ with θ.
Indeed, a slight increase of ρ with θ had been seen in all previous simulations. The best previous estimates were ρ = 0.4075 − const/ √ θ [11] (unfortunately the constant multiplying θ −1/2 was estimated wrongly in [11] due to a simple mistake) and ρ = (0.4084 ± 0.0005) − const/θ 0.47 [14,16]. The results of our present simulations are shown in figure 7. For θ < 10 4 we see a perfect agreement with previous results, but for larger θ there is dramatic disagreement: our measured values are higher than predicted by extrapolation from small θ by up to 100 standard deviations. It is not clear why this was missed in [16]. But we might mention that no data for θ > 10000 are shown in figure 1 of [16], although the authors claimed to have made high statistics measurements up to θ = 32 768.
When plotted against log θ, our values of ρ follow roughly a straight line for θ ≥ 4000. We should, of course, not take this as the asymptotic behaviour, since ρ can never increase beyond p c = 0.5927 . . . , which is the critical value for site percolation in two dimensions (if at any time the density is larger than p c , this would lead to very compact large fires and thus to a very low density afterwards). But we can use it to obtain a very crude estimate of the order of magnitude when ρ ≈ p c should be reached. It is θ ≈ 10 40 . Although any prediction based on such a large extrapolation should be taken with great care, we conjecture that, indeed, ρ converges to p c for θ → ∞. The main reason is that we see no other plausible scenario compatible with our present numerics. It is not clear whether fires in this limit correspond to critical percolation clusters, since we must expect that weak correlations in the tree densities survive in this limit, sufficiently so as to spoil any agreement with uncorrelated percolation on large scales. In any case, ρ → p c is only possible if large fires are very fuzzy and burn only small fractions of trees in the affected areas, because only in this way can strong correlations be avoided which lead to the patchy tree distributions seen at presently reachable θ and to ρ < p c .
A last hint in favour of our claim that large fires asymptotically are dominantly associated with regions of critical tree densities (i.e. ρ ≈ p c ), even though they might not form critical percolation clusters, is obtained by studying the mean sizes of fires which started in regions of given tree density. If a lightning hits a tree at site i, we define the local tree density ρ loc at this site as the number of other trees in the surrounding square of size 9 × 9 divided by 80. This size is, of course, rather arbitrary, but we obtained qualitatively similar results for squares of sizes 7 × 7 and 5 × 5. In figure 8 we plot the average fire size s(ρ loc ), divided by the overall mean s , against ρ loc . For small θ we see a monotonic increase which is easily understood: since even the largest patches of uniform density are not much larger than the square, large fires can only result from regions with high local density. This is no longer true for large θ. There, patches with very high density are probably small, otherwise they would already have been burnt down earlier. At values of θ reached in this work the largest s(ρ loc ) are still for large ρ loc , but a pronounced peak at ρ ≈ p c develops where s(ρ loc ) has a local maximum. Since also the number of sites with ρ loc ≈ p c is much larger than those with ρ loc p c (see the next paragraph), we see that it is fires in regions with critical density which play an increasingly important role as θ → ∞. In any case, the very strong dependence of s(ρ loc ) on θ shows that we are still far from the asymptotic region where we expect this dependence to disappear. Finally, we show in figure 9 the distribution of local densities itself. It depends rather weakly on θ (less than s(ρ loc ), at least), but the precise way it does depend on θ is rather surprising and not yet fully understood. First of all, it develops an increasingly sharp maximum which slowly shifts from ρ loc ≈ 0.4 to ρ loc ≈ 0.6 as θ increases. This is to be expected after our previous observations. What is unexpected and hard to explain is a shoulder at ρ loc ≈ 0.8 which develops for the largest θ values. It is not very large but statistically highly significant (it was seen in all runs with θ ≥ 32 000). There is also a shoulder at small ρ loc (≈0.1) which also increases with θ. Presumably it is related to the shoulder at ρ loc ≈ 0.8: if patches with density >0.7 burn, they leave behind extremely strongly depleted patches. One possible reason why such high density patches can survive at all is that many large fires are fractal and leave behind small disconnected regions of fairly high but subcritical density. These regions then are too small to have a chance to be hit by lightning until their density has grown far beyond p c .

Discussion
The simulations reported in this paper leave little doubt that all scenarios proposed so far for the Drossel-Schwabl forest fire model do not describe the asymptotic behaviour for θ → ∞, where the model should show SOC according to the standard folklore. Indeed, we do not see many indications for any power laws in this model, as all proposed scaling laws seem to be just transient. There are a number of observables which do show straight lines on log-log plots (such as figure 1(a) in the central region of s or the envelope to figure 1(a)), but it seems more likely that these are also spurious.
This situation is, of course, not altogether new. There are a number of models which were supposed to show anomalous scaling, but closer scrutiny proved otherwise. A good example is the Bak-Chen-Tang forest fire model [19] which at least in d = 2 is not critical [20]. Other examples include the Newman-Sneppen model [21] where one can prove semi-numerically that the observed power laws are transient [24] and maybe even the 'classical' Abelian sand pile model in d = 2. While power laws for waves were proven rigorously in that model, it might well be that the observed deviations from finite size scaling [22,23] do not herald multifractality but just simply no scaling at all. One indication for the latter is the fact that some scaling laws show violations which do not seem to vanish for increasing system sizes [8]. Also, some other quantities in the sand pile model which involve superpositions of many waves depend qualitatively on the geometry of the lattice (square versus strip) [24]. For a system with true scaling one would not expect this.
The situation becomes even worse when going to real life phenomena. It does not seem unlikely that many of the observed scaling laws are just artefacts or transients. Problems of spurious scaling in models which can be simulated with very high precision, such as the present one, should be warnings that not every power law supposedly seen in nature is a real power law.