# 𝑑_c=4 is the upper critical dimension for the Bak-Sneppen model

## Abstract

Numerical results are presented indicating $`d_\mathrm{c}=4`$ as the upper critical dimension for the Bak-Sneppen evolution model. This finding agrees with previous theoretical arguments, but contradicts a recent Letter \[Phys. Rev. Lett. 80, 5746-5749 (1998)\] that placed $`d_\mathrm{c}`$ as high as $`d=8`$. In particular, we find that avalanches are compact for all dimensions $`d4`$, and are fractal for $`d>4`$. Under those conditions, scaling arguments predict a $`d_\mathrm{c}=4`$, where hyperscaling relations hold for $`d4`$. Other properties of avalanches, studied for $`1d6`$, corroborate this result. To this end, an improved numerical algorithm is presented that is based on the equivalent branching process. PACS number(s): 64.60.Ak, 05.40.-a, 05.65.+b.

The Bak-Sneppen evolution model B+S  has received considerable attention as an archetype of self-organized criticality BTW , which has been put forward as a general mechanism leading to many non-equilibrium scaling phenomena observed in Nature Bakbook . Although the Bak-Sneppen model was originally developed as an attempt to interpret paleontological data indicating co-evolutionary activity in biological evolution B+S ; PNAS , it has also been used to interpret power law distributions in quiescent periods between earth-quakes Ito . Its generic mechanism of extremal dynamics scaling  has even inspired an algorithm to approximate combinatorial optimization problems EOperc .

In a recent Letter highd , numerical results for several critical exponents of the Bak-Sneppen model in high dimensions were presented. On the basis of that study, and certain scaling arguments, Ref. highd  concluded that the upper critical dimension, where those exponents obtain mean field values, was $`d=8`$. Avalanches, cascading through the system via nearest-neighbor activation events, form domains that were found to be fractal for $`d>2`$. These two findings are surprising: On one hand, avalanches would proceed on a filamentary domain where surface sites outnumber bulk sites. On the other hand, activity would have to return mostly — by pure chance — to sites already within that domain to keep it from growing linearly with time, as required in the mean-field limit.

The claims for Ref. highd  are in stark contrast to earlier derived scaling relations (see Tab. 1 in Ref. scaling ), which would predict an upper critical dimension of $`d=4`$. For instance, the scaling relation for the avalanche cut-off is $`\sigma =1\tau +d/D`$. Its mean-field value is $`\sigma =1/2`$ meanfield ; BoPa1 , and the avalanche distribution exponent attains $`\tau =3/2`$ meanfield ; BoPa1 , while Ref. BoPa1  explicitly derived that $`D=4`$ is the mean field value for the avalanche dimension exponent, implicating $`d=4`$ as the upper critical dimension. But the scaling relations in scaling  were derived under the clearly stated assumption that avalanches would always be compact for $`d<d_\mathrm{c}`$, and activity within them would proceed homogeneously over their domain, with many returns to each site. Ref. highd  argues that avalanches are ramified already for $`d>2`$, allowing for a $`d_\mathrm{c}>4`$.

In this letter we test the claims of Ref. highd  by independent means. To this end, we have developed an alternative algorithm, which has the added benefit of improving temporal cut-offs by more than two decades while eliminating finite lattice size effects entirely. As a result, we find the assumptions underlying the scaling theory in scaling  intact. In particular, our numerical simulation results show that avalanches in $`d4`$ are compact, and that scaling exponents, when asymptotic behavior is extracted, take on mean-field values for $`d=5`$ and 6, very much in contradiction to Ref. highd . Thus we show that the Bak-Sneppen model possesses an upper critical dimension, $`d_c=4`$, below which hyperscaling relations hold, similar to equilibrium critical phenomena. In fact, we have studied a whole range of different properties of the Bak-Sneppen model for $`1d6`$ bsd.pre . Here, we focus directly on those properties relevant with regard to the upper critical dimension.

The Bak-Sneppen model is very easy to state B+S . It consists of random numbers $`\lambda _𝐫`$ between 0 and 1, occupying sites $`𝐫`$ on a $`d`$-dimensional lattice. At each update step, $`s`$, the smallest random number $`\lambda _{min}(s)`$ is located. That site as well as its $`2d`$ nearest neighbors each receive new random numbers drawn independently from a flat distribution on the unit interval. This update step is repeated at the site with the next $`\lambda _{min}(s+1)`$, and so on. The process inevitably evolves toward a self-organized critical state BTW  in which almost all $`\lambda _𝐫`$ are larger than a critical threshold $`\lambda _\mathrm{c}`$ (see Tab. 1), while those $`\lambda _𝐫`$ below are part of an avalanche of activity. The current minimum $`\lambda _{min}(s)`$ is the “active” site while those $`\lambda <\lambda _\mathrm{c}`$ are “unstable” and potentially active; all $`\lambda >\lambda _\mathrm{c}`$ are “stable”. These avalanches are critical and their properties have been described in terms of scaling relations scaling .

