Table Of ContentBroken scaling in the Forest Fire Model
Gunnar Pruessner and Henrik Jeldtoft Jensen
Department of Mathematics, Imperial College, 180 Queen’s Gate, London SW7 2BZ, UK
[email protected] and [email protected]
(Dated: February 1, 2008)
2
0 We investigate the scaling behavior of the cluster size distribution in theDrossel-Schwabl Forest
0 Firemodel(DS-FFM)bymeansoflarge scale numericalsimulations, partlyon(massively) parallel
2 machines. Itturnsoutthatsimplescalingisclearlyviolated,asalreadypointedoutbyGrassberger
n [P. Grassberger, J. Phys. A: Math. Gen. 26, 2081 (1993)], but largely ignored in the literature.
a Most surprisingly the statistics not seems to be described by a universal scaling function, and the
J scaleofthephysicallyrelevantregion seemstobeaconstant. Ourresultsstrongly suggest thatthe
7 DS-FFMis not critical in the sense of being free of characteristic scales.
1
] 1. INTRODUCTION the other. If its state is “empty” turn it into “oc-
h
cupied”.
c
e The Drossel-SchwablForest Fire Model [1] (DS-FFM)
• Relaxation: Choose one site at random. If it is
m is one of the paradigms of non-conservative SOC [2]. Its
emptycontinuewiththefirststep. Otherwisemake
- importancecomesprimarilyfromthefactthatthemodel “empty” all the sites in the cluster the site chosen
t has non-conservative microdynamics. It therefore an-
a belongsto. Inthis casethe update is consideredto
t swered the question whether conservation is necessary
s be successful. Continue with the first step.
. for criticality in driven systems.
t
a The claimthat the DS-FFM is criticalcomesfromthe Here a cluster is defined in the usual fashion as the set
m fact that it shows powerlaw-like behavior for several ge- of occupied sites which are connected via nearest neigh-
ometrical properties of the dissipation events. The aver- borinteractions,i.e. twositesbelongto the samecluster
-
d age size of these is divergent in the so-called SOC limit, if they are nearest neighbors or there is a path between
n where alltimescalesgetseparatedsothat the rate ofthe them along sites, which belong to the same cluster. We
o
external drive becomes infinitely slow and the total in- have applied periodic boundaries in all our simulations
c
[ flow diverges. If one assumes stationarity this is trivial and restrict ourselves to the two dimensional square lat-
to prove [1, 3, 4]. However, as usual in numerical simu- tice.
1 lations, it is not possible to investigate the model in the The cluster removed in the second step is called the
v
limit of divergent drive (θ → 0 in the notation below), “burnt cluster”. To measure the overall distribution of
6
as finite size limits the correlation length and therefore clusterswithinthe system,oneusuallymeasuresthe size
0
3 destroysanypossiblecriticality[4]. Itisremarkablethat of the burnt cluster [4], the distribution of which is bi-
1 most of the literature available for this model is mainly ased by a factor s. To see that, we define n(s) to be the
0 concerned with finding critical exponents and identify- ensemble-averaged,site-normalized density of clusters of
2 ing supposedly critical quantities. It seems that no au- size s in the system. Then, the probability that a ran-
0 thors question whether the model is critical at all and domly chosen site is connected to a cluster of size s is
/
t if so in which sense. In this paper we carefully investi- sn(s),asinstandardpercolation[5]. Thisdistributionis
a
gatethe“scalingfunction”oftheclustersizedistribution probably the most important in the model. Other quan-
m
and show that it is indeed an open question whether the tities are the distribution of the burning time, which is
d- modelis truly critical: Not onlyis there no wayto prove defined as the maximum Manhattan distance (shortest
n its criticality, there is also numerical evidence that the path on the square lattice) from the initially chosen site
o model does not become scale free. of a burnt cluster to all other sites in the same cluster,
c andthe correlationfunctions as defined and discussedin
:
v [6]. In this paper we solely concentrate on the distribu-
i 2. DEFINITION OF THE MODEL AND tion n(s).
X
METHODS Using a new implementation of the model [7] we are
r abletosimulatethesystemonverylargescalesandatthe
a
The model has been described several times and in same time keep track of the entire distribution n(s), in-
greatdetailelsewhere[1,3,4]. Therefore,thedescription steadofmeasuringthebiaseddistributionsn(s),asdone
presented here is rather succinct. The model is defined usually[4]. Betweentwoupdates thechangesinn(s)are
on a d-dimensional lattice of linear size L, where each only of the order (1/θ), so it is a highly correlatedquan-
lattice site can be in one of two states, “occupied” or tity. However, by using standard cluster labeling tech-
“empty”. The lattice is then updated according to the niques [8], it is possible to calculate the full histogram
following rules: n(s) essentially without increasing the computing time,
which depends almost exclusively on θ and is essentially
• Driving: Choose randomly (1/θ) sites, one after independent of the system size. Compared to the stan-
2
dard simulation, we gain up to two orders of magnitude Thatthisparticularchoiceofthe normalizationdoesnot
in performance [19]. A similar method was recently in- affect the overall results can be seen in Tab. I, where
troduced for standard percolation [9]. Using the same the absolute value of n(1,θ) is listed for different values
amount of computing time the results are significantly of θ. Also shown in this table is the first moment of
lessnoisythanthoseofthestandardimplementation(for the distribution or the average density of trees, which is
example [10], which we have used as reference to check defined as
the validity of our results). Large system sizes enable
ρ= sn(s;θ) (3)
us to rule out any finite size effects as described below.
X
s
Theresultshavebeenpartlycross-checkedusingadiffer-
ent randomnumber generator(all results presented here and the second moment of the distribution,
make use of ran2 from [11], and for checking ran1 from s2n(s;θ)
[11] has been used). hsi= s , (4)
P sn(s;θ)
Finite size effects have been ruled out by the follow- s
P
ingdirectmethod: Foreachvalueof(1/θ)asignificantly which is the average size of the cluster connected to a
largersystem wassimulatedwith exactly the same value randomly chosen occupied site.
for (1/θ). The linear size L was typically increased by a Before presenting the actual results, we first discuss
factor 2. The smallest systems we used were L = 1000, the numerical quality of the results.
the largest L = 32000. By comparing the histograms of
different system sizes in conjunction with the standard
deviation calculated for them, it is possible to decide 3.1. Avoiding finite size effects
whether a system size is affected by finite size effects or
not. Compared to other simulations published [12, 13],
which also claimnot to suffer from finite size effects, our 1.05
a L=4000, L=8000
system sizes are huge. This simple method of “redoing” 1 2
1/θ=2000
all simulations and using lattices which are much larger L)2
thanactuallyneededhastheobviousdisadvantageofbe- L,1 1.00
θ,
inginefficient,butthereisprobablynomoreadirectway s;
r(
of identifying finite size effects [14]. This waste of com-
puting poweris overcompensatedby the efficiency ofthe 0.95
algorithm and self-averaging[15].
b L=1000, L=8000
1.0 1 2
L)2 1/θ=2000
3. RESULTS L,1
θ, 0.8
s;
The focus of this paper is the presumably universal r(
scaling function of the distribution n(s). Similar to the
0.6
correlation function one expects 100 101 102 103 104 105
s
n(s;θ)=s−τG(s/s (θ)) (1)
0
if“simplescaling”applies,whichisalreadyknownnotto FIG. 1: Ratio r(s;θ,L1,L2) = n(s;θ,L1)/n(s;θ,L2) with
be the case in the presence of finite size effects [13]. The (1/θ) = 2000 for two pairs L1,L2 with error-bars (one stan-
e e
dard deviation; the error-bars as well as the data shown are
L-dependenceofthisquantityisomittedinthefollowing
exponentiallybinned). Thedataarefromshortruns(106 up-
wherever the context allows it. It is worthwhile to note
dates for statistics). Finite size effects have been considered
thatthis is usually thedefinition ofthe exponentτ. The
negligibleundertheconditionthat(almost all) error-barsfor
functionG isthe(presumably)universalscalingfunction, this ratio have covered 1 (marked by a dashed line) in the
which depends only on the ratio s/s0(θ), where s0(θ) is relevant range. a) L1 = 4000 and L2 = 8000: Almost no
theonlyscaleofthedistribution. Thisscaledependsonly finite size effects, the deviation from 1 is probably due to
ontheexternalparameters,inourcaseθ. Simplescaling noise. Note the fine scale of the ordinate. b) L1 =1000 and
allows another scale, namely the lower cutoff, but this is L2 =8000: Systematic, strong finite size effects for s ' 104.
fixed or at least bounded. Thescaleoftheordinateisfivetimeslargerthanina). Data
of this quality havebeen dismissed.
The function G(x) is usually smooth for small values
of x, therefore it does not make a big difference whether
we investigate n(s) or Throughout this paper we initially performed 5·106
successfulupdates(asdefinedinSec.2)astransient(and
n(s;θ)
n(s;θ) = therefore rejected them) and the same number for pro-
n(1;θ)
ducing statistics, apart from runs for calculating error-
e = G(s/s0(θ))s−τ . (2) bars, where only 106 updates has been used for statis-
G(1/s (θ)) tics, see below. It is known that the transient can be
0
3
100
TABLE I: Static quantities for different choices of L and L=4000,8000
(1/θ). The estimation of the standard deviation of the tree 10−2 1/θ=2000
density ρ, σ2(ρ) = hρ2i−hρi2, where the average runs over
the ensemble, is unfortunately based only on a small subset 10−4
of the configurations produced, and in the case of the large
systems (L ≥ 16000) only on a fraction of the lattice. How- L) 10−6
1;
ever, it is apparent that it behaves like 1/L, as expected for n(
a system without finite size effects. The density of clusters L)/ 10−8
of size 1, n(1), servesas thenormalization of n. The average s;
cluster size is denoted by hsi, for definition see (4), but due n(10−10
e
to a truncation in the histogram for some of the simulations
in the range 2000 ≤ (1/θ) ≤ 16000, the number presented 10−12 L=1000,8000
1/θ=2000
is actually the average size of the burnt cluster. In the sta-
tionary state it is - apart from statistical fluctuations - also 10−14
given by (1−hρi)/(θhρi) [4]. Values of (1/θ) and L printed 100 101 102 103 104 105 106
in bold indicate results shown in Fig. 3, the otherresults are
s
onlyforcomparison. Alldataarebasedon5·106 (successful)
updates(s. Sec.2)forthetransientandstatistics, apartfrom
thoseprinted in italics which are based on short runs(5·106 FIG. 2: The binned histogram n(s;θ,L) for two different
updatesfor thetransient and 1·106 updates for statistics). valuesofLandfixedθ asinFig.1a. Inthisplotthetwohis-
e
(1/θ) L hρi σ2(ρ) n(1) hsi 1−hρi tograms are virtually indistinguishable. However, note that
θhρi the deviations shown on Fig. 1b would also hardly be visible
p
125 1000 0.3797 0.0060 0.04553 204.07 204.18 in this typeof plot, as shown in the inset.
125 1000 0.3798 0.0058 0.04552 203.81 204.15
125 4000 0.3798 0.0014 0.04553 203.88 204.10
125 4000 0.3798 0.0015 0.04552 203.77 204.10
very long [6] (note that the time unit in [6] is expressed
250 1000 0.3876 0.0083 0.04451 395.03 395.06
in our units by multiplying it with (1/θ)/(ρL2)), but in
250 1000 0.3875 0.0082 0.04452 394.08 395.15
all cases presented the number of initial steps seemed to
250 4000 0.3877 0.0022 0.04454 394.97 394.89 be more than sufficient. Numerical checks indicate that
250 4000 0.3877 0.0021 0.04454 395.29 394.91 theclustersizedistributionisverystableagainstthesize
500 1000 0.3932 0.0117 0.04380 764.73 771.75 of the transient, i.e. even a transient, which is presum-
500 1000 0.3932 0.0119 0.04380 764.81 771.77 ably too short, still produces reasonable results for n(s).
500 4000 0.3934 0.0031 0.04382 771.12 770.88 All systems have been initialized by a random inde-
pendent distribution of trees with density 0.41.
500 4000 0.3934 0.0030 0.04382 771.90 770.87
The standarddeviationofthe binned histogramis not
1000 1000 0.3972 0.0169 0.04328 1495.36 1517.91
completely trivial to calculate. In particular, its compu-
1000 1000 0.3971 0.0168 0.04328 1490.05 1518.00
tation requires a significant amount of CPU time, and
1000 4000 0.3976 0.0043 0.04331 1510.85 1515.00
wasthereforeonlycalculatedforthesmallersystemsizes
1000 4000 0.3976 0.0043 0.04331 1513.13 1514.81 (up to L = 8000) and in shorter runs (only 106 updates
1000 8000 0.3976 0.0021 0.04332 1510.10 1514.91 for statistics, but 5·106 for transient). We resorted to
2000 4000 0.4005 0.0060 0.04296 2976.34 2993.35 visual examination for the larger systems when compar-
2000 4000 0.4005 0.0062 0.04297 2990.50 2993.15 ing n(s;θ,L) for different system sizes. Fig. 1a and b
2000 8000 0.4006 0.0030 0.04297 2995.67 2992.56 showthe ratioofn(s;θ,L)fortwodifferentsystemsizes.
e
4000 4000 0.4026 0.0089 0.04273 5929.24 5935.91 A deviation of this ratio from 1 indicates a difference
e
in the statistics and therefore the presence of finite size
4000 4000 0.4025 0.0089 0.04273 5930.97 5938.03
effects. Fig. 1a shows a typical case we accepted as rea-
4000 8000 0.4026 0.0048 0.04274 5931.32 5935.15
sonable agreement. Here L = 4000 and L = 8000 do
1 2
4000 8000 0.4026 0.0046 0.04273 5935.36 5936.47
notseem to differ for (1/θ)=2000. Fig. 1bshows a case
8000 4000 0.4040 0.0135 0.04255 11786.97 11799.72 of finite size corrections we have dismissed (note the dif-
8000 4000 0.4041 0.0135 0.04255 11788.90 11799.07 ferent scales in the two graphs). It differs from Fig. 1a
8000 8000 0.4041 0.0068 0.04257 11801.31 11795.98 only by L =1000.
1
8000 8000 0.4041 0.0068 0.04257 11792.82 11795.38 Fig.2illustratesthe strongagreementofn(s;θ)atthe
16000 4000 0.4052 0.0199 0.04244 23430.01 23481.82 same value of θ for the same two different sizes L as in
e
16000 8000 0.4054 0.0096 0.04243 23466.93 23467.22 Fig. 1a. The two sets of data are virtually indistinguish-
able, but in this kind of plot it is also almost impossible
16000 8000 0.4054 0.0098 0.04243 23446.10 23465.64
to see a difference between the data of L = 1000 and
16000 16000 0.4054 0.0052 0.04245 23449.31 23466.57 1
L = 8000, as shown in the inset of Fig. 2. This is also
2
32000 16000 0.4066 0.0075 0.04232 46443.83 46701.82
the case with the rescaled data below, and the use of
32000 32000 0.4066 0.0032 0.04233 46731.44 46698.51
very large systems throughout this paper might there-
64000 32000 0.4078 0.0042 0.04220 91148.64 92952.40
4
fore be “overcautious” in avoiding finite size effects, al- 2
thoughsuchlargesizesareobviouslyrequiredforanaccu-
1/θ=250 1000 4000 16000 64000
rate quantitative analysis of this model. However, when 1/θ=125 500 2000 8000 32000
it comes only to qualitative analysis, such a judgment 2.22 smax(θ)
seems to be justified. On the other hand, an increase
in system size hardly increases the computing time and
affects“only”thememoryrequirements,whichforcedus *τs smin(θ)
to implement the algorithm for parallel machines. The 1) 2.19
side effect of using multiple CPUs at the same time is s)/n(
a significant reduction of the simulation time especially n(
for large values of (1/θ), a fact which compensates the
complications of parallel coding. 1 2.14(3)
Another indicator for the absence of finite size effects
isthescalingofthestandarddeviationofρ: Ifthelattice 223/91
can be split into independent parts, i.e. if subsets of the 0.8
lattice can be considered as independent, the standard 100 101 102 103 104 105 106
deviationofρshouldscalelike1/Lfordifferentvaluesof s
L at given (1/θ). Such a behavior can be seen in Tab. I,
althoughthe standarddeviationofρ couldbe calculated FIG.3: Therescaledandbinnedhistogramn(s;θ)sτ∗,where
onlyroughly. Thismightexplaintheslightmismatchfor
τ∗=2.10for(1/θ)=125,250,500,··· ,32000,64000(asindi-
(1/θ)=32000,L=16000,32000. e
cated) in a doublelogarithmic plot. The linear size L is cho-
For the highest values of (1/θ) we could not yet do
sen according to the bold printed entries in Tab. I and large
the comparison to another system, so the curve for the enoughtoensureabsenceoffinitesizeeffects. Theerror-bars
largestvalue of (1/θ) in Fig. 3 is dotted, as their quality areestimatedfromshorterruns. Therightmosthistogram(in
is not known. However, it is reasonable to assume that all figures dotted, (1/θ) = 64000) could not be crosschecked
it is not affected by finite size scaling. byanotherrun,seetext. Maximaaremarkedbyarrowspoint-
ing upwards, minima are marked by arrows pointing down-
wards. Thedashedlinesbelongtodifferentexponents,whose
value is specified as the sum of the slope in the diagram and
3.2. The scaling function τ∗,i.e. ahorizontallinewouldcorrespondtoanexponent2.1.
Theshortly dashedlineareestimated exponentsfor different
Comparing the different histograms n(s;θ) for differ- regions of the histogram (2.22 within approx. [20,200] and
ent values of (1/θ) in a plot enables us not only to find 2.19 within [200,2000]), the other exponents are from litera-
the exponent τ, but also to find the ueniversal function ture,namely2.14(3)in[3,4]and223/91≈2.45in[12]. Since
itwasimpossibletorelatetheseexponentstoanypropertyof
G as defined in Eqn. 1. A rough, naive estimate of τ is
given by n(s;θ) fitted against s−τ, which gives a value thedata,theexactposition ofthelinesassociated withthem
of τ∗ ≈ 2.1 in our case. Plotting now n(s;θ)sτ∗ double was chosen arbitrarily.
e
logarithmically should allow us to find the “true” value
e ∗
of τ by performing a data collapse, i.e. choosing τ in
such a way that horizontal shifts (corresponding to the minimacannotbe afeatureofthesameuniversalscaling
choiceofthe scales0(θ) in the scalingfunction) make all function. Otherwise (1) would hold and the quantity
∗
curves collapse. This is shown in Fig. 3, where τ = 2.1
waschosensothatthemaximaforthesecondbumpsare n(s (θ);θ)sτ∗ (θ) , (6)
min min
almost equally high: denoting their position on the ab-
scissa for each value of θ by s (θ), we have chosen τ∗ where s (θ) deneote the position of the minima, would
max min
such that assume the same value for all θ, because they are local
minima of G, which is supposed to be the same for all θ.
n(s (θ);θ)sτ∗ (θ)≈const. . (5) Since these minima cannot be included in the simple
max max
scaling defined in (1), they must be explicitly excluded
e
According to (1) the constant is simply the maximum byintroducing alowercutoff, sothat simple scalingsup-
value of G, namely G(s (θ)/s (θ)), where the value of posedly sets in only above these cutoffs, excluding espe-
max 0
the argument is therefore the same for all θ. cially the minima. However, such a lower cutoff would
The value of τ∗ is close to (but not within the error apparently have to diverge for (1/θ) → ∞ – something
of) the exponent found in the literature, τ = 2.14(3) that is certainly beyond any established concept of scal-
[3, 4] (τ = 2.15(2) in [16], τ = 2.159(6) in [6]), which is ing. Evenwhen accepting this peculiar scaling behavior,
shown in the same figure for comparison. However, it is a data collapse for the second bump still seems to be
impossible to force the minima (see the down pointing unsatisfactory, as shown in Fig. 4.
marks in Fig. 3) to the same height while maintaining If one accepts a divergent lower cutoff of the scaling
theconstraintthatthemaximaremainaligned,i.e. these function,onehastofacethefactthatthiswoulddescribe
5
2 Thisconceptappearstoberathernaive–ontheother
hand, it is hard to assume that (1) can still hold: it
Increasing 1/θ would correspondto (7) with f(s) replacedby sτ, which
is a straight line in a double logarithmic plot. Therefore
(1) can apply only to a region in Fig. 3 where the data
that fall ontop of eachother forma straightline. Those
*τs featuresnotalreadycollapsingwouldthencollapsewhen
1) properlytilted(choosingthe rightτ) andshifted(choos-
s)/n( ing the right s0). Introducing a lower cutoff at s = 10
n( and discarding the data for (1/θ)≤2000 then leads to a
data collapse in a narrow range as shown in Fig. 5. It is
1 worthwhile mentioning that even for some 10<s<200,
namely for values ofs between the squaresand the filled
circles, none of the data collapse. The exponent used in
0.8 this “collapse” is τstp. =2.19, as mentioned above.
10−1 100
s/s
max
6
FIG. 4: The rescaled and binned histogram n(s;θ)sτ∗, 5 s=10
versus s/smax(θ), where τ∗ = 2.10 for (1/θ) = s=200
e
125,250,500,··· ,32000,64000 in a double logarithmic plot.
Thescaless (θ)bywhichthehistogramshavebeenshifted 4
max
are the maxima marked in Fig. 3, so that a data collapse
p.
would bepossible. The arrow indicates theorder of thedata τsst
in increasing (1/θ). 1) 3
n(
s)/
n(
the behavior of n in a region, which becomes physically
less and less interesting in the limit (1/θ)→∞, because 2
e
the vast majority of events are situated at small s and
as the second bump moves out to infinity, the scaling Increasing 1/θ
function hence covers a smaller and smaller part of n.
10−3 10−2 10−1 100 101
However,onlyaregionofnwhichcoversanon-vanishing
fraction of events can be physically relevant. e s/smin(θ)
e
Concentratingnowonthebehaviorofnuptothemin-
imum (see arrows pointing downwards in Fig. 3), one FIG. 5: The rescaled and binned histogram n(s;θ)sτstp.,
finds that this region is also badly descreibed by a func- versus s/smin(θ), where τstp. = 2.19 for (1/θ) =
e
tion like (1). First of all, the question of which region 4000,8000,16000,32000,64000 in a double logarithmic plot.
Thescaless (θ)bywhichthehistogramshavebeenshifted
is supposedly described by the function needs to be an- min
are slightly different from the minima marked in Fig. 3, to
swered. A unique lower cutoff and a θ dependent upper
make the collapse as good as possible. The squares and the
cutoffneeds to be found. At firstviewit looksappealing
filledcircles mark s=10ands=200, respectively,fororien-
to choose these two marks such that they cover the set
tation and relation to other figures. The arrows indicate the
of data, where the curves fall on top of each other. In
order of thedata in increasing (1/θ).
this case the lower cutoff would be 1 and the upper cut-
off, s , would have a value smaller than the minima
naive
marked by downwards pointing arrows in Fig. 3. How-
ever, this would be described by a function like
By considering the function f(s) it becomes apparent
n(s;θ)=f(s)G(s/s ) (7)
naive
that n, and therefore the model, cannot be scale free:
ratherthan(1).eNotetheparameterindependentfunction it depends on the fixed, microscopic scale s = 1. This
f(s)describingthe shape ofthe curve,whileG(s/snaive) entailes that it is always possible to tell (1/θ) by looking
is a sharp cutoff function. Eqn. 7 does not allow for an onlyattheshapeofn;adiagramshowingonlythisshape,
exponent, f(s) is an arbitrary function. Writing it as withoutanyscalesontheaxes,reveals(1/θ),sinceascale
e
f(s)=s−τ(a +higher order corrections) (8) is intrinsically given by the features of f(s). One would
0 only needto rescaleandtilt it until itfits the plotFig. 3
defines τ to be the steepest descent of this part of the andonecouldidentify(1/θ). Onlyiff(s)werescalefree,
curve and gives a value between τ = 2.22 and τ = i.e. astraightlineinadoublelogarithmicplot,wouldthis
stp. stp.
2.19 (see Fig. 3). not be possible.
6
their scaling in (1/θ) does not depend strongly on this
choice. In particular x and x are different for all
min max
∗
choices of τ . A plot of s (θ) versus (1/θ) for differ-
104 ∗ min
ent values of τ is shown in Fig. 6. For small values of
(1/θ) the minimum is not pronouncedenoughto survive
1.20
∗
for large values of τ , so these curves do not give a data
θ) point. Using a linear fit of logs (θ) versus log(1/θ) of
s(min103 Increasing τ∗ the minimum as found in the rmesincaled (τ∗) and binned
histogram,givesan“exponent”between x =0.93and
0.94 min
x = 0.98. The same procedure applied to the max-
min
ima gives an “exponent” in the range x = 1.18 and
max
102 xmax =1.22, shown in Fig. 7. One may expect that xmin
∗
tends towards x for decreasing τ , as s increases
max min
and might enter the scaling region of s , but neither
102 103 104 105 max
“exponent”exhibitsasystematicvariation,andthequal-
1/θ
ity of the fit certainly suffers from the rough procedure
that searches for the extrema in the binned histogram.
FIG. 6: The position of the minimum in the binned Thisisunfortunatelynecessarybecauseofstatisticalfluc-
and rescaled histogram for different values of τ∗ =
tuations,inconjunctionwiththeabsenceoferror-barsfor
2.04,2.08,2.10,2.12,2.16. The exponents shown in the plot
all data points.
are for comparison only.
The scale of the clusters, s is related to the cor-
min/max
relation length ξ by the fractal dimension µ, i.e. (see
106 [4])
s ∝ξµmin/max . (11)
min/max
105 Since ξ ∝ (1/θ)ν, one should expect ν =
x /µ . The minima are supposed to be dom-
1.20 min/max min/max
inated by smaller, fractal events (see [12]), so µ =
min
θs()max104 0.94 1im.9a6(a1r)e[m4]oraenldikethlyerteofobreedνomminin∈at[e0d.4b7y,0c.o5m0]p.acTthfierems,asxo-
ν ∈ [0.59,0.61]. It is unclear how the two exponents
max
ν are related exactly to the exponents of the two
min/max
correlationlengthsfoundbyHoneckerandPeschel[6]for
103 the connected correlation function ν = 0.576(3) and for
Increasing τ∗ the tree-tree correlation function ν =0.541(4).
102 103 104 105
1/θ 4. DISCUSSION
FIG. 7: The position of the maximum in the binned Prima facie the DS-FFM looks like a percolation pro-
and rescaled histogram for different values of τ∗ = cess, and one might naively think that it is indeed a
2.04,2.08,2.10,2.12,2.16. The exponents shown in the plot
percolation process which organizes itself to the criti-
are for orientation only.
cal density: sites are occupied randomly and indepen-
dently and (at least in the thermodynamic limit) there
is only one cluster which is removed with non-vanishing
3.3. Two length scales probability, namely the largest. In this way the den-
sity of occupied sites is automatically reduced below the
percolationthreshold whenever the thresholdis reached.
That n contains features to define at least two scales,
It is puzzling how remarkably close the tree density in
which apparently diverge in (1/θ) with different expo-
e the DS-FFM is to the density of empty sites in critical
nents, becomes clear when analyzing the scaling of the
site percolation on a square lattice (ρ ≈ 0.4078 and
minima and maxima as marked in Fig. 3, using the defi- FFM
1−ρ =0.40725379(13)[9]respectively). However,the
nitions perc
removalprocessinvolvedinthe DS-FFMintroducesspa-
s (θ) ∝ (1/θ)xmin (9) tial correlations which are not present in standard per-
min
s (θ) ∝ (1/θ)xmax . (10) colation. These correlations are expressed, for example,
max
in the form of a patchy tree density distribution [12].
Ofcourse,theexactpositionoftheextremaofn(s;θ)sτ∗ The purpose of this paper is not to add yet another
∗
depends on its tilt, i.e. on the choice of τ . However, model to the enormous zoo of SOC models. However,in
e
7
5 5
1/θ=2000 ρ=0.592746
4000
4 1/θ=1000 4 ρ=0.4005
ρ=0.3975 ρ=0.4025
3 3
*τn(s)/n(1) s2 *τn(s)/n(1) s 2
1
1100 101 102 103 100 101 102 103 104 105
s
s
FIG. 8: The rescaled and binned histogram n(s;θ)sτ∗ FIG.9: Therescaledandbinnedhistogramn(s;ρ)sτ∗ (again
τ∗ = 2.10), for a modified model, where the largest cluster
(again τ∗ = 2.10), for a modified model, where the largest e
e in the system is removed in each relaxation step and the
cluster in the system is removed after each driving step,
corresponding number of trees is filled back into the sys-
for (1/θ) = 1000,2000,4000 (as indicated) with linear sizes
tem afterwards. The three small values of ρ chosen, ρ =
L = 2000,2000,4000. The inset shows the same data on
0.3975,0.4005,0.4025 correspond (uptothelast digit) tothe
the scale of Fig. 3 for comparison. The data for (1/θ) =
valuesof thetreedensityfor(1/θ)=1000,2000,4000 respec-
1000,2000,4000 of the original model as shown in Fig. 3 are
tively, see Tab. I. The linear size was L = 1000,2000,4000.
dotted. The peculiar behavior of the different height-scaling
Thecorresponding dataoftheoriginal modelareshown dot-
of the minimum and themaximum is again visible.
ted. The peculiar behavior of the different height-scaling of
the minimum and the maximum is again visible (a correct
tiltτ∗ wouldmakeitevenmorepronounced),butdisappears
order to investigate certain features of the given model obviously for ρ=ρperc – for these data it is relevant to men-
tion that n(s) was measured after the relaxation. The filled
and identify underlying mechanisms, it makes sense to
circlesshowtheexactresultsforthelatticeanimals[5,17,18]
modify it slightly. The outcome for the histogram of e
at ρ=ρperc.
the DS-FFM modified such that the largest cluster is re-
moved after each driving step is shown for a few values
of (1/θ) in Fig. 8. The distinctive feature of a mini-
mumwhichscalesdifferentlyfromthemaximumisagain vary much if a much smaller system size is simulated at
present, as the peaks of the maxima have approximately thisdensity,soweexpectitessentiallytobefreeoffinite
the same height, while the height of the local minima size corrections. Since it represents a correlated percola-
variesamongdifferentvaluesofθ. Theinsetofthisfigure tionprocess,itisjustconsistentthatthisbumpdoesnot
showsthe histogramonthesamescaleasFig.3together covertheexactresultsforthelatticeanimalsofstandard
withthedataoftheoriginalmodel(dotted)withthecor- percolation [5, 17, 18] at ρ=ρ shown as filled circles
perc
respondingvaluesofθ. Onecanunderstandthattheydo inFig.9. Thedottedgraphsinthefigureshowthecorre-
not fall on top of each other, because the relaxationrule sponding data of the original model. Again they do not
in the modified model erases much larger clusters than match apart from the region of very small s. Unfortu-
in the original model. nately the simulations of the so-modified model are very
Fig.9showsasecondmodificationofthemodel,where expensive in CPU time, because the mass of the largest
againthelargestclusterisremovedduringrelaxationand cluster needs to be refilled after each relaxation, so that
in addition the driving is changed such that the density only 50.000 updates for transient and statistics could be
of trees, ρ, is the same before each relaxation; the trees done.
removedduringtherelaxationarejustfilledinrandomly Since the feature ofdifferentscalingsurvivesthe mod-
afterwards. This model differs from standard percola- ificationsdescribedabove,itseemsreasonableto assume
tion only by its updating scheme [20]. In order to com- that any relaxation rule that favors the largest cluster
pare the outcome with the original model, the values of leads to the peculiar behavior. Its disappearanceathigh
ρ have been chosen close to the values given in Tab. I. densities can be explained by the extremely small cut-
Indeed, the feature of different scaling of the extrema is off in the distribution, which leads to a domination of
still present, but it disappears completely if the density the statistics by very small clusters, while a single, enor-
is increased to ρ = 0.592746 [9], which is shown in mouslylargeonedominatestheburning(theaveragesize
perc
the same figure as the large bump. This curve does not of the burnt cluster for ρ=ρ was 355811). However,
perc
8
much more careful and detailed investigations of mod- This gives rise to an exponent τ = 223/91≈ 2.45, how-
els like the modification described above are required to ever, this is definitely not supported by numerics (see
gain a full understanding of the underlying mechanisms. Fig. 3).
In particular, this should include a modification of the It remains completely unclear how to characterize the
rules such that the feature disappears. scalingofthe DS-FFMin twodimensions. Apparently it
Honecker and Peschel [6] have calculated the correla- is not a mere superposition of two simple scalings,as re-
tionlengthnotonlyfortheprobabilitythattwositesbe- centlyspeculated[12]. Moreoverthemodeldoesnotseem
longtothesamecluster,butalsofortheprobabilitythat to be scale free as described above and it does not seem
twositesareoccupiedatall. Thecorrelationfunctionfor to be possible to identify a unique powerlaw behaviorof
thelatterisofcourseaδ peakinordinarypercolation,as theclustersizedistribution. Neverthelesseffective power
therearenospatialcorrelationsforthedistributionofoc- lawbehavioroverrestrictedregionshasclearlybeenpro-
cupied sites by construction. However, in the DS-FFM duced by the model, making it potentially relevant to
the correlation length for this quantity, ξ, is finite and observation.
seems to diverges when approaching the critical point.
Allwecanconcludeis thatthe DS-FFM is not critical
It is highly remarkable that Honecker and Peschel con-
inthe senseofsimplescaling. Itremindsusthatadiver-
clude from their simulations that this correlation length
gent moment (here hsi, the second moment) can be re-
diverges slightly slower than the correlation length of
gardedasauniquesignofemergentscaleinvarianceonly
the probability for two sites to belong to the same clus-
ifwearecertainthatonesinglescaleissufficienttochar-
ter, ξ . This seems to indicate that for sufficiently large
s acterize the system. If there is more than one relevant
scales the spatial correlation of the occupation proba-
scale,differentpropertiesofthe systemmightdependon
bility becomes arbitrarily small, so that on sufficiently
different scales which may or may not diverge.
large scales the DS-FFM occupation is uncorrelated and
therefore tends to standard percolation. In other words,
it seems to be possible to rescale or “renormalize” the
DS-FFM to make the occupation correlation arbitrarily Acknowledgments
small. Thiswouldintroducehigherorderinteractions,as
known from standard real space renormalization group
The large scale simulations represent the backbone of
and would explain the difference in critical density be-
this paper. Andy Thomas at the Department of Mathe-
tween the rescaled DS-FFM and standard percolation.
matics at Imperial College was eagerto support us tech-
However, if this “mapping” is valid, one should find the
nically, making machines available and upgrading the
exponent for the divergence of ξ /ξ to be the same as in
s
systems in our department such that parallel comput-
standard percolation, but this is precluded by numerics.
ingbecame apleasure. A greatdealofthe resultsinthis
Ithasbeensuggestedatleasttwice[6,12],thattheDS-
paper was possible only because of his work. We thank
FFM is a superposition of cluster distributions n
perc(s,p) him very much.
of standard percolation for a whole range of concentra-
tionsp,weightedbyacertaindistributionfunctionw(p), This work partly relies on resources provided by the
1 Imperial College Parallel Computing Centre. We want
i.e. dpw(p)n(s,p). Obviouslysuch anassumption ne-
0 to especially thank K. M. Sephton for his support.
glectRs spatial correlations. We recall the following result
from standard percolation theory [5], Anotherpartofthisworkwaspossibleonlybecauseof
the generous donation made by “I-D Media AG, Appli-
n(s,p)∝s−τC(−s/(p−p )−1/σ) , (12) cationServers&DistributedApplicationsArchitectures,
c
Berlin”. We especially thank M. Kaulke and O. Kilian
where C denotes the cutofffunction andthe exponents σ for their support.
and τ have their standard definitions. Assuming that G.P. wishes to thank A. Honecker, I. Peschel and K.
perc
theweightingfunctionw(p)isanalyticaroundthecritical Schenk for very helpful communication, as well as N. R.
concentration in standard percolation, pc, (12) leads to Moloney for providing the lattice animal data and for
proofreading.
1
dpw(p)n(s,p)∝s−(τperc+σ) . (13) TheauthorsgratefullyacknowledgethesupportofEP-
Z
0 SRC.
[1] B. Drossel and F. Schwabl, Phys. Rev. Lett. 69, 1629 Matter 8, 6803 (1996).
(1992). [4] S. Clar, B. Drossel, and F. Schwabl, Phys. Rev. E 50,
[2] H.J.Jensen,Self-Organized Criticality (CambridgeUni- 1009 (1994).
versity Press, New York,NY,1998). [5] D. Stauffer and A. Aharony, Introduction to Percolation
[3] S.Clar,B.Drossel,andF.Schwabl,J.Phys.C:Condens. Theory (Taylor & Francis, London, 1994).
9
[6] A.Honeckerand I.Peschel, Physica A 239, 509 (1997). [15] Ferrenberg, Landau, and Binder, J. Stat. Phys. 63, 867
[7] G. Pruessner and H. J. Jensen, in preparation. (1991).
[8] J. Hoshen and R. Kopelman, Phys. Rev. B 14, 3438 [16] P.Grassberger, J.Phys.A:Math.Gen. 26,2081(1993).
(1976). [17] M. F. Sykes and M. Glen, J. Phys. A: Gen. Phys. 9, 87
[9] M. E. J. Newman and R. M. Ziff, Phys. Rev. Lett. 85, (1976).
4104 (2000). [18] S. Mertens, J. Stat. Phys.58, 1095 (1990).
[10] A. Honecker, Program for simulating the two- [19] As thestandard deviation (actually theestimator of the
dimensional forest-fire model, http://www-public.tu- standard deviation of the estimated mean; this quantity
bs.de:8080/˜honecker/software/forest2d.html (1997). includes the correlation time) decreases with the square
[11] W.H.Press,S.A.Teukolsky,W.T.Vetterling,andB.P. rootofthecomputingtime,onehastocomparetheprod-
Flannery,NumericalRecipesinC (CambridgeUniversity ucts of the square root of the computing time and the
Press, NewYork, NY,1992), 2nd ed. standard deviation. Inthealgorithm presented here,the
[12] K. Schenk, B. Drossel, and F. Schwabl (2001), preprint ratio of thecomputing timevaries between 1.4 and 2.1
cond-mat/0105121. [20] Actually,italsodiffersfromstandardpercolationbecause
[13] K. Schenk, B. Drossel, S. Clar, and F. Schwabl, Eur. it fixes the number of occupied sites rather than simply
Phys.J. B 15, 177 (2000). the probability of being occupied. However, this differ-
[14] G.Pruessner,D.Loison,andK.-D.Schotte,Phys.Rev.B ence becomes irrelevant for sufficiently large systems.
64, 134414 (2001).
