text:
# 1/𝑓^𝛼 fluctuations in a ricepile model

## Abstract

The temporal fluctuation of the average slope of a ricepile model is investigated. It is found that the power spectrum $`S(f)`$ scales as $`1/f^\alpha `$ with $`\alpha 1.3`$ when grains of rice are added only to one end of the pile. If grains are randomly added to the pile, the power spectrum exhibits $`1/f^2`$ behavior. The profile fluctuations of the pile under different driving mechanisms are also discussed.

PACS numbers: 05.40.-a, 05.65.+b, 45.70.-n, 05.70.Ln

The term flick noise refers to the phenomenon that a signal $`s(t)`$ fluctuates with a power spectrum $`S(f)1/f^\alpha `$ at low frequency. Since the exponent $`\alpha `$ is often close to $`1`$, flick noise is also called $`1/f`$ noise. The power spectrum of a signal is defined as the Fourier transform of its auto-correlation function $`C(t_0,t)=s(t_0)s(t)`$, where $`\mathrm{}`$ denotes ensemble average. Usually the signal $`s(t)`$ under study is stationary and its auto-correlation function only depends on $`tt_0`$. The auto-correlation function can be alternatively calculated as

$$C(\tau )=_{\mathrm{}}^{\mathrm{}}s(t)s(t+\tau )𝑑t.$$
(1)
If the signal is real-valued, its power spectrum is just

$$S(f)=\left|\widehat{s}(f)\right|^2$$
(2)
according to Wiener-Khinchin Theorem. Here $`\widehat{s}(f)`$ is the Fourier transform of the signal

$$\widehat{s}(f)=_{\mathrm{}}^{\mathrm{}}s(t)e^{i2\pi ft}𝑑t.$$
(3)
In nature and laboratory, many physical systems show flick noise. For example, flick noise appears in a variety of systems ranging from the light of quasars  to water flows in rivers , music and speech , and electrical measurement . Despite its ubiquity, a universal explanation of the flick noise is still lacking. Bak, Tang, and Wiesenfeld have proposed that self-organized criticality (SOC)  may be the mechanism underlying flick noise. They also demonstrated their idea of SOC with a cellular automaton model, the BTW sandpile model. However, the model they proposed does not exhibit flick noise in its 1D and 2D versions . Especially the 1D BTW sandpile model does not even exhibit SOC behavior. Experiments  on piles of rice grains were done to investigate if real granular systems display SOC behaviors. It was found that ricepiles with elongated rice grains exhibit SOC behaviors. The key ingredient here is the heterogeneity of the system. Because elongated rice grains could be packed in different ways , the ricepile thus has many different metastable states. A 1D model of ricepile was proposed by the same group of authors; avalanche distribution and transit-time distributions were discussed . The model was later investigated by several other authors including us . Attentions were paid to the avalanche dynamics and transit-time statistics of the model. In this paper, we are concerned with temporal fluctuation of the average slope of the system. We want to investigate if it has flick noise or flick-like behavior.

The ricepile model is defined as follows: Consider a one-dimensional array of sites $`1`$, $`2`$, $`\mathrm{}`$, $`L`$. Each site contains an integer number $`h_i`$ of rice grains. Here $`h_i`$ is called the local height of the ricepile at site $`i`$. The system is initialized by setting $`h_i=0`$ for all sites. This means the ricepile is to be built up from scratch. The system is driven by dropping rice grains onto the pile. If one grain of rice is dropped at site $`i`$, then the height of column $`i`$ will increase by $`1`$, $`h_ih_i+1`$. With the dropping of rice grains, a ricepile is built up. The local slope of the pile is defined as $`z_i=h_ih_{i+1}`$. Whenever the local slope $`z_i`$ exceeds certain threshold $`z_i^c`$ (specified in the following), site $`i`$ will topple and one grain of rice will be transferred to its neighbor site on the right. That is, $`h_ih_i1`$ and $`h_{i+1}h_{i+1}+1`$. In this model, rice grains are allowed to leave the pile only from the right end of the pile, while the left end of the pile is closed. So the boundary condition is kept as $`h_0=h_1`$ and $`h_{L+1}=0`$. The ricepile is said to be stable if no local slope exceeds threshold value, that is, $`z_iz_i^c`$ for all $`i`$. Rice grains are dropped to the pile only when the pile is stable. In an unstable state, all sites $`i`$ with $`z_i>z_i^c`$ topple in parallel. The topplings of one or more sites is called an avalanche event, and the size of an avalanche is defined as the number of topplings involved in the avalanche event.

The values of $`z_i^c`$’s are essential to the definition of the model. As in Ref. , every time site $`i`$ topples, the threshold slope $`z_i^c`$ will take a new value randomly chosen from $`1`$ to $`r`$ with equal probability. Here $`r`$ is an integer no less than $`1`$. The parameter $`r`$ reflects the level of medium disorder. The larger the $`r`$, the higher the level of medium disorder of the system. If $`r=1`$ the model becomes the $`1`$D BTW sandpile model. When $`r=2`$, the model reduces to the model studied in reference . In ref. , we have investigated the effect of disorder on the universality of the avalanche size distribution and transit time distribution.

In the present work, we performed extensive numerical simulations on the system evolution. Let us first study the case where rice grains are added only to the left end of the pile, i.e., only to the site $`i=1`$, which is in accordance with the experiment setup . We shall refer to this driving mechanism as the fixed driving. When the stationary state is reached, we record the average slope $`z(t)=h_1(t)/L`$ of the pile after every avalanche. Here time $`t`$ is measured in the number of grained added to the system. Thus we got a time series $`z(t)`$. Typical results about the fluctuation of the slope $`z(t)`$ can be seen in Fig.1.a. We calculate the power spectrum of this time series according to Eq.(2). We find that for not too small systems the slope fluctuation displays $`1/f`$-like behavior. In fact, we got a power spectrum $`S(f)`$ scaling as $`1/f^\alpha `$ at low frequency, with the exponent $`\alpha =1.3\pm 0.1`$.

We note the following points:

(a) the temporal behavior of the ricepile model is dramatically different from the 1D BTW sandpile model. In the 1D BTW model, the system has a stable state , which, once reached, cannot be altered by subsequent droppings of sand grains. Thus the 1D BTW model does not display SOC behaviors. The slope of the 1D sandpile model assumes constant values. The power spectrum is thus of the form $`S(f)=\delta (f)`$. In words, there is only dc component in the power spectrum for the 1D BTW model.

(b) With the introduction of varying threshold $`z_i^c`$ into the ricepile model the behavior of the system becomes much more rich. The avalanche distribution follows power law , which is an important signature of SOC. Besides that, the temporal fluctuation of $`z(t)`$ has a power spectrum of flick type at low frequency.

We notice that the power spectrum is flattened at low frequency for small system. This is a size effect. It is believed that flick noise fluctuation is closely related to the long range spatial correlation in the system. For small system, the spatial correlation that can be built up is limited by the system size, thus the long range temporal correlation required by flick noise is truncated off at low frequency. Our numerical results verified the above discussion. In Fig. 2, we show the power spectrum of $`z(t)`$ for different system sizes. It is clear that when the system size increases the $`1/f^\alpha `$ behaviors extends to lower and lower frequency.

We also investigated the effect of the level of disorder on the power spectrum. This was done by simulating the system with different values of $`r`$. In Fig. 3 we show the power spectra for the system with different $`r`$. It can be seen from this figure, that the power-law (straight in the log-log plot) part extends to lower frequency for higher value of $`r`$. The effect can be understood by the following discussion. When the level of disorder is increased, the amplitude of the slope fluctuation also increases. Larger amplitude fluctuation of $`z(t)`$ requires more grains to be added to the system. Hence the period of this fluctuation increase, which gives rise to the increase of low frequency components in the power spectrum. We have made statistics on deviation of the ricepile slope from its average value in the stationary state. We found that the greater the $`r`$, the greater the deviation. Note that the value of $`r>1`$ does not affect the exponent $`\alpha `$, as it does not alter the universal avalanche exponent $`\tau `$ for the avalanche distribution .

It is also interesting to study the effect of driving mechanism on the temporal behavior of the system. Now, instead of dropping rice grains to the left end of the pile, we drop rice grains to randomly chosen sites of the system. This way of dropping rice grains represents a different driving mechanism, which we shall refer to as the random driving. Numerical simulations with this random driving were made and the time series $`z(t)`$ was recorded. For this way of driving, typical result about the fluctuation of $`z(t)`$ can be seen in Fig. 1.b. It seems that $`z(t)`$ now fluctuates more regularly than it does for the fixed driving. We found that under random driving the temporal fluctuation has a trivial $`1/f^2`$ behavior as in the case of various previously studied sandpiles . In Fig. 3, we compare the power spectra of the two cases with different driving mechanisms. Two groups of curves are shown in this figure, each with a different exponent $`\alpha `$. For the fixed driving, we have $`\alpha 1.25`$, while for the random driving we have $`\alpha 2.0`$. The distinction of the two behaviors is quite clear. This shows that the driving mechanism, as an integrated part of the model, has an important role in the temporal behaviors of the system. Recent work on the continuous version of BTW sandpile model, and that on quasi-1D BTW model also showed that certain driving mechanisms are necessary conditions for these models to display $`1/f`$ fluctuation for the total amount of sand (or say, energy).

Here we present an heuristic discussion on why the exponent $`\alpha `$ is smaller for fixed driving than for random driving. For fixed driving, the height $`h_1`$ at site $`i=1`$ changes almost every dropping, and so is the average slope $`z(t)`$. Thus the high frequency component of the $`z(t)`$ fluctuation has a heavier weight in its power spectrum, and this will make the exponent $`\alpha `$ smaller. While for the random driving, the rice grains are dropped to the pile at random sites, the spatial correlations previously built up can be easily destroyed, making the $`z(t)`$ behave more or less as an random walk. So the exponent $`\alpha `$ for this case shall be very close to $`2`$. For random driving, each site has the same probability in receiving a grain in every drop. For the fixed driving, however, only the left end site receive grains, so there is in some sense breaking of symmetry, which would lead to different scalings in the slope fluctuations.

