Analysis of return distributions in the coherent noise model
The return distributions of the coherent noise model are studied for the system size independent case. It is shown that, in this case, these distributions are in the shape of -Gaussians, which are the standard distributions obtained in nonextensive statistical mechanics. Moreover, an exact relation connecting the exponent of avalanche size distribution and the value of appropriate -Gaussian has been obtained as . Making use of this relation one can easily determine the parameter values of the appropriate -Gaussians a priori from one of the well-known exponents of the system. Since the coherent noise model has the advantage of producing different values by varying a model parameter , clear numerical evidences on the validity of the proposed relation have been achieved for different cases. Finally, the effect of the system size has also been analyzed and an analytical expression has been proposed, which is corroborated by the numerical results.
Throughout the last two decades the interest in extended dynamical systems has experienced a steady increase. These systems exhibit avalanches of activity whose size distributions are of power-law type. Although there is not a unique nor unified theory which totally explains all the features of these complex systems, there exist several known mechanisms producing power-law behavior. One of the most popular and well-studied mechanisms is that of self-organized criticality (SOC) introduced by Bak, Tang and Wiesenfeld . Many physical systems and models have shown to exhibit SOC . The most important feature of all these systems is that the entire system is under the influence of a small local driving force, which makes the system evolve towards a critical stationary state having no characteristic spatiotemporal scale, without invoking a fine-tuning of any parameter. On the other hand, SOC is not the only mechanism causing power-law correlations that appear in a nonequilibrium steady state. Another simple and robust mechanism exhibiting the same feature in the absence of criticality is the coherent noise model (CNM) [3, 4]. The CNM is based on the notion of an external stress acting coherently onto all agents of the system without having any direct interaction with agents. Therefore, the model does not exhibit criticality, but it still gives a power-law distribution of event sizes (avalanches).
Recently, it was presented an analysis method to interpret SOC behavior in the limited number of earthquakes from the World and California catalogs by making use of the return distributions (i.e., distributions of the avalanche size differences at subsequent time steps) . In their work Caruso et al obtained the first evidence that the return distributions seem to have the form of -Gaussians, standard distributions appearing naturally in the context of nonextensive statistical mechanics [6, 7]. Based on the assumption that there is no correlation between the size of two events, they were also able to propose a relation between the exponent of the avalanche size distribution and the value of the appropriate -Gaussian as
which is rather important since it makes the parameter determined a priori and therefore it acquits of becoming a fitting parameter. The only little drawback of their work was that the number of data taken from the catalogs is not sufficiently large to obtain a very precise exponent and also clear return distributions with well-defined tails (which is important in order to verify how good the distribution approaches a -Gaussian). Consequently, Eq. (1) could not be rigorously tested until a very recent effort by Bakar and Tirnakli in , where the same analysis was made using a simple SOC model known as the Ehrenfest dog-flea model in the literature  (see also [10, 11]). Thanks to the simplicity of the dog-flea model, it was possible to achieve extensive simulations with very large system sizes (up to ) and also very large number of data elements (up to ). Accordingly, from these extensive simulations, it was obtained a value of , which is in accordance with the “mean-field” exponent determined in several problems [12, 13, 14]. Thence the value of return distributions was deduced a priori from Eq. (1).
In this work, we plod along this way by setting forth the following points: (i) first, we will obtain an exact relation between exponent of the avalanche size distribution and the value of the appropriate -Gaussian without resorting to any assumption and compare it to Caruso et al relation given in Eq. (1), (ii) since the CNM has the advantage of producing different values by varying a model parameter 111In the dog-flea model there is only one available value of since the only parameter is the number of fleas., we now have the opportunity to test the validity of our exact relation (and also the Caruso et al relation) not only for one case but for various cases, (iii) since the corresponding return distributions are expected to converge to the -Gaussian as the system size goes to infinity, the effect of finite system size is also important and we shall try to analyze this effect proposing an analytical expression, (iv) and finally since this model is not a SOC model, our results also give us the possibility of checking the generality of this behavior observed so far in SOC models.
2 The coherent noise model
Let us start by introducing the CNM. It is a system of agents, each one having a threshold against an external stress . The threshold levels and the external stress are randomly chosen from probability distributions and , respectively. Throughout our simulations we use the exponential distribution for the external stress, namely, and the uniform distribution () for . The dynamics of the model is very simple: (i) generate a random stress from and replace all agents with by new agents with new threshold drawn from , (ii) choose a small fraction of agents and assign them new thresholds drawn again from , (iii) repeat the first step for the next time step. The model can be described in the form of a two step-master equation that we present in the appendix. The number of agents replaced in the first step of the dynamics determines the event size for this model. Although the CNM has been introduced for analyzing biological extinctions , it has then been adopted as a very simple mean field model for earthquakes even though no geometric configuration space is introduced in the model . It is shown that the model obeys the Omori law for the temporal decay pattern of aftershocks , exhibits aging phenomena  and power-law sensitivity to initial conditions .
3 Avalanche Size and Return Distributions
3.1 Size independent case
As pointed out in , there is advantage in choosing the uniform distribution () for the thresholds of the CNM agents seeing that the model can be simulated in the limit using a fast algorithm which acts directly on the threshold distribution instead of acting on the agents of the system. This enables us to obtain the avalanche size distribution of the model as being independent of the system size. The distribution is expected to be a power-law over many decades until the values reach a particular point , thereafter it falls off exponentially. From our point of view this is rather important since it means that if we measure the avalanche size exponent using the region , then we must use this value to predict a priori the value of the -Gaussian that the return distribution is expected to converge in the entire region without any deterioration (not only in the central part but also in the tails). The results obtained for the avalanche size distributions of three representative cases with , and are given in the left column of Fig. 1.
Since each case with different values has a different size exponent , this allows us to check the validity of Caruso et al relation given in Eq. (1) or any other equation relating values to the values of the appropriate -Gaussians. From the master-equation of the CNM is theoretically possible to compute the probability of and bringing to bear standard techniques  to obtain the return distribution. However, its level of complexity turns out the solution almost analytically impossible or its (asymptotic) behavior deeply unclear as it happens in several other problems of this class . Regardless, we are in the position where we can propose an exact relation for the return distribution bringing into play no other assumption than the distribution of avalanche sizes, where is the difference between two consecutive event sizes, i.e., . Let us mathematically define the avalanche size distribution,
with being a constant value describing the asymptotic limit . The process of avalanches is completely Markovian (independent) and therefore the probability of the difference of sizes is
where is the Heaviside step function and denotes the previous avalanche size. Making use of  and attending to the symmetric nature of we can explicit the negative branch,
where is a coefficient only depending on and related to the convolution of very large values of with very large values of yielding a dependence 222This can be flatly checked out performing the calculation with .. Thus, the distribution is mainly described by the product of the isolated factor by the incomplete Beta function . Applying the asymptotic behavior of  we finally get,
Taking note of the -Gaussian distribution,
we straightforwardly obtain
This relation is slightly different from the approximate relation presented in  as can be seen in Fig. 2. For , both relations approach and they are almost identical except in the region where values are smaller than . Moreover, only the relation (5) correctly achieves value when . These values define the limits of the domain of each parameter so that the distributions (2) and (4) are normalizable. The approximate relation Eq. (1) does not fulfill this condition as . It should be noted that, since this discrepancy is only meaningful for , the approximate relation predicts the values with a difference from the exact one, which are also acceptable for all the cases we present.
We are now ready to proceed analyzing the return distributions. The centered returns are given in terms of variable
where represents the mean value of a given data set. As can be seen from the right column of Fig. 1, in our simulations we generated the return distributions of the three representative cases of the CNM in order to check the validity of the relation (5). In each case, an extremely large number of events () has been used to build the numerical distribution, namely, the central part and tails. It is clear that the return distribution (green dots) can by no means be approached by a Gaussian. They actually exhibit fat tails which agree with -Gaussians Eq. (4) where characterizes the width of the distribution and is the parameter which should be determined directly from Eq. (5) a priori and therefore is no longer a fitting parameter. In each panel on the right column of Fig. 1, the dashed black lines represent the appropriate -Gaussian with the value obtained from Eq. (5). Perfect agreement with the data can be easily appreciated not only for the tails but also for the intermediate and the very central part as it is demonstrated in the insets.
3.2 Size dependent case
Although we might think that the size independent (i.e., infinite size) case would be enough for such an analysis, we believe that it is still instructive to look also at the size dependent case at least from two different perspectives: (i) we can check how the size of the system affects the shape of the return distributions and whether the tendency is consistent with the infinite size case as the size of the system increases, (ii) unlike the CNM, generic size independent cases cannot be achieved for such systems and thus the only possibility is to always analyze the size dependent case.
It is very easy to implement the size dependent algorithm for the CNM. We just need to apply the previously described steps of the dynamics to a system of agents. As increases, this algorithm clearly slows down and for the same number of events () the larger value of that we can simulate in a reasonable time is . In Fig. 3(a) the behavior of the avalanche size distribution is given for case for various values. It is clearly seen that the power-law regime is always followed by an exponential decay of all the curves and this decay is postponed to larger sizes as increases. For each case, we estimate the value using the standard regression method in the region before the exponential decay (we determine the size of this interval by looking at the regression coefficient to become always more than 0.9997 in each case). Therefore, we should expect that the exponential decay part would tamper with the -Gaussian behavior of the return distributions and this meddling must diminish as gets larger and larger, which is in fact observed in Fig. 3(b) for the return distributions of four representative values. When values are very small, avalanche size distribution has a very short power-law region and the exponential decay part dominates, which simply causes the return distributions to deviate immediately from the -Gaussian shape. As increases, return distributions start approaching the thermodynamic limit (dotted black line), which is a full -Gaussian with , yielding better and better from the central part to the tails, i.e., as the expected scale-free regime sets in.
In order to explain this gradual approach to -Gaussians when finite-size effects are present, let us try to develop a simple mathematical model by considering the differential equation
This equation has very interesting and different solutions depending on the choice of and values (see refs.[7, 21, 22]), but for our purpose, let us concentrate on case and , whose solution is given by
If , then the solution coincides with the -Gaussian, whereas if (which means that ), the solution turns out to be the Gaussian. On the other hand, between these two extremes, namely if and , we obtain a crossover between them. Specifically, for , Eq. (8) approaches a -Gaussian, . Our results, which are depicted in Fig. 3, show that the small values of imply that the -Gaussian form is valid up to rather large values of . On the other hand, for , the exponential outnumbers the remaining terms leading to the Gaussian behaviour,
The approximate dependence of Eq. (8) can thus be split into different regions defined by three values of . Namely the first value is
( is the Lambert function ) whence the curve assumes a power-law dependence described by the exponent that persists up to
when the it starts being perturbed by the Gaussian dependence. Last, there is the final convergence to the Gaussian functional form which occurs at
This crossover seems to coincide with the behavior of the return distributions of the dependent cases as plotted with dashed black lines on top of each curve in Fig. 3(b). This behavior simply reveals that the longer the power-law regime persists for avalanche size distribution, the better the appropriate -Gaussian dominates in the return distribution. Finally, as , the power-law regime prevails for the avalanche size distribution giving forth a return distribution following the appropriate -Gaussian for the entire region.
In this work, we have studied the behavior of the return distributions for the CNM by directly simulating the size independent case. By means of extensive simulations, it is clearly shown that these distributions converge to -Gaussians with appropriate values which are deduced a priori from the exact relation (5) that we developed here. It is worth noting that although the -Gaussian description is actually an analytical approximation the result provides for an understandable depiction of the distribution, which hardly occurs when we keep a special functions representation, with no fundamental accuracy lost. This relation makes the parameter be related to one of the well-known exponents (avalanche size exponent ) of such complex systems and therefore it rescues from being a fitting parameter in this analysis. Moreover, since the model parameter allows us to obtain different values, we were able to check this behavior for various cases. These results clearly imply that the observed behavior is not restricted to self-organized critical models, but instead it seems to be a rather generic feature presented by many complex systems which exhibit asymptotic power-law distribution of avalanche sizes.
We have also investigated the finite-size effect by simulating directly the model dynamics and found that the convergence to appropriate -Gaussian starts from the central part and gradually evolves towards the tails as the system size increases. This is in complete agreement with the gradual extension of the power-law regime in the avalanche size distribution before the appearance of the exponential decay due to finite-size of the system. These results corroborate the analysis of size independent case since it is clearly seen that, as , curves of return distributions for size dependent case converge to the one comes from the size independent case.
Finally it should be noted that, since it is generically extremely difficult (if not impossible) to achieve the size independent case for such complex systems, the size dependent case has its particular importance. Therefore, although the return distributions appear to be -Gaussians for the entire region in the thermodynamic limit, we have tried to propose a mathematical model in order to explain the behavior of return distributions for the size dependent case.
We are indebted to M E J Newman for providing us his fast (size independent) code for the coherent noise model and C Anteneodo for interesting discussions about passage problems and related references. This work has been supported by TUBITAK (Turkish Agency) under the Research Project number 104T148 and by Ege University under the Research Project number 2009FEN027.
Appendix A The CNM master equation
The dynamics of the CNM can be described according to the probability of having agents in the system that at time present a critical value up to , . In conformity with step 1 we can write the master-equation,
with the probability transitions given by
where the first term on the rhs comes from the case and the second one otherwise. The inverse cumulative probability 333For our case, i.e., implies . and means the number of agents with critical value below and . The following elements are
This corresponds to a matrix with vanishing elements below the diagonal. From these relations is then possible to spell out the occurrence of an avalanche of size
which is numerically well described by the power-law (2) with a small value of .
Regarding step 2 the master equation is abstractly pretty much the same,
with the probability transition matrix is given by
where is used to define the subfraction of agents, , whose critical value before updating was less than .
-  Bak P, Tang C and Wiesenfeld K 1987 Phys. Rev. Lett. 59 381
-  Jensen H J 1988 Self-Organized Criticality: Emergent Complex Behavior in Physical and Biological Systems (Cambridge: Cambridge University Press); Bak P, How Nature Works: The Science of Self-organized Criticality (New York: Copernicus)
-  Newman M E J 1996 Proc. R. Soc. London, Ser. B 263 1605
-  Newman M E J and Sneppen K 1996 Phys. Rev. E 54 6226; Sneppen K and Newman M E J 1997 Physica D 110 209
-  Caruso F, Pluchino A, Latora V, Vinciguerra S and Rapisarda A 2007 Phys. Rev. E 75 055101(R)
-  Tsallis C 1988 J. Stat. Phys. 52 479; Curado E M F and Tsallis C 1991 J. Phys. A 24 L69; Corrigenda: 1991 24 3187 and 1992 25 1019; Tsallis C, Mendes R S and Plastino A R 1998 Physica A 261 534
-  Tsallis C 2009 Introduction to Nonextensive Statistical Mechanics - Approaching a Complex World (New York:Springer)
-  Bakar B and Tirnakli U 2009 Phys. Rev. E 79 040103(R)
-  Ehrenfest P and Ehrenfest T 1907 Phys. Z. 8 311
-  Nagler J, Hauert C and Schuster H G 1999 Phys. Rev. E 60 2706
-  Hauert C, Nagler J and Schuster H G 2004 J. Stat. Phys. 116 1453
-  Kac M 1947 Amer. Math. Month. 54 369
-  Redner S 2001 A guide to first-passage processes (Cambridge: Cambridge University Press); Fisher M E 1984 J. Stat. Phys. 34 667
-  Anteneodo C 2009 Phys. Rev. E 80 041131
-  Wilke C, Altmeyer S and Martinetz T 1998 Physica D 120 401
-  Tirnakli U and Abe S 2004 Phys. Rev. E 70 056120
-  Ergun E and Tirnakli U 2005 Eur. Phys. J. B 46 377
-  Anteneodo C and Duarte Queirós S M 2009 unpublished
-  Gradshteyn I S and Ryzhik I M 1980 Table of Integrals, Series, and Products (New York: Academic Press)
-  http://functions.wolfram.com
-  Tsallis C, Bemski G and Mendes R S 1999 Phys. Lett. A 257 93
-  Tsallis C and Tirnakli U 2010 J. Phys.: Conf. Ser. 201 012001