Boundary effects in a random neighbor model of earthquakes Stefano Lise1,a and Attilio L. Stella2,b 8 9 9 1 INFM – Istituto Nazionale per la Fisica della Materia, and 1 SISSA/ISAS – International School for Advanced Studies, n a Via Beirut 2–4, I-34014 Trieste, Italy J 2 INFM–Dipartimento di Fisica e Sezione INFN, 7 Universita` di Padova, 35131 Padova, Italy ] h c e m Abstract - t a t Weintroducespatialinhomogeneities (boundaries)inarandomneighborver- s . sion of the Olami, Feder and Christensen model [Phys. Rev. Lett. 68, 1244 t a (1992)] and study the distributions of avalanches starting both from the bulk m and from the boundaries of the system. Because of their clear geophysical in- - d terpretation, two different boundary conditions have been considered (named n free and open, respectively). In both cases the bulk distribution is described o c by the exponent τ ≃ 3. Boundary distributions are instead characterized by 2 [ two different exponents τ′ ≃ 3 and τ′ ≃ 7, for free and open boundary condi- 2 4 1 tions, respectively. These exponents indicate that the mean–field behavior of v 4 this model is correctly described by a recently proposed inhomogeneous form 5 of critical branching process. 0 1 0 8 9 / t PACS numbers: 64.60.Lx, 64.60.Ht, 5.40.+j a m - d n o c : v i X r a Typeset using REVTEX 1 Some ten years ago Bak,Tang and Wiesenfeld (BTW) introduced the concept of self– organized criticality (SOC) [1]. The term refers to the tendency of a large class of extended dynamical systems to spontaneously organize into a dynamical critical state characterized by long range correlations, in both space and time. Many models displaying SOC have been introduced and used to describe, for instance, sandpiles [1], earthquakes [2] and biological evolution [3]. The commonfeatureof allthese models isthe simplicity ofthe localdynamical rules, to be contrasted to the global complex structures (e.g., fractal) that emerge as a result of the continued local interactions. A model which has recently attracted much attention in the literature is that proposed by Olami, Feder and Christensen (OFC) [2], to mimic earthquake dynamics. This model has revealed to be particularly challenging to understand and has raised many fundamental questions on the very nature of SOC [2,4–6]. The OFC model is a coupled–map lattice model, where to each site (i,j) of a square lattice, is associated a continuous “energy”, E , initially fixed at a random value in the ij interval (0,E ). All the energies are increased uniformly and simultaneously at the same c speed, until one of them reaches the threshold value E and becomes unstable (E ≥ E ). c ij c The uniform driving is then stopped and an “earthquake” starts: E → 0 E ≥ E ⇒ i,j (1) i,j c E → E +αE ( nn nn i,j where nn denotes the set of nearest neighbor sites of (i,j) and α ∈ [0, 1] is a parameter 4 1 controlling the level of conservation of the dynamics (α = corresponds to the conservative 4 case). The “toppling” rule (1) can possibly produce a chain reaction, which ends when all sites are stable again (E < E , ∀kl). The uniform growth then starts again. Boundary kl c conditions imply that sites close to boundaries of the system have a smaller coordination number (except in the periodic case) and, in principle, also a different value for α (α ). BC Open boundary conditions are generally set. In this case α = α. BC The absence of a characteristic length scale in the system is reflected by the behavior of the probability P(s) that an earthquake involves s sites. Indeed, in the stationary state, P(s) ∝ s−τ, where τ is a critical exponent. Contrary to the BTW sandpile model [7], the OFC model is believed to be self–organized 1 critical even when the dynamics is non conserving (α < ) [2,4–6]. In a sense this is a rather 4 surprising, not yet fully understood result, which possibly implies a peculiar mechanism leading to SOC, different from that of conserved models. Moreover, the avalanche distribu- tion exponent τ is non–universal: it seems to depend on the conservation parameter α and even on the type of boundary conditions [2]. One of the peculiarities of this model is the role played by boundaries. Besides being “sinks” for the energy in excess in the system (which, moreover, can be dissipated even in the bulk, due to the non-conserving dynamics), they act as inhomogeneities which frustrate the natu- ral tendency of the model to synchronize. This is believed to be the fundamental mechanism producing SOC in this system [5,6]. Indeed, it has been shown that, with periodic boundary conditions and for sufficiently small values of the conservation parameter α (α ≃ 0.18), the system reaches a strictly periodic state, in which each avalanche involves just one site [4–6]. For larger values of α the situation is slightly more complicated with multiple topplings involved in a single avalanche, but the avalanches are still localized and criticality is not observed. Open boundaries create a discontinuity in an homogeneous system and break 2 the periodic state it would otherwise reach. In ref. [6] it was suggested that sites close to the boundaries start to organize themselves first, building up long range correlations. The critical region grows with time, until, in the stationary state, it invades the whole lattice. Despite the simplicity of local dynamical rules in SOC systems, analytical approaches haverarelybeensuccessful(foranexception, seee.g.[8]). Thisappliesinparticulartomodels of the OFC type. Most of the results achieved so far have been obtained through computer simulations. To overcome this limitation it has sometimes been useful to consider random neighbor (RN) models [9–11], where each site interacts with randomly chosen sites instead of its nearest neighbors on the lattice. This considerably simplifies the problem, though retaining some essential features of the original model. Moreover, RN models can be seen as mean field descriptions of their nearest neighbor counterparts, since spatial correlations are neglected. TheRNversions oftheBTW sandpilemodel[1]andtheBakandSneppen(BS) evolutionary model [3] have been solved exactly in ref. [9] and ref. [10], respectively. The solutions show that these models become equivalent to a critical branching process. Accordingly, the 3 avalanche size distribution exponents is τ = . 2 The RN version of the OFC model [11] is described by the following toppling rule E → 0 E ≥ E ⇒ i,j i,j c E → E +αE ( rn rn i,j wherernstandsfor4siteschosenrandomlyinafinitelatticebox. Openboundaryconditions are implemented by requiring that sites at the box boundary can collect an infinite amount of energy without toppling. In this model all bulk sites are equivalent and, therefore, no geometrical meaning can be attached to boundaries. A numerical investigation of this version of the RN OFC model gave evidence that in the stationary state avalanches are power law distributed in a whole range of α values (α ≤ α ≤ c 1/4; α ≃ 0.225), with an exponent τ ≃ 3/2. These results have been questioned by some c recent works aimed at an exact analytical control of the RN model [12,13]. By studying a continuum–time equation it was deduced that avalanches are localized as soon as α < 1/4, although the mean avalanche size grows exponentially fast as dissipation tends to zero. Anyway, as far as the description of boundary effects is concerned, this approach is even less satisfactory. In fact no “boundary” dissipation is explicitly introduced in ref. [12,13], which, strictly, would make the system “explode” in the conservative limit (α → 1/4). ThemainpurposeofthepresentworkisthatofprovidingaformulationinwhichtheOFC model is still of a RN nature, while allowing a meaningful distinction between boundary and bulk of the system. This formulation, which can be handled numerically, is worth analyzing in view of the important role boundary inhomogeneities are expected to play in the nearest neighborOFCmodel. Moreover, ourmodel, introducingthenotionofpositioninthesystem, constitutes a substantial improvement of the standard RN one, just like a Landau–Ginzburg approach compared to the Weiss theory of ferromagnetism. In order to introduce a proper inhomogeneity due to boundaryeffects, we have proceeded as follows. We have considered the 2D–system as divided into columns, say from 1 to L (each column has two “boundary” sites, which are never allowed to discharge). When a site with energy E > E in column i topples, it distributes an energy αE to two randomly c chosen sites, one in column i−1 and another one in column i+1. Moreover, an energy αE 3 is also received by two randomly chosen sites in column i. In this way the notion of position along the direction perpendicular to the columns acquires a meaning and the effect of the boundaries (located at columns 0 and L+1) can propagate into the system. At the same time, the random choice of sites within columns should guarantee that we are dealing with an inhomogeneous mean–field model. The set of all columns behaves as a finite chain with boundaries, from which also bulk behavior can be extracted in the L → ∞ limit. Boundary conditions are determined by setting the level of conservation of the dynamics α at sites close to the boundaries (α ). Two different possibilities are suggested by the BC Burridge–Knopoff model of earthquakes [14]: a) free boundary conditions, in which α = BC α and b) open boundary conditions, in which α = α. The Burridge–Knopoff model 1−α BC is a driven spring–block model, which can be directly mapped into the OFC model [2]. It can be schematized as a 2–dimensional network of blocks interconnected by springs. In addition, all the blocks are subject to an external driving force, which pulls them, and to a static friction, which opposes their motion. The case of free boundary conditions corresponds to boundary blocks connected only to blocks belonging to the system. The case of open boundary conditions, instead, corresponds to blocks at the boundary coupled also to an imaginary external block (see ref. [2] for further details). Although both of them are probably not realistic, it is believed that more adequate boundary conditions should somehow interpolate between these two extreme limits [15]. Note moreover that free boundary conditions are more conservative than open ones and become strictly conservative 1 for α = . Below we will show that this distinction leads to a radical difference in the 4 distributions of avalanches propagating from the borders of the chain. We performed extensive numerical simulations in order to sample bulk and boundary avalanche size distributions. The bulk avalanche distribution of the system appear not to affected by the choice of boundary conditions. Fig. 1(a) and 2(a) report the results in the conservative case (α = 0.25) for free and open boundary conditions, respectively. Only avalanches starting from the deep interior of the lattice have been taken into account in the statistics. As a matter of fact, bulk behavior is also easily extracted upon, e.g., averaging the avalanche distribution over all possible starting columns in the lattice. The estimated exponents are τ = 1.45 ± 0.1 and τ = 1.5 ± 0.1 for free and open boundary conditions, 3 respectively. In both cases they are consistent with the usual mean–field one, i.e. τ = . 2 On the contrary, avalanches starting from the borders, i.e. from columns 1 or L, are strongly influenced by the choice of boundary conditions. Boundary avalanches are distributed as P(s) ∝ s−τ′, where τ′ = 1.4 ± 0.1 for free boundaries (fig. 1(b)) and τ′ = 1.75 ± 0.1 for open boundaries (fig. 2(b)). We have also checked that boundary conditions interpolating between free and open behave as the open case, i.e. τ′ ≃ 1.75. In this respect τ′ ≃ 1.75 reflects a more robust behavior of the model, appropriate to any non zero level of dissipation realized at the border. Finally we have verified that the introduction of bulk dissipation (α < 1/4) has no apparent effect on bulk and boundary exponents, as long as α > α , as determined in ref. [11]. c However, it is controversial [12,13] whether for α < 1/4 the model should be considered critical. An apparent criticality could result numerically from finite size effects. Our results indicate that the mean–field OFC model with open boundary conditions is correctly described by an inhomogeneous branching process [16]. This generalization of the standardbranching process hasbeenrecently proposedandstudied asa paradigmfor amore 4 complete description of SOC models in the mean–field limit. The inhomogeneous branching process takes place in a situation in which translation invariance is broken. An example is given by the geometry of a semi–infinite chain. Each site of the chain involved in the process at a given stage of the process can activate its neighbors (and/or reactivate itself) at the subsequent stage with well defined probabilities. The probability of generating a tree (or avalanche) with s individuals becomes a function of the position n = 1,...,∞ where the tree starts along a chain. It was shown exactly in ref. [16] that critical trees (avalanches) starting from the boundary are distributed as P1(s) ∝ s−τ′, with τ′ = 47, different from the 3 bulk exponent τ = . Boundary sites were there identified as sites with an average number 2 of branchings smaller than 1 [17], reflecting in this way a sort of “dissipative” behavior. In this respect open boundary conditions in our model correspond to the existence of these “dissipative” boundaries of the inhomogeneous branching process. On the other hand, if the average number of branchings is maintained equal to 1 for all sites of the chain (including the end), it is possible to verify numerically that, in the inhomogeneous branching process, the exponent τ′ coincides with the bulk one (i.e. τ′ ≃ 3) [18]. This case seems to be realized 2 by our model with free boundary conditions. Thus, free conditions here seem adequate to keep at its critical bulk value (i.e. 1) the average number ofdescendants for each generations occurring at the border, in the branching process underlying earthquake dynamics. In conclusion we have shown how it is possible to introduce proper boundary effects at the level of a RN OFC model. The model we have considered can also be interpreted as a 1D nearest–neighbor model with “∞–components”. Indeed, the energies in a given column can be thought of as the many “components” of an energy vector associated to a single site of a 1D chain. Since we have always considered very large L, our approach effectively corresponds to an “∞–components” limit of the 1D model. We have verified that, with open boundary conditions, boundary avalanches in the OFC model are distributed as P(s) ∝ s−47, as predicted by the analysis of the inhomogeneous branching process. Free boundary conditions, instead, correspond to “conservative” bound- ary conditions in the inhomogeneous branching process and imply τ′ ≃ τ ≃ 3/2. As a prospective for future work, it would be interesting to investigate whether the mecha- nism of invasion from the boundaries of the “self–organized” region (named “phase locking” in ref. [6]) is actually present also in the random neighbor OFC model with boundaries. We expect that also in the nearest neighbor OFC model boundary scalings different from the bulk ones should be observed with appropriate conditions at the borders. We acknowledge Claudio Tebaldi for many useful conversations on the subject and for a critical reading of the manuscript. a Email: [email protected] b Email [email protected] 5 REFERENCES [1] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59,381 (1987); Phys. Rev. A 38, 364 (1988) [2] Z. Olami, H.J.S. Feder, and K. Christensen, Phys. Rev. Lett. 68, 1244 (1992). K. Christensen and Z. Olami, Phys. Rev. A 46, 1829 (1992) [3] P. Bak and K. Sneppen, Phys. Rev. Lett. 71, 4083 (1993) [4] J.E.S. Socolar, G. Grinstein, and C. Jayaprakash, Phys. Rev E, 47, 2366 (1993). [5] P. Grassberger, Phys. Rev. E 49, 2436 (1994). [6] A.A. Middleton and C. Tang, Phys. Rev. Lett. 74, 742 (1995). [7] S.S. Manna, L. B. Kiss, and J. Kerte´sz, J. Stat. Phys. 61, 923 (1990). [8] D. Dhar, Phys. Rev. Lett. 64, 1613 (1990). See also E.V. Ivashkevich, D.V. Ktitarev, and V.B. Priezzhev, J. Phys. A 27 L585 (1994). [9] K. Christensen and Z. Olami, Phys. Rev. E 48, 3361 (1993). [10] H. Flyvbjerg, K. Sneppen and P. Bak, Phys. Rev. Lett. 71, 4087 (1993). J. de Boer, B.Derrida, H.Flyvberg, A.D.JacksonandT.Wetting, Phys. Rev.Lett.73, 906(1994). [11] S. Lise and H.J. Jensen, Phys. Rev. Lett. 76, 2326 (1996). [12] M. L. Chabanol and V. Hakim, Phys. Rev. E 56, R2343 (1997). [13] H. M. Broker and P. Grassberger, preprint adap–org/9707002. [14] R. Burridge and L. Knopoff, Bull. Seismol. Soc. Am. 57, 341 (1967). [15] K. Christensen and Z. Olami, J. Geophys. Res. 97, 8729 (1992). [16] G. Caldarelli, C. Tebaldi and A. L. Stella, Phys. Rev. Lett. 76, 4983 (1996). [17] In the critical branching process bulk sites give rise on average to exactly one branch within an avalanche. See T. E. Harris, The Theory of Branching Processes (Dover, New York, 1989) [18] E. Montevecchi and A. Stella, in preparation. FIGURE CAPTIONS Fig. 1: Simulation results for (a) bulk and (b) boundary avalanches in a system with free boundary conditions. The system size is L = 100 and the conservation parameter is α = 0.25. The estimated exponents are (a) τ = 1.45±0.1 and (b) τ′ = 1.4±0.1. Fig 2: Results for (a) bulk and (b) boundary avalanches in a system with open boundary conditions. System size are L = 100, 200 and conservation parameter is α = 0.25. The estimated exponents are (a) τ = 1.5±0.1 and (b) τ′ = 1.75±0.1. 6 Fig. 1 -2 10 (a) P(s) -4 10 -6 10 -2 10 (b) P(s) -4 10 -6 10 1 2 3 4 10 10 10 10 s Fig. 2 -2 10 (a) P(s) -4 10 -6 10 -2 10 (b) P(s) -4 10 -6 10 1 2 3 4 5 10 10 10 10 10 s