1/f Noise and Long Configuration Memory in Bak-Tang-Wiesenfeld Models on Narrow Stripes

We report our findings of an 1/f power spectrum for the total amount of sand in directed and undirected Bak-Tang-Wiesenfeld models confined on narrow stripes and driven locally. The underlying mechanism for the 1/f noise in these systems is an exponentially long configuration memory giving rise to a very broad distribution of time scales. Both models are solved analytically with the help of an operator algebra to explicitly show the appearance of the long configuration memory.


I. INTRODUCTION
The ubiquitous 1/f noise fascinated physicists for generations [1,2]. There are many examples from wildly different systems in which the power spectra S(f ) ∼ 1/f α with α close to one. This phenomenon usually indicates the presense of a broad distribution of time scales in the system. The common case where α = 1 is particularly intriguing, in that it implies a kind of "equal partition" of power among every decade of frequency range, i.e. the integral is independent of f . One mechanism to generate such a distribution of time scales in systems at thermal equilibrium is through thermal activation events over a sufficiently broad and flat distribution of energy barriers [2]. The "local" power spectrum generated by any single barrier has a characteristic frequency which decreases exponentially with increasing barrier height. But the superposition of the power spectra from all the barriers gives rise to an 1/f spectrum. This mechanism is often employed to explain, for example, the 1/f spectrum of low frequency voltage fluctuations in semi-conductors [2]. In search for a more general answer, applicable to nonequilibrium and dynamic systems, Bak, Tang, and Wiesenfeld (BTW) introduced the notion of Self-Organized Criticality (SOC) [3]. In particular, they proposed a simple "sandpile" (BTW) model which shows the emergent scale free behavior in both space and time. However, the original BTW model did not exhibit the 1/f noise [4]. In this paper we report the observation of 1/f noise for directed and standard (undirected) BTW models confined on narrow stripes (quasi-one-dimensional geometries). In these models, sand flows in the long direction, with periodic or closed boundary conditions in the other direction. The system is driven locally by randomly adding sand to a unique set of sites that have the same coordinate along the long axis. The total amount of sand in the sandpile as a function of time, measured by the number of added grains, exhibits a clean 1/f power spectrum with an exponentially small lower cutoff. Surprisingly, the mechanism for the 1/f noise in this athermal nonequilibrium model is rather similar to the above mentioned thermal mechanism. In our model the local characteristic frequency also falls off exponentially as a function of some parameter, which in this case is the distance from the driving point.