The obvious way to implement the model is to specify a lattice, place a random number $`\lambda _𝐫`$ on each site $`𝐫`$, keep an efficient list of those numbers, ordered by size, and repeatedly replace the smallest $`\lambda _{min}(s)`$ and its neighbors. In this method, knowledge of $`\lambda _\mathrm{c}`$ is a priori not necessary to determine the critical properties of the avalanches. This is, in fact, the algorithm employed in Ref. highd  and some other, low-dimensional studies. But in higher dimensions feasible length scales in the pre-defined lattice, with a fixed number of sites $`N=L^d`$, exponentially diminish with dimension. For instance, Ref. highd  used $`N=2^{16}`$, giving $`L=16`$ in $`d=4`$ (only $`L=4`$ in $`d=8`$), or a maximal spatial separation of $`r_{\mathrm{max}}=dL/2=32`$, and a temporal cut-off of maximally $`s_{\mathrm{co}}r_{\mathrm{max}}^D10^6`$, since $`D4`$.

A more efficient way to describe an avalanche utilizes a generalized branching process PMB ; scaling  that has been shown to be exactly equivalent to the Bak-Sneppen model. In the update procedure, the values of the new random numbers that enter the system are unrelated to each other, their predecessors, or any other number in the system. Thus, any number that has no chance of becoming active in itself (i. e. those with $`\lambda >\lambda _\mathrm{c}`$), will not affect the dynamics of the critical avalanche through its particular value. Now, assume that we are at the beginning of a critical avalanche, i. e. all numbers in the system are above $`\lambda _\mathrm{c}`$. All we need to know to describe the ensuing critical avalanche is the existence of a smallest number in the system to initiate the avalanche. In updating that number and its neighbors, either no number below $`\lambda _\mathrm{c}`$ will be created and that critical avalanche happens to be empty, or it will produce numbers below $`\lambda _\mathrm{c}`$ and will progress through more updates until no more numbers below $`\lambda _\mathrm{c}`$ remain at some later time step. Thereafter, a new critical avalanche will start at a new, unrelated location. To describe an individual avalanche, we can assume that it started at the origin, and at all times we need only keep track of those numbers below $`\lambda _\mathrm{c}`$. Lattice addresses appear only as parameters characterizing exclusively those numbers currently unstable. While the number of sites covered by an avalanche $`n_{\mathrm{cov}}`$ can grow as much as linear in time $`s`$, currently unstable sites at most grow like $`s^{1/2}`$ scaling . Thus, our temporal cut-off due to memory constraints $`N`$ behaves like $`s_{\mathrm{co}}N^2`$.

Since our results depend so strongly on results from numerical simulations, we describe our improved method in great detail. To simulate this process we create a list $``$ of currently unstable numbers, $`=\{(\lambda _𝐫,𝐫)|\lambda _𝐫<\lambda _\mathrm{c}\}`$. We initialize $``$ with a single entry $`(\lambda _𝐫=0,𝐫=\mathrm{𝟎})`$. At each update we first remove the smallest $`\lambda _𝐫`$ in $``$ and any of its neighbors in the lattice, if those happen to be in $``$. Then, we draw a new number for each of those sites, but only sort into $``$ those numbers below $`\lambda _\mathrm{c}`$ by storing $`(\lambda _𝐫,𝐫)`$. Addresses $`𝐫`$ could extend arbitrarily far from the origin, thus eliminating any spatial cut-off. It is very easy to keep $``$ compact and sorted according to $`\lambda `$. In fact, it is sufficient to sort in a new $`\lambda _𝐫`$ linearly from the bottom (instead of using a heap, say): the dynamics quickly leads to a list in which almost all numbers are very densely packed just below $`\lambda _\mathrm{c}`$, and almost every number inserted at the bottom moves only a few steps through $``$.

Without explicit lattice reference, it is not easy for this algorithm to check the minimum’s neighbors in the lattice which, if unstable, would need to be removed from $``$. To be sure, we would have to search $``$ at each update to eliminate those neighbors, which would be unreasonably time consuming. Instead, we utilize a procedure similar to hash tables hash : When a lattice address $`𝐫`$ with $`\lambda _𝐫<\lambda _\mathrm{c}`$ is stored in some entry $`_k`$, $`𝐫`$ is mapped into an index $`i=i(𝐫)`$ for a large, sparse array $`𝒜`$, $`|𝒜|||`$. $`𝒜_i`$ in turn stores the index $`k`$ of $`_k`$, or is empty otherwise. When activity returns to the neighborhood of $`𝐫`$, $`i(𝐫)`$ is calculated and $`𝒜_i`$ can be checked at once to track an unstable number in $``$. It is crucial that $`i(𝐫)`$ is unique, of course, but since $`|𝒜|`$ is finite there exist lattice addresses $`𝐫𝐫^{}`$ such that $`i(𝐫)=i(𝐫^{})`$. Those conflicts are rare between any two currently unstable sites, if $`|𝒜|||`$. They can be resolved by moving to the next available index $`i(𝐫)+I`$, where $`I`$ is the offset to the nearest free entry in $`𝒜`$ (for important details, see Ref. bsd.pre ).