To see more about the effect of driving mechanism, it is helpful to investigate the profile fluctuations of the pile under different drivings. Because avalanches change the surface of the ricepile, the pile fluctuates around an average profile, and the size of the fluctuations characterize the active zone of the ricepile surface. Let the standard deviation of height at site $`i`$ be $`\sigma _h(i,L)=\sqrt{h_i^2h_i^2}`$. Here $`\mathrm{}`$ represents average over time. Then we calculate the profile width of the ricepile, $`w=\frac{1}{L}_i\sigma _h(i,L)`$, which is a function of the system size. In Fig. 4, we show the profile width of the ricepile for different $`r`$ and different driving mechanisms. It can be seen that $`w`$ scales with $`L`$ as $`wL^\chi `$, and that $`\chi =0.25\pm 0.01`$ for fixed driving, $`\chi =0.09\pm 0.01`$ for random driving. For given driving mechanism, $`w`$ increases with increasing $`r`$. As we stated in Ref. , the parameter $`r`$ reflects the level of medium disorder in the rice pile. Although the level of disorder does not change the scaling exponent $`\chi `$, it does affect the amplitude of fluctuations. For greater $`r`$, the profile width is larger. It can be seen in the figure, that the data points for $`r=4`$ are above that for $`r=2`$, for given driving mechanism. A recent experiment showed that piles of polished rice grains has a smaller profile width than unpolished ones, which have higher level of medium disorder.

In conclusion, we have studied the temporal fluctuation of the slope of ricepile model in its critical stationary state. The power spectrum of this fluctuation is closely related to the driving mechanism. When the rice grains are dropped to the left end of the pile, the slope fluctuates with a flick-type power spectrum, with the exponent $`\alpha =1.3\pm 0.1`$. When the driving mechanism is changed to random driving, the model displays $`1/f^2`$ behaviors. Greater system size and higher level of disorder will extend the frequency range where the power spectrum has the form $`1/f^\alpha `$.

The author thanks Prof. Vespignani for helpful discussions. This work was supported by the National Natural Science Foundation of China under Grant No. 19705002, and the Research Fund for the Doctoral Program of Higher Education (RFDP).


format: text

text:
# Boundary-induced phase transitions in traffic flow

\[

## Abstract

Boundary-induced phase transitions are one of the surprising phenomena appearing in nonequilibrium systems. These transitions have been found in driven systems, especially the asymmetric simple exclusion process. However, so far no direct observations of this phenomenon in real systems exists. Here we present evidence for the appearance of such a nonequilibrium phase transition in traffic flow occurring on highways in the vicinity of on- and off-ramps. Measurements on a German motorway close to Cologne show a first-order nonequilibrium phase transition between a free-flow phase and a congested phase. It is induced by the interplay of density waves (caused by an on-ramp) and a shock wave moving on the motorway. The full phase diagram, including the effect of off-ramps, is explored using computer simulations and suggests means to optimize the capacity of a traffic network.

\]

One-dimensional physical systems with short-ranged interactions in thermal equilibrium do not exhibit phase transitions. This is no longer true if the action of external forces sets up a steady mass transport and drives the system out of equilibrium. Then boundary conditions, usually insignificant for an equilibrium system, can induce nonequilibrium phase transitions. Moreover, such phase transitions may occur in a rather wide class of driven complex systems, including biological and sociological mechanisms involving many interacting agents. In spite of the importance of this phenomenon and a number of theoretical studies , it has never been observed directly. So far only indirect experimental evidence for a boundary-induced phase transition exists in older studies of the kinetics of biopolymerization on nucleic acid templates .

In the present work we report the first direct observation of a boundary-induced phase transition in traffic flow. Analysis of traffic data sets taken on a motorway near Cologne exhibits transitions from free-flow traffic to congested traffic, caused by a boundary effect, viz. the presence of an on-ramp. These transitions are characterized by a discontinuous reduction of the average speed .

Vehicular traffic on a motorway is controlled by a mixture of bulk and boundary effects caused by on- and off-ramps, varying number of lanes, speed limits, leaving aside temporary effects of changing weather conditions and other non-permanent factors. The fundamental characteristic of the bulk motion is the stationary flow-density diagram, i.e. the fundamental diagram, which incorporates the collective effects of individual drivers behavior such as trying to maintain an optimal high speed while keeping safety precautions. The qualitative shape of the flow-density diagram $`j(\rho )`$ is largely independent of the precise details of the road and hence amenable to numerical analysis using either stochastic lattice gas models or partial differential equations .

A by now well-established lattice gas model for traffic flow, the cellular automaton model by Nagel and Schreckenberg , reproduces empirical traffic data  rather well. Fig. 1 shows simulation data for the fundamental diagram obtained from the Nagel-Schreckenberg (NaSch) model. This has to be compared to measurements of the flow $`j`$ taken with the help of detectors on the motorway A1 near Cologne which show a maximum of about 2000 vehicles/hour at a density of about $`\rho ^{}=20`$ vehicles per lane and km . At densities below $`\rho ^{}`$ one observes free flow, while for larger densities one observes congested traffic.

In addition to the density dependence of the flow two important characteristics are derived directly from the fundamental diagram: the shock velocity of a ‘domain wall’  between two stationary regions of densities $`\rho ^{},\rho ^+`$

$$v_{shock}=\frac{j(\rho ^+)j(\rho ^{})}{\rho ^+\rho ^{}},$$
(1)
obtained from mass conservation, and the collective velocity

$$v_c=\frac{j(\rho )}{\rho }$$
(2)
which is the velocity of the center of mass of a local perturbation in a homogeneous, stationary background of density $`\rho `$. Both velocities are readily observed in real traffic. The collective velocity $`v_c`$ describes the upstream movement of a local, compact jam. In the density range $`25\mathrm{}130`$ cars/km $`v_c`$ ranges from approximately $`10`$ km/hour to $`20`$ km/hour (Fig. 1) which has to be compared with the empirically observed value $`v15`$ km/hour . The shock velocity is the velocity of the upstream front of an extended, stable traffic jam. The formation of a stable shock is usually a boundary-driven process, caused by a ‘bottleneck’ on a road. Bottlenecks on a highway arise from a reduction in the number of lanes and from on-ramps where additional cars enter a road .

The experimental data considered here (see Fig. 2 for the relevant part of the highway) show boundary effects caused by the presence of an on-ramp. Far upstream from the on-ramp, free flow of vehicles with density $`\rho ^{}`$ and flow $`j^{}j(\rho ^{})`$ is maintained. Just before the on-ramp the vehicle density is $`\rho ^+`$ with corresponding flow $`j^+j(\rho ^+)`$. Note that no experimental data are available for $`\rho ^{}`$, $`j^{}`$ and $`\rho ^+`$, $`j^+`$ as well as the activity of the ramp. The only data come from a detector located upstream from the on-ramp  which measures a traffic density $`\widehat{\rho }`$ and the corresponding flow $`\widehat{ȷ}`$.

Next the effects of the on-ramp are considered. Cars entering the motorway cause the mainstream of vehicles to slow down locally. Therefore, the vehicle density just before the on-ramp increases to $`\rho ^+>\rho ^{}`$. Then a shock, formed at the on-ramp, will propagate with mean velocity $`v_{shock}`$ (see (1)). Depending on the sign of $`v_{shock}`$, two scenarios are possible:
1) $`v_{shock}>0`$ (i.e. $`j^+>j^{}`$): In this case the shock propagates (on average) downstream towards the on-ramp. Only by fluctuations a brief upstream motion is possible. Therefore the detector will measure a traffic density $`\widehat{\rho }=\rho ^{}`$ and flow $`\widehat{ȷ}=j^{}`$.
2) $`v_{shock}<0`$ (i.e. $`j^+<j^{}`$): The shock wave starts propagating with the mean velocity $`v_{shock}`$ upstream, thus expanding the congested traffic region with density $`\rho ^+`$. The detector will now measure $`\widehat{\rho }=\rho ^+`$ and flow $`\widehat{ȷ}=j^+`$.

Let us now discuss the transition between these two scenarios. Suppose one starts with a situation where $`j^+>j^{}`$ is realized. If now the far-upstream-density $`\rho ^{}`$ increases it will reach a critical point $`\rho _{crit}<\rho ^{}`$ above which $`j^{}>j^+`$, i.e., the free flow upstream $`j^{}`$ prevails over the flow $`j^+`$ which the ‘bottleneck’, i.e. the on-ramp, is able to support. At this point shock wave velocity $`v_{shock}`$ will change sign (see (1)) and the shock starts traveling upstream. As a result, the stationary bulk density $`\widehat{\rho }`$ measured by the detector upstream from the on-ramp will change discontinuously from the critical value $`\rho _{crit}`$ to $`\rho ^+`$. This marks a nonequilibrium phase transition of first order with respect to order parameter $`\widehat{\rho }`$. The discontinuous change of $`\widehat{\rho }`$ leads also an abrupt reduction of the local velocity. Notice that the flow $`\widehat{ȷ}=j^+`$ through the on-ramp (then also measured by the detector) will stay independent of the free flow upstream from the congested region $`j^{}`$ as long as the condition $`j^{}>j^+`$ holds.

Empirically this phenomenon can be seen in the traffic data taken from measurements at the detector D1 on the motorway A1 close to Cologne . Fig. 3 shows a typical time series of the one-minute velocity averages. One can clearly see the sharp drop of the velocity at about 8 a.m.

Also the measurements of the flow versus local density, i.e. the fundamental diagram (Fig. 4), support our interpretation. Two branches can be distinguished. The increasing part corresponds to an almost linear rise of the flow with density in the free-flow regime . In accordance with our considerations this part of the flow diagram is not affected by the presence of the on-ramp at all and one measures $`\widehat{ȷ}=j^{}`$ which is the actual upstream flow. The second branch are measurements taken during congested traffic hours, the transition period being omitted for better statistics. The transition from free flow to congested traffic is characterized by a discontinuous reduction of the local velocity. However, as predicted above the flow does not change significantly in the congested regime. In contrast, in local measurements large density fluctuations can be observed. Therefore in this regime the density does not take the constant value $`\rho ^+`$ as suggested by the argument given above, but varies from 20 veh/km/lane to 80 veh/km/lane (see Fig. 4).