II. DIRECTED MODELS
Let us first consider the simpler directed model, defined as follows. An integer variable z(x, y) is assigned on every site (x, y) of a two-dimensional lattice of size L x × L y (1 ≤ x ≤ L x , 1 ≤ y ≤ L y ). Throughout the paper, we refer to z(x, y) as the number of grains of sand (or height) at the site (x, y). The dynamics consists of the following steps: (i) Add a grain of sand to a randomly selected site in the first column, (1, y): z(1, y) → z(1, y) + 1.
(ii) If as a result of the process the height z(x, y) exceeds a critical value z c = 2, the site topples and three grains of sand are redistributed from this site to three of its nearest neighbors up, down, and to the right, that is (iii) Repeat step (ii) until all sites are stable, i.e. z(x, y) ≤ 2 everywhere. This chain reaction of updates is referred to as an avalanche.
(iv) When the avalanche is over, measure the total amount of sand in the system Z(t) = z(x, y). Then go to step (i).
Notice that the flow of sand is directed to the right along the x-axis. Also note the separation of time scales -the duration of individual avalanches is taken to be much faster than the unit time interval which is defined by the addtion of sand grains. The boundary condition in the x direction is always set to open: z(L x +1, y) = 0. While in the y direction we either choose the periodic boundary condition (which we refer to as Model 1) or the closed boundary condition. In the latter case we restrict ourselves to L y = 2 and refer to it as Model 1A. In Model 1A, we set z c = 1 and the redistribution rule (ii) prescribes to move two grains of sand from the toppling site: one to the right along the x direction and the other to its nearest neighbor (up or down) in the y direction. After some transient period the above dynamics brings the system to a stationary state, where the total amount of sand in the system Z saturates and fluctuates about its average value. At this point we start recording Z(t) and measure its power spectrum S(f ) = |Ẑ(f )| 2 , wherê Z(f ) is the Fourier transform of Z(t). A time window of Z(t) is shown in Fig. 1 together with its running averages. Notice the fluctuations on many time scales. In Fig. 2, we show the power spectra for Models 1 and 1A. Even for small systems one observes a very broad 1/f region. In fact, as we will demonstrate later, the lower cutoff of the 1/f region falls off with L x exponentially. Our simulations indicate that as the width of the stripe in Model 1, L y is increased the 1/f region shrinks, as shown in Fig. 2(c). Direct observations of configurational changes at each time step clearly indicate that the rate of configurational changes at x decreases drastically with increasing x [5], suggesting that there are many time scales and some kind of long memory in the system. To understand this, we proceed with solving Model 1A using the group of operators introduced by Dhar [6,7] to treat sandpile models. To simplify the notation let us denote the configuration at the pair of sites z(x, 1) and z(x, 2) by the column z(x, 2) z(x, 1) , and let L x = L. Any pair configuration with both z(x, y) ≤ z c = 1 is stable. However, the recurrent pair configurations, present in the stationary SOC state, are 0 1 , 1 0 , and 1 1 , while 0 0 is never realized in the system after some transient period. As usual in directed models [8], there are no additional restrictions on pair configurations at different columns. The total number of recurrent (SOC) states is thus 3 L . Let us define operators U x and D x acting on recurrent configurations. The action of these operators consists of adding one grain of sand at sites (x, 2) and (x, 1) correspondingly, and, if necessary, relaxing the resulting configuration according to the avalanche rules of the dynamics. The final stable configuration is the result of the operator acting on the initial configuration. One can demonstrate that these operators commute with each other, i.e.
= 0, and the model is therefore an Abelian model [6]. The following operator identities [6] result directly from the relaxation rules of the model These identities simply state that the addition of two grains of sand to any site (of a recurrent state) will certainly make it unstable and, therefore, two grains will be transferred to the neighbors according to the relaxation rule of Model 1A. Open boundary at x = L + 1 corresponds to U L+1 = D L+1 = I, where I is the identity operator. From Eqs. (1) and (2) it immediately follows that In other words, addition of one grain of sand to both upper and lower sites at any column x triggers a downstream avalanche, in which two grains fall off the right edge, but the underlying configuration of sand columns remains unchanged.
Repeated application of U 1 3 L times makes the system visit every one of the 3 L recurrent states exactly once and return to its original configuration. If U 1 and D 1 are applied in random order (that is how we drive the system in Model 1A), eventually all recurrent states will still be visited. But, since U m , for a random sequence of U 1 's and D 1 's the average time required to visit all 3 L states is given by √ T = 3 L , or T = 9 L . Now let us study the "microscopic" details of how operators U x and D x change the configurations of the sandpile. Note that It is clear that the only cases where the operator U x or D x can "propagate" to the right and change the state of the next neighbor are U So, in order for the operator U 1 to be able to change the configuration of a pair at x, all x − 1 pairs to the left of this pair must be in the state 1 0 . Thus, the probability that a pair at x changes its configuration as a result of the addition of a grain of sand at the left end of the system is 1/3 x−1 . In other words, the system has an exponentially long configuration memory. If the system is driven by a random sequence of U 1 and D 1 , on average it will take τ x ∼ (3 x−1 ) 2 = 9 x−1 grains of sand to change the configuration of a pair at x. This exponential increase of the characteristic time τ x with x is manifested in the local autocorrelation func- is the number of grains at column x and at time t. In Fig. 3 we plot C(x, t) vs. t/9 x−1 for several x's. One sees that C(x, t) = F (t/9 x−1 ). This form of the local autocorrelation function implies a scaling form for the local power spectrum: S loc (f, propagates through a column it leaves the number of grains on that column Z(x) unchanged. It follows that the addition of a grain at the left end of the system can change at most the number of grains at one column. If we assume that the local events of changing Z(x) are independent for different x's (or the correlations of which is not too strong), which is a reasonable approximation when we drive Model 1A with a random sequence of U 1 's and D 1 's, then the global power spectrum of the total number of grains in the system is the superposition of the local power spectrum. The exponential fall off of the local characteristic frequencies of configuration changes would give rise to a global 1/f power spectrum, similar to the case of thermal activations in equilibrium systems mentioned earlier [2].

III. UNDIRECTED MODELS
We now turn our attention to the undirected sandpile model on a stripe L x × L y . In this model an unstable site with z(x, y) > z c = 3 redistributes one grain of sand to each of its four neighbors. In our simulations at each time step we randomly select a site on the central column (x = (L x +1)/2, or, if L x is an even number we randomly select one of the 2L y sites on the two central columns) and add one grain of sand to that site. We choose to have open boundaries at x = 0 and x = L x + 1. In the first version of the model the boundary condition along the y direction is periodic. We refer to this version as Model 2. In addition to Model 2 we consider a simpler model which we refer to as Model 2A. Model 2A is defined on an L × 2 stripe with closed boundary condition in the ydirection. In Model 2A z c = 2 and a site with z(x, y) > 2 moves one grain of sand to each of its three neighbors. The advantage of studying Model 2A is that it has fewer recurrent configurations which are easier to classify. In Fig. 4, we show the power spectra of the total amount of sand in Models 2 and 2A. Similar to the case of the directed Model 1, in Model 2 the 1/f region shrinks for increasing L y . The dynamics of the undirected models is apparently more complex than that of the directed ones. However, much of the apparent complexity is due to the motion of "troughs" [5] -columns in which all z ≤ z c − 1, so that avalanches can not propagate beyond them [9]. Let us first understand the trough dynamics.

