Power spectrum of rare events in two-dimensional BTW model, violation of $1/f$ noise hypothesis
Z. Moghadam M. N. Najafi A. Saber Z. Ebadi
cond-mat.stat-mech






Department of Physics, University of Mohaghegh Ardabili, P.O. Box 179, Ardabil, Iran
zahramoghadam.physics@gmail.com 


Department of Physics, University of Mohaghegh Ardabili, P.O. Box 179, Ardabil, Iran
morteza.nattagh@gmail.com


Department of Physics, University of Mohaghegh Ardabili, P.O. Box 179, Ardabil, Iran
ahad.saber@gmail.com 


Department of Physics, University of Mohaghegh Ardabili, P.O. Box 179, Ardabil, Iran
z.ebadii@uma.ac.ir





One of the primitive aims of the two-dimensional BTW model had been to explain the 1/f^α noise which is widely seen in the natural systems. In this paper we study some time signals, namely the activity inside an avalanche (x(t)), the avalanches sizes (s(T)) and the rare events waiting time (REWT) (τ(n) as a new type of noise). The latter is expected to be important in predicting the period and also the vastness of the upcoming large scale events in a sequence of self-organized natural events. Especially we report some exponential anti-correlation behaviors for s(T) and τ(n) which are finite-size effects. Two characteristic time scales δ T_s and δ T_τ emerge in our analysis. It is proposed that the power spectrum of s(T) and τ(n) behave like ( b_s,τ(L)^2+ω^2)^-1/2, in which b_s and b_τ are some L-dependent parameters and ω is the angular frequency. The 1/f noise is therefore obtained in the limit ω≫ b_s,τ. b_s and b_τ decrease also in a power-law fashion with the system size L, which signals the fact that in the thermodynamic limit the power spectrum tends to the Dirac delta function.



05., 05.20.-y, 05.10.Ln, 05.65.+b, 05.45.Df


Power spectrum of rare events in two-dimensional BTW model,

violation of 1/f noise hypothesis
    Z. Ebadi
    December 30, 2023
==============================================================================================






§ INTRODUCTION

Many natural processes are explainable in terms of some local simple rules and show degrees of self-organized criticality (SOC) <cit.>, most of which show 1/f^α power-law decay in the tails of their power spectra of time signals. The term flick noise refers to the phenomenon that a signal s(t) fluctuates with a power spectrum PS_s(f)=1/f^α at very low frequencies. Since the exponent α is often close to 1, flick noise is also called 1/f noise <cit.>. Examples are electrical noise <cit.>, solar flares <cit.>, stock market price variations, rain <cit.> and earthquake <cit.>. The destructive features of the large scale events in the natural disasters motivates one to model the harmony of the occurrence of them which helps to predict their behaviors <cit.>. Due to the random nature of the problem, statistics and probability analysis are extensively deployed as basic mathematic tools to analyze rare events based on historical data. However, it is difficult to acquire accurate and adequate data in practice as very limited information can be recorded. In this regard the theoretical modeling and the statistical analysis on the corresponding rare events yields some valuable information concerning their sequence pattern and also their intensities. Despite of random nature of the SOC systems, their power-law harmony restricts strongly their spatial and temporal correlations and other behaviors, giving the chance to make some predictions on the future of the sequence <cit.>.
 
The sandpile model as a prototype of self-organized critical systems is a dissipative system which is slowly driven and displays irregular bursts of scale free large activity named as avalanche. This model was designed firstly by Bak et al. to explain the 1/f^α noise in natural phenomena <cit.>. These models show critical behavior without fine tuning of any external parameter. The examples of the systems that this idea have had impact are ecology <cit.>, cognitive processes <cit.>, magnetic systems <cit.>, superconductors <cit.>, mechanics <cit.> and magnetosphere <cit.>. Studying the BTW sandpile model as the first example of such systems is very helpful to find the temporal and spatial structure of correlations, which is vital for predicting the future of a sequence of avalanches. 

Soon after the work of BTW, many authors investigated the correlations and also the power spectrum of the sandpile models and many estimations of the decay exponents were obtained, ranging from Lorentzian α=2 <cit.> to more exact exponent α=1.59 <cit.>. Two types of correlations should be distinguished: within a single avalanche, and over time scales greater than the interval between successive particle additions, which measures correlations between avalanches. The α exponent has been obtained for the relaxation events inside an avalanche in the mentioned references, and less attention has been paid to the statistical correlations between distinct avalanches <cit.>. The importance of the latter case can be understood in rare events statistics which is expected to have connections with natural self-organized rare events such as earthquake and making predictions. The precise determination of correlation/anti-correlation and their decay times is also important, e.g. in estimating the time required for rare events to become nearly independent.