One should stress here that congested traffic data are usually not easy to interpret, because the traffic conditions (mean inflow and outflow of cars on the on- and off-ramps, and so the bulk mean flow) are changing in time. According to our arguments, in a congested regime the detector measures $`j^+`$, solely due to the on-ramp activity. Therefore, $`j^+(t)=j(\rho ^+(t))<j^{}(t)`$ must be satisfied. During times of very dense traffic one expects always cars ready to enter the motorway at the on-ramp, thus guaranteeing a sufficient and approximately constant on-ramp activity. The measured flow is constant over long periods of time which is in agreement with the notion that the transition is due a stable traffic jam. Spontaneously emerging and decaying jams would lead to the observation of a non-constant flow.

The use of our approach is not limited to a qualitative explanation of the traffic data. Beyond that it can also be used to calculate the phase diagrams of systems with open boundary conditions for a large class of traffic models. We modeled a section of a road with on-ramp on the left and off-ramp (on-ramp) on the right using the NaSch cellular automaton . We modify the basic model by using open boundary conditions with injection of cars at the left boundary (corresponding to in-flow into the road segment) and removal of cars at the right boundary (corresponding to outflow). Therefore it can also be regarded as a generalization of the asymmetric simple exclusion process  to particles with higher velocity.

During the simulations local measurements of the velocity have been performed analogous to the experimental setup. For comparison the results of the computer simulations have been included in Fig. 3. Note that even the quantitive agreement with the empirical data is very good. This has been achieved by using a finer discretization of the model, i.e. the length of the cell is considered as $`l=2.5m`$. The results were obtained for $`L=960`$, $`p=0.25`$ and $`v_{max}=13`$. We kept the input probability $`\alpha =0.65`$ constant. Then the free-flow part is obtained using $`\beta =1.`$ and the congested part using $`\beta =0.55`$. The transition was observed at ten minutes after we reduced the output probability. The “detector” was located at the link from site $`480`$ to $`481`$.

Fig. 5 shows the full phase diagram of the NaSch model with open boundary conditions. It describes the stationary bulk density $`\widehat{\rho }`$ as a function of the far-upstream in-flow boundary density $`\rho ^{}`$ and the effective right boundary density $`\rho ^+`$. For the case of an on-ramp (or shrinking road width etc.) at the right boundary corresponds to the situation discussed above. Here, the density is increased locally to $`\rho ^+>\rho ^{}`$. In agreement with the empirical observation we find a line of first order transitions from a free flow (FF) phase with bulk density $`\widehat{\rho }=\rho ^{}`$ to a congested (CT) phase with $`\widehat{\rho }=\rho ^+`$. On this line $`v_{shock}`$ changes sign.

The case of an off-ramp (or expansion of road space etc.) leads to a local decrease $`\rho ^+<\rho ^{}`$ of the density. Here the collective velocity $`v_c`$ (2) plays a prominent role. As long as $`v_c`$ is positive (i.e. in the free-flow regime $`\rho ^{}<\rho ^{}`$, see Fig. 1), perturbations caused by a small increase of the upstream boundary density $`\rho ^{}`$ gradually spread into the bulk, rendering $`\widehat{\rho }=\rho ^{}`$ (FF regime). At $`\rho ^{}=\rho ^{}`$, $`v_c`$ changes sign  and an overfeeding effect occurs: a perturbation from the upstream boundary does not spread into the bulk  and therefore further increase of the upstream boundary density does not increase the bulk density. The system enters the maximal flow (MF) phase with constant bulk density $`\widehat{\rho }=\rho ^{}`$ and flow $`j(\rho ^{})=j_{max}`$. The transition to the MF phase is of second order, because $`\widehat{\rho }`$ changes continuously across the phase transition point.

The existence of a maximal flow phase was not emphasized in the context of traffic flow up to now. At the same time, it is the most desirable phase, carrying the maximal possible throughput of vehicles $`j_{max}`$. For practical purposes our observations may directly be used in order to operate a highway in the optimal regime. E.g. the flow near a lane reduction could be increased significantly if the traffic state at the entry would allow to attain the maximal possible flow of the bottleneck. This could be achieved by controlling the density far upstream, e.g. by closing temporally an on-ramp, such that the cars still enter the bottleneck with high velocity.

We stress that the stationary phase diagram Fig. 5 is generic in the sense that it is determined solely by the macroscopic flow-density relation. The number of lanes of the road, the distribution of individual optimal velocities, speed limits, and other details enter only in so far as they determine the exact values characterizing the flow-density relation for that particular road. We also note that throughout the paper we assumed the external conditions to vary slowly, so that the system has enough time to readjust to its new stationary state. Experimenting with different cellular traffic models in a real time scale shows that the typical time to reach a stationary state in a road segment of about 1.2 km is of the order of 3-5 min, which is reasonably small. Close to phase transitions lines, however, where the shock velocity vanishes, this time diverges and intrinsically non-stationary dynamic phenomena  take the lead.

In conclusion, we have shown that traffic data collected on German motorways provide evidence for a boundary-induced nonequilibrium phase transition of first order from the free flowing to the congested phase. The features of this phenomenon are readily understood in terms of the flow-density diagram. The dynamical mechanism leading to this transition is an interplay of shocks and local fluctuations caused by an on-ramp. Full investigation of a cellular automaton model for traffic flow reproduces this phase transition, but also exhibits a richer phase diagram with an interesting maximal flow phase. These results are not only important from the point of view of nonequilibrium physics, but also suggest new mechanisms of traffic control.

Acknowledgments: We thank Lutz Neubert for useful discussions and help in producing Figs. 3 and 4. L. S. acknowledges support from the Deutsche Forschungsgemeinschaft under Grant No. SA864/1-1.


text:
# Synchronization of Chaotic Systems by Common Random Forcing

## I Introduction

The issue of whether chaotic systems can be synchronized by common random noise sources has attracted much attention recentlymar94 ; mal96 ; san97 ; min98 ; per99 ; ali97 . It has been reported that for some chaotic maps, the introduction of the same (additive) noise in independent copies of the same map could lead (for large enough noise intensity) to a collapse onto the same trajectory, independently of the initial condition assigned to each of the copiesmar94 . This synchronization of chaotic systems by the addition of random terms is a remarkable and counterintuitive effect of noise. Nowadays, some contradictory results exist for the existence of this phenomenon of noise–induced synchronization. It is purpose of this paper to give explicit examples in which it is shown that one can indeed obtain such a synchronization. Moreover, the examples open the possibility to obtain such a synchronization in electronic circuits, hence suggesting that noise-induced synchronization of chaotic circuits can indeed be used for encryption purposes.

Although the issue of which is the effect of noise in chaotic systems was considered at the beginning of the 80’smt83 , to our knowledge the first atempt to synchronize two chaotic systems by using the same noise signal was considered by Maritan and Banavarmar94 . These authors analysed the logistic map in the presence of noise:

$$x_{n+1}=4x_n(1x_n)+\xi _n$$
(1)
where $`\xi _n`$ is the noise term, considered to be uniformly distributed in a symmetric interval $`[W,+W]`$. They showed that, if $`W`$ was large enough (i.e. for a large noise intensity) two different trajectories which started with different initial conditions but used otherwise the same sequence of random numbers, would eventually coincide into the same trajectory. This result was heavily critisized by Pikovskypik94  who argued that two systems can synchronize only if the largest Lyapunov exponent is negative. He then shows that the largest Lyapunov exponent of the logistic map in the presence of noise is always positive and concludes that the synchronization is, in fact, a numerical effect of lack of precision of the calculation. Furthermore, Malesciomal96  pointed out that the noise used to simulate Eq.(1) in mar94  was not really symmetric. This is because the requirement $`x_n(0,1),n`$, actually leads to discard those values for the random number $`\xi _n`$, which do not fulfill such condition. The average value of the random numbers which have been accepted is different from zero, hence producing an effective biased noise, i.e. one which does not have zero mean. The introduction of a non-zero mean noise means that we are altering essentially the properties of the deterministic map.

Noise induced synchronization has been since studied for other chaotic systems such as the Lorenz modelmar94 ; mal96  and the Chua circuitsan97 ; per99 . Synchronization of trajectories starting with different initial conditions but using otherwise the same sequence of random numbers was observed in the numerical integration of a Lorenz system in the presence of a noise distributed uniformly in the interval $`[0,W_L]`$, i.e. again a noise which does not have a mean of zeromar94 . Other detailed studiesmal96  also conclude that it is not possible to synchronize trajectories in a Lorenz system by adding an unbiased noise. Similarly, the studies of the Chua circuit always conclude that a biased noise is needed for synchronizationper99 . Therefore a widespread belief exists that it is not possible to synchronize two chaotic systems by injecting the same noisy signal to both of them. However, in this paper we give numerical evidence that it is possible to synchronize two chaotic systems by the addition of a common noise which is Gaussian distributed and not biased. We analyse specifically a 1-d map and the Lorenz system, both in the chaotic region. The necessary criterion introduced in ref.pik94  is fully confirmed and some heuristic arguments are given about the general validity of our results. Finally, we conclude with some open questions relating the general validity of our results.

## II Results

The first example is that of the map:

$$x_{n+1}=f(x_n)+ϵ\xi _n$$
(2)
where $`\xi _n`$ is a set of uncorrelated Gaussian variables of zero mean and variance 1. We use explicitely

$$f(x)=\mathrm{exp}\left[\left(\frac{x0.5}{\omega }\right)^2\right]$$
(3)
We plot in Fig.(1a) the bifurcation diagram of this map. In the noiseless case, we can see the typical windows in which the system behaves chaotically. The associated Lyapunov exponent, $`\lambda `$, in these regions is positive. For instance, for $`\omega =0.3`$ (the case we will be considereing throughout the paper) it is $`\lambda 0.53`$. In Fig.(1b) we observe that the Lyapunov exponent becomes negative for most values of $`\omega `$ for large enough noise level $`ϵ`$. Again for $`\omega =0.3`$ and now for $`ϵ=0.2`$ it is $`\lambda =0.17`$.