We have first used this simulated branching process to extrapolate for the values of $`\lambda _\mathrm{c}`$ bsd.pre  given in Tab. 1 by simulating avalanches at various $`\lambda `$ below $`\lambda _\mathrm{c}`$ (see Ref. scaling ). This extrapolation uses the domain covered by an avalanche, $`n_{\mathrm{cov}}`$, requiring $``$ to record all currently and previously unstable sites. At $`\lambda _\mathrm{c}`$, $`n_{\mathrm{cov}}`$ can increase as much as linear with $`s`$, quickly exhausting memory (see below). Fortunately, for $`\lambda `$ well below $`\lambda _\mathrm{c}`$, avalanche duration and coverage are quickly cut off.

We have run up to $`2\times 10^6`$ critical avalanches in each $`d`$, $`1d6`$, to determine their statistical properties. We used $`||_{\mathrm{max}}=2^{16}`$ and $`|𝒜|=2^{20}`$, easily sufficient to run avalanches up to $`s_{\mathrm{max}}=10^8||_{\mathrm{max}}^2`$ update steps. \[Avalanches were terminated at $`s_{\mathrm{max}}`$ due to time constraints and due to the error in $`\lambda _\mathrm{c}`$, see Tab. 1.\] $`𝒜`$ is already much larger than the lattice used in Ref. highd , but we point out that even on a fixed lattice with $`N=2^{24}`$ sites in $`d=6`$, the largest distance is $`r_{\mathrm{max}}=48`$, i. e. $`s_{\mathrm{co}}<10^6`$. In contrast, our algorithm has $`s_{\mathrm{co}}||_{\mathrm{max}}^25\times 10^9`$ in any dimension, yielding significant spatial correlations for $`r>100`$ even in $`d=6`$ (see Fig. 1).

In a long critical avalanche, almost all unstable numbers cluster densely in $``$ with values just below $`\lambda _\mathrm{c}`$, which puts high demands on a random number generator. Typically, in an ongoing avalanche that has grown to $`10^4`$ unstable numbers at some time $`s`$, $`10^3`$ of those are placed within $`\mathrm{\Delta }\lambda /\lambda _\mathrm{c}10^4`$ just below $`\lambda _\mathrm{c}`$. Using $`\sigma =1/2`$, those numbers can expect to survive up to $`s(\mathrm{\Delta }\lambda /\lambda _\mathrm{c})^{1/\sigma }10^8`$ update steps. To record temporal correlations accurately, it is crucial to resolve these numbers accurately. This requires numbers that are sufficiently random for more than $`8`$ digits! Our data has been obtained with a sophisticated 64-bit random number generator provided to us by the authors of Ref. RNG .

First, we discuss the data for the distribution of avalanche activity, $`P(r,s)`$, which leads to the avalanche dimension exponent $`D`$. To find $`P(r,s)`$, we record the instances of having activity at a (positive) distance $`r`$ relative to the origin at update step $`s`$. Similar to a random walk, the moments of the distribution define the exponent $`D`$ via $`r^q_ss^{q/D}`$. In Fig. 1 we have plotted the data as $`s/r^D`$ vs. $`r`$ for $`d4`$. For $`d<4`$ we have used the value of our best fit for the exponent $`D`$ scaling ; bsd.pre  as given for each dimension in Tab. 1. For $`d4`$, the fitted value corresponded to the mean-field result BoPa1 , $`D=4`$. For $`d=4`$, we conjecture a simple asymptotic form for the scaling behavior,

$`s{\displaystyle \frac{r^4}{\mathrm{ln}r}}(d=4).`$ (1)   
Thus, in Fig. 1, unlike for the other dimensions $`d`$, for $`d=4`$ we have plotted $`s\mathrm{ln}r/r^4`$. For all dimensions the data levels out horizontally for increasing $`r`$. Logarithmic factors, as seen in Eq. (1) for $`d=4`$, are common for scaling behavior at the upper critical dimension logscale .