Many attempts have been made to find the structure of temporal and spatial correlations/anticorrelations in sanpile models <cit.>. Anticorrelations in sandpiles are caused by discharge effects, meaning that after a large scale event occurs, the system loses dramatically its energy content and needs some time to re-obtain the necessary energy to trigger another large event. In this paper we consider 2D BTW sandpile model and analyze three time scales: x(t) which is the number of unstable sites in the internal time t, s(T) which is Tth avalanche and τ(n) which is the waiting time between nth and (n+1)th avalanches. We observe that a time scale emerge for s(T) and τ(n) which limits the domain of the temporal anti-correlations, leaving the only possibility that the corresponding power spectrum behaves like (b_s,τ^2+ω)^-1/2 in which b_s≡ 1/δ T_s, b_τ≡ 1/δ T_τ and 1/δ T_s and 1/δ T_τ are the time scales for s(T) and τ(n) respectively. The x(t) noise however show scale-free behaviors with two distinct intervals, having its roots in the multi-fractality of the avalanches. The nearly 1/f^1.6 behavior is observed for this noise in accordance with the previous results.


The paper has been organized as follows: In the following section, we explain the method and define the functions. In SEC <ref> we present the results for the three noises. We end the paper by a conclusion.






§ THE CONSTRUCTION OF THE PROBLEM


The main ingredients of avalanche propagation in natural phenomena (as well as the sandpile models) are a slow driving, a local threshold (or non-linearity) for the dynamics and a dissipation mechanism.

Before going into the details, let us first briefly introduce the standard BTW model on a regular d-dimensional hypercubic lattice <cit.>. Let each site of the lattice i has an integer height (energy) E_i ≥ 1. At initial state, one can set randomly the height of each site in which E_i ≤ E_c. E_c is the threshold height equal to the number of nearest neighbors of each site (e.g, for hypercubic lattice E_c=2d). At each time step a grain is added on a randomly chosen site (E_i → E_i+1). If the height of this site exceeds E_c, a toppling occurs: E_i→ E_i+Δ_i,j in which Δ_i,j=-E_c if i=j, Δ_i,j=1 if i and j are neighbors and zero otherwise. A toppling may cause the nearest-neighbor sites to become unstable (have hight higher than E_c) and topple in their own turn and so on, until the entire lattice sites are below the critical threshold (stable state). The total process which starts by a local perturbation (making the first site unstable) until reaching another unstable configuration is called an avalanche. The model is conservative and the energy is dissipated only from the boundary sites. The properties of the model in d=2 has been investigated extensively and well understood in the literature <cit.>, as well as d=3 <cit.>.