For the noiseless case, it is $`\lambda >0`$ and trajectories starting with differential initial conditions, obviously, remain different for all the iteration steps, see Fig.(2a). However, when moderated levels of noise ($`ϵ0.1`$) are used, $`\lambda `$ becomes negative and trajectories starting with different initial conditions, but using the same sequence of random numbers, synchronize perfectly, see Fig.(2b).

According to pik94 , convergence of trajectories to the same one, or loss of memory of the initial condition, can be stated as negativity of the Lyapunov exponent. The Lyapunov exponent of a map $`x_{n+1}=F(x_n)`$ is defined as

$$\lambda =\underset{N\mathrm{}}{lim}\frac{1}{N}\underset{i=1}{\overset{N}{}}\mathrm{ln}|F^{}(x_i)|$$
(4)
It is the average of (the logarithm of the absolute value of) the successive slopes $`F^{}`$ found by the trajectory. Slopes in $`[1,1]`$ contribute to $`\lambda `$ with negative values, indicating trajectory convergence. Larger or smaller slopes contribute with positive values, indicating trajectory divergence. Since the deterministic and noisy maps satisfy $`F^{}=f^{}`$ one is tempted to conclude that the Lyapunov exponent is not modified by the presence of noise. However, there is noise-dependence through the trajectory values $`x_i`$, $`i=1,\mathrm{},N`$. In the absence of noise, $`\lambda `$ is positive, indicating trajectory separation. When synchronization is observed, the Lyapunov exponent becomes negative, as required by the argument in pik94 .

By using the definition of the invariant measure on the attractor, or stationary probability distribution $`P_{st}(x)`$, the Lyapunov exponent can be calculated also as

$$\lambda =\mathrm{log}|F^{}(x)|=\mathrm{log}|f^{}(x)|P_{st}(x)\mathrm{log}|f^{}(x)|dx$$
(5)
Here we see clearly the two contributions to the Lyapunov exponent: although the derivative $`f^{}(x)`$ does not change when including noise in the trajectory, the stationary probability does change (see Fig.3), thus producing the observed change in the Lyapunov exponents. Synchronization, then, can be a general feature in maps which have a large region in which the derivative $`|f^{}(x)|`$ is smaller than one. Noise will be able, then , to explore that region and yield, on the average, a negative Lyapunov exponent.

The second system we have studied is the well known Lorenz model with random terms addedlorenz ; mar94 :

$`\dot{x}`$ $`=`$ $`p(yx)`$
$`\dot{y}`$ $`=`$ $`xz+rxy+ϵ\xi `$ (6)
$`\dot{z}`$ $`=`$ $`xybz`$
$`\xi `$ now is white noise: a Gaussian random process of mean zero and delta correlated, $`\xi (t)\xi (t^{})=\delta (tt^{})`$. We have used $`p=10`$, $`b=8/3`$ and $`r=28`$ which, in the deterministic case, $`ϵ=0`$ are known to lead to a chaotic behavior (the largest Lyapunov exponent is $`\lambda 0.9>0`$). We have integrated numerically the above equations using the Euler method with a time step $`\mathrm{\Delta }t=0.001`$. For the deterministic case, trajectories starting with different initial conditions are completely uncorrelated, see Fig. (4a). This is also the situations for small values of $`ϵ`$. However, when using a noise intensity $`ϵ=40`$ the noise is strong enough to induce synchronization of the trajectories. Again the presence of the noise terms makes the largest Lyapunov exponent become negative (for $`ϵ=40`$ it is $`\lambda 0.2`$). As in the example of the map, after some transient time, two different evolutions which have started in completely different initial conditions synchronize towards the same value of the three variables (see Fig. (4b) for the $`z`$ coordinate). One could argue that the intensity of the noise is very large. However, the basic structure of the “butterfly” Lorenz attractor remains unchanged as shown in Fig. (5).

In conclusion, we have shown that it is possible for noise to synchronize trajectories of a system which, deterministically, is chaotic. The novelty of our results is that the noise used in the two examples, a 1-d map and the Lorenz system, is unbiased, i.e. has always zero mean.

There still remain many open questions in this field. They involve the development of a general theory, probably based in the invariant measure, that could give us a criterion to determine the range of parameters (including noise levels) for which the Lyapunov exponent becomes negative, thus allowing synchronization. Another important question relates the structural stability of this phenomenon. Any practical realization of this system can not produce two identical samples. If one wants to use sthocastic synchronization of electronic emitters and receivers (as a means of encryption) one should be able to determine which is the allowed discrepancy between circuits before the lack of synchronization becomes unacceptable.

Acknowledgements We thank financial support from DGESIC (Spain) projects numbers PB94-1167 and PB97-0141-C02-01.

## Figure captions

FIGURE 1: (a) Bifurcation diagram of the map given by Eqs.(2) and (3) in the absence of noise terms. (b) Lyapunov exponent for the noiseless map ($`ϵ=0`$, continuous line) and the map with a noise intensity $`ϵ=0.1`$ (dotted line) and $`ϵ=0.2`$ (dot-dashed line).

FIGURE 2: Plot of two realizations $`x^{(1)}`$, $`x^{(2)}`$ of of the map given by Eqs. (2) and (3). Each realization consists of 10,000 points which have been obtained by iteration of the map starting in each case from a different initial condition (100,000 initial iterations have been discarded and are not shown). In figure (a) there is no noise, $`ϵ=0`$ and the trajectories are independent of each other. In figure (b) we have use a level of noise $`ϵ=0.2`$ producing a perfect synchronization (after discarding some initial iterations).

FIGURE 3: Plot of the stationary distribution for the map given by Eqs.(2) and (3) in the (a) deterministic case $`ϵ=0`$, and (b) the case with noise along the trajectory, $`ϵ=0.2`$.

FIGURE 4: Same than figure (1) for the $`z`$ variable of the Lorenz system, Eqs.(6) in the (a) deterministic case $`ϵ=0`$ and (b) $`ϵ=40`$. Notice the perfect synchronization in case (b).

“Butterfly” attractor of the Lorenz system in the cases (a) of no noise $`ϵ=0`$ and (b) $`ϵ=40`$.


text:
# QSO Absorption Line Constraints on Intragroup High–Velocity Clouds

## 1. Introduction

The origin(s) of high–velocity clouds (HVCs), gaseous material that departs from the Galactic rotation law by more than 100 km s-1, is a topic under debate. Undoubtedly, some HVCs arise from tidal streams (e.g. the Magellanic Stream), and from fountain processes local to the Galaxy (Wakker & van Woerden (1997)). Recently, however, the hypothesis that most HVCs are distributed ubiquitously throughout the Local Group and are relics of group formation has returned to favor (Blitz et al. (1999); Braun & Burton (1999)).

In the intragroup HVC hypothesis (1) the cloud kinematics follow the Local Group standard of rest (LGSR), not the Galactic standard of rest (GSR), with the exception of some HVCs related to tidal stripping or Galactic fountains (Blitz et al. (1999); Braun & Burton (1999)); (2) the cloud Galactocentric distances are typically 1 Mpc; (3) the extended HVC cloud complexes are presently accreting onto the Milky Way; (4) the clouds are local analogs of the Lyman limit absorbers observed in quasar spectra; (5) the clouds have masses of $`10^7`$ M and greater; and (6) the metallicities are lower than expected if the material originated from blowout or fountains from the Milky Way (Wakker et al. (1999); Bowen & Blades (1993)).

Blitz et al. (1999; hereafter BSTHB) suggest that there are $`300`$ clouds above the $`21`$–cm $`N(\text{H}\text{i})`$ detection threshold of $`2\times 10^{18}`$ cm-2. These clouds have radii $`15`$ kpc and are ubiquitous throughout the group.

Braun and Burton (1999, hereafter BB) cataloged $`65`$ Local Group CHVCs, which represent a homogeneous subset of the HVC population discussed by BSTHB. High resolution $`21`$–cm observations (Braun & Burton (2000)) show that the CHVCs, have compact, $`N(\text{H}\text{i})>10^{19}`$ cm-2, cold cores with few km s-1 FWHM surrounded by extended “halos” with FWHM $`25`$ km s-1. The typical radius is $`5`$$`8`$ kpc at the estimated distance of $`700`$ kpc. The BB sample is homogeneous but is not complete; they estimate that there could be as many as $`200`$ Local Group CHVCs.

Recently, Zwaan and Briggs (2000) reported evidence in contradiction of the intragroup hypothesis. In a blind Hi $`21`$–cm survey of extragalactic groups, sensitive to $`N(\text{H}\text{i})10^{18}`$ cm-2 (capable of detecting $`10^7`$ M Hi clouds), they failed to locate any extragalactic counterparts of the Local Group HVCs. This is in remarkable contrast to the numbers predicted. If intragroup HVCs exist around all galaxies or galaxy groups, and the Hi mass function is the same in extragalactic groups as measured locally, then Zwaan & Briggs should have detected $`70`$ in groups and $`250`$ around galaxies ($`10`$ and $`40`$ for the CHVCs, respectively).

Thus, the Zwaan and Briggs result is in conflict with the intragroup HVC hypothesis. Since the hypothesized intragroup clouds are remnants of galaxy formation and are shown to be stable against destruction mechanisms (BSTHB), they are predicted to form at very high redshifts and to be ubiquitous in galaxy groups to the present epoch. In this Letter, we argue that the version of the intragroup HVC hypothesis presented by BSTHB is also in conflict with the observed redshift number density of moderate redshift $`(z0.5)`$ Mgii and Lyman limit (LLS) quasar absorption line systems. We also find that the properties of the extragalactic analogs of the BB CHVCs are severely constrained.

In general, the redshift number density of a non–evolving population of objects, to be interpreted as the number per unit redshift, is written