A. Trough dynamics
Let us focus on Model 2A. Like in Model 1A, one can define operators U x and D x . The operator algebra now satisfies the following operator relations: The open boundaries at two ends imply U 0 = D 0 = U L+1 = D L+1 = I, where, as before, I is the identity  To understand the physical nature of this subgroup let us take a closer look at the set of recurrent configurations in Model 2A. In a stable configuration each z(x, y) can take values 0, 1, 2. Each pair thus has 9 stable configurations. However, the number of recurrent configurations is much smaller than 9 L . To check which stable configurations are forbidden in the recurrent set one applies the rule developed in [6]. According to this rule, a subconfiguration at a subset of sites F is forbidden if for every site (x, y) ∈ F z(x, y) is strictly smaller than the number of its neighbors in the subset F . It is clear that the pair 0 0 is forbidden. Let us refer to pairs Having understood that operators O k = U k D k , do not destroy the long term configuration memory of the system but only move, create, and annihilate troughs, we proceed with describing how individual operators U k = D −1 k O k change the configuration. We separate the trivial trough dynamics from others by defining the equivalence relation of operators: if A = BO k , we say that A is equivalent to B and denote it by A ∼ = B. Thus U k ∼ = D −1 k . One can rewrite the basic operator identity (8) as We now derive the relation between operators U k at different sites. Since U 0 = I, Eq. (11) implies U 4 , which gives the recursion relation N (k + 1) = 4N (k) − N (k − 1) with initial conditions N (0) = 1, N (1) = 4. It is easy to show that In a system of size L one has U In other words, any recurrent configuration can be obtained from a given one by the action of some power of U 1 (there are N(L) choices of this power before configurations start repeating themselves) and, if necessary, creation, annihilation, or change of the position of the trough achieved by the the action of L operators O k .
[12]. Asymptotically, only 2 + √ 3 ≃ 3.732 pair configurations per site are allowed in a recurrent state, compared to 9 stable pair configurations. The above operator relations can be easily modified for Model 2 on an L × 2 lattice. For this model the subgroup (10) of operators O k remains unchanged, while (11) becomes U 6 The number of recurrent configurations in this model is given by N

C. Configuration memory
Now we are in the position to address the question of long memory in Model 2A. (For Model 2 on an L×2 stripe these arguments can be repeated step by step with some minor changes.) Let us restrict ourselves to the operator U k acting on a state that has no trough. Indeed, since we are interested in general properties of the equivalency class produced by the action of O k 's, one can alway select from this class a representative state that has no trough. The action of the operator U k on the pair at column k is given by of the rules of avalanche dynamics shows that in this case all pairs associated with this FSC and the pairs next to it will be updated. The other way of creating the above mentioned FSC is for U k to act on a string of pairs of 2 1 ended with 2 0 , 2 2 · · · 2 2 2 1 1 · · · 1 1 0 .
Again, in this case all pairs associated with this FSC and the pairs next to it will be changed. Both of the scenarios require that the starting configuration con-

IV. CONCLUSION
In spite of the apparent differences between directed and undirected models, the mechanism for a long term memory and the 1/f spectrum in these two systems is the same. We summarize our analysis as follows: changes in the configuration. In the directed models these operators do not change the configuration at all, while in the undirected models their ability to change the configuration is restricted to creation, annihilation and motion of the trough.
2) In order to produce irreversible changes in a configuration at a distance k from the place of sand addition, all k pairs in between have to be in a unique peculiar configuration. Since this particular subconfiguration is just one among N SOC (k) ∼ A k possible recurrent subconfigurations, the characteristic frequency of irreversible updates falls off with distance exponentially.
3) Such exponential dependence of the characteristic frequency of updates leads to the 1/f spectrum of the total amount of sand in the sandpile.
It is straightforward to generate the case to higher dimensions in which sand flows in one (say the x) direction with closed or periodic boundaries in other directions. One would still observe the 1/f spectrum. Recently, De Los Rios and Zhang [13] observed an 1/f spectrum in a nonconserved sandpile-like model in which certain fraction of sand is lost in each toppling process. Due to the absence of conservation avalanches themselves are exponentially unlikely to reach a distant site, giving rise naturally to an exponential distribution of time scales. In contrast, in our model avalanches constantly pass through the system but they produce only small changes of the configuration. An 1/f spectrum was also observed previously for a continuously boundary-driven BTW model [14]. Its origin was attributed to a (linear) diffusion of z(x, y) with a noisy boundary condition [14,15], which gives a powerlaw lower cutoff f c ∼ 1/L 2 x for the 1/f spectrum -a mechanism very different than ours.
[8] D. Dhar , y), where a(x, y) is the operator of adding a grain of sand at (x, y). It is easy to see that the Lx operators defined this way and the identity operator O0 = I form the cyclic subgroup.
[11] There is one minor detail one has to take into account in order for this statement to be exactly true: the action of O −1 k is a bit more tricky than simply decreasing the height of both sites in the k-th pair by one. can cause even more changes in the configuration. The system however does not forget its original configuration and operators Om acting on other sites restore the original pattern of z(x, y) up to the creation of "true" troughs 1 1 , 1 0 , and 0 1 .
[12] One can also calculate N