We have two kinds of time, namely the internal time (denoted by t) and the external time (denoted by T). For defining these times, we have partitioned the total number of topplings in an avalanche into some internal time steps. One unit of internal time is defined as a L^2 search for unstable sites. Therefore one can define the sequence { x(t) }_t=1^t_max in which x(t) is the number of topplings at time t. The sum of activities in an avalanche ∑_t=1^t_max(T)x(t) ( in which t_max(T)≡ total internal time needed for the Tth avalanche to be ended) equals the total size of that avalanche s(T) (≡# toppling events in the Tth avalanche) which can be considered as another noise. In fact x(t) records the number of topplings (local relaxation events) taking place in the sandpile during each parallel update of the whole lattice, one such update defining the unit of internal time, whereas s(T) records the avalanche size in Tth avalanche.

For the internal time signal (x(t)) of avalanche one can present the following scaling argument to obtain the α exponent. The power spectrum is defined as 

    PS_x(f)≡⟨| 1/t_0∫_0^t_0 dt x(t,t_0)exp[-2iπ ft] |^2 ⟩

in which ⟨⟩ show ensemble averaging and t_0 is the time duration of the avalanche and x(t,t_0) is time series conditioned to be terminated in t=t_0. Note that in this ensemble averaging, the duration of t_0 is also sample-dependent. To facilitate the notation we use the quantity PS(f,t_0) which is the power spectrum of avalanches with fixed duration t_0, for which:

    PS_x(f|t_0)≡   1/t_0^2∬_0^t_0dtdt' exp[-2iπ f(t-t')] 
       ×⟨ x(t,t_0)x^*(t',t_0)⟩

The total power spectrum is readily obtained by a second averaging over the avalanche sizes, i.e. PS_x(f)=∫ P(s)PS_x(f|s)ds in which P(s) is the probability density function of s. Knowing that the avalanche size scales with the its time duration via the relation s∼ t_0^γ_st_0, and using the facts that PS_x(f|s)=s^2F_1(f^γ_st_0s) and x(t,t_0)=t_0^γ_st_0-1F(t/t_0) <cit.> and also P(s)∼ s^-τ_s, one obtains:

    PS_x(f)= f^-γ_st_0(3-τ_s)∫_0^s^*f^γ_st_0dxx^2-τ_sF_1(x).

in which s^* is the size cut-off of avalanches which is determined by finite size of the system. If we let s^*→∞ the integral equals a constant and one obtains α=γ_st_0(3-τ_s).

The similar arguments also hold for the time series s(T) as a signal. Here the power spectrum is obtained without ensemble averaging, by the following relation:

    PS_s(f)≡lim_T_0→∞| 1/T_0∫_0^T_0 dT s(T)exp[-2iπ fT] |^2

In calculating this quantity, we can begin by one initial configuration and let it to reach the steady state, after which the above quantity is obtained up to some large time T_0.


Another quantity which can be seen as a time series is the rare events waiting time (REWT). Apparently we have freedom to define a rare (large scale) event. Fixing a proper definition of the rare events needs some statistical analysis to minimize the sensibility of the results to the definition. In this paper we have defined a rare avalanche as an avalanche whose size is more than twice the lattice size, i.e. s^threshold=2L^2. In fact we have examined various rates for this threshold value (for defining the rare event) and have done simulations for s^threshold=L^2,1.5L^2,2L^2,3L^2 and 5L^5 and the results have not been changed. In Fig. <ref> we have shown n(s^>)/N in terms of s^threshold/L^2 in which n(s^>) is the number of avalanches larger than threshold and N is the total number of avalanches. It is seen that for s^threshold≳ 10^0 n(s^>)/N becomes slow-varying. The more general quantity is the waiting time between avalanches with sizes larger than s, i.e. t_w^>s <cit.>. To avoid unnecessary complications, we have fixed the threshold. After fixing the definition of large events, we define REWT, denoted by τ(n), as the time interval between two successive rare events s(T_n) and s(T_n+1). In other words if the (n)th rare event occurs at time T_n (i.e. s(T_n)>s^threshold) and the next rare event occurs at time T_n=1, then τ(n)≡ T_n+1-T_n. The REWT has been shown in Fig. <ref> for a system with size L=128 in which the threshold has been shown by a broken line and some τ's have been shown. Therefore a time series of rare events is formed, which is denoted by {τ(n)}_n=1^n_max in which n_max is the maximum of n in the simulations and T_0≡0. Now the corresponding power spectrum is calculated by the following relation:

    PS_n(f)≡lim_n_max→∞| 1/n_max∑_n=0^n_maxτ(n)exp[-2iπ fn] |^2.

The statistics of REWT enables us to make some statistical predictions on the large scale events, which is important in natural disasters as explained in the introduction. Especially the autocorrelation of s(T) REWT is very important in the rare event analysis. These correlation demonstrate the tendency of a large scale event to have a more/less vast event in the future, i.e. answers the question whether the rare events amplify or weaken each other in a time sequence. The answer whatever it is, seems to be of vital importance in the natural self-organized critical systems, provided that the connection to the underlying theoretical model is correct. Regardless of the exact form of the autocorrelation functions, the positivity or negativity of the time-correlations of the rare events is expected to have a deep impact in the analysis of the corresponding natural phenomena. In this paper we concentrate mainly on the time signals s(T) and τ(n), and also regenerate the results for x(t) which has been investigated vastly in the literature (see for example <cit.>).




§ RESULTS


We consider two-dimensional BTW model on square lattice of linear sizes L=32,64,128,256,512 and 1024. For calculating our desired quantities, more than 10^7 avalanches have been taken into account for each lattice size. We start with a random height distribution h_i∈ [1,4] and inject the sand grains randomly through the sample. Once the system reached the steady state, the statistical observables are analyzed. 




Figure <ref> shows the autocorrelation of the noise x(t) (A_x(t,t')≡⟨ x(t)x(t')⟩) and its power spectrum (PS_x(ω)) for L=512. For x(t) the autocorrelation function A_x(t,t') depends separately on both t and t' and not t-t' (remember that it is not invariant under the transformation t→ t+a and t'→ t'+a). This function shows multi-fractal behavior. For t=1, there are two distinct power-law regimes, one with t'^0.58± 0.0.02 and another with t'^-0.17± 0.02. In the inset the power-law behavior for the power spectrum is seen with the exponent ≈ 1.6 in agreement with the previous results <cit.>. To illustrate the other features of this noise, we have investigate its other statistical properties. x̅(t)≡⟨ x(t)⟩ and δ x(t)≡√(⟨ x(t)^2⟩-⟨ x(t)⟩^2 ) have been reported in Fig. <ref>. For long enough times x̅(t)∼ t^-5 and δ x(t)∼ t^-1.0± 0.1, whereas for the small times the behaviors are different. The point at which the behavior changes is compatible with the corresponding point in Fig <ref>, and results from the multi-fractality of the avalanches. In the Fig. <ref> we have shown the distribution function P(x,t) (≡ the number of events that at internal time t, there are x unstable sites) for various times t in terms of x, in which x̅(t) is the position of their peaks and also δ x(t) is their variance. 




Now let us turn to the calculation of the correlation between the avalanches s(T). The autocorrelation function is defined by:

    f_s(T_0)≡⟨  s(T)s(T+T_0)⟩_T -( ⟨  s(T)⟩_T)^2

in which the T-average of an arbitrary statistical observable is defined by ⟨ O⟩_T≡1/T_max∑_T=0^T_max O(s(T)). Figure <ref> shows the re-scaled autocorrelation function of s(T) defined by A_s(T)≡( ⟨  s(T)⟩_T)^-2f_s(T). An interesting feature of this graph is the exponential anti-correlation of the noise s(T). It is seen that A_s(T) is negative for small times and grows in an exponential fashion with time approaching zero for long times. This effect is magnified for small lattice sizes. This function is properly fitted to the analytic expression -A_s0exp(-b_sT) for all lattice sizes considered in this paper. To show this dependence, we have shown A_s(T)/A_s0 in terms of b_sT in the inset of this figure for various lattice sizes. The dependence of b and A_s0 to the inverse of the lattice size 1/L has been shown in <ref> in a log-log plot, from which we see that A_s0 tends to zero for L→∞ with the exponent γ_A_s0,L≃5/4 defined by A_s0∼ L^-γ_A_s0,L. This shows that, although the correlations become more long range for larger systems, the amplitude of the anti-correlation vanishes in the thermodynamic limit and is therefore a finite size effect. The anti-correlation behavior means that when a rare event takes place at time T, it is not favorable for the system to have such a large scale avalanche in the approximate time interval (T+1,T+1/b).

Having the fitted analytic form of the autocorrelation functions, one can calculate the power spectrum simply by taking the Furrier transform, which leads to the relation |PS_s(ω)|∼( b_s(L)^2+ω^2)^-1/2. Therefore the large frequency (small time) limit of power spectrum is PS_s(ω)∼ω^-1. The power spectrum of the s(T) noise has been shown in Figs <ref> and <ref>. The linear dependence on 1/√(b_s^2+ω^2) is evident in Fig <ref> with L-dependent slope which is not of primary importance. The (nearly) -1 exponent for large ω values is seen in Fig <ref>. This is the flick noise for 2D BTW model.




Empirically the analysis of REWT seems to be more important. The time signal τ(n) has direct connection to the statistics of the time signal s(T). The statistical analysis of this quantity helps to predict vastness of the upcoming large scale events. A complete analysis of this quantity has been presented in Fig. <ref> in which all quantities presented for S(T) have been regenerated. In Fig. <ref> A_τ(n) has been shown, which is fitted by -A_0τexp[-b_τT] (A_0τ and b_τ are fitting parameters). The exponential behavior signals the presence of a characteristic time scale in the problem, i.e. δ T_τ=1/b_τ above which the correlation between the signals becomes negligibly small. The dependence of b_τ as well as A_0τ on the system size has been shown in Fig. <ref>. We see that as the lattice size becomes larger, the correlations become more long range, just like s(T). Having this fitting form of autocorrelations, one easily obtains that the corresponding power spectrum should be proportional to PS_τ(ω)∼(ω^2+b_τ(L)^2 )^-1/2. This has been shown in the Figs. <ref> and <ref> which display this dependence for various rates of lattice sizes. The linear dependence on (b_τ^2 +ω^2)^-1/2 confirms that 1/b_τ (as 1/b_s in the s(T) signal) is a new time scale in the problem. It is worthy to note here that δ T_τ=1/b_τ is the number of the rare event above which the REWT signals become nearly independent, and the simple time word should not be confused with the real times such as T. Here for s(T) and τ(n) we see that the 1/f noise systematically arises for ω≫ b_τ and ω≫ b_s respectively. Despite this 1/f noise behavior, it is very important to note that for small enough frequencies, there is a time scale which destroy the power-law behavior of the power spectrum of these functions in the BTW model.




§ DISCUSSION AND CONCLUSION



Many features of the sandpile models on the various systems are known. One of the primitive aims of the BTW model had been to explain the 1/f^α noise which is widely seen in the natural systems. Our work is motivated by the existence of two well-separated time scales in self-organized sandpile models, one related to the spreading of avalanches and the other imposed by the external driving. In this paper we have analyzed many features of the some noises in the BTW system. The activity inside an avalanche (x(t)), the avalanche sizes (s(T)) and the rare events waiting time (REWT) (τ(n)) have been studied and the autocorrelation and the power spectrum of these noises were reported. For x(t) a power-law (1/f^α) decay was seen with α≃ 1.6 in accordance with the previous results. The multi-fractal behavior of avalanches was observed in the time dependence of A_x(t,t')=⟨ x(t)x(t')⟩ for which two different power-law dependence were obtained. The analysis of the probability density P(x,t) shows also these distinct behaviors. For s(T) and τ(n) however some exponential anti-correlation behaviors were observed. Both signals become more long-range as the lattice size increases. Despite of this fact, the amplitudes of these anti-correlations weakens as the system size increases.

The anti-correlation behavior of s(T) shows that when a large scale event takes place, it is favorable for the system to generate smaller avalanches in the close next steps and vice versa. This behavior demonstrates that sandpiles have discharging effects. The interpretation of the anti-correlation behaviors of τ's is more complicated. This behavior means that the largeness of the time period between two successive rare events, induce a small time period between two successive rare invents in the next step and vice versa. These temporal anti-correlations have also been observed previously in the BTW model <cit.> and one-dimensional sandpile model with log-normal form autocorrelations <cit.>. The spatial correlation/anticorrelation behaviors of the sandpile models can be deduced from the Hurst exponent H (defined by C_2(r)≡⟨ |x(r_0+r)-x(r_0)|^2⟩∼ r^2H in which x(r) is the toppling number of the sand column at t) in which H≃ 0.66>0.5 for the BTW model shows that the BTW model is smooth and less fluctuating and corresponds to a correlated surface <cit.>, whereas for the model on the Bethe lattice, they are weakly anti-correlated <cit.>. This correlation/anti-correlation behavior is model dependent in sandpile models and can be regarded as the measure of cross over between them <cit.>.

Besides the anti-correlation, the exact form of the correlation is also of central importance. The exponential dependence of the autocorrelation functions of s(T) and τ(n) time signals is the result of the violation of temporal scale-invariance in the BTW model. The emergent time scales (δ T_s≡ 1/b_s and δ T_τ≡ 1/b_τ) also cast the overall behavior of power spectrum to two disjoint intervals, e.g. ω≪ b_s and ω≫ b_s (and the same relations for b_τ). For the latter case we have observed the famous 1/f noise, whereas for the former case the power spectrum is dominated by b_s. The overall dependence of the power spectrum of s(T) and τ(n) are PS_s,τ(ω)∼(ω^2+b_s,τ(L)^2)^-1/2. If our extrapolation to the large systems L→∞ is true, and b_s and b_τ go to zero in a power-law fashion, then the power-spectrum PS_s and PS_τ tend to the Dirac delta function in the thermodynamic limit. Therefore our results support the fact that the 1/f noise in the BTW model is a finite size effect.