$$\frac{dN}{dz}=C_f\frac{n\sigma c}{H_0}\left(1+z\right)\left(1+2q_0z\right)^{1/2},$$
(1)
where $`C_f`$ is the covering factor, $`n`$ is the number density of absorbing structures, and $`\sigma `$ and $`C_f`$ are the surface area presented by each structure and its covering factor for detectable absorption. Throughout, we use $`H_0=100`$ km s-1 Mpc-1 and $`q_0=0.5`$, which gives $`dN/dz(1+z)^{1/2}`$.

## 2. Mg II Systems

The statistics of Mgii absorbers are well–established at $`0.3z2.2`$. For rest–frame equivalent widths of $`W(\text{Mg}\text{ii})>0.3`$ Å (“strong” Mgii absorption) Steidel & Sargent (1992) found $`dN/dz=0.8\pm 0.2`$ for $`z0.5`$ with a redshift dependence consistent with no evolution expectations. Normal, bright ($`L0.1L^{}`$) galaxies are almost always found within $`40`$ kpc of strong Mgii absorbers (Bergeron & Boissé (1991); Bergeron et al. (1992); Le Brun et al. (1993); Steidel, Dickinson, & Persson (1994); Steidel (1995); Steidel et al. (1997)). From the Steidel, Dickinson, and Persson survey, all but $`3`$ of $`58`$ strong MgII absorbers, detected toward $`51`$ quasars, have identified galaxies with a coincident redshift within that impact parameter (sky projected separation from the quasar line of sight) (see Charlton & Churchill (1996)). Also, it is rare to observe a galaxy with an impact parameter less than $`40h^1`$ kpc that does not give rise to Mgii absorption with $`W(\text{Mg}\text{ii})>0.3`$ Å (Steidel (1995)). In $`25`$ “control fields” of quasars, without observed strong Mgii absorption in their spectra, only two galaxies had impact parameters less than $`40h^1`$ kpc (see also Charlton & Churchill (1996)). As such, the regions within $`40h^1`$ kpc of typical galaxies account for the vast majority of Mgii absorbers above this equivalent width threshold; there is nearly a “one–to–one” correspondence. If we accept these results, it would imply that there is little room for a contribution to $`dN/dz`$ from a population of intragroup clouds which have impact parameters much greater than $`40h^1`$ kpc.

However, the predicted cross section for Mgii absorption from the extragalactic intragroup clouds analogous to HVCs would be substantial. We quantify the overprediction of the redshift path density by computing the ratio of $`dN/dz`$ of the intragroup clouds to that of Mgii absorbing galaxies,

$$\frac{(dN/dz)_{cl}}{(dN/dz)_{gal}}=F\left(\frac{f_{cl}}{f_{gal}}\right)\left(\frac{N_{cl}}{N_{gal}}\right)\left(\frac{R_{cl}}{R_{gal}}\right)^2,$$
(2)
where $`F`$ is the fraction of Mgii absorbing galaxies that reside in groups having intragroup HVC–like clouds, $`f_{cl}`$ is the fraction of the area of the clouds and $`f_{gal}`$ is the fraction of the area of the galaxies that would produce $`W(\text{Mg}\text{ii})>0.3`$ Å along the line of sight, and $`N_{cl}`$ and $`N_{gal}`$ are the number of clouds and galaxies per group, respectively. The cross section of the group times the intragroup cloud covering factor, $`C_f\pi R_{gr}^2`$, is equal to $`N_{cl}\pi R_{cl}^2`$. The total predicted $`dN/dz`$ for Mgii absorbers with $`W(\text{Mg}\text{ii})>0.3`$ Å is then,

$$\left(\frac{dN}{dz}\right)_{tot}=\left(\frac{dN}{dz}\right)_{gal}\left[1+\frac{(dN/dz)_{cl}}{(dN/dz)_{gal}}\right].$$
(3)
If virtually all Mgii absorbers are accounted for by galaxies, it is required that $`(dN/dz)_{tot}(dN/dz)_{gal}`$; the left hand side of Equation 2 must be very close to zero.

In the BSTHB version of the intragroup HVC model, the “best” expected values are $`N_{cl}=300`$ and $`R_{cl}=15`$ kpc (BSTHB; Blitz 2000, private communication ); if we take $`R_{gal}=40`$ kpc and $`f_{gal}=1`$ (Steidel (1995)), and assuming $`N_{gal}=4`$, we find that the covering factor for Mgii absorption from extragalactic analogs to the Local Group HVCs would exceed that from galaxies by a factor of $`10`$ for $`F=1`$ and $`f_{cl}=1`$, giving $`(dN/dz)_{tot}9`$. More recently, Blitz and Robinshaw (2000) have suggested that sizes may be smaller ($`R_{cl}=8`$ kpc) when beam–smearing is considered. Considering this as an indication of uncertainties in the BSTHB parameters, and considering $`2<N_{gal}<6`$ for the typical number of group galaxies, we find, for $`F=1`$ and $`f_{cl}=1`$, a range $`2<(dN/dz)_{cl}/(dN/dz)_{gal}<21`$. This corresponds to $`2.6<(dN/dz)_{tot}<17.6`$. It is unlikely that $`F`$ is significantly less than unity; the majority of galaxies reside in groups like the Local Group that would have HVC analogs. In order that $`(dN/dz)_{tot}(dN/dz)_{gal}`$, $`f_{cl}0.2`$ is required.

It is not clear what fraction $`f_{cl}`$ of HVCs with $`N(\text{H}\text{i})`$ above the $`10^{18}`$ cm-2 detection threshold will give rise to $`W(\text{Mg}\text{ii})0.3`$ Å because the equivalent width is sensitive to the metallicity and internal velocity dispersion of the clouds. Based upon Cloudy (Ferland (1996)) photoionization equilibrium models, a cloud with $`N(\text{H}\text{i})10^{18}`$ cm-2, subject to the ionizing metagalactic background (Haardt & Madau (1996)), would give rise to Mgii absorption with $`N(\text{Mg}\text{ii})10^{14}N_{18}(Z/Z_{})`$ cm-2, where $`N_{18}`$ is the Hi column density in units of $`10^{18}`$ cm-2 and $`Z/Z_{}`$ is the metallicity in solar units. For optically thick clouds, those with $`N(\text{H}\text{i})10^{17.5}`$ cm-2, this result is insensitive to the assumed ionization parameter111The ionization parameter is the ratio of the number density of hydrogen ionizing photons to the number density of electrons, $`n_\gamma /n_e`$. over the range $`10^{4.5}`$ to $`10^{1.5}`$.

BSTHB expect HVCs to have $`Z/Z_{}0.1`$. For $`N_{18}=2`$ and $`Z/Z_{}=0.1`$, clouds with internal velocity dispersions of $`\sigma _{cl}20`$ km s-1 (Doppler $`b28`$ km s-1) give rise to $`W(\text{Mg}\text{ii})0.5`$ Å. For $`\sigma _{cl}=10`$ km s-1 ($`b=14`$ km s-1), $`W(\text{Mg}\text{ii})=0.3`$ Å. The CHVC “halos” typically have FWHM of $`29`$$`34`$ km s-1, which corresponds to $`\sigma _{cl}12`$$`14`$ km s-1 (Braun & Burton (1999)). Thus it appears that most lines of sight through the BSTHB extragalactic analogs will produce strong Mgii absorption. Certainly $`f_{cl}>0.2`$, so there is a serious discrepancy between the predicted $`(dN/dz)_{tot}`$ and the observed value.

However, if the intragroup clouds have lower metallicities, this would result in smaller $`W(\text{Mg}\text{ii})`$. Unfortunately, there has only been one metallicity estimate published for an HVC, which may or may not be related to the Galaxy. Braun and Burton (2000) estimate that CHVC 125+41-207, with $`W(\text{Mg}\text{ii})=0.15`$ Å, has a metallicity of $`0.04<Z/Z_{}<0.07`$, however this is quite uncertain because of the effects of beam smearing on measuring the $`N(\text{H}\text{i})`$ value. Because of the uncertainties, we simply state that a population of low metallicity clouds could reduce the discrepancy between the predicted redshift density for intragroup clouds, $`(dN/dz)_{cl}`$, and the observed value of $`(dN/dz)_{tot}`$. However, then the expected number of smaller $`W(\text{Mg}\text{ii})`$ systems to arise from intragroup clouds would be increased.

The observed Mgii equivalent width distribution rises rapidly below $`0.3`$ Å (“weak” Mgii absorbers), such that $`dN/dz=2.2\pm 0.3`$ for $`W(\text{Mg}\text{ii})>0.02`$ Å at $`z=0.5`$ (Churchill et al. (1999)). To this equivalent width limit, Mgii absorption could be observed from intragroup HVCs with $`N_{18}=2`$ and metallicities as low as $`Z/Z_{}=0.0025`$ \[for $`N(\text{Mg}\text{ii})=10^{11.7}`$ cm-2, $`W(\text{Mg}\text{ii})`$ is independent of $`\sigma _{cl}`$\]. However, almost all ($`9`$ out of a sample of $`10`$) Mgii absorbers with $`W(\text{Mg}\text{ii})<0.3`$ Å do not have associated Lyman limit breaks (Churchill et al. 2000a ); that is, their $`N(\text{H}\text{i})`$ is more than a decade below the sensitivity of $`21`$–cm surveys. Thus, based upon available data, roughly $`90`$% of the “weak” Mgii absorbers do not have the properties of HVCs, and therefore are not analogous to the clouds invoked for the intragroup HVC scenario. If $`10`$% of the weak Mgii absorbers are analogs to the intragroup HVCs, they would contribute an additional $`0.20`$ to $`(dN/dz)_{cl}`$.

Since the BB CHVC extragalactic analogs have smaller cross sections, we should separately consider whether they would produce a discrepancy with the observed Mgii absorption statistics. BB observed $`N_{cl}=65`$ and inferred a typical $`R_{cl}=5`$$`8`$ kpc for the CHVCs, however a complete sample might have $`N_{cl}=200`$. Assuming $`N_{gal}=2`$$`6`$, $`R_{gal}=40`$ kpc, $`f_{gal}=1`$, and $`F=1`$, for the BB subsample of the HVC population, we obtain $`0.17f_{cl}<(dN/dz)_{cl}/(dN/dz)_{gal}<4.0f_{cl}`$.