Similar evidence for $`d_\mathrm{c}=4`$ is provided by the data for the backward avalanche exponent $`\tau _\mathrm{b}^{\mathrm{all}}`$, which is related to the avalanche distribution exponent through $`\tau =3\tau _\mathrm{b}^{\mathrm{all}}`$ scaling . This particular scaling relation does not rely on avalanches being compact and is valid both above and below the upper critical dimension. The data exhibits mean-field behavior, $`\tau _\mathrm{b}^{\mathrm{all}}=3/2(=\tau )`$, already for $`d=5`$ and 6. In Fig. 2 we show the data for the backward-avalanche distribution reduced by the mean field scaling behavior, $`P_\mathrm{b}^{\mathrm{all}}(s)/s^{3/2}`$ as a function of $`1/\mathrm{ln}(s)`$, for $`d=1`$ through $`d=6`$, and for the multi-trait model BoPa1 , which analytic calculations show exhibits mean field behavior. The curves for $`d<4`$ clearly approach zero rapidly for increasing times $`s`$ on this scale, indicating that $`\tau _\mathrm{b}^{\mathrm{all}}>3/2`$ in those dimensions. The corresponding data for $`d>4`$ is clearly bounded away from zero, comparable to the data from the multi-trait model, where $`\tau _\mathrm{b}^{\mathrm{all}}=3/2`$ exactly. This data also shows that there are significant corrections to scaling and cut-off effects, even for the multi-trait model. In this case, simply fitting a power law over the largest region available from the simulation data can give misleading results due to systematic curvature on increasing scales. (Fitted values for $`\tau _\mathrm{b}^{\mathrm{all}}`$ are discussed in Ref. bsd.pre .) In $`d=4`$, the data is again consistent with some logarithmic correction to scaling behavior.

These findings receive strong theoretical support. The only assumption in the scaling theory of Ref. scaling  is that the avalanches are compact below the upper critical dimension where each site that is visited in an avalanche is typically visited many times. Our numerical results indicate that avalanches indeed remain compact for $`d4`$. To test this explicitly, we have determined the probability distribution $`P(n_{\mathrm{cov}},R)`$ of having a finished avalanche that has covered a domain $`n_{\mathrm{cov}}`$ with a radius of gyration $`R`$. As mentioned above, $`n_{\mathrm{cov}}`$ may grow linear with $`s`$, rapidly exhausting memory for $``$. Thus, we used $`s_{\mathrm{max}}=||_{\mathrm{max}}=2^{20}10^6`$ and $`|𝒜|=2^{22}`$. The exponent $`d_{cov}`$ is defined through the moments of the distribution $`P(n_{\mathrm{cov}},R)`$ via $`n_{\mathrm{cov}}R^{d_{cov}}`$ for large $`R`$, providing a measure of the fractal structure of avalanches if $`d_{cov}<d`$. In turn, avalanches are compact if $`d_{cov}=d`$. Using our numerical results, we plot in Fig. 3 the quantity $`n_{\mathrm{cov}}/R^d`$ as a function of $`1/\mathrm{ln}(R^d)`$. If avalanches are fractal ($`d_{cov}<d`$), we would find that asymptotically this quantity should approach zero (like $`R^{(dd_{cov})}`$). In fact, we find that this quantity clearly remains finite for large $`R`$ for all dimensions $`d<4`$. For $`d>4`$, we find a rapid approach to zero on this logarithmic scale, indicating fractal behavior. In $`d=4`$ this quantity appears to vanish as well, but as a linear function of $`1/\mathrm{ln}(R^d)`$. Similar to Eq. (1), we find that the numerical simulations indicate marginal behavior in $`d=4`$ with logarithmic corrections. Thus, in $`d=4`$ avalanches are marginally compact (a “fat fractal”), and the conditions underlying the scaling theory in Ref. scaling  are upheld for all $`d4`$, implying $`d_\mathrm{c}=4`$.

We demonstrate the consistency of our argument by also measuring the scaling of the coverage with the lifetime of an avalanche, $`n_{cov}s^\mu `$ highd . Since there should be only one characteristic length for a compact avalanche, we expect that $`rR`$, and thus, $`\mu =d/D`$ for $`dD`$, and $`\mu =1`$ in the mean-field limit. In Fig. 4, we plot the coverage reduced by its mean-field scaling, $`n_{cov}/s`$, as a function of duration $`s`$. Clearly, the data for $`d=5`$ and 6 is in perfect agreement with mean-field behavior, while for $`d=4`$ the deviation from mean-field behavior are minute, $`\mu _{\mathrm{meas}}0.99`$, even without considering the effect of logarithmic factors. For $`d=3`$ we find $`\mu 0.905`$, quite consistent with the value predicted from $`\mu =d/D0.896`$ using our measured value $`D=3.35`$ in $`d=3`$ (see Tab. 1).

SB acknowledges helpful discussions with M. J. Creutz and T. T. Warnock.