The cores of the CHVCs have $`N(\text{H}\text{i})>10^{19}`$ cm-2 and they occupy only about $`15`$% of the detected extent. For $`Z/Z_{}>0.01`$ and $`\sigma _{cl}=10`$ km s-1, these cores can produce $`W(\text{Mg}\text{ii})0.3`$ Å over their full area. It follows that $`f_{cl}=0.15`$, which yields $`0.025<(dN/dz)_{cl}/(dN/dz)_{gal}<0.6`$. Depending on the specific parameters, there may or may not be a conflict with the observed $`(dN/dz)_{tot}`$ for strong Mgii absorption.

The “halos” of the CHVCs have $`N(\text{H}\text{i})>10^{18}`$ cm-2 and, as discussed above, would produce weak Mgii absorption for $`Z/Z_{}>0.005`$ over most of the cloud area. This implies contribution to $`(dN/dz)`$ from BB CHVC analogs that is in the range $`0.14<(dN/dz)_{cl}<3.2`$. If the number were at the high end of this range, the cross section would be comparable to the observed $`(dN/dz)`$ for weak Mgii absorbers at $`z=0.5`$. However, as noted above when considering the BSTHB scenario, there is a serious discrepancy. Only $`10`$% of the weak Mgii absorbers show a Lyman limit break, so extragalactic analogs of the BB CHVCs can only be a fraction of the weak Mgii population. Regions of CHVCs at larger radii, with $`N(\text{H}\text{i})`$ below the threshold of present $`21`$–cm observations, are constrained to have $`Z/Z_{}0.01`$ in order that they do not produce a much larger population of weak Mgii absorbers with Lyman limit break than is observed.

## 3. Lyman Limit Systems

The redshift number density of LLS also places strong constraints on intragroup environments that give rise to Lyman breaks in quasar spectra. This argument is not sensitive to the assumed cloud velocity dispersion and/or metallicity.

Statistically, $`dN/dz`$ for Mgii systems is consistent (1 $`\sigma `$) with $`dN/dz`$ for LLS. At $`z0.5`$, LLS have $`dN/dz=0.5\pm 0.3`$ (Stengler–Larrea et al. (1995)) and Mgii systems have $`dN/dz=0.8\pm 0.2`$ (Steidel & Sargent (1992)). Churchill et al. (2000a) found a Lyman limit break \[i.e. $`N(\text{H}\text{i})10^{16.8}`$ cm-2\] for each system in a sample of ten having $`W(\text{Mg}\text{ii})>0.3`$ Å. LLS and Mgii absorbers have roughly the same redshift number density and therefore Mgii–LLS absorption must almost always arise within $`40`$ kpc of galaxies (Steidel (1993)). As such, there is little latitude for a substantial contribution to $`dN/dz`$ from intragroup Lyman limit clouds.

Using Equation 1, we could estimate this contribution by considering the volume density of galaxy groups and the cross section for HVC Lyman limit absorption in each. However, the volume density of groups is not well measured, particularly out to $`z=0.5`$.

Instead, we make a restrictive argument based upon a comparison between the cross sections for Lyman limit absorption of $`L^{}`$ galaxies and for HVCs in a typical group (similar to the discussion of Mgii absorbers in § 2). Again, we simply compare the values of $`C_f`$ for the different populations of objects in a typical group. The covering factor for HVCs within the group is

$$C_f=N_{cl}\frac{\left(R_{cl}\right)^2}{\left(R_{gr}\right)^2}.$$
(4)
The best estimate for the BSTHB version of the intragroup HVC model, with $`N_{cl}=300`$, $`R_{cl}=15`$ kpc, and a group radius $`R_{gr}=1.5`$ Mpc, gives $`C_f=0.03`$ for $`N(\text{H}\text{i})2\times 10^{18}`$ cm-2. If instead we use the BB number of observed CHVCs, $`N_{cl}=65`$, and $`R_{cl}=5`$$`8`$ kpc, we obtain a much smaller number, $`0.0007<C_f<0.0018`$. However, if the BB sample is corrected for incompleteness such that $`N_{cl}=200`$, these numbers increase so that $`0.002<C_f<0.006`$.

In comparison, a typical group with $`4`$ $`L^{}`$ galaxies, each with a Lyman limit absorption cross section of $`R_{cl}40`$ kpc, would have $`C_f=0.002`$. If they existed with the properties discussed, the extragalactic analogs to the BSTHB HVCs would dominate the contribution of $`L^{}`$ galaxies to the $`dN/dz`$ of LLS by at least a factor of $`15`$, and this is only considering HVC regions with $`N(\text{H}\text{i})>2\times 10^{18}`$ cm-2 that are detected in the $`21`$–cm surveys. Any extensions in area below this threshold value \[down to $`N(\text{H}\text{i})5\times 10^{16}`$ cm-2\] would worsen the discrepancy. As such, the BSTHB hypothesis is definitively ruled out.

For regions of BB CHVCs with $`N(\text{H}\text{i})>2\times 10^{18}`$ cm-2, the covering factor ranges from $`C_f=0.0007`$ to $`C_f=0.006`$ depending on assumed sizes and corrections for incompleteness. This ranges from $`35`$$`300`$% of the cross section for the $`L^{}`$ galaxies. The total observed $`dN/dz`$ for Lyman limit absorption (down to $`\mathrm{log}N(\text{H}\text{i})=17`$ cm-2) is only $`0.5`$, even a $`35`$% contribution to the Lyman limit cross section from HVCs that are separate from galaxies creates a discrepancy. This would imply that the result that most lines of sight within $`40`$ kpc of a typical $`L^{}`$ galaxy produce Lyman limit absorption is incorrect. This would further imply that there is a substantial population of strong Mgii absorbers without Lyman limit breaks (to account for $`dN/dz=0.8`$ for strong Mgii absorption) or of strong Mgii absorbers not associated with galaxies. Both types of objects are rarely observed (Churchill et al. 2000a ; Bergeron & Boissé (1991); Bergeron et al. (1992); Le Brun et al. (1993); Steidel, Dickinson, & Persson (1994); Steidel (1995); Steidel et al. (1997)). Furthermore, $`C_f=0.002`$ for BB CHVCs only takes into account the fraction of the BB CHVC areas with $`N(\text{H}\text{i})>2\times 10^{18}`$ cm-2. Therefore, the extended “halos” around the CHVCs are also constrained not to contribute substantial cross section for Lyman limit absorption along extragalactic lines of sight.

## 4. Summary

We have made straight–forward estimates of the predicted redshift number density at $`z0.5`$ of Mgii and LLS absorption from hypothetical extragalactic analogs to intragroup HVCs as expected by extrapolating from the BSTHB and BB Local Group samples. We find that it is difficult to reconcile the intragroup hypothesis for HVCs with the observed $`dN/dz`$ of Mgii and LLS systems.

The discrepancy between the $`dN/dz`$ of Mgii–LLS absorbers and the observed covering factor of “intragroup” HVCs could be reduced if the HVCs have a clumpy structure. Such structure would result in Mgii–LLS absorption observable only in some fraction, $`f_{los}`$, of the lines of sight through the cloud. Effectively, this reduces the covering factor for Lyman limit absorption, or the value of $`f_{cl}`$ in equation (2) for Mgii absorbers. Considering beam smearing in $`21`$–cm surveys, substructures would be detected above a $`N(\text{H}\text{i})>2\times 10^{18}`$ cm-2 $`21`$–cm detection threshold if their column densities were $`N(\text{H}\text{i})_{sub}>2\times 10^{18}/f_{los}`$ cm-2. The predicted $`dN/dz`$ for HVC–like clouds could be reduced by a factor of ten if $`f_{los}0.1`$, giving $`N(\text{H}\text{i})_{sub}>2\times 10^{19}`$ cm-2. All the gas outside these higher Hi column density substructures would need to be below the Lyman limit or the arguments in § 3 would hold. It is difficult to reconcile such a density distribution with the high resolution observations of BB CHVCs which show diffuse halos around the core concentrations (Braun & Burton (2000)), but these ideas merit further consideration.

### 4.1. The BSTHB Scenario

We conclude that the predicted $`dN/dz`$ from the hypothetical population of intragroup HVCs along extragalactic sight lines to quasars from the BSTHB scenario would exceed:

1) the $`dN/dz`$ of Mgii absorbers with $`W(\text{Mg}\text{ii})0.3`$ Å. This class of absorber is already known to arise within $`40h^1`$ kpc of normal, bright galaxies (Bergeron & Boissé (1991); Bergeron et al. (1992); Le Brun et al. (1993); Steidel, Dickinson, & Persson (1994); Steidel (1995); Steidel et al. (1997)).

2) the $`dN/dz`$ of “weak” Mgii absorbers with $`0.02<W(\text{Mg}\text{ii})<0.3`$ Å absorption. In principle, weak Mgii absorption could arise from low metallicity, $`0.005Z/Z_{}<0.1`$, intragroup HVCs. However, the majority of observed weak systems are already known to be higher metallicity, $`Z/Z_{}0.1`$, sub–Lyman limit systems (Churchill et al. (1999); Rigby et al. (2000)).

3) the $`dN/dz`$ of Lyman limit systems. These would be produced by all extragalactic BSTHB HVC analogs regardless of metallicity. However, most Lyman limit systems are seen to arise within $`40h^1`$ kpc of luminous galaxies (Steidel (1993); Churchill et al. 2000b ).

These points do not preclude a population of infalling intragroup clouds which do not present a significant cross section for $`21`$–cm absorption, as predicted by CDM models (Klypin et al. (1999); Moore et al. (1999)). In fact, such intragroup objects could be related to sub–Lyman limit weak Mgii absorbers (Rigby et al. (2000)).

### 4.2. The BB Scenario

The properties of the BB CHVC population are also significantly constrained by Mgii and Lyman limit absorber statistics:

1) They could produce $`W(\text{Mg}\text{ii})0.3`$ Å in excess of what is observed if a large incompleteness correction is applied (i.e. so that $`N_{cl}=200`$), or if relatively large sizes ($`R_{cl}8`$ kpc) are assumed.

2) They would be expected to contribute to the $`dN/dz`$ of weak \[$`W(\text{Mg}\text{ii})>0.02`$ Å\] Mgii absorption. However, based upon observations (Churchill et al. 2000a ), only $`10`$% of the population of weak Mgii absorbers have Lyman–limit breaks. Therefore, only a small fraction of weak Mgii absorption could arise in extragalactic BB CHVC analogs.

3) The $`dN/dz`$ for Lyman limit absorption from the hypothesized BB CHVC population could be a significant fraction, or comparable to that expected from the local environments of $`L^{}`$ galaxies (within $`40`$ kpc); the observed value is already consistent with that produced by the galaxies.

4) The CHVCs are observed to have a cool core with $`N(\text{H}\text{i})>10^{19}`$ cm-2, surrounded by a halo which typically extends to $`R_{cl}5`$ kpc. It is natural to expect that the Hi extends out to larger radii at smaller $`N(\text{H}\text{i})`$ and should produce a Lyman limit break out to the radius at which $`N(\text{H}\text{i})<10^{16.8}`$ cm-2. Although there is expected to be a sharp edge to the Hi disk at $`N(\text{H}\text{i})10^{17.5}`$ or $`10^{18}`$ cm-2 (Maloney (1993); Corbelli & Salpeter (1994); Dove & Shull (1994)), physically we would expect that this edge would level off at $`10^{17.5}`$ cm-2, such that a significant cross section would be presented at $`10^{16.8}<N(\text{H}\text{i})<10^{17.5}`$ cm-2. Another possibility is that there is a sharp cutoff of the structure at $`N(\text{H}\text{i})2\times 10^{18}`$ cm-2, but this is contrived.

## 5. Conclusion

We are forced to the conclusion that there can only be a limited number of extragalactic infalling group HVC analogs at $`z0.5`$. Future data could force a reevaluation of the relationships between galaxies, Lyman limit systems, and Mgii absorbers, but it seems unlikely that the more serious inconsistencies we have identified could be reconciled in this way. A clumpy distribution of Hi could be constructed that would reduce the discrepancy, but would require very diffuse material (below the Lyman limit) around dense cores. Evolution in the population of HVCs is another possibility. If the extragalactic background radiation declined from $`z=0.5`$ to the present, the clouds would have been more ionized in the past, and therefore would have had smaller cross sections at a given $`N(\text{H}\text{i})`$. However, this does not explain why Zwaan and Briggs (2000) do not see the $`z=0`$ extragalactic analogs to the HVCS or CHVCs. Our results are entirely consistent with theirs, and the implications are the same: the discrepancies between the Local Group HVC population and the statistics of Mgii and Lyman limit absorbers can only be reconciled if most of the extragalactic HVC analogs are within $`100`$$`200`$ kpc of galaxies, and not at large throughout the groups.

We thank L. Blitz, J. Bregman, J. Mulchaey, B. Savage, K. Sembach, T. Tripp, and especially B. Wakker and Buell Jannuzi, and our referees for stimulating discussions and comments. Support for this work was provided by NSF grant AST–9617185 (J. R. R. was supported by an REU supplement) and by NASA grant NAG 5–6399.


text:
# From waves to avalanches: two different mechanisms of sandpile dynamics

## Abstract

Time series resulting from wave decomposition show the existence of different correlation patterns for avalanche dynamics. For the $`d=2`$ Bak-Tang-Wiesenfeld model, long range correlations determine a modification of the wave size distribution under coarse graining in time, and multifractal scaling for avalanches. In the Manna model, the distribution of avalanches coincides with that of waves, which are uncorrelated and obey finite size scaling, a result expected also for the $`d=3`$ Bak et al. model.

PACS numbers:05.65.+b,05.40.-a,45.70.Ht,64.60.Ak

Full information on a self–organized critical (SOC) process  is contained in the time series, if the time step is the most microscopic conceivable. The self-similarity of such a process, due to its intermittent, avalanche character, should be revealed by the scaling of the time autocorrelation or its power spectrum. In spite of this, time series analyses have been seldom performed on sandpiles or similar systems, and mostly concerned the response to finite, random external disturbances, i.e. the problem of $`1/f`$ noise . Most efforts in the characterization of SOC scaling concentrated on probability distribution functions (PDF’s) of global properties of the avalanche events, which occupy large intervals of the microscopic evolution time . The numerical analysis of such PDF’s is often difficult, and universality issues can not be easily solved. The situation is even more problematic if, like for the two-dimensional (2D) Bak-Tang-Wiesenfeld sandpile (BTW) , the usually assumed finite size scaling (FSS) form reveals inadequate for the PDF’s, and needs to be replaced by a multifractal one . These findings rise the additional issue of why the 2D BTW displays such multifractal scaling, while apparently very similar sandpiles, like the Manna model (M) , do not .

In recent theoretical approaches to the BTW and other Abelian sandpiles , a major role has been played by the waves of toppling into which avalanches can be decomposed . For the BTW, the PDF’s of waves, as sampled from a large collection of successive avalanches, obey FSS with exactly known exponents . By analyzing the succession of waves as a stationary time series, one could hope to determine the statistical properties of avalanches. However, so far, no precise information on the correlations of such series was obtained.

In this Letter we generalize the wave description to the M in 2D. The study of the respective wave time series reveals that BTW and M are prototypes of two very different scenarios for avalanche dynamics. In the M case, successive wave sizes are totally uncorrelated. As a consequence, avalanche and wave PDF’s have identical scaling properties, consistent with FSS. For the BTW, to the contrary, wave sizes show long range correlations and persistency in time. PDF’s of “block variables” given by sums of $`n`$ successive wave sizes show that these correlations are responsible for the fact that avalanche scaling differs from the wave one and has multifractal features. For the 3D BTW, on the other hand, our results suggest validity of a scenario identical to the M one, and lead to conjecture exact avalanche exponents coinciding with those of the wave PDF.

The 2D BTW  is defined on a square $`L\times L`$ lattice. An integer $`z_i`$, the number of grains, is assigned to site $`i`$. Starting from a stable configuration ($`z_iz_c=4`$, $`i`$), a grain is added at a randomly chosen site; after each addition, all sites exceeding the stability threshold, $`z_k>z_c`$, undergo toppling, distributing one grain to each one of the nearest neighbors. The topplings, which dissipate grains when occurring at the edges, continue until all sites are stable, and a new grain is added. The $`s`$ topplings between two consecutive additions form an avalanche. After many additions, the system organizes into a stationary critical state.

Manna  studied a two–state version, M, of the sandpile. The sites can be either empty or occupied; grains are added randomly, and when one of them drops onto an occupied site, a “hard core” repulsion pushes two particles out to randomly chosen nearest neighbors. Compared to the BTW, in which toppling is deterministic, this model has an extra stochastic ingredient in the microscopic evolution, and $`z_c=1`$.

One can define a special ordering in the topplings of the BTW by introducing the wave decomposition of avalanches . After the site of addition, $`O`$, topples, any other unstable site is allowed to topple except possibly $`O`$. This is the first wave of toppling. If $`O`$ is still unstable, it is allowed to topple once again, leading to the second wave. This continues until $`O`$ becomes stable. Thus, an avalanche is broken into a sequence of waves. During a wave, sites topple only once, and for an $`m`$-wave avalanche $`s=_{k=1}^ms_k`$, where $`s_k`$ are the topplings of the $`k`$-th wave. Following the definitions for the BTW, we implement here a wave decomposition of avalanches also for the M. Unlike for the BTW, sites involved in a wave may topple more than once. Furthermore, the toppling order chosen implies now a peculiar sequence of stable configurations visited upon addition of grains, but the realization probabilities of possible configurations are independent of this order .

The BTW wave size PDF has FSS form $`P_w(s)s^{\tau _w}p_w(s/L^{D_w})`$, with $`\tau _w=1`$, $`D_w=2`$ in 2D, and $`p_w`$ a suitable scaling function . For the M waves we also found a PDF obeying FSS, with $`\tau _w=1.31\pm 0.02`$ and $`D_w=2.75\pm 0.01`$. These exponents are remarkably consistent with those estimated for avalanches . Such a coincidence certainly does not apply to the BTW. Attempts to derive exact BTW avalanche exponents within FSS were based on the observation that often waves within an avalanche show a final, long contraction phase, and on a scaling assumption for the corresponding $`s_ks_{k+1}`$ . A more adequate ansatz for the conditional PDF of $`s_{k+1}`$ given $`s_k`$ , and Markovianity assumptions, did not help in better characterizing avalanche scaling along these lines . In fact the 2D BTW obeys a multifractal form of scaling, which is not catched by such simplified schemes .

The wave time series $`\{s_k\}`$ provide coarse grained dynamical descriptions. In the $`L\mathrm{}`$ limit, these descriptions are infinitely rescaled with respect to those at microscopic time scale, but still infinitely finer than the mere records of successive avalanche sizes. This intermediate time scale reveals essential in order to understand the dynamics indside avalanches, whose size sequence we found to be uncorrelated in the sense discussed below for waves. The microscopic scale gave no significant results in comparing the two models, since, at that level, similar strong correlations exists in both of them, due to the parallel updating algorithm.

We determined for BTW and M the autocorrelation

$$C(t,L)=\frac{s_{k+t}s_k_L\mu ^2}{s_k^2_L\mu ^2},$$
(1)
with $`t=1,2,\mathrm{}`$, and $`\mu =s_k_L`$, the time averages being taken over up to $`10^7`$ waves for $`L=128,256,512,1024`$ and $`2048`$. A first, striking result is that, as soon as $`t>0`$ ($`C(0,L)=1`$ by normalization), waves are uncorrelated in the M. Indeed, as $`L`$ grows, $`C`$ manifestly approaches $`0`$ as soon as $`t>0`$ (Fig. 1). To the contrary, $`C`$ is long range for the BTW, because it approaches $`0`$ only for $`t`$ exceeding the maximum number of waves in an avalanche (Fig. 1), which we found to scale $`L`$, for $`L\mathrm{}`$. For the BTW we further tested the FSS form

$$C(t,L)=t^{\tau _c}g(t/L^{D_c})\text{(}\text{t,L≫1}\text{)},$$
(2)
on the basis of the $`L\mathrm{}`$ scaling of the moments

$$t^q_L=\underset{t}{}C(t,L)t^qL^{\sigma _c(q)}.$$
(3)
FSS would imply the piecewise linear form

$$\sigma _c(q)=\{\begin{array}{cc}D_c(q\tau _c+1)\hfill & \text{if }q\tau _c1,\hfill \\ 0\hfill & \text{if }q<\tau _c1.\hfill \end{array}$$
(4)
Fig. 2 shows the extrapolated $`\sigma _c(q)`$, which has an approximately linear part for $`1q4`$, consistent with $`D_c=1.02\pm 0.05`$ and $`\tau _c=0.40\pm 0.05`$. The curvature for $`q<1`$ is due to the fact that, for finite $`L`$, a logarithm ($`\tau _c=0`$) can not be easily distinguished from a power law with $`\tau _c0`$. The inset of Fig. 2 shows this logarithmic growth of $`t^q_L`$. We conclude that, for the BTW, a simple power law tail $`C(t,\mathrm{})t^{\tau _c}`$, is a first, rough, approximation. The increment $`y(t)=_{k=1}^ts_k`$ is comparable to the trail of a fractional Brownian motion with Hurst exponent $`H=(2\tau _c)/2`$ , such that $`H=0.80\pm 0.03`$ for the BTW. $`H`$ can be measured directly from the fluctuation

$$F(t,L)=\left[\mathrm{\Delta }y(t)^2_L\mathrm{\Delta }y(t)_{L}^{}{}_{}{}^{2}\right]^{1/2},$$
(5)
with $`\mathrm{\Delta }y(t)=y(k+t)y(k)`$ , which should scale as $`F(t,\mathrm{})t^H`$. Fig. 3 reports tests of this scaling for $`L=2048`$. For the BTW, $`H0.85`$ at low $`t`$, in agreement with $`\tau _c=0.40\pm 0.05`$. A crossover to $`H=1/2`$ is observed for large $`t`$. The crossover time of course increases with $`L`$. $`H=1/2`$ corresponds to a process with $`C`$ exponentially decaying with $`t`$, or with $`C(t,\mathrm{})=\delta _{t,0}`$, as is the case for the M. Thus, the crossover is due to the fact that, beyond the maximal time duration of avalanches, waves are uncorrelated. $`H1/2`$ implies long range correlations, as we find for the BTW. Furthermore, $`H>1/2`$ corresponds to persistency: an increasing or decreasing trend of $`y`$ in the past mostly implies a similar tendency in the future. This accounts for the observed expansion and contraction phases in avalanche growth . The above features of $`C`$ should be responsible for the peculiar long range on–site correlations of the noise expected when mapping the BTW onto a discrete interface growth equation .

The BTW has been shown to display a non constant gap in the high $`q`$ moments of some avalanche PDF’s . Thus, it should not surprise if, unlike assumed in Eq. (4), also $`C`$ would manifest similar multifractal properties. Indeed, a more accurate analysis shows that, for high $`q`$, the gap $`d\sigma _c(q)/dq`$ grows slowly with $`q`$, beyond the above estimate of $`D_c`$. This multifractal character can be embodied in the more general scaling ansatz

$$C(t=L^\alpha ,L)=L^{f_c(\alpha )\alpha },\text{(}\text{L≫1}\text{)},$$
(6)
with a nonlinear singularity spectrum $`f_c>\mathrm{}`$ in an $`\alpha `$-interval covering the whole range of possible gaps and linked to $`\sigma _c`$ by Legendre transform. $`C`$ would satisfy the FFS ansatz (2) only if $`f_c`$ were a linear function, i.e. $`f_c(\alpha )=(\tau _c1)\alpha `$ if $`\alpha [0,D_c]`$, $`f_c(\alpha )=\mathrm{}`$ otherwise. The FSS picture given above is in fact only an approximation. This explains also the slight curvature of the $`F`$ plot for low $`t`$, which makes the direct measurement of $`H`$ ambiguous (Fig. 3). Rather than attempting a more precise determination of $`\sigma _c`$ and $`f_c`$, below we clarify the difference between BTW and M, and the origin of multiscaling in the former, in the light of probabilistic concepts.

Waves have a relatively simple behavior. So, in a renormalization group (RG) spirit, we can coarse grain the time, by looking at the PDF, $`P^{(n)}(s,L)`$, of the sum of the sizes of $`n`$ consecutive waves, regardless of the avalanche they belong to. Since avalanches are constituted by an infinite number  of waves for $`L\mathrm{}`$, by sending also $`n\mathrm{}`$, we expect $`P^{(n)}`$ to approach the PDF of avalanche sizes as a RG fixed point. This approach can be monitored on the basis of the effective moment scaling exponents defined by

$$s^q_{n,L}=𝑑ss^qP^{(n)}(s,L)L^{\sigma _{s,n}(q)}.$$
(7)
For the M, $`\sigma _{s,n}(q)`$ does not depend on $`n`$, and, within our accuracy, is equal from the start to $`\sigma _s(q)`$, such that $`𝑑ss^qP(s,L)L^{\sigma _s(q)}`$, with $`P(s,L)`$ representing the avalanche size PDF. To the contrary, in the BTW, as $`n`$ increases, $`\sigma _{s,n}(q)`$ varies and moves gradually towards the appropriate $`\sigma _s(q)`$ (Fig. 4).

The result for the M can be explained on the basis of the fact that PDF’s of independent variables, satisfying a scaling ansatz of the type (6), have a spectrum which does not change under convolution . For example, if we put $`P_w(s=L^\alpha ,L)L^{f_w(\alpha )\alpha }`$, where $`f_w`$ is the spectrum of the wave size PDF, one can verify that

$`P^{(2)}(L^\alpha ,L)=`$ (9)
$`{\displaystyle d(L^\beta )P_w(L^\alpha L^\beta ,L)P_w(L^\beta ,L)}L^{f_w(\alpha )\alpha }.`$
Thus, also the spectrum associated with $`P^{(n)}(s,L)`$, which is the convolution of $`n`$ $`P_w`$’s, does not depend on $`n`$. This implies the $`n`$–independence of $`\sigma _{s,n}(q)`$, which is determined once $`f_w`$ is given. Now, let $`P(s,m,L)`$ be the probability of having an avalanche with $`s`$ topplings and $`m`$ waves: due to the uncorrelated character of different wave sizes, for the M $`P(s,m,L)`$ is also the convolution of $`m`$ $`P_w`$’s, i.e $`P(s,m,L)=P^{(m)}(s,L)`$. Therefore, if $`P(m,L)`$ is the PDF of the total number, $`m`$, of waves in an avalanche, one has

$$P(s,L)=\underset{m}{}P^{(m)}(s,L)P(m,L),$$
(10)
such that

$$L^{\sigma _s(q)}\underset{m}{}L^{\sigma _{s,m}(q)}P(m,L),$$
(11)
which, together with the above results, imply $`\sigma _{s,n}(q)=\sigma _s(q)`$. The coarse graining of waves in the M does not modify the block PDF spectrum, which is from the start at its fixed point, representing also the scaling of avalanches.

In the BTW case, the nontrivial composition of correlated waves is responsible for the change with $`n`$ of the effective singularity spectrum of $`P^{(n)}`$ (Fig. 4). For $`n=1`$, the simple linear FSS form $`f_w(\alpha )=(\tau _w1)\alpha `$ ($`\alpha [0,D_w]`$) applies, while for $`n\mathrm{}`$ one should recover the nonlinear form needed to describe the complex scaling of $`P(s,L)`$ . In practice, sampling limitations prevent us from reaching very high $`n`$. However, even if convergence is relatively slow, the tendency of $`\sigma _{s,n}(q)`$ to move towards $`\sigma _s(q)`$ is very manifest. We could clearly detect the increase with $`n`$ of the gap $`d\sigma _{s,n}(q)/dq`$, at fixed high $`q`$: for example, we determined $`d\sigma _{s,n}(4)/dq2.00,2.20,2.25,2.38`$, for $`n=1,8,12,24`$, respectively. Furthermore, while the asymtoptic gap for $`P_w(s,L)`$ gets readily to the maximum value, i.e. $`d\sigma _{s,1}(q)/dq=2`$ for $`q1`$, as soon as $`n>1`$ a constant gap could not be detected for $`\sigma _{s,n}(q)`$; this confirms the progressive appearance of effective multifractal scaling for the block variable.

In summary, a comparative analysis of wave time series shows that the different forms of universal scaling in the BTW and M are determined by distinct dynamical correlation patterns. For the M case, the wave level of description is in fact coarse grained enough to fully account, without further modifications, also for the avalanche level. Waves are uncorrelated and their PDF, satisfying FSS, coincides with the avalanche PDF, as far as exponents are concerned. To the contrary, in the BTW, under coarse graining, long time correlations substantially modify the scaling properties of waves, determining also multiscaling features. We regard these as prototype mechanisms to be generally expected in sandpile and similar SOC models. The 2D BTW behavior is probably less generic than the M one. For example, in the 3D BTW, the fact that avalanches are constituted by an average number of waves which remains finite for $`L\mathrm{}`$, strongly suggests a M type mechanism and FSS for the avalanche PDF. Indeed, even if subsequent waves are not strictly uncorrelated, like in the 2D Manna model, $`C`$ should decay extremely fast, with an $`L`$ independent time cutoff. Thus, coarse graining in time can not substantially modify the PDF of block variables, with respect to $`P_w`$. We expect the exactly known wave exponents to apply also to the avalanche size PDF. Indeed, numerically determined avalanche exponents for 3D BTW turn out to be strikingly close to the wave values .

We thank C. Tebaldi and C. Vanderzande for useful discussions. Partial support from the European Network Contract No. ERBFMRXCT980183 is acknowledged.